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

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

index 99ef653..13560ca 100644 (file)
@@ -8268,7 +8268,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        int MNSUA,NEQNS,\r
             TNSUA;\r
        final double DELTA = 0.2;\r
-       double LGTOL,SUPER,THET0;\r
+       double LGTOL,THET0;\r
        double CT[] = new double[2];\r
        double theta;\r
        double mag;\r
@@ -8288,7 +8288,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        MNEQN=ISNPH[5];\r
        NJIND=NARCS+1;\r
        TNGQP=NJIND*NQPTS;\r
-       SUPER=RGEOM[0];\r
+       //SUPER=RGEOM[0];\r
        LGTOL=RGEOM[1];\r
        \r
        IQUPH[1]=TNSUA;\r
@@ -9623,7 +9623,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
             \r
        //     LOCAL VARAIBLES\r
        \r
-       int MNSUA,TNSUA;\r
+       int TNSUA;\r
        // CHARACTER*6 IERTXT\r
        \r
        // EXTERNAL PHTCA1,IERTXT\r
@@ -9631,7 +9631,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        NARCS=ISNPH[0];\r
        NQPTS=ISNPH[1];\r
        TNSUA=ISNPH[2];\r
-       MNSUA=ISNPH[4];\r
+       //MNSUA=ISNPH[4];\r
        MNEQN=ISNPH[5];\r
        MQUPH=IQUPH[3];\r
        \r
@@ -10632,10 +10632,10 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                // OTHER SCALAR VARIABLES\r
                \r
-               int I,MNCOF,MNEQN,MNSUA,MNSUC,NEQNS,NJCOG,PT,TNSUA,TNSUC;\r
+               int I,MNCOF,MNSUC,NEQNS,NJCOG,PT,TNSUA,TNSUC;\r
                double COLSC[] = new double[]{-1.0};\r
                double INNRAD = 0.0;\r
-               double LGTOL,PI,SUPER,THET0,NEWTL;\r
+               double LGTOL,PI,THET0,NEWTL;\r
                double DCAP0[] = new double[2];\r
                double FACTR[] = new double[2];\r
                double mult;\r
@@ -10653,11 +10653,11 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                NARCS=IGEOM[0];\r
                NQPTS=IGEOM[1];\r
-               MNSUA=IGEOM[3];\r
+               //MNSUA=IGEOM[3];\r
                TNSUA=ISNPH[2];\r
                NEQNS=ISNPH[3];\r
                MNEQN=ISNPH[5];\r
-               SUPER=RGEOM[0];\r
+               //SUPER=RGEOM[0];\r
                LGTOL=RGEOM[1];\r
                MNCOF=IBNDS[1];\r
                MNSUC=IBNDS[0];\r
@@ -11325,7 +11325,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           final int MNITS = 10;\r
           final int NBSCT = 3;\r
           int NITS,STEPS;\r
-          double DFT,EPS,FL,FT,FU,JACSUM,RDIFF,TUPI;\r
+          double DFT,EPS,FL,FT,FU,RDIFF,TUPI;\r
           double CO[];\r
           // EXTERNAL JACSUM\r
           int i;\r
@@ -12048,7 +12048,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
           // LOCAL VARIABLES\r
                \r
-          int AJT,DG,I,IC,L,LOM,LOD,PSA,PT;\r
+          int AJT,DG,I,IC,LOM,LOD,PSA,PT;\r
           double MOD;\r
           double COF[] = new double[2];\r
           // COMPLEX COF\r
@@ -12192,7 +12192,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
           final int NINTS = 5;\r
           int TNPQP[] = new int[1];\r
-          int MNCOF,MNSUC,TNSUC;\r
+          int MNSUC,TNSUC;\r
           final double DELTA = 0.2;\r
           double TOLOU[] = new double[1];\r
           double LGTOL;\r
@@ -12207,7 +12207,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           NQPTS=ISNCA[1];\r
           TNSUC=ISNCA[2];\r
           MNSUC=ISNCA[4];\r
-          MNCOF=ISNCA[5];\r
+          //MNCOF=ISNCA[5];\r
           LGTOL=RSNCA[0];\r
           ZQUCA[0]=ZSNCA[0];\r
           ZQUCA[1]=ZSNCA[1];\r
@@ -12277,7 +12277,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                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
+       //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
@@ -12678,8 +12678,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        int NPRVS[] = new int[1];\r
        int NZERD[] = new int[1];\r
        int IMNLA = 0;\r
-       int I,J,L,MNSUA,TNSUA;\r
-       final double BIG = 4.4;\r
+       int I,J,TNSUA;\r
+       //final double BIG = 4.4;\r
        double CCAPH[] = new double[1];\r
        double COCAP[] = new double[1];\r
        double COPHC[] = new double[1];\r
@@ -12688,7 +12688,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        double EXPHC[] = new double[1];\r
        double TOTLN[] = new double[1];\r
        double ANGSP,CR,LA,\r
-            OFLOW,PI,R1MACH,MCHEP;\r
+            OFLOW,PI,MCHEP;\r
        String CHPC, CHCP;\r
        //CHARACTER OFLC*6,OFP0*6,OFP1*6,JBNM*4,CHPC*2,CHCP*2\r
        \r
@@ -12727,7 +12727,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        NARCS=ISNPH[0];\r
        NQPTS=ISNPH[1];\r
        TNSUA=ISNPH[2];\r
-       MNSUA=ISNPH[4];\r
+       //MNSUA=ISNPH[4];\r
        MNEQN=ISNPH[5];\r
        \r
        NJIND=NARCS+1;\r
@@ -12971,8 +12971,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        final int MXCOF = 32;\r
        final int QP = 4;\r
        int AJT,DG,I,I1,IA,JT,K,LOD,LOM,NINTS;\r
-       double AL,ATOL,BT,CC,COF,D,DSDT,H0,HH,JACSUM,MCHEP,MPT,\r
-            PHI,R1MACH,RTOL,SEND,SINC,SJT,SS,SUM,TERM,TINC,TT,TUPI,X,XX,YMAX,\r
+       double AL,BT,CC,COF,D,H0,HH,MCHEP,MPT,\r
+            PHI,SEND,SINC,SJT,SS,SUM,TERM,TINC,TT,TUPI,X,XX,YMAX,\r
             YMIN,YY;\r
        double T1[] = new double[2];\r
        double T2[] = new double[2];\r
@@ -12989,8 +12989,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
        TUPI=2.0*Math.PI;\r
        MCHEP=EPS;\r
-       RTOL=10.0*MCHEP;\r
-       ATOL=100.0*MCHEP;\r
+       //RTOL=10.0*MCHEP;\r
+       //ATOL=100.0*MCHEP;\r
        NCRVS[0]=0;\r
        NPRVS[0]=0;\r
        CCAPH[0]=0.0;\r
@@ -13003,7 +13003,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        //     TOTAL LENGTH (TOTLN) OF THE BOUNDARY\r
        \r
        TOTLN[0]=0.0;\r
-       /*for (IA=1; IA <= TNSUA; IA++) {\r
+       for (IA=1; IA <= TNSUA; IA++) {\r
            PT=PARNT[IA-1];\r
            MD=MIDPT[IA-1];\r
            HL=HALEN[IA-1];\r
@@ -13145,254 +13145,254 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        //       ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS FOR INFINITE\r
        //       DERIVATIVE CASES.\r
        \r
-               if (BT < 0.0) {\r
-                       A = new double[DG-1];\r
-                       B = new double[DG-1];\r
-                       for (N = 0; N < DG-1; N++) {\r
-                               A[N] = A1COF[LOD+N-1];\r
-                               B[N] = B1COF[LOD+N-1];\r
-                       }\r
-                       CO = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                               CO[N] = BCFSN[LOM+N];\r
-                       }\r
-                   PHI=JACSUM(-1E+0,DG-1,A,B,H1VAL[AJT-1],CO);\r
-                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
-                   COF=Math.abs(PHI)/Math.pow(D,(BT+1.0));\r
-                   TERM=Math.pow(MCHEP,BT)*COF;\r
-                   if (TERM > CPHCA[0]) {\r
-                       CPHCA[0]=TERM;\r
-                       COPHC[0]=COF;\r
-                       EXPHC[0]=BT;\r
-                   }\r
-               } // if (BT < 0.0)\r
-               if (BT > 0.0) {\r
-                       A = new double[DG-1];\r
-                       B = new double[DG-1];\r
-                       for (N = 0; N < DG-1; N++) {\r
-                               A[N] = A1COF[LOD+N-1];\r
-                               B[N] = B1COF[LOD+N-1];\r
-                       }\r
-                       CO = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                               CO[N] = BCFSN[LOM+N];\r
-                       }\r
-                   PHI=JACSUM(-1.0,DG-1,A,B,H1VAL[AJT-1],CO);\r
-                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
-                   if (Math.abs(PHI) == 0.0) {\r
-                       CCAPH[0]=YMAX;\r
-                       COCAP[0]=YMAX;\r
+               if (BT < 0.0) {\r
+                A = new double[DG-1];\r
+                B = new double[DG-1];\r
+               for (N = 0; N < DG-1; N++) {\r
+                       A[N] = A1COF[LOD+N-1];\r
+                       B[N] = B1COF[LOD+N-1];\r
+               }\r
+               CO = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                       CO[N] = BCFSN[LOM+N];\r
+               }\r
+                   PHI=JACSUM(-1E+0,DG-1,A,B,H1VAL[AJT-1],CO);\r
+                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
+                   COF=Math.abs(PHI)/Math.pow(D,(BT+1.0));\r
+                   TERM=Math.pow(MCHEP,BT)*COF;\r
+                   if (TERM > CPHCA[0]) {\r
+                       CPHCA[0]=TERM;\r
+                       COPHC[0]=COF;\r
+                       EXPHC[0]=BT;\r
+                   }\r
+               } // if (BT < 0.0)\r
+               if (BT > 0.0) {\r
+                       A = new double[DG-1];\r
+               B = new double[DG-1];\r
+               for (N = 0; N < DG-1; N++) {\r
+                       A[N] = A1COF[LOD+N-1];\r
+                       B[N] = B1COF[LOD+N-1];\r
+               }\r
+               CO = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                       CO[N] = BCFSN[LOM+N];\r
+               }\r
+                   PHI=JACSUM(-1.0,DG-1,A,B,H1VAL[AJT-1],CO);\r
+                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
+                   if (Math.abs(PHI) == 0.0) {\r
+                       CCAPH[0]=YMAX;\r
+                       COCAP[0]=YMAX;\r
+                       EXCAP[0]=AL-1.0;\r
+                   }\r
+                   else {\r
+                   COF=D/Math.pow(Math.abs(PHI),AL);\r
+                   TERM=Math.pow(MCHEP,(AL-1.0))*COF;\r
+                   if (TERM > CCAPH[0]) {\r
+                       CCAPH[0]=TERM;\r
+                       COCAP[0]=COF;\r
                        EXCAP[0]=AL-1.0;\r
-                   }\r
-                   else {\r
-                           COF=D/Math.pow(Math.abs(PHI),AL);\r
-                           TERM=Math.pow(MCHEP,(AL-1.0))*COF;\r
-                           if (TERM > CCAPH[0]) {\r
-                               CCAPH[0]=TERM;\r
-                               COCAP[0]=COF;\r
-                               EXCAP[0]=AL-1.0;\r
-                           } // if (TERM > CCAPH[0])\r
-                   } // else\r
-               // if (BT > 0.0)\r
-       \r
-                // "DO 50" LOOP FOR POINTS INTERIOR TO ARC NUMBER IA\r
-       \r
-               for (I=1; I <= NINTS-1; I++) {\r
-                   TT=TT+TINC;\r
-       \r
-                   // ARC LENGTH INCREASE BY GAUSS-LEGENDRE\r
-       \r
-                   SUM=0.0;\r
-                   for (K=1; K <= NQPTS; K++) {\r
-                       X=TT+0.5*TINC*(QUPTS[K-1]-1.0);\r
-                       SUM=SUM+QUWTS[K-1]*DSDT(X);\r
-                   } //  for (K=1; K <= NQPTS; K++)\r
-                   SINC=0.5*TINC*SUM;\r
-                   SS=SS+SINC;\r
-                   XX=SS/TOTLN[0];\r
-       \r
-                   // EVALUATE DIMENSIONLESS BCF DERIVATIVE *YY*\r
-                   A = new double[DG];\r
-                       B = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                               A[N] = ACOEF[LOD+N-1];\r
-                               B[N] = BCOEF[LOD+N-1];\r
-                       }\r
-                   PHI=JACSUM(SJT*TT,DG,A,B,H0,JACOF);\r
-                   D=DSDT(TT);\r
-                   if (D == 0.0) {\r
-                       IER[0]=51;\r
-                       return;\r
-                   } // if (D == 0.0)\r
-                   YY=TOTLN[0]*Math.pow((1.0+SJT*TT),BT)*PHI/D;\r
-                   Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
-                   YMIN=Math.min(YY,YMIN);\r
-       \r
-                   // ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS FOR FINITE\r
-                   // DERIVATIVE CASES.\r
-       \r
-                   if (NINFD[0] == 0.0) {\r
-                       CPHCA[0]=Math.max(CPHCA[0],TUPI*Math.abs(YY)/TOTLN[0]);\r
-                   }\r
-                   if (NZERD[0] == 0.0) {\r
-                       if (YY == 0.0) {\r
-                           CCAPH[0]=YMAX;\r
-                       }\r
-                       else {\r
-                           CCAPH[0]=Math.max(CCAPH[0],TOTLN[0]/TUPI/Math.abs(YY));\r
-                       }\r
-                   } // if (NZERD[0] == 0.0)\r
-       \r
-                   // EVALUATE DIMENSIONLESS BCF *YY*\r
-                   A = new double[DG-1];\r
-                       B = new double[DG-1];\r
-                       for (N = 0; N < DG-1; N++) {\r
-                               A[N] = A1COF[LOD+N-1];\r
-                               B[N] = B1COF[LOD+N-1];\r
-                       }\r
-                       CO = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                               CO[N] = BCFSN[LOM+N];\r
-                       }\r
-                   PHI=JACSUM(SJT*TT,DG-1,A,B,H1VAL[AJT-1],CO);\r
-                   PHI=BCFSN[LOM-1]-(1.0-SJT*TT)*PHI;\r
-                   YY=(CC+SJT*Math.pow((1.0+SJT*TT),(1.0+BT))*PHI)/TUPI;\r
-                   Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
-               } // for (I=1; I <= NINTS-1; I++)\r
-       \r
-               // NEXT TAKE END POINT OF ARC NUMBER IA\r
-       \r
-               TT=1.0;\r
-               D=DSDT(TT);\r
-               SEND=SEND+ARCLN[IA-1];\r
-               SS=SEND;\r
-               XX=SS/TOTLN[0];\r
-       \r
-               // EVALUATE DIMENSIONLESS BCF DERIVATIVE *YY*\r
-       \r
-               if (JT < 0.0) {\r
-                   if (BT < 0.0) {\r
-                       YY=YMAX;\r
-                       NASYM[0]=NASYM[0]+1;\r
-                       ASYMP[NASYM[0]-1]=XX;\r
-                   } // if (BT < 0.)\r
-                   else if (BT > 0.0) {\r
-                       YY=0.0;\r
-                   }\r
-                   else {\r
-                       A = new double[DG];\r
-                       B = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                                A[N] = ACOEF[LOD+N-1];\r
-                                B[N] = BCOEF[LOD+N-1];\r
-                       }\r
-                       PHI=JACSUM(SJT*TT,DG,A,B,H0,JACOF);\r
-                       if (D == 0.0) {\r
-                           IER[0]=51;\r
-                           return;\r
-                       }\r
-                       YY=TOTLN[0]*PHI/D;\r
-                   } // else\r
-               } // if (JT < 0.0)\r
-               else {\r
-                       A = new double[DG];\r
+                   } // if (TERM > CCAPH[0])\r
+                   } // else\r
+              } // if (BT > 0.0) \r
+       \r
+            // "DO 50" LOOP FOR POINTS INTERIOR TO ARC NUMBER IA\r
+       \r
+               for (I=1; I <= NINTS-1; I++) {\r
+                   TT=TT+TINC;\r
+       \r
+                   // ARC LENGTH INCREASE BY GAUSS-LEGENDRE\r
+       \r
+                   SUM=0.0;\r
+                   for (K=1; K <= NQPTS; K++) {\r
+                       X=TT+0.5*TINC*(QUPTS[K-1]-1.0);\r
+                       SUM=SUM+QUWTS[K-1]*DSDT(X);\r
+                   } //  for (K=1; K <= NQPTS; K++)\r
+                   SINC=0.5*TINC*SUM;\r
+                   SS=SS+SINC;\r
+                   XX=SS/TOTLN[0];\r
+       \r
+                   // EVALUATE DIMENSIONLESS BCF DERIVATIVE *YY*\r
+                   A = new double[DG];\r
+               B = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                       A[N] = ACOEF[LOD+N-1];\r
+                       B[N] = BCOEF[LOD+N-1];\r
+               }\r
+                   PHI=JACSUM(SJT*TT,DG,A,B,H0,JACOF);\r
+                   D=DSDT(TT);\r
+                   if (D == 0.0) {\r
+                       IER[0]=51;\r
+                       return;\r
+                   } // if (D == 0.0)\r
+                   YY=TOTLN[0]*Math.pow((1.0+SJT*TT),BT)*PHI/D;\r
+                   Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
+                   YMIN=Math.min(YY,YMIN);\r
+       \r
+                   // ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS FOR FINITE\r
+                   // DERIVATIVE CASES.\r
+       \r
+                   if (NINFD[0] == 0.0) {\r
+                       CPHCA[0]=Math.max(CPHCA[0],TUPI*Math.abs(YY)/TOTLN[0]);\r
+                   }\r
+                   if (NZERD[0] == 0.0) {\r
+                       if (YY == 0.0) {\r
+                           CCAPH[0]=YMAX;\r
+                       }\r
+                       else {\r
+                           CCAPH[0]=Math.max(CCAPH[0],TOTLN[0]/TUPI/Math.abs(YY));\r
+                       }\r
+                   } // if (NZERD[0] == 0.0)\r
+       \r
+                   // EVALUATE DIMENSIONLESS BCF *YY*\r
+                   A = new double[DG-1];\r
+               B = new double[DG-1];\r
+               for (N = 0; N < DG-1; N++) {\r
+                       A[N] = A1COF[LOD+N-1];\r
+                       B[N] = B1COF[LOD+N-1];\r
+               }\r
+               CO = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                       CO[N] = BCFSN[LOM+N];\r
+               }\r
+                   PHI=JACSUM(SJT*TT,DG-1,A,B,H1VAL[AJT-1],CO);\r
+                   PHI=BCFSN[LOM-1]-(1.0-SJT*TT)*PHI;\r
+                   YY=(CC+SJT*Math.pow((1.0+SJT*TT),(1.0+BT))*PHI)/TUPI;\r
+                   Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
+               } // for (I=1; I <= NINTS-1; I++)\r
+       \r
+               // NEXT TAKE END POINT OF ARC NUMBER IA\r
+       \r
+               TT=1.0;\r
+               D=DSDT(TT);\r
+               SEND=SEND+ARCLN[IA-1];\r
+               SS=SEND;\r
+               XX=SS/TOTLN[0];\r
+       \r
+               // EVALUATE DIMENSIONLESS BCF DERIVATIVE *YY*\r
+       \r
+               if (JT < 0.0) {\r
+                   if (BT < 0.0) {\r
+                       YY=YMAX;\r
+                       NASYM[0]=NASYM[0]+1;\r
+                       ASYMP[NASYM[0]-1]=XX;\r
+                   } // if (BT < 0.)\r
+                   else if (BT > 0.0) {\r
+                       YY=0.0;\r
+                   }\r
+                   else {\r
+                       A = new double[DG];\r
                        B = new double[DG];\r
                        for (N = 0; N < DG; N++) {\r
                                 A[N] = ACOEF[LOD+N-1];\r
                                 B[N] = BCOEF[LOD+N-1];\r
                        }\r
-                   PHI=JACSUM(TT,DG,A,B,H0,JACOF);\r
-                   if (D == 0.0) {\r
-                       IER[0]=51;\r
-                       return;\r
-                   }\r
-                   YY=TOTLN[0]*Math.pow(2.0,BT)*PHI/D;\r
-               } // else\r
-               Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
-               YMIN=Math.min(YY,YMIN);\r
-               if (YMIN < 0.0 && (VTARG[IA] >= VTARG[IA-1])) {\r
-                   NPRVS[0]=NPRVS[0]+1;\r
-                   IPRVS[NPRVS[0]-1]=IA;\r
-                   BCDMN[NPRVS[0]-1]=YMIN;\r
-                   MAP11[0]= false;\r
-               } // if (YMIN < 0.0 && (VTARG[IA] >= VTARG[IA-1]))\r
-       \r
-               // ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS\r
-       \r
-               if (NINFD[0] == 0.0) {\r
-                   CPHCA[0]=Math.max(CPHCA[0],TUPI*Math.abs(YY)/TOTLN[0]);\r
-               }\r
-               if (NZERD[0] == 0.0) {\r
-                   if (YY == 0.0) {\r
-                       CCAPH[0]=YMAX;\r
-                   }\r
-                   else {\r
-                       CCAPH[0]=Math.max(CCAPH[0],TOTLN[0]/TUPI/Math.abs(YY));\r
-                   }\r
-               } // if (NZERD[0] == 0.0)\r
-               if (BT < 0.0) {\r
-                       A = new double[DG-1];\r
-                       B = new double[DG-1];\r
-                       for (N = 0; N < DG-1; N++) {\r
-                               A[N] = A1COF[LOD+N-1];\r
-                               B[N] = B1COF[LOD+N-1];\r
-                       }\r
-                       CO = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                               CO[N] = BCFSN[LOM+N];\r
-                       }\r
-                   PHI=JACSUM(-1.0,DG-1,A,B,H1VAL[AJT-1],CO);\r
-                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
-                   COF=Math.abs(PHI)/Math.pow(D,(BT+1.0));\r
-                   TERM=Math.pow(MCHEP,BT)*COF;\r
-                   if (TERM > CPHCA[0]) {\r
-                       CPHCA[0]=TERM;\r
-                       COPHC[0]=COF;\r
-                       EXPHC[0]=BT;\r
-                   } // if (TERM > CPHCA[0])\r
-               } // if (BT < 0.0)\r
-               if (BT > 0.0) {\r
-                       A = new double[DG-1];\r
-                       B = new double[DG-1];\r
-                       for (N = 0; N < DG-1; N++) {\r
-                               A[N] = A1COF[LOD+N-1];\r
-                               B[N] = B1COF[LOD+N-1];\r
-                       }\r
-                       CO = new double[DG];\r
-                       for (N = 0; N < DG; N++) {\r
-                               CO[N] = BCFSN[LOM+N];\r
-                       }\r
-                   PHI=JACSUM(-1.0,DG-1,A,B,H1VAL[AJT-1],CO);\r
-                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
-                   if (Math.abs(PHI) == 0.0) {\r
-                       CCAPH[0]=YMAX;\r
-                       COCAP[0]=YMAX;\r
-                       EXCAP[0]=AL-1.0;\r
-                   }\r
-                   else {\r
-                       COF=D/Math.pow(Math.abs(PHI),AL);\r
-                       TERM=Math.pow(MCHEP,(AL-1.0))*COF;\r
-                       if (TERM > CCAPH[0]) {\r
-                           CCAPH[0]=TERM;\r
-                           COCAP[0]=COF;\r
-                           EXCAP[0]=AL-1.0;\r
-                       } // if (TERM > CCAPH[0])\r
-                   } // else\r
-               } // if (BT > 0.0)\r
-       \r
-               // EVALUATE DIMENSIONLESS BCF *YY*\r
-       \r
-               YY=(VTARG[IA]-VTARG[0])/TUPI;\r
-               Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
-               if (JT < 0) {\r
-                   CORXX[PT]=XX;\r
-               }\r
+                       PHI=JACSUM(SJT*TT,DG,A,B,H0,JACOF);\r
+                       if (D == 0.0) {\r
+                           IER[0]=51;\r
+                           return;\r
+                       }\r
+                       YY=TOTLN[0]*PHI/D;\r
+                   } // else\r
+               } // if (JT < 0.0)\r
+               else {\r
+                       A = new double[DG];\r
+               B = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                        A[N] = ACOEF[LOD+N-1];\r
+                        B[N] = BCOEF[LOD+N-1];\r
+               }\r
+                   PHI=JACSUM(TT,DG,A,B,H0,JACOF);\r
+                   if (D == 0.0) {\r
+                       IER[0]=51;\r
+                       return;\r
+                   }\r
+                   YY=TOTLN[0]*Math.pow(2.0,BT)*PHI/D;\r
+               } // else\r
+               Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
+               YMIN=Math.min(YY,YMIN);\r
+               if (YMIN < 0.0 && (VTARG[IA] >= VTARG[IA-1])) {\r
+                   NPRVS[0]=NPRVS[0]+1;\r
+                   IPRVS[NPRVS[0]-1]=IA;\r
+                   BCDMN[NPRVS[0]-1]=YMIN;\r
+                   MAP11[0]= false;\r
+               } // if (YMIN < 0.0 && (VTARG[IA] >= VTARG[IA-1]))\r
+       \r
+               // ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS\r
+       \r
+               if (NINFD[0] == 0.0) {\r
+                   CPHCA[0]=Math.max(CPHCA[0],TUPI*Math.abs(YY)/TOTLN[0]);\r
+               }\r
+               if (NZERD[0] == 0.0) {\r
+                   if (YY == 0.0) {\r
+                       CCAPH[0]=YMAX;\r
+                   }\r
+                   else {\r
+                       CCAPH[0]=Math.max(CCAPH[0],TOTLN[0]/TUPI/Math.abs(YY));\r
+                   }\r
+               } // if (NZERD[0] == 0.0)\r
+               if (BT < 0.0) {\r
+                       A = new double[DG-1];\r
+               B = new double[DG-1];\r
+               for (N = 0; N < DG-1; N++) {\r
+                       A[N] = A1COF[LOD+N-1];\r
+                       B[N] = B1COF[LOD+N-1];\r
+               }\r
+               CO = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                       CO[N] = BCFSN[LOM+N];\r
+               }\r
+                   PHI=JACSUM(-1.0,DG-1,A,B,H1VAL[AJT-1],CO);\r
+                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
+                   COF=Math.abs(PHI)/Math.pow(D,(BT+1.0));\r
+                   TERM=Math.pow(MCHEP,BT)*COF;\r
+                   if (TERM > CPHCA[0]) {\r
+                       CPHCA[0]=TERM;\r
+                       COPHC[0]=COF;\r
+                       EXPHC[0]=BT;\r
+                   } // if (TERM > CPHCA[0])\r
+               } // if (BT < 0.0)\r
+               if (BT > 0.0) {\r
+                       A = new double[DG-1];\r
+               B = new double[DG-1];\r
+               for (N = 0; N < DG-1; N++) {\r
+                       A[N] = A1COF[LOD+N-1];\r
+                       B[N] = B1COF[LOD+N-1];\r
+               }\r
+               CO = new double[DG];\r
+               for (N = 0; N < DG; N++) {\r
+                       CO[N] = BCFSN[LOM+N];\r
+               }\r
+                   PHI=JACSUM(-1.0,DG-1,A,B,H1VAL[AJT-1],CO);\r
+                   PHI=BCFSN[LOM-1]-2.0*PHI;\r
+                   if (Math.abs(PHI) == 0.0) {\r
+                       CCAPH[0]=YMAX;\r
+                       COCAP[0]=YMAX;\r
+                       EXCAP[0]=AL-1.0;\r
+                   }\r
+                   else {\r
+                       COF=D/Math.pow(Math.abs(PHI),AL);\r
+                       TERM=Math.pow(MCHEP,(AL-1.0))*COF;\r
+                       if (TERM > CCAPH[0]) {\r
+                           CCAPH[0]=TERM;\r
+                           COCAP[0]=COF;\r
+                           EXCAP[0]=AL-1.0;\r
+                       } // if (TERM > CCAPH[0])\r
+                   } // else\r
+               } // if (BT > 0.0)\r
+       \r
+               // EVALUATE DIMENSIONLESS BCF *YY*\r
+       \r
+               YY=(VTARG[IA]-VTARG[0])/TUPI;\r
+               Preferences.debug("XX == " + XX + " YY = " + YY + "\n", Preferences.DEBUG_ALGORITHM);\r
+               if (JT < 0) {\r
+                   CORXX[PT]=XX;\r
+               } \r
        \r
        } // for (IA=1; IA <= TNSUA; IA++)\r
        \r
        //     NORMAL EXIT\r
        \r
-       IER[0]=0;*/\r
+       IER[0]=0;\r
         \r
     } // private void DIAGN4\r
     \r
@@ -13416,6 +13416,499 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         return result;\r
 \r
     } // private double DSDT\r
+    \r
+    private void DMCANP(int NPTS, double PHYPT[][], double CANPT[][],\r
+        boolean WANTM, int IER[]) {\r
+       \r
+        //INTEGER NPTS,IER\r
+       //INTEGER IGEOM(*),ISNCA(*),IQUCA(*)\r
+       //REAL RGEOM(*),RSNCA(*)\r
+       //COMPLEX CENTR\r
+       //COMPLEX PHYPT(*),CANPT(*),ZSNCA(*),ZQUCA(*)\r
+       //LOGICAL INTER,WANTM\r
+       \r
+       // ......................................................................\r
+       \r
+       // 1.     DMCANP\r
+       //           DOMAIN MAPPING FOR THE CANONICAL --> PHYSICAL MAP.\r
+       \r
+       // 2.     PURPOSE\r
+       //           GIVEN A VECTOR OF ARBITRARY POINTS IN THE CANONICAL DOMAIN, \r
+       //           THIS ROUTINE COMPUTES THE VECTOR OF APPROXIMATE IMAGE POINTS\r
+       //           IN THE PHYSICAL DOMAIN.\r
+       \r
+       // 3.     CALLING SEQUENCE\r
+       //           CALL DMCANP(NPTS,PHYPT,CANPT,INTER,CENTR,IGEOM,RGEOM,ISNCA,\r
+       //                       RSNCA,ZSNCA,IQUCA,ZQUCA,WANTM,IER)\r
+       \r
+       //        PARAMETERS\r
+       //         ON ENTRY\r
+       //            NPTS   - INTEGER\r
+       //                     THE NUMBER OF POINTS TO BE MAPPED.\r
+       \r
+       //            CANPT  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST NPTS.  THIS IS\r
+       //                     THE VECTOR OF GIVEN POINTS IN THE CANONICAL \r
+       //                     DOMAIN.\r
+       \r
+       //            INTER  - LOGICAL\r
+       //                     TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE \r
+       //                     OTHERWISE. (AS PREVIOUSLY USED IN JAPHYC, GQPHYC)\r
+       \r
+       //            CENTR  - COMPLEX\r
+       //                     THE POINT IN THE PHYSICAL PLANE THAT IS TO BE\r
+       //                     MAPPED TO THE CENTRE OF THE UNIT DISC.  FOR\r
+       //                     EXTERIOR DOMAINS CENTR MUST BE SOME POINT IN THE\r
+       //                     COMPLEMENTARY INTERIOR  PHYSICAL DOMAIN. (AS PREV-\r
+       //                     IOUSLY USED IN JAPHYC, GQPHYC)\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
+       \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
+       //            IQUCA  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR PREVIOUSLY SET UP BY GQCANP.\r
+       \r
+       //            RQUCA  - REAL ARRAY\r
+       //                     THE REAL VECTOR PREVIOUSLY SET UP BY GQCANP.\r
+       \r
+       //            ZQUCA  - COMPLEX ARRAY\r
+       //                     THE COMPLEX VECTOR PREVIOUSLY SET UP BY GQCANP.\r
+       \r
+       //            WANTM  - LOGICAL\r
+       //                     IF WANTM IS TRUE THEN, ON AN ABNORMAL EXIT, AN\r
+       //                     ERROR MESSAGE IS WRITTEN ON THE STANDARD OUTPUT\r
+        //                     CHANNEL.  IF WANTM IS FALSE THEN NO MESSAGE IS\r
+       //                     WRITTEN.\r
+       //         ON EXIT\r
+       //            PHYPT  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST NPTS.  PHYPT(K)\r
+       //                     IS THE COMPUTED IMAGE IN THE PHYSICAL DOMAIN OF \r
+       //                     THE GIVEN CANONICAL POINT CANPT(K), K=1,...,NPTS.\r
+       \r
+       //            IER    - INTEGER\r
+       //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+       //                     A MESSAGE TO DESCRIBE THE ERROR IS WRITTEN ON \r
+       //                     THE STANDARD OUTPUT CHANNEL IF WANTM IS TRUE.\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
+       //             ROUTINES JACANP AND GQCANP HAVE SUCCESSFULLY EXECUTED, \r
+       //             AND THAT MANY INPUT ARGUMENTS FOR DMCANP ARE OUTPUT VALUES\r
+       //             FROM JACANP AND GQCANP.\r
+       //           - THIS ROUTINE MAY BE USED FOR MAPPING POINTS ON THE UNIT\r
+       //             CIRCLE, BUT THE ROUTINE BMCANP WILL BE SOMEWHAT MORE \r
+       //             EFFICIENT FOR THIS CASE.\r
+       \r
+       // ......................................................................\r
+       //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+       //     LAST UPDATE: 3 JULY 1990\r
+       // ......................................................................\r
+            \r
+       //     LOCAL VARAIBLES\r
+       \r
+       int MNCOF,MNSUA,MNSUC,TNSUC;\r
+       //CHARACTER*6 IERTXT\r
+       \r
+       //EXTERNAL CATPH4,IERTXT\r
+       \r
+       NARCS=ISNCA[0];\r
+       NQPTS=ISNCA[1];\r
+       TNSUC=ISNCA[2];\r
+       MNSUC=ISNCA[4];\r
+       MNCOF=ISNCA[5];\r
+        MQUCA=IQUCA[3];\r
+       MNSUA=IGEOM[3];\r
+       \r
+       NJIND=NARCS+1;\r
+       TNGQP=NJIND*NQPTS;\r
+       \r
+       // SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC\r
+       \r
+       // SET UP POINTERS TO ISNCA, RSNCA AND ZSNCA, AS IN JACANP\r
+       \r
+       // SET UP POINTERS TO IQUCA AND ZQUCA, AS IN GQCANP\r
+       \r
+       // GET THE REQUIRED PHYSICAL POINTS\r
+       \r
+       CATPH4(NPTS,PHYPT,CANPT,NARCS,NQPTS,TNSUC,DGPOC,\r
+            JTYPC,LSUBC,LQSBG,NPPQG,PARNT,\r
+            PRNSA,AICOC,ACOFC,BICOC,BCOFC,\r
+            H0VLC,HIVLC,HALEN,JAINC,RGEOM[1],\r
+            MIDPT,PHPAS,QUPTC,QUWTC,VTARG[0],\r
+            VARGC,BFSNC,CENTR,ZSNCA,SOLNC,\r
+            WPPQG,ZPPQG,INTER,IER);\r
+       \r
+       // SEND ERROR MESSAGE TO STANDARD OUTPUT OF NECESSARY\r
+       \r
+       if (IER[0] > 0 && WANTM) System.out.println(IERTXT(IER[0]));\r
+       \r
+    } // private void DMCANP\r
+\r
+    private void CATPH4(int NPTS, double PHYPT[][], double CANPT[][], int NARCS, int NQPTS,\r
+        int TNSUC, int DGPOC[], int JTYPC[], int LSUBC[], int LQSBG[], int NPPQG[],\r
+        int PARNT[], int PRNSA[], double A1COC[], double ACOFC[], double B1COC[], double BCOFC[],\r
+        double H0VLC[], double H1VLC[], double HALEN[], double JAINC[], double LGTOL,\r
+        double MIDPT[], double PHPAS[], double QUPTC[], double QUWTC[], double THET0,\r
+       double VARGC[], double BFSNC[][], double CENTR[], double FACTR[], double SOLNC[][],\r
+       double WPPQG[][], double ZPPQG[][], boolean INTER, int IER[]) {\r
+       //INTEGER IER,NPTS,NARCS,NQPTS,TNSUC\r
+       //INTEGER DGPOC(*),JTYPC(*),LSUBC(*),LQSBG(*),NPPQG(*),PARNT(*),\r
+       //+PRNSA(*)\r
+       //REAL LGTOL,THET0\r
+       //REAL A1COC(*),ACOFC(*),B1COC(*),BCOFC(*),H0VLC(*),H1VLC(*),\r
+       //+HALEN(*),JAINC(*),MIDPT(*),PHPAS(*),QUPTC(*),QUWTC(*),VARGC(*)\r
+       //COMPLEX CENTR,FACTR\r
+       //COMPLEX BFSNC(*),CANPT(*),PHYPT(*),SOLNC(*),WPPQG(*),ZPPQG(*)\r
+       //LOGICAL INTER\r
+       \r
+       //     GIVEN THE ARRAY CANPT OF NPTS POINTS IN THE CANONICAL PLANE, THIS\r
+       //     ROUTINE COMPUTES THE ARRAY PHYPT OF IMAGES IN THE PHYSICAL\r
+       //     PLANE.\r
+       \r
+       //     IER=0  - NORMAL EXIT\r
+       //     IER=47 - LOCAL PARAMETER MXNQD NEEDS INCREASING\r
+       //     IER=48 - LOCAL PARAMETER MNCOF NEEDS INCREASING\r
+       //     IER=49 - LOCAL PARAMETER MQIN1 NEEDS INCREASING\r
+       \r
+       //.......................................................................\r
+       //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+       //     LAST UPDATE: 7 JULY 90\r
+       //.......................................................................\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MNCOF = 32;\r
+       final int MQIN1 = 21;\r
+       final int MXNQD = 144;\r
+       final int NBSCT = 3;\r
+       int AJTC,DGC,I,IA,IP,K,J,J1,JQ,JTC,LODC,LOL,\r
+            LOMC,NQ,NQUAD,PSA,QINTS;\r
+       final double PTHTL = 1.0E-3;\r
+       final double DELTA = 2.0E-1;\r
+       double AA,ARG,ARGW,AWW,BB,BETAC,DIST,EFPTL,EPS_CATPH4,HA,ILW,IMXI,MD,\r
+            MEAN,PI,REXI,RLW,RR,RRB,S2C,SCO,SJTC,SUM1,TOLOU,TT,\r
+            TUPI,UPPER;\r
+       double CT0[] = new double[2];\r
+       double CT2[] = new double[2];\r
+       double CTA[] = new double[2];\r
+       double CTB[] = new double[2];\r
+       double FP[] = new double[2];\r
+       double CSUM[] = new double[2];\r
+       double CT[] = new double[2];\r
+       double PHI[] = new double[2];\r
+       double TERM[] = new double[2];\r
+       double XI[] = new double[2];\r
+       double WW[] = new double[2];\r
+       //COMPLEX CT0,CT2,CTA,CTB,FP,CCJACS,CSUM,CT,JACSUC,PARFUN,PHI,TERM,\r
+       //+XI,WW\r
+        double XENPT[] = new double[MQIN1];\r
+        double JCOFC[][] = new double[MNCOF][2];\r
+        double WSPEC[][] = new double[MXNQD][2];\r
+        double ZSPEC[][] = new double[MXNQD][2];\r
+       //COMPLEX JCOFC(MNCOF),WSPEC(MXNQD),ZSPEC(MXNQD)\r
+       //EXTERNAL CCJACS,JACSUC,PARFUN,PPSBI1,R1MACH\r
+        double cr[] = new double[1];\r
+        double ci[] = new double[1];\r
+        double cr2[] = new double[1];\r
+        double ci2[] = new double[1];\r
+       \r
+       EPS_CATPH4=10.0 * EPS;\r
+       PI=Math.PI;\r
+       TUPI=2.0*PI;\r
+       LOL=NARCS*NQPTS;\r
+       /*loopIP: for (IP=1; IP <= NPTS; IP++) {\r
+           WW[0]=CANPT[IP-1][0];\r
+           WW[1]=CANPT[IP-1][1];\r
+           AWW=zabs(WW[0],WW[1]);\r
+           if (AWW <= EPS_CATPH4) {\r
+               PHYPT[IP-1][0]=CENTR[0];\r
+               PHYPT[IP-1][1]=CENTR[1];\r
+               continue;\r
+           }\r
+           RLW=Math.log(AWW);\r
+           ILW=Math.atan2(WW[1],WW[0]);\r
+           while (ILW > THET0+TUPI) {\r
+               ILW=ILW-TUPI;\r
+           }\r
+           while (ILW < THET0) {\r
+               ILW=ILW+TUPI;\r
+           }\r
+           CSUM[0] = 0.0;\r
+           CSUM[1] = 0.0;\r
+           for (IA=1; IA <= TNSUC; IA++) {\r
+       \r
+               // PRELIMINARIES FOR ARC IA\r
+       \r
+               HA=(VARGC[IA]-VARGC[IA-1])*0.5;\r
+               MD=(VARGC[IA]+VARGC[IA-1])*0.5;\r
+               EFPTL=Math.max(PTHTL,EPS_CATPH4/HA);\r
+               IMXI=-RLW/HA;\r
+       \r
+               if (ILW > (MD+PI)) {\r
+                   ARGW=ILW-TUPI;\r
+               }\r
+               else if (ILW < (MD-PI)) {\r
+                   ARGW=ILW+TUPI;\r
+               }\r
+               else {\r
+                   ARGW=ILW;\r
+               }\r
+       \r
+               REXI=(ARGW-MD)/HA; \r
+               if (REXI > 1.0) {\r
+                   DIST=Math.sqrt(IMXI*IMXI+(REXI-1.0)*(REXI-1.0));\r
+               }\r
+               else if (REXI < -1E+0) {\r
+                   DIST=Math.sqrt(IMXI*IMXI+(REXI+1.0)*(REXI+1.0));\r
+               }\r
+               else {\r
+                   DIST=Math.abs(IMXI);\r
+               }\r
+       \r
+               if (DIST >= DELTA) {\r
+       \r
+                   // USE THE STANDARD PRECOMPUTED COMPOSITE GAUSSIAN RULE\r
+       \r
+                   NQ=NPPQG[IA-1];\r
+                   K=LQSBG[IA-1]-1;\r
+                   for (JQ=1;JQ <= NQ; JQ++) {\r
+                       K=K+1;\r
+                       zdiv(WW[0],WW[1],ZPPQG[K-1][0],ZPPQG[K-1][1],cr,ci);\r
+                       CT[0] = cr[0];\r
+                       CT[1] = ci[0];\r
+                       if (! INTER) {\r
+                               zdiv(1.0,0.0,CT[0],CT[1],cr,ci);\r
+                               CT[0] = cr[0];\r
+                               CT[1] = ci[0];\r
+                       }\r
+                       double LOGR = Math.log(zabs(1.0-CT[0],-CT[1]));\r
+                       double LOGI = Math.atan2(-CT[1], 1.0 - CT[0]);\r
+                       zmlt(WPPQG[K-1][0],WPPQG[K-1][1],LOGR,LOGI,cr,ci);\r
+                       CSUM[0] = CSUM[0] + cr[0];\r
+                       CSUM[1] = CSUM[1] + ci[0];\r
+                   } // for (JQ=1;JQ <= NQ; JQ++)\r
+               } // if (DIST >= DELTA)\r
+               else if (DIST < EFPTL) {\r
+       \r
+                   // WW IS PATHOLOGICALLY CLOSE TO ARC IA (OR MAYBE JUST CLOSE TO\r
+                   // ARC IA AND APPARENTLY OUTSIDE THE CANONICAL DOMAIN) AND WE \r
+                   // USE THE CONTINUATION OF THE BOUNDARY CORRESPONDENCE FUNCTION\r
+                   // TO ESTIMATE PHYPT.\r
+       \r
+                   // INITIALISE SOME DATA\r
+       \r
+                   XI[0] = REXI;\r
+                   XI[1] = IMXI;\r
+                   DGC=DGPOC[IA-1];\r
+                   if (DGC+1 > MNCOF) {\r
+                       IER[0]=48;\r
+                       return;\r
+                   }\r
+                   JTC=JTYPC[IA-1];\r
+                   AJTC=Math.abs(JTC);\r
+                   BETAC=JAINC[AJTC-1];\r
+                   LOMC=LSUBC[IA-1];\r
+                   LODC=(AJTC-1)*NQPTS+1;\r
+                   if (JTC >= 0) {\r
+                       SJTC = 1.0;\r
+                   }\r
+                   else {\r
+                       SJTC = -1.0;\r
+                   }\r
+                   PSA=PRNSA[IA-1];\r
+       \r
+                   if (PHPAS[IA] <= PHPAS[IA-1]) {\r
+                       BB=1.0;\r
+                   }\r
+                   else {\r
+                       BB=PHPAS[IA];\r
+                   }\r
+                   AA=PHPAS[IA-1];\r
+       \r
+                   if (INTER) {\r
+                       S2C=SJTC;\r
+                   }\r
+                   else {\r
+                       S2C=-SJTC;\r
+                   }\r
+       \r
+                   CTA[0]=MIDPT[PSA-1]+AA*HALEN[PSA-1];\r
+                   CTA[1] = 0.0;\r
+                   CTB[0]=MIDPT[PSA-1]+BB*HALEN[PSA-1];\r
+                   CTB[1] = 0.0;\r
+                   if (SJTC > 0) {\r
+                     CT0[0]=CTA[0];\r
+                     CT0[1]=CTA[1];\r
+                     CT2[0]=CTB[0];\r
+                     CT2[1]=CTB[1];\r
+                   }\r
+                   else {\r
+                     CT0[0]=CTB[0];\r
+                     CT0[1]=CTB[1];\r
+                     CT2[0]=CTA[0];\r
+                     CT2[1]=CTA[1];\r
+                   }\r
+       \r
+                   TERM[0]=1.+0+SJTC*XI[0];\r
+                   TERM[1] = SJTC*XI[1];\r
+                   if ((TERM[0] == 0.0) && (TERM[1] == 0.0)) {\r
+                       PHYPT[IP-1]=PARFUN(PARNT[PSA-1],CT0);\r
+                       continue IPloop;\r
+                   }\r
+    \r
+                   IF (TERM.EQ.(2E+0,0E+0)) THEN\r
+                     PHYPT(IP)=PARFUN(PARNT(PSA),CT2)\r
+                     GOTO 300\r
+                   ENDIF\r
+       C\r
+                   ARG=ATAN2(AIMAG(TERM),REAL(TERM))\r
+                   UPPER=(1E+0+S2C*5E-1)*PI\r
+                   IF (ARG.GT.UPPER) THEN\r
+                     ARG=ARG-TUPI\r
+                   ELSE IF (ARG.LE.(UPPER-TUPI)) THEN\r
+                     ARG=ARG+TUPI\r
+                   ENDIF\r
+       C\r
+                   FP=ABS(TERM)**(BETAC+1E+0)*CMPLX(COS((BETAC+1E+0)*ARG),\r
+            +                                       SIN((BETAC+1E+0)*ARG))\r
+                   IF (INTER) THEN\r
+                     FP=-FP\r
+                   ENDIF\r
+                   PHI=CCJACS(SJTC*XI,DGC-1,A1COC(LODC),B1COC(LODC),\r
+            +                 H1VLC(AJTC),BFSNC(LOMC+1))\r
+                   PHI=(BFSNC(LOMC)-(1E+0-SJTC*XI)*PHI)*SJTC*CMPLX(0E+0,1E+0)\r
+                   PHYPT(IP)=CENTR+(PARFUN(PARNT(PSA),CT0)-CENTR)*EXP(FP*PHI)\r
+                   GOTO 300\r
+               } // else if (DIST < EFPTL)\r
+               else {\r
+       C\r
+       C           SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS\r
+       C           PARTICULAR POINT WW.\r
+       C\r
+       C           INITIALISE SOME DATA\r
+       C\r
+                   DGC=DGPOC(IA)\r
+                   IF (DGC+1 .GT. MNCOF) THEN\r
+                     IER=48\r
+                     RETURN\r
+                   ENDIF\r
+                   JTC=JTYPC(IA)\r
+                   AJTC=ABS(JTC)\r
+                   BETAC=JAINC(AJTC)\r
+                   LOMC=LSUBC(IA)\r
+                   LODC=(AJTC-1)*NQPTS+1\r
+                   SJTC=SIGN(1E+0,REAL(JTC))\r
+                   SCO=SJTC\r
+       C\r
+                   DO 100 J=1,DGC+1\r
+                     J1=LOMC+J-1\r
+                     SCO=SCO*SJTC\r
+                     JCOFC(J)=SOLNC(J1)*SCO\r
+       100         CONTINUE\r
+       C\r
+                   XI=SJTC*CMPLX(REXI,IMXI)            \r
+                   CALL PPSBI1(XI,BETAC,NQPTS,DGC,ACOFC(LODC),BCOFC(LODC),\r
+            +                  H0VLC(AJTC),JCOFC,LGTOL,TOLOU,XENPT,QINTS,\r
+            +                  MQIN1,IER)\r
+                   IF (IER .GT. 0) THEN\r
+                     IF (IER .EQ. 29) THEN\r
+                       IER=49\r
+                     ENDIF\r
+                     RETURN\r
+                   ENDIF\r
+                   NQUAD=QINTS*NQPTS\r
+                   IF (NQUAD .GT. MXNQD) THEN\r
+                     IER=47\r
+                     RETURN\r
+                   ENDIF\r
+                   K=0\r
+                   SUM1=BETAC+1E+0\r
+                   DO 130 I=1,QINTS\r
+                     RR=(XENPT(I+1)-XENPT(I))*5E-1\r
+                     MEAN=(XENPT(I+1)+XENPT(I))*5E-1\r
+                     IF (I .EQ. 1) THEN\r
+                       RRB=RR**SUM1\r
+                       DO 110 J=1,NQPTS\r
+                         J1=LODC+J-1\r
+                         K=K+1\r
+                         TT=MEAN+RR*QUPTC(J1)\r
+                         WSPEC(K)=RRB*QUWTC(J1)*JACSUC(TT,DGC,ACOFC(LODC),\r
+            +                     BCOFC(LODC),H0VLC(AJTC),JCOFC)\r
+                         TT=MD+TT*SJTC*HA\r
+                         ZSPEC(K)=CMPLX(COS(TT),SIN(TT))\r
+       110             CONTINUE\r
+                     ELSE\r
+                       DO 120 J=1,NQPTS\r
+                         J1=LOL+J\r
+                         K=K+1\r
+                         TT=MEAN+RR*QUPTC(J1)\r
+                         WSPEC(K)=RR*QUWTC(J1)*(1E+0+TT)**BETAC*JACSUC(TT,\r
+            +                     DGC,ACOFC(LODC),BCOFC(LODC),H0VLC(AJTC),\r
+            +                     JCOFC)\r
+                         TT=MD+TT*SJTC*HA\r
+                         ZSPEC(K)=CMPLX(COS(TT),SIN(TT))\r
+       120             CONTINUE\r
+                     ENDIF\r
+       130         CONTINUE\r
+       C\r
+       C           THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS\r
+       C           AND POINTS WSPEC AND ZSPEC.  NOW ESTIMATE THE INTEGRAL.\r
+       C\r
+                   DO 140 K=1,NQUAD\r
+                     CT=WW/ZSPEC(K)\r
+                     IF (.NOT. INTER) THEN\r
+                       CT=1E+0/CT\r
+                     ENDIF\r
+                     CSUM=CSUM+WSPEC(K)*CLOG(1E+0-CT)\r
+       140         CONTINUE\r
+       \r
+                    // END OF ELSE BLOCK RELATING TO SPECIAL QUADRATURE RULE FOR\r
+                    // WW NEAR ARC IA\r
+       \r
+               } // else             \r
+       \r
+                // END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA\r
+       \r
+           } // for (IA=1; IA <= TNSUC; IA++)\r
+           double EXPBASE = Math.exp(CSUM[0]);\r
+           double EXPR = EXPBASE * Math.cos(CSUM[1]);\r
+           double EXPI = EXPBASE * Math.sin(CSUM[1]);\r
+           zmlt(FACTR[0],FACTR[1],WW[0],WW[1],cr,ci);\r
+           zmlt(cr[0],ci[0],EXPR,EXPI,cr2,ci2);\r
+           PHYPT[IP-1][0] = CENTR[0] + cr2[0];\r
+           PHYPT[IP-1][1] = CENTR[1] + ci2[0];\r
+       \r
+           // END OF MAP CALCULATION FOR FIELD POINT NUMBER IP\r
+       \r
+       } // loopIP: for (IP=1; IP <= NPTS; IP++) */\r
+       \r
+       IER[0]=0;\r
+       \r
+    } // private void CATPH4\r
 \r
 \r
       /**\r