Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 19 Jan 2018 23:14:05 +0000 (23:14 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 19 Jan 2018 23:14:05 +0000 (23:14 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15347 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/SymmsIntegralMapping.java

index ba4b1e1..ad7775e 100644 (file)
@@ -244,14 +244,44 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     private double COARG[] = new double[NJIND];\r
     private double PHPAS[] = new double[IBNDS[0]];\r
     private double VARGC[] = new double[IBNDS[0]];\r
-       //            CSNCA  - COMPLEX ARRAY\r
+       //            ZSNCA  - COMPLEX ARRAY\r
+    //            Also called CSNCA elsewhere\r
        //                     A COMPLEX VECTOR OF SIZE AT LEAST 2*IBNDS(2) + 1;\r
        //                     STORES THE JACOBI COEFFICIENTS FOR THE COMPLEX\r
        //                     (INVERSE) BOUNDARY CORRESPONDENCE FUNCTION AND\r
        //                     ITS DERIVATIVE.\r
-    private double CSNCA[] = new double[2]; // At first location of CSNCA\r
+    private double ZSNCA[] = new double[2]; // At first location of ZSNCA\r
     private double BFSNC[][] = new double[IBNDS[1]][2];\r
     private double SOLNC[][] = new double[IBNDS[1]][2];\r
+    // From GQCANP:\r
+    //            MQIN1  - INTEGER\r
+       //                     DEFINES THE NUMBER OF PANELS ALLOWED IN A \r
+       //                     COMPOSITE RULE.  SPECIFICALLY, MQIN1 = 1 + (THE\r
+       //                     MAXIMUM NUMBER OF PANELS IN A COMPOSITE RULE FOR\r
+       //                     A SINGLE SUB-ARC ON THE BOUNDARY)\r
+       private int MQIN1;\r
+       //            MQUCA  - INTEGER\r
+       //                     THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
+       //                     IN THE FINAL GLOBAL RULE.  THE VALUE OF THIS\r
+       //                     ARGUMENT IS LINKED TO THOSE OF ARGUMENTS NQPTS\r
+       //                     AND IBNDS(1) PREVIOUSLY SUPPLIED TO JACANP VIA\r
+       //                     MQUCA <= (MQIN1-1)*NQPTS*IBNDS(1).  (NOTE THAT\r
+    //                     NQPTS = ISNCA(2) 'JACANP'IBNDS(1) =ISNCA(5) )\r
+       private int MQUCA = (MQIN1-1)*NQPTS*IBNDS[0];\r
+    //            IQUCA  - INTEGER ARRAY\r
+       //                     AN INTEGER VECTOR OF SIZE AT LEAST 2*IBNDS(1) + 4,\r
+       //                     WHERE IBNDS(1) (=ISNCA(5)) IS THE VALUE PREVIOUSLY\r
+       //                     SUPPLIED TO JACANP;  IQUCA MAINLY STORES POINTERS\r
+       //                     TO ACCESS ZQUCA.\r
+    private int IQUCA[] = new int[4];\r
+    private int LQSBG[] = new int[IBNDS[0]];\r
+    private int NPPQG[] = new int[IBNDS[0]];\r
+       //            ZQUCA  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST 2*MQUCA+1;\r
+       //                     STORES THE QUADRATURE POINTS AND WEIGHTS.\r
+    private double ZQUCA[] = new double[2]; // At first location of ZQUCA\r
+    private double WPPQG[][] = new double[MQUCA][2];\r
+    private double ZPPQG[][] = new double[MQUCA][2];\r
        public SymmsIntegralMapping() {\r
                \r
        }\r
@@ -10545,7 +10575,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                //                     AND THE ARGUMENTS OF SUB-ARC ENDPOINTS ON THE UNIT\r
                //                     DISC.\r
                \r
-               //            CSNCA  - COMPLEX ARRAY\r
+               //            ZSNCA  - COMPLEX ARRAY\r
                //                     A COMPLEX VECTOR OF SIZE AT LEAST 2*IBNDS(2) + 1;\r
                //                     STORES THE JACOBI COEFFICIENTS FOR THE COMPLEX\r
                //                     (INVERSE) BOUNDARY CORRESPONDENCE FUNCTION AND\r
@@ -10751,8 +10781,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    DCAP0[0] = DCAP0[0] * mult;\r
                    DCAP0[1] = DCAP0[1] * mult;\r
                }\r
-               CSNCA[0]=DCAP0[0];\r
-               CSNCA[1]=DCAP0[1];\r
+               ZSNCA[0]=DCAP0[0];\r
+               ZSNCA[1]=DCAP0[1];\r
                \r
                // GET THE COEFFICIENTS BFSNC FOR THE COMPLEX BOUNDARY CORRESPONDENCE\r
                // FUNCTION FOR THE INVERSE MAP.\r
@@ -10767,6 +10797,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                // OUTPUT ALL RESULTS REQUIRED FOR SUBSEQUENT PROCESSING\r
                \r
                // CALL OUPTCA(ISNCA,RSNCA,ZSNCA,CHNL)\r
+               NJIND=ISNCA[0]+1;\r
+           TNGQP=ISNCA[1]*NJIND;\r
+\r
                \r
                WRTAIL(3,0,IER[0],null);\r
            return;\r
@@ -12025,34 +12058,871 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                   Preferences.debug("INNER RADIUS = " + INNRAD + "\n", Preferences.DEBUG_ALGORITHM);\r
           }\r
                        \r
-       /*            WRITE(OC,'(/,A,/)') 'JACOBI COEFFICIENTS FOR INVERSE DENSITY FUNCT\r
-                    +IONS'\r
-                     DO 20 IC=1,TNSUC\r
-                       DG=DGPOC(IC)\r
-                       PSA=PRNSA(IC)\r
-                       PT=PARNT(PSA)\r
-                       WRITE(OC,*)\r
-                       WRITE(OC,101) IC,PSA,PT\r
-                       LOM=LSUBC(IC)\r
-                       AJT=ABS(JTYPC(IC))\r
-                       LOD=(AJT-1)*NQPTS+1\r
-                       DO 10 I=0,DG\r
-                         COF=SOLNC(LOM+I)\r
-                         MOD=ABS(COF)\r
-                         WRITE(OC,102) I,COF,MOD,LGTOL/RIGLL(LOD+I)\r
-               10      CONTINUE\r
-               20    CONTINUE\r
-               C\r
-               101   FORMAT('SUB ARC = ',I3,'; PHYSICAL PARENTAL SUB ARC = ',I3,' ON GL\r
-                    +OBAL ARC ',I2/,' N',T6,'REAL PART',T24,'IMAGINARY PART',T42,'MODUL\r
-                    +US',T56,'IGNORE LVL')\r
-               102   FORMAT(I2,T4,E16.8,T22,E16.8,T40,E11.3,T54,E11.3)\r
-               C\r
-                     CALL WRTAIL(3,OC,IER)\r
-                     CLOSE(OC)\r
-               C */\r
+          Preferences.debug("JACOBI COEFFICIENTS FOR INVERSE DENSITY FUNCTIONS\n", Preferences.DEBUG_ALGORITHM);\r
+          for (IC=1; IC <= TNSUC; IC++) {\r
+              DG=DGPOC[IC-1];\r
+                  PSA=PRNSA[IC-1];\r
+                  PT=PARNT[PSA-1];\r
+                  Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+                  Preferences.debug("SUB ARC = " + IC +"; PHYSICAL PARENTAL SUN ARC = " + PSA + " ON GLOBAL ARC " + PT + "\n",\r
+                                  Preferences.DEBUG_ALGORITHM);\r
+                  Preferences.debug(" N   REAL PART         IMAGINARY PART    MODULUS       IGNORE LVL\n", Preferences.DEBUG_ALGORITHM);\r
+                  LOM=LSUBC[IC-1];\r
+                  AJT=Math.abs(JTYPC[IC-1]);\r
+                  LOD=(AJT-1)*NQPTS+1;\r
+                  for (I=0; I <= DG; I++) {\r
+                      COF[0]=SOLNC[LOM+I-1][0];\r
+                      COF[1]=SOLNC[LOM+I-1][1];\r
+                      MOD=zabs(COF[0],COF[1]);\r
+                      Preferences.debug(I + " " + COF[0] + " " + COF[1] + " " + MOD + " " + (LGTOL/RIGLL[LOD+I-1]) + "\n",\r
+                                  Preferences.DEBUG_ALGORITHM);\r
+                  } // for (I=0; I <= DG; I++)\r
+          } // for (IC=1; IC <= TNSUC; IC++)\r
+               \r
+               WRTAIL(3, 0,IER, null);\r
+                    \r
    } // private void OUPTCL\r
 \r
+   private void GQCANP(double RWORK[], int IER[]) {\r
+               \r
+       // INTEGER MQIN1,MQUCA,CHNL,IER\r
+          // INTEGER IQUCA(*),ISNCA(*)\r
+          // REAL RSNCA(*),RWORK(*)\r
+          // COMPLEX ZQUCA(*),ZSNCA(*)\r
+               \r
+          //  ......................................................................\r
+       \r
+          // 1.     GQCANP\r
+          //           COMPUTES A GLOBAL QUADRATURE RULE FOR APPROXIMATING THE\r
+          //           BOUNDARY INTEGRAL REPRESENTATION OF THE MAP : CANONICAL -->\r
+          //           PHYSICAL.\r
+               \r
+          // 2.     PURPOSE\r
+          //           THE ROUTINE SETS UP THE BOUNDARY QUADRATURE POINTS AND\r
+          //           CORRESPONDING WEIGHTS FOR A COMPOSITE GAUSS-JACOBI/GAUSS-\r
+          //           LEGENDRE RULE FOR ESTIMATING THE BOUNDARY INTEGRAL THAT\r
+          //           APPEARS IN THE RESPRESENTATION FOR THE CONFORMAL MAP OF THE\r
+          //           CANONICAL DOMAIN ONTO THE PHYSICAL DOMAIN.  THIS QUADRATURE\r
+          //           RULE IS USED IN THE STANDARD NON-SINGULAR CASE WHEN THE\r
+          //           FIELD POINT IN THE CANONICAL DOMAIN DOES NOT LIE CLOSE TO \r
+          //           THE UNIT CIRCLE.\r
+               \r
+          // 3.     CALLING SEQUENCE\r
+          //           CALL GQCANP(MQIN1,MQUCA,ISNCA,RSNCA,ZSNCA,RWORK,CHNL,IQUCA,\r
+          //                       ZQUCA,IER)\r
+               \r
+          //        PARAMETERS\r
+          //         ON ENTRY\r
+          //            MQIN1  - INTEGER\r
+          //                     DEFINES THE NUMBER OF PANELS ALLOWED IN A \r
+          //                     COMPOSITE RULE.  SPECIFICALLY, MQIN1 = 1 + (THE\r
+          //                     MAXIMUM NUMBER OF PANELS IN A COMPOSITE RULE FOR\r
+          //                     A SINGLE SUB-ARC ON THE BOUNDARY)\r
+               \r
+          //            MQUCA  - INTEGER\r
+          //                     THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
+          //                     IN THE FINAL GLOBAL RULE.  THE VALUE OF THIS\r
+          //                     ARGUMENT IS LINKED TO THOSE OF ARGUMENTS NQPTS\r
+          //                     AND IBNDS(1) PREVIOUSLY SUPPLIED TO JACANP VIA\r
+          //                     MQUCA <= (MQIN1-1)*NQPTS*IBNDS(1).  (NOTE THAT\r
+          //                     NQPTS = ISNCA(2) 'JACANP'IBNDS(1) =ISNCA(5) )\r
+               \r
+          //            ISNCA  - INTEGER ARRAY\r
+          //                     THE INTEGER VECTOR PREVIOUSLY SET UP BY JACANP.\r
+               \r
+          //            RSNCA  - REAL ARRAY\r
+          //                     THE REAL VECTOR PREVIOUSLY SET UP BY JACANP.\r
+               \r
+          //            ZSNCA  - COMPLEX ARRAY\r
+          //                     THE COMPLEX VECTOR PREVIOUSLY SET UP BY JACANP.\r
+               \r
+          //            RWORK  - REAL ARRAY\r
+          //                     A WORKING VECTOR OF SIZE AT LEAST MQIN1.\r
+               \r
+          //             CHNL  - INTEGER\r
+          //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+          //                     WRITING THE FILE <JBNM>cq.\r
+               \r
+          //         ON EXIT\r
+          //            IQUCA  - INTEGER ARRAY\r
+          //                     AN INTEGER VECTOR OF SIZE AT LEAST 2*IBNDS(1) + 4,\r
+          //                     WHERE IBNDS(1) (=ISNCA(5)) IS THE VALUE PREVIOUSLY\r
+          //                     SUPPLIED TO JACANP;  IQUCA MAINLY STORES POINTERS\r
+          //                     TO ACCESS ZQUCA.\r
+               \r
+          //            ZQUCA  - COMPLEX ARRAY\r
+          //                     A COMPLEX VECTOR OF SIZE AT LEAST 2*MQUCA+1;\r
+          //                     STORES THE QUADRATURE POINTS AND WEIGHTS.\r
+               \r
+          //            IER    - INTEGER\r
+          //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+          //                     A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
+          //                     WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+          //                     IER=0 - NORMAL EXIT.\r
+          //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+          //                             BE SELF EXPLANATORY.\r
+               \r
+               \r
+          // 4.     SUBROUTINES OR FUNCTIONS NEEDED\r
+          //              - THE CONFPACK LIBRARY.\r
+          //              - THE REAL FUNCTION R1MACH.\r
+          //              - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN.\r
+               \r
+               \r
+          // 5.     FURTHER COMMENTS\r
+          //           - NOTE THAT THIS ROUTINE CAN ONLY BE USED  A F T E R  THE  \r
+          //             ROUTINE JACANP HAS SUCCESSFULLY EXECUTED, AND THAT SOME  \r
+          //             INPUT ARGUMENTS FOR GQCANP ARE OUTPUT VALUES FROM JACANP.\r
+          //           - THE GLOBAL QUADRATURE DATA ARE AUTOMATICALLY OUTPUT ON THE\r
+          //             FILE <JBNM>cq, WHERE <JBNM> IS COLLECTED FROM THE FILE \r
+          //             jbnm PREVIOUSLY CREATED BY JAPHYC.\r
+          //           - A SUMMARY LISTING OF ACTIONS TAKEN IS AUTOMATICALLY\r
+          //             WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+               \r
+          // ......................................................................\r
+          //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+          //     LAST UPDATE: 3 JULY 1990\r
+          // ......................................................................\r
+                    \r
+          //     LOCAL VARIABLES\r
+               \r
+          final int NINTS = 5;\r
+          int TNPQP[] = new int[1];\r
+          int MNCOF,MNSUC,TNSUC;\r
+          final double DELTA = 0.2;\r
+          double TOLOU[] = new double[1];\r
+          double LGTOL;\r
+               \r
+          // EXTERNAL POPQG1,OUPTCQ\r
+               \r
+          // **** WRITE HEADING TO STANDARD OUTPUT CHANNEL\r
+               \r
+          WRHEAD(4,0,null);\r
+               \r
+          NARCS=ISNCA[0];\r
+          NQPTS=ISNCA[1];\r
+          TNSUC=ISNCA[2];\r
+          MNSUC=ISNCA[4];\r
+          MNCOF=ISNCA[5];\r
+          LGTOL=RSNCA[0];\r
+          ZQUCA[0]=ZSNCA[0];\r
+          ZQUCA[1]=ZSNCA[1];\r
+               \r
+          NJIND=NARCS+1;\r
+          TNGQP=NJIND*NQPTS;\r
+               \r
+          IQUCA[1]=TNSUC;\r
+          IQUCA[2]=MNSUC;\r
+          IQUCA[3]=MQUCA;\r
+               \r
+          // COPY RELEVANT POINTERS TO ISNCA, RSNCA AND ZSNCA FROM JACANP\r
+               \r
+          // SET UP POINTERS FOR QUADRATURE ARRAYS IQUCA AND ZQUCA\r
+               \r
+          System.out.println("QUADRATURE RULES STARTED:");\r
+          POPQG1(NPPQG,LQSBG,TNPQP,TOLOU,WPPQG,\r
+                    ZPPQG,MQUCA,MQIN1,NARCS,NINTS,NQPTS,TNSUC,DGPOC,\r
+                    JTYPC,LSUBC,DELTA,LGTOL,ACOFC,BCOFC,\r
+                    H0VLC,JAINC,QUPTC,QUWTC,SOLNC,\r
+                    VARGC,RWORK,IER);\r
+          IQUCA[0]=TNPQP[0];\r
+          System.out.println("QUADRATURE RULES DONE:");\r
+               \r
+          // WRITE CLOSING MESSAGE TO STANDARD OUTPUT CHANNEL\r
+               \r
+          WRTAIL(4,0,IER[0], null);\r
+               \r
+    } // private void GQCANP\r
+\r
+    private void POPQG1(int NPPQG[], int LPQSB[], int TNPQP[], double TOLOU[], double WPPQG[][],\r
+        double ZPPQG[][], int MNQUA, int MQIN1, int NARCS, int NINTS, int NQPTS, int TNSUC,\r
+        int DGPOC[], int JTYPC[], int LSUBC[], double DELTA, double LGTOL, double ACOFC[],\r
+               double BCOFC[], double H0VLC[], double JAINC[], double QUPTC[], double QUWTC[],\r
+               double SOLNC[][], double VARGC[], double XENPT[], int IER[]) {\r
+               // INTEGER IER,MQIN1,MNQUA,NARCS,NINTS,NQPTS,TNPQP,TNSUC\r
+               // INTEGER DGPOC(*),JTYPC(*),LSUBC(*),LPQSB(*),NPPQG(*)\r
+               // REAL DELTA,LGTOL,TOLOU\r
+               // REAL ACOFC(*),BCOFC(*),H0VLC(*),JAINC(*),QUPTC(*),QUWTC(*),\r
+               // +VARGC(*),XENPT(*)\r
+               // COMPLEX SOLNC(*),WPPQG(*),ZPPQG(*)\r
+               \r
+           //      THE MAIN PURPOSE OF THIS ROUTINE IS TO SET UP THE COMPLEX ARRAY\r
+               //      WPPQG OF QUADRATURE WEIGHTS AND THE COMPLEX ARRAY ZPPQG OF\r
+               //      QUADRATURE POINTS ON THE PHYSICAL BOUNDARY; ALL THESE\r
+               //      DATA DEFINE THE POST-PROCESSING COMPOSITE GAUSSIAN\r
+               //      QUADRATURE RULES FOR THE ESTIMATE OF THE MAP G (CANONICAL ONTO\r
+               //      PHYSICAL) WHEN Z IS NOT TOO CLOSE TO THE BOUNDARY.\r
+               \r
+               //     THE ARRAY ELEMENT NPPQG(I) RECORDS THE NUMBER OF QUADRATURE\r
+               //     POINTS FOR THE SUBARC NUMBER I, I=1,...,TNSUC.  THE WEIGHTS AND\r
+               //     QUADRATURE POINTS FOR ARC NUMBER I ARE STORED IN WPPQG AND ZPPQG\r
+               //     STARTING AT THE ELEMENT WITH INDEX LPQSB(I).\r
+               \r
+               //     IER=0  - NORMAL EXIT \r
+               //     IER=41 - THE TOTAL NUMBER OF QUADRATURE POINTS REQUESTED FOR \r
+               //              THE WHOLE BOUNDARY EXCEEDS THE LIMIT DEFINED BY THE\r
+               //              INPUT PARAMETER MNQUA; MNQUA MUST BE INCREASED IN A\r
+               //              HIGHER LEVEL ROUTINE.\r
+               //     IER=42 - THE LOCAL PARAMETER MNCOF SHOULD BE INCREASED\r
+               //     IER=43 - THE REQUIRED NUMBER OF QUADRATURE INTERVALS EXCEEDS THAT\r
+               //              SPECIFIED BY THE GLOBAL PARAMETER MQIN1; MQIN1 MUST BE\r
+               //              INCREASED.\r
+               \r
+               //     LOCAL VARIABLES\r
+               \r
+               final int MNCOF = 32;\r
+               int QINTS[] = new int[1];\r
+       int AJT,DEG,HI,HI1,I,I1,J,J1,J2,JT,K,LIM,LOD,LOL,LOM,M;\r
+       final double RHO = 0.13;\r
+               double BETA,HA,MD,MEAN,RR,RRB,SS,SUM1,TT,SCO;\r
+               double CTT[] = new double[2];\r
+               // COMPLEX CTT,DPARFN,JACSUC,PARFUN\r
+               double JACOF[][] = new double[MNCOF][2];\r
+               // COMPLEX JACOF(MNCOF)\r
+               // EXTERNAL DPARFN,JACSUC,PARFUN,PPSBI7\r
+               double JOUT[];\r
+               \r
+               HI=0;\r
+               LOL=NARCS*NQPTS;\r
+               for (I1=1; I1 <= TNSUC; I1++) {\r
+                   JT=JTYPC[I1-1];\r
+                   if (JT > 0) {\r
+                       SS=1.0;\r
+                   }\r
+                   else {\r
+                       SS=-1.0;\r
+                   }\r
+                   AJT=Math.abs(JT);\r
+                   BETA=JAINC[AJT-1];\r
+                   DEG=DGPOC[I1-1];\r
+                   if (DEG+1 > MNCOF) {\r
+                       IER[0]=42;\r
+                       return;\r
+                   }\r
+                   LOM=LSUBC[I1-1];\r
+                   LOD=(AJT-1)*NQPTS+1;\r
+                   HA=(VARGC[I1]-VARGC[I1-1])*0.5;\r
+                   MD=(VARGC[I1]+VARGC[I1-1])*0.5;\r
+               \r
+                   SCO=SS;\r
+                   for (J=1; J <= DEG+1; J++) {\r
+                       J1=LOM+J-1;\r
+                       SCO=SCO*SS;\r
+                       JACOF[J-1][0]=SOLNC[J1-1][0]*SCO;\r
+                       JACOF[J-1][1]=SOLNC[J1-1][1]*SCO;\r
+                   } // for (J=1; J <= DEG+1; J++)\r
+               \r
+                   double ACOEF[] = new double[DEG];\r
+                   double BCOEF[] = new double[DEG];\r
+                   for (I = 0; I < DEG; I++) {\r
+                       ACOEF[I] = ACOFC[LOD+I-1];\r
+                       BCOEF[I] = BCOFC[LOD+I-1];\r
+                   }\r
+                   PPSBI7(DELTA,NINTS,BETA,NQPTS,DEG,ACOEF,\r
+                   BCOEF,H0VLC[AJT-1],JACOF,LGTOL,TOLOU,XENPT,\r
+                      QINTS,MQIN1,IER);\r
+                   if (IER[0] > 0) {\r
+                       if (IER[0] == 24) {\r
+                           IER[0]=43;\r
+                       }\r
+                       return;\r
+                   } // if (IER[0] > 0)\r
+                   NPPQG[I1-1]=QINTS[0]*NQPTS;\r
+                   LPQSB[I1-1]=HI+1;\r
+                   HI1=HI+NPPQG[I1-1];\r
+                   if (HI1 > MNQUA) {\r
+                       IER[0]=41;\r
+                       return;\r
+                   }\r
+                   K=HI;\r
+                   SUM1=BETA+1.0;\r
+                   for (I=1; I <= QINTS[0]; I++) {\r
+                       RR=(XENPT[I]-XENPT[I-1])*0.5;\r
+                       MEAN=(XENPT[I]+XENPT[I-1])*0.5;\r
+                       if (I == 1) {\r
+                           RRB=Math.pow(RR,SUM1);\r
+                           for (J=1; J <= NQPTS; J++) {\r
+                               J1=LOD+J-1;\r
+                               K=K+1;\r
+                               TT=(MEAN+RR*QUPTC[J1-1]);\r
+                               double A[] = new double[DEG];\r
+                               double B[] = new double[DEG];\r
+                               for (M = 0; M < DEG; M++) {\r
+                                       A[M] = ACOFC[LOD+M-1];\r
+                                       B[M] = BCOFC[LOD+M-1];\r
+                               }\r
+                               JOUT = JACSUC(TT,DEG,A,B,H0VLC[AJT-1],JACOF);\r
+                               WPPQG[K-1][0]=RRB*QUWTC[J1-1]*JOUT[0];\r
+                               WPPQG[K-1][1]=RRB*QUWTC[J1-1]*JOUT[1];\r
+                               TT=TT*SS;\r
+                               TT=MD+TT*HA;\r
+                               ZPPQG[K-1][0]=Math.cos(TT);\r
+                               ZPPQG[K-1][1]=Math.sin(TT);\r
+                           } // for (J=1; J <= NQPTS; J++)\r
+                       } // if (I == 1)\r
+                       else { // I != 1\r
+                               for (J=1; J <= NQPTS; J++) {\r
+                               J1=LOL+J;\r
+                               K=K+1;\r
+                               TT=(MEAN+RR*QUPTC[J1-1]);\r
+                               double A[] = new double[DEG];\r
+                               double B[] = new double[DEG];\r
+                               for (M = 0; M < DEG; M++) {\r
+                                       A[M] = ACOFC[LOD+M-1];\r
+                                       B[M] = BCOFC[LOD+M-1];\r
+                               }\r
+                               JOUT = JACSUC(TT,DEG,A,B,H0VLC[AJT-1],JACOF);\r
+                               WPPQG[K-1][0]=RR*QUWTC[J1-1]*Math.pow((1.0+TT),BETA)*JOUT[0];\r
+                               WPPQG[K-1][1]=RR*QUWTC[J1-1]*Math.pow((1.0+TT),BETA)*JOUT[1];\r
+                               TT=TT*SS;\r
+                               TT=MD+TT*HA;\r
+                               ZPPQG[K-1][0]=Math.cos(TT);\r
+                               ZPPQG[K-1][1]=Math.sin(TT);\r
+                               } // for (J=1; J <= NQPTS; J++)\r
+                       } // else I != 1\r
+                   } // for (I=1; I <= QINTS[0]; I++)\r
+                   if (SS < 0.0) {\r
+                       LIM=NPPQG[I1-1];\r
+                       if ((LIM%2) == 0) {\r
+                           LIM=LIM/2;\r
+                       }\r
+                       else {\r
+                           LIM=(LIM-1)/2;\r
+                       }\r
+                       J1=LPQSB[I1-1]-1;\r
+                       J2=HI1+1;\r
+                       for (J=1; J <= LIM; J++) {\r
+                           J1=J1+1;\r
+                           J2=J2-1;\r
+                           CTT[0]=WPPQG[J1-1][0];\r
+                           CTT[1]=WPPQG[J1-1][1];\r
+                           WPPQG[J1-1][0]=WPPQG[J2-1][0];\r
+                           WPPQG[J1-1][1]=WPPQG[J2-1][1];\r
+                           WPPQG[J2-1][0]=CTT[0];\r
+                           WPPQG[J2-1][1]=CTT[1];\r
+                           CTT[0]=ZPPQG[J1-1][0];\r
+                           CTT[1]=ZPPQG[J1-1][1];\r
+                           ZPPQG[J1-1][0]=ZPPQG[J2-1][0];\r
+                           ZPPQG[J1-1][1]=ZPPQG[J2-1][1];\r
+                           ZPPQG[J2-1][0]=CTT[0];\r
+                           ZPPQG[J2-1][1]=CTT[1];\r
+                       } // for (J=1; J <= LIM; J++)\r
+                   } // if (SS < 0.0)\r
+                   HI=HI1;\r
+               } // for (I1=1; I1 <= TNSUC; I1++)\r
+               \r
+               TNPQP[0]=HI;\r
+               \r
+               IER[0]=0;\r
+               \r
+    } // private void POPQG1\r
+    \r
+    // COMPLEX FUNCATION\r
+    private double[] JACSUC(double X, int N, double A[], double B[],\r
+        double H, double CO[][]) {\r
+        // INTEGER N\r
+        // REAL A(*),B(*),H,X\r
+        // COMPLEX CO(*)\r
+\r
+        // ..TO CALCULATE SUMMATION{CO(K+1)*P(K,X)},K=0(1)N, WHERE P(K,X)\r
+        // ..DENOTES THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE K\r
+        // ..EVALUATED AT X, ARRAY CO STORES A GIVEN SET OF COEFFICIENTS,\r
+        // ..ARRAYS A,B STORE THE COEFFICIENTS IN THE THREE-TERM\r
+        // ..RECURRENCE FORMULA FOR THE JACOBI POLYNOMIALS (SEE ASONJ7)\r
+        // ..AND H IS THE SQUARED 2-NORM OF UNITY.\r
+\r
+        double PREV[] = new double[2];\r
+        double CURR[] = new double[2];\r
+        double NEXT[] = new double[2];\r
+       // COMPLEX PREV,CURR,NEXT\r
+        int K;\r
+        double result[] = new double[2];\r
+\r
+        if (N == 0) {\r
+            result[0]=CO[0][0]/Math.sqrt(H);\r
+            result[1]=CO[0][1]/Math.sqrt(H);\r
+        }\r
+        else if (N > 0) {\r
+            PREV[0]=CO[N][0];\r
+            PREV[1]=CO[N][1];\r
+            CURR[0]=CO[N-1][0]+(X-B[N-1])*PREV[0]/A[N-1];\r
+            CURR[1]=CO[N-1][1]+(X-B[N-1])*PREV[1]/A[N-1];\r
+            for (K=N-2; K >= 0; K--) {\r
+                NEXT[0]=CO[K][0]+(X-B[K])*CURR[0]/A[K]-A[K]*PREV[0]/A[K+1];\r
+                NEXT[1]=CO[K][1]+(X-B[K])*CURR[1]/A[K]-A[K]*PREV[1]/A[K+1];\r
+                PREV[0]=CURR[0];\r
+                PREV[1]=CURR[1];\r
+                CURR[0]=NEXT[0];\r
+                CURR[1]=NEXT[1];\r
+            } //for (K=N-2; K >= 0; K--) \r
+            result[0]=CURR[0]/Math.sqrt(H);\r
+            result[1]=CURR[1]/Math.sqrt(H);\r
+        } // else if (N > 0)\r
+        else {\r
+            result[0] = 0.0;\r
+            result[1] = 0.0;\r
+        }\r
+        return result; \r
+    } // private double[] JACSUC\r
+\r
+    private void CNDPLT(boolean MAP11[], double RESMN[], double UPHYC[], double UCANP[], double CRRES,\r
+        String DASH[], String NEWD[], int IER[]) {\r
+       \r
+       // INTEGER CH0,CH1,IER\r
+       // INTEGER IGEOM(*),ISNPH(*)\r
+       // REAL RESMN,UPHYC,UCANP,CRRES\r
+       // REAL RGEOM(*),RSNPH(*)\r
+       // LOGICAL MAP11\r
+       // CHARACTER DASH*(*),NEWD*(*)\r
+       \r
+       // ......................................................................\r
+       \r
+       // 1.     CNDPLT\r
+       //           REPORTS ON THE CONDITION OF THE PROBLEM OF EVALUATING THE\r
+       //           MAPPING FUNCTIONS AND ALSO OUTPUTS DATA FOR GRAPH PLOTTING.\r
+       \r
+       // 2.     PURPOSE\r
+       //           THE ROUTINE COMPUTES CONDITION NUMBERS FOR THE PROBLEMS OF\r
+       //           EVALUATING THE TWO MAPS PHYSICAL --> CANONICAL AND\r
+       //           CANONICAL --> PHYSICAL AND COMPUTES THE ERROR THAT MAY BE\r
+       //           EXPECTED (IN THE WORST CASE) IN THE RANGE OF EACH APPROX-\r
+       //           IMATE MAP FROM A MACHINE PRECISION LEVEL ROUNDING ERROR IN \r
+       //           THE DOMAIN OF EACH MAP.\r
+       //           THE ROUTINE ALSO COMPUTES THE LEAST RESOLUTION OF THE\r
+       //           COMPUTED MAP : PHYSICAL --> CANONICAL OVER ALL SUB-ARCS ON\r
+       //           THE PHYSICAL BOUNDARY.  THE RESOLUTION OF THE MAP FOR ANY\r
+       //           PHYSICAL SUB-ARC IS DEFINED AS THE COMPUTED ANGULAR WIDTH OF\r
+       //           THE IMAGE SUB-ARC ON THE UNIT DISC DIVIDED BY THE ESTIMATED\r
+       //           MAXIMUM ERROR IN THE MODULUS OF THE MAP ON THE GIVEN SUB-\r
+       //           ARC.  A LEAST RESOLUTION OF LESS THAN, SAY, 10 INDICATES\r
+       //           THAT THERE ARE REGIONS OF SEVERE CROWDING AND THAT IT WILL\r
+       //           BE PRACTICALLY IMPOSSIBLE TO COMPUTE THE INVERSE MAP\r
+       //           EVERYWHERE ON THE CANONICAL DOMAIN.\r
+       //           THE ROUTINE ALSO SEARCHES (NOT VERY EXHAUSTIVELY) FOR\r
+       //           CHANGES OF SIGN IN THE COMPUTED BOUNDARY CORRESPONDENCE\r
+       //           DERIVATIVE FOR THE MAP : PHYSICAL --> CANONICAL.  SUCH\r
+       //           SIGN CHANGES MEAN THAT THE COMPUTED MAP IS NOT ONE-TO-ONE\r
+       //           AND HENCE ONE SHOULD EXPECT DIFFICULTIES IN TRYING TO\r
+       //           COMPUTE THE INVERSE MAP : CANONICAL --> PHYSICAL.\r
+       //           FINALLY THREE OUTPUT FILES\r
+       //                <JBNM>cn, <JBNM>p0, <JBNM>p1\r
+       //           ARE WRITTEN.  THE FIRST OF THESE IS A SUMMARY OF THE ABOVE\r
+       //           RESULTS INTENDED TO BE READ BY THE USER.   THE TWO FILES\r
+       //           <JBNM>p0 AND <JBNM>p1 ARE NOT INTENDED TO BE READ BY THE \r
+       //           USER, BUT COULD BE USED TO CREATE PLOTS OF THE BOUNDARY \r
+       //           CORRESPONDENCE FUNCTION AND ITS DERIVATIVE; SEE FURTHER \r
+       //           COMMENTS BELOW.\r
+       \r
+       // 3.     CALLING SEQUENCE\r
+       //           CALL CNDPLT(MAP11,RESMN,UPHYC,UCANP,CRRES,IGEOM,RGEOM,ISNPH,\r
+       //                       RSNPH,CH0,CH1,DASH,NEWD,IER)\r
+       \r
+       //        PARAMETERS\r
+       //         ON ENTRY\r
+       //            CRRES  - REAL\r
+       //                     THE CRITICAL RESOLUTION.  IF THE COMPUTED RESOL-\r
+       //                     UTION OF THE PHYSICAL-->CANONICAL MAP ON ANY ARC\r
+       //                     FALLS BELOW CRRES THAN A WARNING MESSAGE IS\r
+       //                     OUTPUT.  \r
+       \r
+       //            IGEOM  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR IGEOM PREVIOUSLY SET UP BY \r
+       //                     JAPHYC.\r
+       \r
+       //            RGEOM  - REAL ARRAY\r
+       //                     THE REAL VECTOR RGEOM PREVIOUSLY SET UP BY JAPHYC.\r
+       \r
+       //            ISNPH  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR ISNPH PREVIOUSLY SET UP BY \r
+       //                     JAPHYC.\r
+       \r
+       //            RSNPH  - REAL ARRAY\r
+       //                     THE REAL VECTOR RSNPH PREVIOUSLY SET UP BY JAPHYC.\r
+       \r
+       //            CH0    - INTEGER\r
+       //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+       //                     WRITING THE FILES <JBNM>cn AND <JBNM>p0.\r
+       \r
+       //            CH1    - INTEGER\r
+       //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+       //                     WRITING THE FILE <JBNM>p1; MUST HAVE CH0.NE.CH1. \r
+       \r
+       //            DASH   - CHARACTER\r
+       //                     A CHARACTER VARIABLE OF USER-DEFINED LENGTH WHICH\r
+       //                     DEFINES THE DASH-PATTERN THAT THE USER MAY REQUIRE\r
+       //                     FOR GRAPH PLOTTING;  SEE FURTHER COMMENTS BELOW\r
+       \r
+       //            NEWD   - CHARACTER\r
+       //                     A CHARACTER VARIABLE OF USER-DEFINED LENGTH WHICH\r
+       //                     DENOTES THE START OF A NEW DATA GROUP THAT THE \r
+       //                     USER MAY REQUIRE FOR GRAPH PLOTTING;  SEE FURTHER \r
+       //                     COMMENTS BELOW\r
+       //         ON EXIT \r
+       //            MAP11  - LOGICAL\r
+       //                     IF BOUNDARY REVERSALS ARE DETECTED THEN MAP11 IS\r
+       //                     SET TO .FALSE. (THE COMPUTED PHYSICAL-->CANONICAL\r
+       //                     MAP ISN'T 1-1) OTHERWISE MAP11 IS SET TO .TRUE.\r
+       \r
+       //            RESMN  - REAL\r
+       //                     THE MINIMUM COMPUTED RESOLUTION OF THE PHYSICAL-->\r
+       //                     CANONICAL MAP OVER ALL SUBARCS ON THE PHYSICAL\r
+       //                     BOUNDARY.  IF RESMN IS LESS THAN CRRES THEN\r
+       //                     A WARNING MESSAGE IS OUTPUT.\r
+       \r
+       //            UPHYC  - REAL\r
+       //                     ESTIMATED MAXIMUM POSSIBLE ERROR IN THE RANGE OF\r
+       //                     THE PHYSICAL-->CANONICAL MAP DUE UNIT ROUNDOFF IN \r
+       //                     THE PHYSICAL DOMAIN.\r
+       \r
+       //            UCANP  - REAL\r
+       //                     ESTIMATED MAXIMUM POSSIBLE ERROR IN THE RANGE OF\r
+       //                     THE CANONICAL-->PHYSICAL MAP DUE UNIT ROUNDOFF IN \r
+       //                     THE CANONICAL DOMAIN.\r
+       \r
+       //            IER    - INTEGER\r
+       //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+       //                     A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
+       //                     WRITTEN ON THE STANDARD OUTPUT CHANNEL AND THE\r
+       //                     LISTING FILE <JBNM>cn.\r
+       //                     IER=0 - NORMAL EXIT.\r
+       //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+       //                             BE SELF EXPLANATORY.\r
+       \r
+       \r
+       // 4.     SUBROUTINES OR FUNCTIONS NEEDED\r
+       //              - THE CONFPACK LIBRARY.\r
+       //              - THE REAL FUNCTION R1MACH.\r
+       //              - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN.\r
+       \r
+       \r
+       // 5.     FURTHER COMMENTS\r
+       //           - NOTE THAT THIS ROUTINE CAN ONLY BE USED  A F T E R  THE  \r
+       //             ROUTINE JAPHYC HAS SUCCESSFULLY EXECUTED, AND THAT MOST  \r
+       //             INPUT ARGUMENTS FOR CNDPLT ARE OUTPUT VALUES FROM JAPHYC.\r
+       //           - A DETAILED LISTING OF RESULTS IS WRITTEN ON THE FILE\r
+       //             <JBNM>cn.\r
+       //           - DATA FOR PLOTTING A GRAPH OF THE DIMENSIONLESS BOUNDARY\r
+       //             CORRESPONDENCE FUNCTION AGAINST DIMENSIONLESS ARC LENGTH\r
+       //             ARE WRITTEN ON THE FILE <JBNM>p0.  THE CONTENTS OF THIS \r
+       //             FILE ARE AS FOLLOWS:\r
+       //             1.  ABOUT 200 COORDINATE PAIRS X Y, ONE PAIR PER LINE,\r
+       //                 WHERE X = DIMENSIONLESS ARC LENGTH (0 <= X <=1) AND\r
+       //                 Y = DIMENSIONLESS BOUNDARY CORRESPONDENCE FUNCTION\r
+       //                 (0 <= Y <=1);  THE NUMBER OF COORDINATE PAIRS IS \r
+       //                 CONTROLLED BY THE LOCAL PARAMETER NXINT.\r
+       //             2.  THE SINGLE LINE\r
+       //                   <DASH>\r
+       //                 WHERE <DASH> DENOTES THE VALUE OF THE ARGUMENT DASH;\r
+       //                 THIS CAN BE USED TO INDICATE A CHANGE OF DASH PATTERN\r
+       //                 TO THE LOCAL GRAPH PLOTTER.\r
+       //             3.  SEVERAL REPETITIONS OF THE FOLLOWING 3-LINE GROUP:\r
+       //                   <NEWD>\r
+       //                   X 0E+0\r
+       //                   X 1E+0\r
+       //                 HERE <NEWD> DENOTES THE VALUE OF THE ARGUMENT NEWD AND\r
+       //                 X (WITH 0 < X < 1) IS THE DIMENSIONLESS ARC \r
+       //                 LENGTH OF A CORNER POINT.  THE ABOVE GROUP MAY THEN\r
+       //                 BE USED TO CONSTRUCT A DASHED LINE FROM (X,0) TO \r
+       //                 (X,1).  THE NUMBER OF REPETITIONS IS EQUAL TO THE\r
+       //                 NUMBER OF CORNERS WITH ARC LENGTH IN THE INTERVAL\r
+       //                 0 < X < 1.\r
+       //           - DATA FOR PLOTTING A GRAPH OF THE DERIVATIVE OF THE DIMEN-\r
+       //             SIONLESS BOUNDARY CORRESPONDENCE FUNCTION WITH RESPECT TO\r
+       //             DIMENSIONLESS ARC LENGTH ARE WRITTEN ON THE FILE <JBNM>p1.\r
+       //             THE CONTENTS OF THIS FILE ARE AS FOLLOWS:\r
+       //             1.  ABOUT 200 COORDINATE PAIRS X Y, ONE PAIR PER LINE,\r
+       //                 WHERE X = DIMENSIONLESS ARC LENGTH (0 <= X <=1) AND\r
+       //                 Y = DIMENSIONLESS BOUNDARY CORRESPONDENCE DERIVATIVE\r
+       //                 (0 <= Y <=1);  THE NUMBER OF COORDINATE PAIRS IS \r
+       //                 CONTROLLED BY THE LOCAL PARAMETER NXINT.\r
+       //             2.  THE SINGLE LINE\r
+       //                   <DASH>\r
+       //                 WHERE <DASH> DENOTES THE VALUE OF THE ARGUMENT DASH;\r
+       //                 THIS CAN BE USED TO INDICATE A CHANGE OF DASH PATTERN\r
+       //                 TO THE LOCAL GRAPH PLOTTER.\r
+       //             3.  SEVERAL REPETITIONS OF THE FOLLOWING 3-LINE GROUP:\r
+       //                   <NEWD>\r
+       //                   X 0E+0\r
+       //                   X 4.4E+0\r
+       //                 HERE <NEWD> DENOTES THE VALUE OF THE ARGUMENT NEWD AND\r
+       //                 X (WITH 0 < X < 1) IS THE DIMENSIONLESS ARC \r
+       //                 LENGTH OF A RE-ENTRANT CORNER POINT.  THE ABOVE GROUP \r
+       //                 MAY THEN BE USED TO CONSTRUCT A DASHED LINE FROM (X,0)\r
+       //                 TO (X,4.4), TO INDICATE THE PRESENCE OF AN ASYMPTOTE.\r
+       //                 SINCE THE AVERAGE VALUE OF THE DIMENSIONLESS BCF DERI-\r
+       //                 VATIVE IS 1,  4.4 IS AN ARBITRARILY CHOSEN BUT RELATI-\r
+       //                 VELY LARGE HEIGHT AT WHICH TO TERMINATE AN ASYMPTOTE; \r
+       //                 THIS HEIGHT IS CONTROLLED BY THE LOCAL PARAMETER BIG. \r
+       //                 THE NUMBER OF REPETITIONS OF THIS GROUP IS EQUAL TO \r
+       //                 THE NUMBER OF RE-ENTRANT CORNERS WITH DIMENSIONLESS \r
+       //                 ARC LENGTH IN THE INTERVAL 0 < X < 1.\r
+       //           - A SUMMARY LISTING OF RESULTS IS AUTOMATICALLY\r
+       //             WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+       \r
+       // ......................................................................\r
+       //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+       //     LAST UPDATE: 17 JULY 1990\r
+       // ......................................................................\r
+            \r
+       //     LOCAL VARAIBLES\r
+       \r
+       final int NIXINT = 200;\r
+       final int MAXSA = 100;\r
+       int I,IMNLA,J,L,LODP,LODW,MNEQN,MNSUA,NARCS,NASYM,NCRVS,\r
+            NINFD,NJIND,NPRVS,NQPTS,NZERD,TNGQP,TNSUA;\r
+       final double BIG = 4.4;\r
+       double ANGSP,CCAPH,COCAP,COPHC,CPHCA,CR,EXCAP,EXPHC,LA,\r
+            OFLOW,PI,R1MACH,TOTLN,MCHEP;\r
+       //CHARACTER OFLC*6,OFP0*6,OFP1*6,JBNM*4,CHPC*2,CHCP*2\r
+       \r
+       \r
+       //**** NXINT  = GLOBAL NUMBER OF INTERVALS ON [0,1] FOR SAMPLING THE\r
+       //****          DIMENSIONLESS DERIVATE OF THE BOUNDARY CORRESPONDENCE\r
+       //****          FUNCTION.\r
+       //**** IER=52 - LOCAL PARAMETER MAXSA MUST BE INCEASED TO AT LEAST THE\r
+       //****          VALUE OF ARGUMENT ISNPH(3)=TOTAL NUMBER OF SUBARCS ON\r
+       //****          PHYSICAL BOUNDARY.\r
+       \r
+       int ICRVS[] = new int[MAXSA];\r
+       int IPRVS[] = new int[MAXSA];\r
+       double ARCLN[] = new double[MAXSA];\r
+       double ASYMP[] = new double[MAXSA];\r
+       double BCDMN[] = new double[MAXSA];\r
+       double CORXX[] = new double[MAXSA];\r
+       // EXTERNAL DIAGN4,R1MACH,WRHEAD,WRTAIL\r
+       \r
+       //**** WRITE HEADING TO STANDARD OUTPUT CHANNEL\r
+       \r
+       WRHEAD(5,0,null);\r
+       \r
+       // GET JOBNAME FROM FILE *jbnm*\r
+       \r
+        // OPEN(CH0,FILE='jbnm')\r
+       // READ(CH0,'(A4)') JBNM\r
+       // CLOSE(CH0)\r
+       // L=INDEX(JBNM,' ')-1\r
+       // IF (L.EQ.-1) L=4\r
+       \r
+       // OFLC=JBNM(1:L)//'cn'\r
+       // OFP1=JBNM(1:L)//'p1'\r
+       // OFP0=JBNM(1:L)//'p0'\r
+       \r
+       NARCS=ISNPH[0];\r
+       NQPTS=ISNPH[1];\r
+       TNSUA=ISNPH[2];\r
+       MNSUA=ISNPH[4];\r
+       MNEQN=ISNPH[5];\r
+       \r
+       NJIND=NARCS+1;\r
+       TNGQP=NJIND*NQPTS;\r
+       \r
+       if (TNSUA > MAXSA) {\r
+           IER[0]=52;\r
+           WRTAIL(5,0,IER[0],null);\r
+           return;\r
+       }\r
+       \r
+       //**** COPY POINTERS FROM JAPHYC\r
+       \r
+            \r
+       /*\r
+             OPEN(CH0,FILE=OFP0)\r
+             OPEN(CH1,FILE=OFP1)\r
+             WRITE(*,5) 'EVALUATION OF BCF STARTED:'\r
+       5     FORMAT(A45)\r
+             LODP=QUPTS+NARCS*NQPTS\r
+             LODW=QUWTS+NARCS*NQPTS\r
+             CALL DIAGN4(CCAPH,COCAP,COPHC,CPHCA,EXCAP,EXPHC,ICRVS,IER,\r
+            +IPRVS,NASYM,NCRVS,NINFD,NPRVS,NZERD,ARCLN,ASYMP,BCDMN,CORXX,\r
+            +TOTLN,RGEOM(VTARG),MAP11,ISNPH(DGPOL),ISNPH(JATYP),ISNPH(LOSUB),\r
+            +NARCS,NQPTS,NXINT,CH0,CH1,IGEOM(PARNT),TNSUA,RSNPH(AICOF),\r
+            +RSNPH(ACOEF),RSNPH(BICOF),RSNPH(BCFSN),RSNPH(BCOEF),RSNPH(H0VAL),\r
+            +RSNPH(HIVAL),RGEOM(HALEN),RSNPH(JACIN),RGEOM(MIDPT),RSNPH(SOLUN),\r
+            +RSNPH(LODP),RSNPH(LODW))\r
+             WRITE(*,5) 'EVALUATION OF BCF DONE:'\r
+       C\r
+             IF (IER .GT. 0) THEN\r
+               GOTO 999\r
+             ENDIF\r
+       C\r
+             IF (NASYM .GT. 0) THEN\r
+               WRITE(CH1,*) DASH\r
+               DO 10 I=1,NASYM\r
+                 WRITE(CH1,*) NEWD\r
+                 WRITE(CH1,20) ASYMP(I),0E+0\r
+                 WRITE(CH1,20) ASYMP(I),BIG\r
+       10      CONTINUE\r
+       20      FORMAT(2E16.8)\r
+             ENDIF\r
+             CLOSE(CH1)\r
+       C\r
+             WRITE(CH0,*) DASH\r
+             DO 30 I=2,NARCS\r
+               WRITE(CH0,*) NEWD\r
+               WRITE(CH0,20) CORXX(I),0E+0\r
+               WRITE(CH0,20) CORXX(I),1E+0\r
+       30    CONTINUE\r
+             CLOSE(CH0)\r
+             WRITE(*,5) 'DATA FOR PLOTS DONE:'\r
+       C\r
+             OFLOW=R1MACH(2)\r
+             MCHEP=R1MACH(4)\r
+             UPHYC=MCHEP*CPHCA\r
+             UCANP=MCHEP*CCAPH\r
+             WRITE(*,*)\r
+             WRITE(*,35) 'PHYSICAL ROUNDOFF MAGNIFIES TO:',UPHYC\r
+             WRITE(*,35) 'CANONICAL ROUNDOFF MAGNIFIES TO:',UCANP\r
+       35    FORMAT(A45,2X,E9.2)\r
+       C\r
+             OPEN(CH0,FILE=OFLC)\r
+       C\r
+       C**** WRITE CONFPACK HEADING ON LISTING FILE\r
+       C\r
+             CALL WRHEAD(5,CH0)\r
+       C\r
+             WRITE(CH0,*)\r
+             WRITE(CH0,40)\r
+       40    FORMAT(T4,'MAP',T18,'ESTIMATED EVALUATION',T42,'ESTIMATED MAXIMUM'\r
+            +,/,T18,'CONDITION NUMBER',T42,'ROUNDOFF ERROR *',/)\r
+       C\r
+             IF (NINFD .GT. 0) THEN\r
+               CHPC='**'\r
+             ELSE\r
+               CHPC='  '\r
+             ENDIF\r
+             IF (NZERD .GT. 0) THEN\r
+               CHCP='**'\r
+             ELSE\r
+               CHCP='  '\r
+             ENDIF\r
+       C\r
+             WRITE(CH0,50) CPHCA,CHPC,UPHYC\r
+             WRITE(CH0,60) CCAPH,CHCP,UCANP\r
+       50    FORMAT('PHY --> CAN',T20,E11.3,A2,T44,E11.3,/)\r
+       60    FORMAT('CAN --> PHY',T20,E11.3,A2,T44,E11.3,/)\r
+             WRITE(CH0,*) '* BASED ON UNIT ROUNDOFF IN DOMAIN OF MAP'\r
+             IF (NINFD.GT.0 .OR. NZERD.GT.0) THEN\r
+               WRITE(CH0,*)'** CONDITION NUMBER DEPENDS ON UNIT ROUNDOFF,U:'\r
+               IF (NINFD .GT. 0) THEN\r
+                 WRITE(CH0,70) COPHC,EXPHC\r
+               ENDIF\r
+               IF (NZERD .GT. 0) THEN\r
+                 WRITE(CH0,80) COCAP,EXCAP\r
+               ENDIF\r
+             ENDIF\r
+       70    FORMAT(T4,'PHY --> CAN : CONDTN NO = ',E11.3,'*U**',E11.3)\r
+       80    FORMAT(T4,'CAN --> PHY : CONDTN NO = ',E11.3,'*U**',E11.3)\r
+       C\r
+             PI=4E+0*ATAN(1E+0)\r
+             WRITE(CH0,90) 'END PT.','PARENT','ARGUMENT/PI'\r
+       90    FORMAT(//,A,T10,A,T18,A)\r
+             DO 100 I=1,TNSUA\r
+               WRITE(CH0,110) I,IGEOM(PARNT+I-1),RGEOM(VTARG+I-1)/PI\r
+       100   CONTINUE\r
+       110   FORMAT(I3,T10,I3,T18,E16.8)\r
+       C\r
+             WRITE(CH0,120) 'SUBARC','% PHYSICAL','% CIRCLE'\r
+       120   FORMAT(/,A,T10,A,T29,A)\r
+             DO 130 I=1,TNSUA\r
+               ANGSP=RGEOM(VTARG+I)-RGEOM(VTARG+I-1)\r
+               WRITE(CH0,140) I,ARCLN(I)/TOTLN,ANGSP/2E+0/PI\r
+       130   CONTINUE\r
+       140   FORMAT(I4,T10,E14.7,T29,E14.7)\r
+       C\r
+             WRITE(CH0,150) 'SUB','ACHIEVED','CROWDING','ARC','RESOLUTION',\r
+            +'FACTOR'\r
+       150   FORMAT(/,A,T7,A,T19,A,/,A,T7,A,T19,A)\r
+             RESMN=OFLOW\r
+             DO 160 I=1,TNSUA\r
+               ANGSP=RGEOM(VTARG+I)-RGEOM(VTARG+I-1)\r
+               IF (ANGSP.EQ.0E+0) THEN\r
+                 CR=OFLOW\r
+                 LA=0E+0\r
+               ELSE\r
+                 CR=2E+0*PI*ARCLN(I)/ABS(ANGSP)/TOTLN\r
+                 IF (RSNPH(ERARC+I-1).EQ.0E+0) THEN\r
+                   LA=OFLOW\r
+                 ELSE\r
+                   LA=ABS(ANGSP)/(2E+0*RSNPH(ERARC+I-1))\r
+                 ENDIF\r
+               ENDIF\r
+               IF (LA .LT. RESMN) THEN\r
+                 RESMN=LA\r
+                 IMNLA=I\r
+               ENDIF\r
+               WRITE(CH0,170) I,LA,CR\r
+       160   CONTINUE\r
+       170   FORMAT(I2,T4,2E12.3)\r
+       C\r
+             WRITE(CH0,180) RESMN,IMNLA\r
+       180   FORMAT(/,'MINIMUM SUBARC RESOLUTION IS ',E11.3,' ON SUBARC ',I2) \r
+             WRITE(*,*)\r
+             WRITE(*,35) 'MINIMUM SUBARC RESOLUTION:',RESMN\r
+       C\r
+             WRITE(CH0,*)\r
+             IF (.NOT.MAP11 .OR. RESMN.LT.CRRES) THEN\r
+       C\r
+       C****   MESSAGE TO STANDARD OUTPUT\r
+       C\r
+               WRITE(*,185) '*** W A R N I N G  ***'\r
+       185     FORMAT(//,T20,A)\r
+               IF (RESMN.LT.CRRES) THEN\r
+                 WRITE(*,5) 'THE ABOVE RESOLUTION IS TOO SMALL:'\r
+               ENDIF\r
+               IF (.NOT.MAP11) THEN\r
+                 WRITE(*,5) 'BCF DERIVATIVE CHANGES SIGN:'\r
+               ENDIF\r
+       C\r
+               WRITE(CH0,*) '        ***  W A R N I N G ***'\r
+               WRITE(CH0,*)\r
+               I=0\r
+       C\r
+               IF (RESMN.LT.1) THEN\r
+                 I=I+1\r
+                 WRITE(CH0,190) I,'.  THE ABOVE SUBARC RESOLUTION MEANS THAT IT\r
+            + WILL BE PRACTICALLY'\r
+                 WRITE(CH0,*) '    IMPOSSIBLE FOR THE INVERSE MAP TO DISCRIMINA\r
+            +TE CORRECTLY'\r
+                 WRITE(CH0,200) '    BETWEEN NEIGHBOURING POINTS NEAR SUB ARC '\r
+            +,IMNLA \r
+               ELSE IF (RESMN.LT.CRRES) THEN\r
+                 I=I+1\r
+                 WRITE(CH0,190) I,'. THE ABOVE SUBARC RESOLUTION MEANS THAT THE\r
+            + INVERSE MAP MAY NOT'\r
+                 WRITE(CH0,*) '    BE ABLE TO RELIABLY DISCRIMINATE CORRECTLY B\r
+            +ETWEEN'\r
+                 WRITE(CH0,200) '    NEIGHBOURING POINTS NEAR ARC ',IMNLA \r
+               ENDIF\r
+       190   FORMAT(/,I1,A)\r
+       200   FORMAT(A,I3)\r
+       C\r
+               IF (NCRVS .GT. 0) THEN\r
+                 I=I+1\r
+                 WRITE(CH0,190) I,'.  THERE IS A COMPLETE REVERSAL OF DIRECTION\r
+            + ON THE FOLLOWING SUB ARCS:'\r
+                 WRITE(CH0,'(T10,I3)') (ICRVS(J),J=1,NCRVS)\r
+               ENDIF\r
+       C\r
+               IF (NPRVS .GT. 0) THEN\r
+                 I=I+1\r
+                 WRITE(CH0,190) I,'.  THERE IS A REVERSAL OF DIRECTION WITHIN T\r
+            +HE FOLLOWING SUB ARCS:'\r
+                 WRITE(CH0,'(T10,I3)') (IPRVS(J),J=1,NPRVS)\r
+                 WRITE(CH0,*) '    THE CORRESPONDING MINIMUM VALUES OF THE BOUN\r
+            +DARY CORRESPONDENCE'\r
+                 WRITE(CH0,*) '    DERIVATIVE ARE:'\r
+                 WRITE(CH0,'(T10,E9.2)') (BCDMN(J),J=1,NPRVS)\r
+               ENDIF\r
+             ENDIF\r
+             CLOSE(CH0)\r
+       999   CONTINUE\r
+       C\r
+       C**** WRITE CLOSING MESSAGE TO STANDARD OUTPUT CHANNEL AND LISTING FILE\r
+       C\r
+             CALL WRTAIL(5,0,IER)\r
+             CALL WRTAIL(5,CH0,IER)\r
+       C */\r
+    } // private void CNDPLT\r
 \r
 \r
       /**\r