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

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

index b2e5827..4c0113d 100644 (file)
@@ -8061,7 +8061,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 \r
     } // private void RESCAL\r
 \r
-    private void GQPHYC(int MQIN1, int IER[]) {\r
+    private void GQPHYC(int MQIN1, double RWORK[], int IER[]) {\r
+       // not the RWORK set up by JAPHYC\r
        \r
        // INTEGER MQIN1,MQUPH,CHNL,IER\r
        // INTEGER IQUPH(*),IGEOM(*),ISNPH(*)\r
@@ -8188,14 +8189,14 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        //     LOCAL VARAIBLES\r
        \r
        final int NINTS = 5;\r
-       int ACOEF,BCOEF,DGPOL,FACTR,H0VAL,HALEN,JACIN,JATYP,LOSUB,\r
-            LQSBF,MIDPT,NPPQF,PARNT,QUPTS,QUWTS,SOLUN,TPPQF,TRRAD,VTARG,WPPQF,\r
-            ZPPQF;\r
-       int MNEQN,MNSUA,NARCS,NEQNS,NJIND,NQPTS,TNGQP,TNPQP,\r
+       int TNPQP[] = new int[1];\r
+       int MNSUA,NEQNS,\r
             TNSUA;\r
-       final double DLETA = 0.2;\r
-       double LGTOL,SUPER,THET0,TOLOU;\r
+       final double DELTA = 0.2;\r
+       double LGTOL,SUPER,THET0;\r
        double CT[] = new double[2];\r
+       double theta;\r
+       double mag;\r
        //COMPLEX CT\r
        \r
         // EXTERNAL POPQF1,OUPTPQ,WRHEAD,WRTAIL\r
@@ -8224,7 +8225,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        // SET UP POINTERS FOR QUADRATURE DATA.\r
        \r
        \r
-       /*RQUPH[0]=SOLUN[NEQNS-1];\r
+       RQUPH[0]=SOLUN[NEQNS-1];\r
        \r
        System.out.println("QUADRATURE RULES STARTED:");\r
        POPQF1(NPPQF,LQSBF,TNPQP,TOLOU,TPPQF,\r
@@ -8233,40 +8234,54 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
             DELTA,LGTOL,ACOEF,BCOEF,H0VAL,HALEN,\r
             JACIN,MIDPT,QUPTS,QUWTS,SOLUN,\r
             RWORK,IER);\r
-             WRITE(*,10) 'QUADRATURE RULES DONE:'\r
-       C\r
-             IF (IER .GT. 0) THEN\r
-               GOTO 999\r
-             ENDIF\r
-       C\r
-             IQUPH(1)=TNPQP\r
-       C\r
-       C**** SET UP THE CONSTANT FACTOR FOR THE MAPPING FORMULA\r
-       C\r
-             THET0=RGEOM(VTARG)\r
-             IF (INTER) THEN\r
-               ZQUPH(FACTR)=CMPLX(COS(THET0),SIN(THET0))\r
-             ELSE\r
-               ZQUPH(FACTR)=(1E+0,0E+0)\r
-               CALL DMPHYC(1,CENTR,CT,INTER,CENTR,IGEOM,RGEOM,ISNPH,RSNPH,\r
-            +              IQUPH,RQUPH,ZQUPH,.TRUE.,IER)\r
-               IF (IER .GT. 0) THEN\r
-                 GOTO 999\r
-               ENDIF\r
-               CT=CMPLX(RQUPH(1),THET0-AIMAG(LOG(CT)))\r
-               ZQUPH(FACTR)=CEXP(CT)\r
-             ENDIF\r
-       C\r
-             CALL OUPTPQ(IQUPH,RQUPH,ZQUPH,CHNL)\r
-       999   CONTINUE\r
-       C\r
-       C**** WRITE CLOSING MESSAGE TO STANDARD OUTPUT CHANNEL\r
-       C\r
-             CALL WRTAIL(2,0,IER)\r
-       C*/\r
+       System.out.println("QUADRATURE RULES DONE:");\r
+       \r
+       if (IER[0] > 0) {\r
+           WRTAIL(2,0,IER[0],null);\r
+           return;\r
+       }\r
+       \r
+       IQUPH[0]=TNPQP[0];\r
+       \r
+       // SET UP THE CONSTANT FACTOR FOR THE MAPPING FORMULA\r
+       \r
+       THET0=VTARG[0];\r
+       if (INTER) {\r
+           FACTR[0]=Math.cos(THET0);\r
+           FACTR[1] = Math.sin(THET0);\r
+       }\r
+       else {\r
+           FACTR[0]=1.0;\r
+           FACTR[1]= 0.0;\r
+           \r
+           double centrArr[][]= new double[1][2];\r
+           centrArr[0][0] = CENTR[0];\r
+           centrArr[0][1] = CENTR[1];\r
+           double ctArr[][] = new double[1][2];\r
+           DMPHYC(1,centrArr,ctArr,true,IER);\r
+           CT[0] = ctArr[0][0];\r
+           CT[1] = ctArr[0][1];\r
+           if (IER[0] > 0) {\r
+               WRTAIL(2,0,IER[0],null);\r
+                   return;\r
+           }\r
+           theta = Math.atan2(CT[1],CT[0]);\r
+           CT[0] = RQUPH[0];\r
+           CT[1] = THET0 - theta;\r
+           mag = Math.exp(CT[0]);\r
+           FACTR[0] = mag * Math.cos(CT[1]);\r
+           FACTR[1] = mag * Math.sin(CT[1]);\r
+       }\r
+       \r
+       // OUPTPQ(IQUPH,RQUPH,ZQUPH,CHNL);\r
+       \r
+        // WRITE CLOSING MESSAGE TO STANDARD OUTPUT CHANNEL\r
+       \r
+       WRTAIL(2,0,IER[0],null);\r
+       \r
     } // private void GQPHYC\r
 \r
-    private void POPQF1(int NPPQF[], int LPQSB[], int TNPQP[], double TOLOU, double TPPQF[],\r
+    private void POPQF1(int NPPQF[], int LPQSB[], int TNPQP[], double TOLOU[], double TPPQF[],\r
         double TRRAD[], double WPPQF[], double ZPPQF[][], int MNQUA, int MQIN1, int NARCS,\r
         int NINTS, int NQPTS,int TNSUA, int DGPOL[], int JATYP[], int LOSUB[], int PARNT[],\r
        double DELTA, double LGTOL, double ACOEF[], double BCOEF[], double H0VAL[], double HALEN[],\r
@@ -8317,20 +8332,22 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        //     LOCAL VARIABLES\r
        \r
        final int MNCOF = 32;\r
-       int AJT,DEG,HI,HI1,I,I1,J,J1,J2,JT,K,LIM,LOD,LOL,LOM,PT,QINTS;\r
+       int QINTS[] = new int[1];\r
+       int AJT,DEG,HI,HI1,I,I1,J,J1,J2,JT,K,LIM,LOD,LOL,LOM,PT;\r
        final double RHO = 0.13;\r
-       double BETA,DIST,HL,JACSUM,MD,MEAN,MXDIS,RR,RRB,SS,SUM1,TT,SCO;\r
+       double BETA,DIST,HL,MD,MEAN,MXDIS,RR,RRB,SS,SUM1,TT,SCO;\r
        double CTT[] = new double[2];\r
        // COMPLEX CTT,DPARFN,PARFUN\r
        double JACOF[] = new double[MNCOF];\r
        double JCOFC[][] = new double[MNCOF][2];\r
+       double DOUT[];\r
        // COMPLEX JCOFC(MNCOF)\r
        // EXTERNAL DPARFN,JACSUM,PARFUN,PPSBI7\r
        \r
        MXDIS=DELTA;\r
        HI=0;\r
        LOL=NARCS*NQPTS;\r
-       /*for (I1=1; I1 <= TNSUA; I1++) {\r
+       for (I1=1; I1 <= TNSUA; I1++) {\r
            JT=JATYP[I1-1];\r
            if (JT > 0) {\r
                SS=1.0;\r
@@ -8360,133 +8377,158 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                JCOFC[J-1][1] = 0.0;\r
            } // for (J=1; J <= DEG+1; J++)\r
        \r
-           PPSBI7(DELTA,NINTS,BETA,NQPTS,DEG,ACOEF(LOD),\r
-            BCOEF(LOD),H0VAL(AJT),JCOFC,LGTOL,TOLOU,XENPT,\r
+           double AIN[] = new double[DEG];\r
+           double BIN[] = new double[DEG];\r
+           for (I = 0; I < DEG; I++) {\r
+               AIN[I] = ACOEF[LOD-1+I];\r
+               BIN[I] = BCOEF[LOD-1+I];\r
+           }\r
+           double HIN = H0VAL[AJT-1];\r
+           PPSBI7(DELTA,NINTS,BETA,NQPTS,DEG,AIN,\r
+            BIN,HIN,JCOFC,LGTOL,TOLOU,XENPT,\r
             QINTS,MQIN1,IER);\r
-               IF (IER .GT. 0) THEN\r
-                 RETURN\r
-               ENDIF\r
-               NPPQF(I1)=QINTS*NQPTS\r
-               LPQSB(I1)=HI+1\r
-               HI1=HI+NPPQF(I1)\r
-               IF (HI1 .GT. MNQUA) THEN\r
-                 IER=22\r
-                 RETURN\r
-               ENDIF\r
-               K=HI\r
-               SUM1=BETA+1E+0\r
-               DO 30 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 10 J=1,NQPTS\r
-                     J1=LOD+J-1\r
-                     K=K+1\r
-                     TT=(MEAN+RR*QUPTS(J1))\r
-                     WPPQF(K)=RRB*QUWTS(J1)*JACSUM(TT,DEG,ACOEF(LOD),\r
-            +                 BCOEF(LOD),H0VAL(AJT),JACOF)\r
-                     TT=TT*SS\r
-                     TPPQF(K)=TT\r
-                     CTT=CMPLX(MD+TT*HL)\r
-                     ZPPQF(K)=PARFUN(PT,CTT)\r
-                     TRRAD(K)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
-       10          CONTINUE\r
-                 ELSE\r
-                   DO 20 J=1,NQPTS\r
-                     J1=LOL+J\r
-                     K=K+1\r
-                     TT=(MEAN+RR*QUPTS(J1))\r
-                     WPPQF(K)=RR*QUWTS(J1)*(1E+0+TT)**BETA*JACSUM(TT,DEG,\r
-            +                 ACOEF(LOD),BCOEF(LOD),H0VAL(AJT),JACOF)\r
-                     TT=TT*SS\r
-                     TPPQF(K)=TT\r
-                     CTT=CMPLX(MD+TT*HL)\r
-                     ZPPQF(K)=PARFUN(PT,CTT)\r
-                     TRRAD(K)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
-       20          CONTINUE\r
-                 ENDIF\r
-       30      CONTINUE\r
-               IF (SS .LT. 0E+0) THEN\r
-                 LIM=NPPQF(I1)\r
-                 IF (MOD(LIM,2) .EQ. 0) THEN\r
-                   LIM=LIM/2\r
-                 ELSE\r
-                   LIM=(LIM-1)/2\r
-                 ENDIF\r
-                 J1=LPQSB(I1)-1\r
-                 J2=HI1+1\r
-                 DO 40 J=1,LIM\r
-                   J1=J1+1\r
-                   J2=J2-1\r
-                   TT=WPPQF(J1)\r
-                   WPPQF(J1)=WPPQF(J2)\r
-                   WPPQF(J2)=TT\r
-                   TT=TRRAD(J1)\r
-                   TRRAD(J1)=TRRAD(J2)\r
-                   TRRAD(J2)=TT\r
-                   TT=TPPQF(J1)\r
-                   TPPQF(J1)=TPPQF(J2)\r
-                   TPPQF(J2)=TT\r
-                   CTT=ZPPQF(J1)\r
-                   ZPPQF(J1)=ZPPQF(J2)\r
-                   ZPPQF(J2)=CTT\r
-       40        CONTINUE\r
-               ENDIF\r
-       C\r
-       C       NEXT WE INSERT ANY NECESSARY FICTICIOUS QUADRATURE POINTS\r
-       C\r
-               J1=LPQSB(I1)\r
-               IF (TPPQF(J1)+1E+0 .GT. RHO*MXDIS) THEN\r
-                 J2=HI1\r
-                 DO 50 I=J2,J1,-1\r
-                   WPPQF(I+1)=WPPQF(I)\r
-                   TRRAD(I+1)=TRRAD(I)\r
-                   TPPQF(I+1)=TPPQF(I)\r
-                   ZPPQF(I+1)=ZPPQF(I)\r
-       50        CONTINUE\r
-                 HI1=J2+1\r
-                 NPPQF(I1)=NPPQF(I1)+1\r
-                 WPPQF(J1)=0E+0\r
-                 TPPQF(J1)=-1E+0\r
-                 CTT=CMPLX(MD-HL)\r
-                 ZPPQF(J1)=PARFUN(PT,CTT)\r
-                 TRRAD(J1)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
-               ENDIF\r
-       C\r
-       60      CONTINUE\r
-               J1=J1+1\r
-               IF (J1 .LT. HI1) THEN\r
-                 DIST=TPPQF(J1+1)-TPPQF(J1)\r
-                 IF (DIST .GT. MXDIS) THEN\r
-                   J2=HI1\r
-                   DO 70 I=J2,J1+1,-1\r
-                     WPPQF(I+1)=WPPQF(I)\r
-                     TRRAD(I+1)=TRRAD(I)\r
-                     TPPQF(I+1)=TPPQF(I)\r
-                     ZPPQF(I+1)=ZPPQF(I)\r
-       70          CONTINUE\r
-                   HI1=J2+1\r
-                   NPPQF(I1)=NPPQF(I1)+1\r
-                   WPPQF(J1+1)=0E+0\r
-                   TPPQF(J1+1)=(TPPQF(J1)+TPPQF(J1+2))*5E-1\r
-                   CTT=CMPLX(MD+HL*TPPQF(J1+1))\r
-                   ZPPQF(J1+1)=PARFUN(PT,CTT)\r
-                   TRRAD(J1+1)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
-                 ENDIF\r
-                 GOTO 60\r
-               ELSE IF (1E+0-TPPQF(J1) .GT. RHO*MXDIS) THEN\r
-                 J1=J1+1\r
-                 HI1=J1\r
-                 NPPQF(I1)=NPPQF(I1)+1\r
-                 WPPQF(J1)=0E+0\r
-                 TPPQF(J1)=1E+0\r
-                 CTT=CMPLX(MD+HL)\r
-                 ZPPQF(J1)=PARFUN(PT,CTT)\r
-                 TRRAD(J1)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
-               ENDIF\r
-               HI=HI1\r
-       } // for (I1=1; I1 <= TNSUA; I1++)*/\r
+           if (IER[0] > 0) {\r
+               return;\r
+           }\r
+           NPPQF[I1-1]=QINTS[0]*NQPTS;\r
+           LPQSB[I1-1]=HI+1;\r
+           HI1=HI+NPPQF[I1-1];\r
+           if (HI1 > MNQUA) {\r
+               IER[0]=22;\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*QUPTS[J1-1]);\r
+                       WPPQF[K-1]=RRB*QUWTS[J1-1]*JACSUM(TT,DEG,AIN,BIN,HIN,JACOF);\r
+                       TT=TT*SS;\r
+                       TPPQF[K-1]=TT;\r
+                       CTT[0]=MD+TT*HL;\r
+                       CTT[1] = 0.0;\r
+                       ZPPQF[K-1]=PARFUN(PT,CTT);\r
+                       DOUT = DPARFN(PT,CTT);\r
+                       TRRAD[K-1]=HL*DELTA*zabs(DOUT[0],DOUT[1]);\r
+                   } // for (J=1; J <= NQPTS; J++)\r
+               } // if (I == 1)\r
+               else {\r
+                   for (J=1; J <= NQPTS; J++) {\r
+                       J1=LOL+J;\r
+                       K=K+1;\r
+                       TT=(MEAN+RR*QUPTS[J1-1]);\r
+                       WPPQF[K-1]=RR*QUWTS[J1-1]*Math.pow((1E+0+TT),BETA)*JACSUM(TT,DEG,\r
+                             AIN,BIN,HIN,JACOF);\r
+                       TT=TT*SS;\r
+                       TPPQF[K-1]=TT;\r
+                       CTT[0]=MD+TT*HL;\r
+                       CTT[1] = 0.0;\r
+                       ZPPQF[K-1]=PARFUN(PT,CTT);\r
+                       DOUT = DPARFN(PT,CTT);\r
+                       TRRAD[K-1]=HL*DELTA*zabs(DOUT[0],DOUT[1]);\r
+                   } // for (J=1; J <= NQPTS; J++)\r
+               } // else\r
+           } // for (I=1; I <= QINTS[0]; I++)\r
+           if (SS < 0.0) {\r
+               LIM=NPPQF[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
+                   TT=WPPQF[J1-1];\r
+                   WPPQF[J1-1]=WPPQF[J2-1];\r
+                   WPPQF[J2-1]=TT;\r
+                   TT=TRRAD[J1-1];\r
+                   TRRAD[J1-1]=TRRAD[J2-1];\r
+                   TRRAD[J2-1]=TT;\r
+                   TT=TPPQF[J1-1];\r
+                   TPPQF[J1-1]=TPPQF[J2-1];\r
+                   TPPQF[J2-1]=TT;\r
+                   CTT[0]=ZPPQF[J1-1][0];\r
+                   CTT[1]=ZPPQF[J1-1][1];\r
+                   ZPPQF[J1-1][0]=ZPPQF[J2-1][0];\r
+                   ZPPQF[J1-1][1]=ZPPQF[J2-1][1];\r
+                   ZPPQF[J2-1][0]=CTT[0];\r
+                   ZPPQF[J2-1][1]=CTT[1];\r
+               } // for (J=1; J <= LIM; J++)\r
+           } // if (SS < 0.0)\r
+       \r
+           // NEXT WE INSERT ANY NECESSARY FICTICIOUS QUADRATURE POINTS\r
+       \r
+           J1=LPQSB[I1-1];\r
+           if (TPPQF[J1-1]+1.0 > RHO*MXDIS) {\r
+               J2=HI1;\r
+               for (I=J2; I >= J1; I--) {\r
+                   WPPQF[I]=WPPQF[I-1];\r
+                   TRRAD[I]=TRRAD[I-1];\r
+                   TPPQF[I]=TPPQF[I-1];\r
+                   ZPPQF[I][0]=ZPPQF[I-1][0];\r
+                   ZPPQF[I][1]=ZPPQF[I-1][1];\r
+               } // for (I=J2; I >= J1; I--)\r
+               HI1=J2+1;\r
+               NPPQF[I1-1]=NPPQF[I1-1]+1;\r
+               WPPQF[J1-1]=0.0;\r
+               TPPQF[J1-1]=-1.0;\r
+               CTT[0]=MD-HL;\r
+               ZPPQF[J1-1]=PARFUN(PT,CTT);\r
+               DOUT = DPARFN(PT,CTT);\r
+               TRRAD[J1-1]=HL*DELTA*zabs(DOUT[0],DOUT[1]);\r
+           } // if (TPPQF[J1-1]+1.0 > RHO*MXDIS)\r
+       \r
+           while (true) {\r
+               J1=J1+1;\r
+               if (J1 < HI1) {\r
+                   DIST=TPPQF[J1]-TPPQF[J1-1];\r
+                   if (DIST > MXDIS) {\r
+                       J2=HI1;\r
+                       for (I=J2; I >= J1+1; I--) {\r
+                           WPPQF[I]=WPPQF[I-1];\r
+                           TRRAD[I]=TRRAD[I-1];\r
+                           TPPQF[I]=TPPQF[I-1];\r
+                           ZPPQF[I][0]=ZPPQF[I-1][0];\r
+                           ZPPQF[I][1]=ZPPQF[I-1][1];\r
+                       } // for (I=J2; I >= J1+1; I--)\r
+                       HI1=J2+1;\r
+                       NPPQF[I1-1]=NPPQF[I1-1]+1;\r
+                       WPPQF[J1]=0.0;\r
+                       TPPQF[J1]=(TPPQF[J1-1]+TPPQF[J1+1])*0.5;\r
+                       CTT[0]=MD+HL*TPPQF[J1];\r
+                       CTT[1] = 0.0;\r
+                       ZPPQF[J1]=PARFUN(PT,CTT);\r
+                       DOUT = DPARFN(PT,CTT);\r
+                       TRRAD[J1]=HL*DELTA*zabs(DOUT[0],DOUT[1]);\r
+                   } // if (DIST > MXDIS)\r
+                   continue;\r
+               }\r
+               else if (1.0-TPPQF[J1-1] > RHO*MXDIS) {\r
+                   J1=J1+1;\r
+                   HI1=J1;\r
+                   NPPQF[I1-1]=NPPQF[I1-1]+1;\r
+                   WPPQF[J1-1]=0.0;\r
+                   TPPQF[J1-1]=1.0;\r
+                   CTT[0]=MD+HL;\r
+                   CTT[1] = 0.0;\r
+                   ZPPQF[J1-1]=PARFUN(PT,CTT);\r
+                   DOUT = DPARFN(PT,CTT);\r
+                   TRRAD[J1-1]=HL*DELTA*zabs(DOUT[0],DOUT[1]);\r
+               } // else if (1.0-TPPQF[J1-1] > RHO*MXDIS)\r
+               break;\r
+           } // while (true)\r
+           HI=HI1;\r
+       } // for (I1=1; I1 <= TNSUA; I1++)\r
        \r
        TNPQP[0]=HI;\r
        \r
@@ -8495,7 +8537,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     } // private void POPQF1\r
     \r
     private void PPSBI7(double DELTA, int NINTS, double BETA, int NQUAD, int DGPOL, double ACOEF[],\r
-        double BCOEF[], double H0VAL, double SOLUN[][], double TOLIN, double TOLOU, double XENPT[],\r
+        double BCOEF[], double H0VAL, double SOLUN[][], double TOLIN, double TOLOU[], double XENPT[],\r
         int QINTS[], int MQIN1, int IER[]) {\r
        // INTEGER DGPOL,MQIN1,NQUAD,QINTS,IER,NINTS\r
        // REAL BETA,H0VAL,TOLIN,TOLOU,DELTA\r
@@ -8538,66 +8580,73 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        int INTS;\r
        double TAU[] = new double[1];\r
        double MAXRM[] = new double[1];\r
-       double TOL,RIGHT;\r
+       double TOL[] = new double[1];\r
+       double RIGHT[] = new double[1];\r
        boolean T1FXD;\r
+       double TAU1[] = new double[1];\r
        // EXTERNAL DEPPJ8,DEPPL8\r
        \r
        QINTS[0]=1;\r
        XENPT[0]=-1.0;\r
-       TOL=TOLIN;\r
+       TOL[0]=TOLIN;\r
        INTS=NINTS;\r
        DEPPJ8(BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL,\r
             SOLUN,TOL,MAXRM,INTS,DELTA,IER);\r
-       /*      IF (IER .GT. 0) THEN\r
-               RETURN\r
-             ENDIF\r
-             TOLOU=MAXRM\r
-             XENPT(2)=TAU\r
-       C\r
-             IF (XENPT(2) .LT. 1E+0) THEN\r
-               QINTS=2\r
-               T1FXD=.FALSE.\r
-               TAU=1E+0\r
-               RIGHT=-1E+0\r
-               INTS=NINTS\r
-               CALL DEPPL8(BETA,RIGHT,TAU,T1FXD,NQUAD,DGPOL,ACOEF,\r
-            +              BCOEF,H0VAL,SOLUN,TOL,MAXRM,INTS,DELTA,IER)\r
-               IF (IER .GT. 0) THEN\r
-                 RETURN\r
-               ENDIF\r
-               TOLOU=TOLOU+MAXRM\r
-               T1FXD=.TRUE.\r
-       C\r
-       100     CONTINUE\r
-       C\r
-               IF (XENPT(QINTS) .GT. RIGHT) THEN\r
-                 XENPT(QINTS)=5E-1*(XENPT(QINTS)+RIGHT)\r
-                 XENPT(QINTS+1)=1E+0\r
-               ELSE\r
-                 TAU=1E+0\r
-                 INTS=NINTS\r
-                 CALL DEPPL8(BETA,XENPT(QINTS),TAU,T1FXD,NQUAD,DGPOL,\r
-            +                ACOEF,BCOEF,H0VAL,SOLUN,TOL,MAXRM,INTS,DELTA,IER) \r
-                 IF (IER .GT. 0) THEN\r
-                   RETURN\r
-                 ENDIF\r
-                 TOLOU=TOLOU+MAXRM\r
-                 QINTS=QINTS+1\r
-                 IF (QINTS .GE. MQIN1) THEN\r
-                   IER=24\r
-                   RETURN\r
-                 ENDIF\r
-                 XENPT(QINTS)=TAU\r
-                 GOTO 100\r
-               ENDIF\r
-             ENDIF\r
-       C\r
-             IER=0\r
-       C */\r
+       if (IER[0] > 0) {\r
+           return;\r
+       }\r
+       TOLOU[0]=MAXRM[0];\r
+       XENPT[1]=TAU[0];\r
+       \r
+       if (XENPT[1] < 1.0) {\r
+           QINTS[0]=2;\r
+           T1FXD=false;\r
+           TAU[0]=1.0;\r
+           RIGHT[0]=-1.0;\r
+           INTS=NINTS;\r
+           DEPPL8(BETA,RIGHT,TAU,T1FXD,NQUAD,DGPOL,ACOEF,\r
+                   BCOEF,H0VAL,SOLUN,TOL,MAXRM,INTS,DELTA,IER);\r
+           if (IER[0] > 0) {\r
+               return;\r
+           }\r
+           TOLOU[0]=TOLOU[0]+MAXRM[0];\r
+           T1FXD=true;\r
+       \r
+           while (true) {\r
+       \r
+               if (XENPT[QINTS[0]-1] > RIGHT[0]) {\r
+                   XENPT[QINTS[0]-1]=0.5*(XENPT[QINTS[0]-1]+RIGHT[0]);\r
+                   XENPT[QINTS[0]]=1.0;\r
+                   break;\r
+               }\r
+               else {\r
+                   TAU[0]=1.0;\r
+                   INTS=NINTS;\r
+                   TAU1[0] = XENPT[QINTS[0]-1];\r
+                   DEPPL8(BETA,TAU1,TAU,T1FXD,NQUAD,DGPOL,\r
+                          ACOEF,BCOEF,H0VAL,SOLUN,TOL,MAXRM,INTS,DELTA,IER);\r
+                   XENPT[QINTS[0]-1] = TAU1[0];\r
+                   if (IER[0] > 0) {\r
+                       return;\r
+                   }\r
+                   TOLOU[0]=TOLOU[0]+MAXRM[0];\r
+                   QINTS[0]=QINTS[0]+1;\r
+                   if (QINTS[0] >= MQIN1) {\r
+                       IER[0]=24;\r
+                       return;\r
+                   }\r
+                   XENPT[QINTS[0]-1]=TAU[0];\r
+                   continue;\r
+               }\r
+           } // while (true)\r
+       } // if (XENPT[1] < 1.0)\r
+       \r
+       IER[0]=0;\r
+    \r
     } // private void PPSBI7\r
     \r
     private void DEPPJ8(double BETA, double TAU[], int NQUAD, int DGPOL, double ACOEF[],\r
-        double BCOEF[], double H0VAL, double SOLUN[][], double TOL, double MAXRM[], int NINTS,\r
+        double BCOEF[], double H0VAL, double SOLUN[][], double TOL[], double MAXRM[], int NINTS,\r
         double DELTA,int IER[]) {\r
         //INTEGER NQUAD,IER,DGPOL,NINTS\r
        //REAL BETA,TAU,TOL,MAXRM,ACOEF(*),BCOEF(*),H0VAL,DELTA\r
@@ -8725,7 +8774,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     \r
     private void DEPPJ9(double ZZ[][], int NZZ, double BETA, double TAU[], int NQUAD,\r
         int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][],\r
-        double TOL, double MAXRM[], int IER[]) {\r
+        double TOL[], double MAXRM[], int IER[]) {\r
         // INTEGER NQUAD,NZZ,IER,DGPOL\r
         // REAL BETA,TAU,TOL,MAXRM,ACOEF(*),BCOEF(*),H0VAL\r
         // COMPLEX SOLUN(*),ZZ(*)\r
@@ -8860,7 +8909,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         while (true) {\r
 \r
             MAXRM[0]=0.0;\r
-            HTOL=0.5*TOL;\r
+            HTOL=0.5*TOL[0];\r
             for (I=1; I <= NZZ; I++) {\r
                XI[0] = (2.0*ZZ[I-1][0]+1.0-TAU[0])/(1.0+TAU[0]);\r
                XI[1] = (2.0*ZZ[I-1][1])/(1.0+TAU[0]);\r
@@ -8894,7 +8943,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                 MAXRM[0]=Math.max(MAXRM[0],TERM);\r
             } // for (I=1; I <= NZZ; I++)\r
 \r
-            if (MAXRM[0] < TOL) {\r
+            if (MAXRM[0] < TOL[0]) {\r
 \r
                 // ACCURACY IS ACHIEVED, BUT MAYBE TAU COULD BE INCREASED.\r
 \r
@@ -8910,13 +8959,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    \r
                 } // if (MAXRM[0] < HTOL)\r
                 break;\r
-            } // if (MAXRM[0] < TOL)\r
+            } // if (MAXRM[0] < TOL[0)\r
             else {\r
 \r
                 // ACCURACY NOT ACHIEVED AND TAU NEEDS DECREASING.\r
 \r
                if (TAU[0] == 1.0) {\r
-                   TOL=HTOL;\r
+                   TOL[0]=HTOL;\r
                }\r
                UPPER=TAU[0];\r
                TAU[0]=0.5*(LOWER+UPPER);\r
@@ -8983,7 +9032,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     } // private double[] CCJACS\r
 \r
     private void DEPPL8(double BETA, double TAU1[], double TAU2[], boolean T1FXD, int NQUAD,\r
-        int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][], double TOL,\r
+        int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][], double TOL[],\r
         double MAXRM[], int NINTS, double DELTA,int IER[]) {\r
        // INTEGER NQUAD,IER,DGPOL,NINTS\r
        // REAL BETA,TAU1,TAU2,TOL,MAXRM,H0VAL,DELTA\r
@@ -9126,7 +9175,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     } // private void DEPPL8\r
     \r
     private void DEPPL9(double ZZ[][], int NZZ, double BETA, double TAU1[], double TAU2[], boolean T1FXD,\r
-        int NQUAD, int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][], double TOL,\r
+        int NQUAD, int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][], double TOL[],\r
         double MAXRM[], boolean FIRST, int IER[]) {\r
        // INTEGER NQUAD,DGPOL,IER,NZZ\r
        // REAL BETA,TAU1,TAU2,H0VAL,TOL,MAXRM\r
@@ -9187,6 +9236,17 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        // COMPLEX CCJACS\r
        double GG[][] = new double[MAXNZ][2];\r
        double HH[][] = new double[MAXNZ][2];\r
+       double CIN[] = new double[2];\r
+       double r;\r
+       double mag;\r
+       double theta;\r
+       double ang;\r
+       double realVar;\r
+       double imagVar;\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
        // COMPLEX GG(MAXNZ),HH(MAXNZ)\r
        // EXTERNAL GAMMA,CCJACS,LGGAM\r
        \r
@@ -9225,87 +9285,666 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        KK=KK/FORK;\r
        \r
        HCO=Math.sin(PI*BETA)*GAMMA(BETA+1.0)/PI/Math.pow(FORK,BETA);\r
-       /*      FPHI1=CCJACS((-1E+0,0E+0),DGPOL,ACOEF,BCOEF,H0VAL,SOLUN)\r
-             DO 125 I=1,NZZ\r
-               FNPHI=CCJACS(ZZ(I),DGPOL,ACOEF,BCOEF,H0VAL,SOLUN)\r
-               GG(I)=-(1E+0+ZZ(I))**BETA*KK*FNPHI\r
-               HH(I)=-CLOG(1E+0+ZZ(I))*HCO*KK*FPHI1\r
-       125   CONTINUE\r
+       CIN[0] = -1.0;\r
+       CIN[1] = 0.0;\r
+       FPHI1=CCJACS(CIN,DGPOL,ACOEF,BCOEF,H0VAL,SOLUN);\r
+       for (I=1; I <= NZZ; I++) {\r
+           FNPHI=CCJACS(ZZ[I-1],DGPOL,ACOEF,BCOEF,H0VAL,SOLUN);\r
+           r = zabs(1.0+ZZ[I-1][0],ZZ[I-1][1]);\r
+           mag = Math.pow(r, BETA);\r
+           theta = Math.atan2(ZZ[I-1][1], 1.0 + ZZ[I-1][0]);\r
+           ang = theta * BETA;\r
+           realVar = -mag * Math.cos(ang);\r
+           imagVar = -mag * Math.sin(ang);\r
+           zmlt(realVar,imagVar,KK*FNPHI[0],KK*FNPHI[1],cr,ci);\r
+           GG[I-1][0] = cr[0];\r
+           GG[I-1][1] = ci[0];\r
+           mag = -Math.log(r);\r
+           zmlt(mag,-theta,HCO*KK*FPHI1[0],HCO*KK*FPHI1[1],cr,ci);\r
+           HH[I-1][0] = cr[0];\r
+           HH[I-1][1] = ci[0];\r
+       } // for (I=1; I <= NZZ; I++)\r
               \r
-       C\r
-       C     NOW COME THE FACTORS DEPENDENT ON TAU1 AND TAU2.\r
-       C\r
-             LOWER=TAU1\r
-             UPPER=TAU2\r
-       C\r
-       250   CONTINUE\r
-       C\r
-             HTOL=5E-1*TOL\r
-             RR=(TAU2-TAU1)*5E-1\r
-             MEAN=(TAU1+TAU2)*5E-1\r
-             BB=(1E+0+MEAN)/RR\r
-             BB1=BB+SQRT(BB*BB-1E+0)\r
-             FFH=(RR*(BB1-1E+0/BB1))**(BETA+1E+0)/BB1**TUK\r
-             MAXRM=0E+0\r
-             DO 325 I=1,NZZ\r
-               XI=(ZZ(I)-MEAN)/RR\r
-               Z1=SQRT(XI*XI-1E+0)\r
-               XI1=XI+Z1\r
-               IF (ABS(XI1) .LT. 1E+0) THEN\r
-                 XI1=XI-Z1\r
-               ENDIF\r
-               FFG=XI1**(-TUK-1E+0)*(XI1*XI1-1E+0)*RR\r
-               REMND=GG(I)*FFG+HH(I)*FFH\r
-               TERM=ABS(REMND)\r
-               MAXRM=MAX(MAXRM,TERM)\r
-       325   CONTINUE\r
-       C\r
-             IF (MAXRM .LT. TOL) THEN\r
-       C\r
-       C       ACCURACY IS ACHIEVED, BUT MAYBE TAU2 COULD BE INCREASED OR\r
-       C       TAU1 DECREASED\r
-       C\r
-               IF (MAXRM .LT. HTOL) THEN\r
-       C\r
-       C         TAU2 NEEDS INCREASING IF T1FXD (BUT THIS IS ONLY POSSIBLE IF \r
-       C         TAU2<TAU2I) OR TAU1 NEED DECREASING OTHERWISE (BUT THIS IS \r
-       C         ONLY POSSIBLE IF TAU1>TAU1I)\r
-       C\r
-                 IF (T1FXD .AND. TAU2 .LT. TAU2I) THEN\r
-                   LOWER=TAU2\r
-                   TAU2=5E-1*(LOWER+UPPER)\r
-                   GOTO 250\r
-                 ELSE IF (.NOT.T1FXD .AND. TAU1 .GT. TAU1I) THEN\r
-                   UPPER=TAU1\r
-                   TAU1=5E-1*(LOWER+UPPER)\r
-                   GOTO 250\r
-                 ENDIF\r
-               ENDIF\r
-             ELSE\r
-       C\r
-       C       ACCURACY NOT ACHIEVED AND TAU2 NEEDS DECREASING OR TAU1 NEEDS\r
-       C       INCREASING.\r
-       C\r
-               IF (FIRST) THEN\r
-                 TOL=HTOL\r
-                 FIRST=.FALSE.\r
-               ENDIF\r
-               IF (T1FXD) THEN\r
-                 UPPER=TAU2\r
-                 TAU2=5E-1*(LOWER+UPPER)\r
-               ELSE\r
-                 LOWER=TAU1\r
-                 TAU1=5E-1*(LOWER+UPPER)\r
-               ENDIF\r
-               GOTO 250\r
-             ENDIF\r
-       C\r
-       C     NORMAL EXIT\r
-       C\r
-             IER=0\r
-       C */\r
+       \r
+       // NOW COME THE FACTORS DEPENDENT ON TAU1 AND TAU2.\r
+       \r
+       LOWER=TAU1[0];\r
+       UPPER=TAU2[0];\r
+       \r
+       while (true) {\r
+       \r
+           HTOL=0.5*TOL[0];\r
+           RR=(TAU2[0]-TAU1[0])*0.5;\r
+           MEAN=(TAU1[0]+TAU2[0])*0.5;\r
+           BB=(1.0+MEAN)/RR;\r
+           BB1=BB+Math.sqrt(BB*BB-1.0);\r
+           FFH=Math.pow((RR*(BB1-1.0/BB1)),(BETA+1.0))/Math.pow(BB1,TUK);\r
+           MAXRM[0]=0.0;\r
+           for (I=1; I <= NZZ; I++) {\r
+               XI[0]=(ZZ[I-1][0]-MEAN)/RR;\r
+               XI[1] = ZZ[I-1][1]/RR;\r
+               zmlt(XI[0],XI[1],XI[0],XI[1],cr,ci);\r
+               r = zabs(cr[0]-1.0,ci[0]);\r
+               theta = Math.atan2(ci[0], cr[0]-1.0);\r
+               mag = Math.sqrt(r);\r
+               ang = 0.5 * theta;\r
+               Z1[0] = mag * Math.cos(ang);\r
+               Z1[1] = mag * Math.sin(ang);\r
+               XI1[0]=XI[0]+Z1[0];\r
+               XI1[1]=XI[1]+Z1[1];\r
+               if (zabs(XI1[0],XI1[1]) < 1.0) {\r
+                 XI1[0]=XI[0]-Z1[0];\r
+                 XI1[1]=XI[1]-Z1[1];\r
+               }\r
+               r = zabs(XI1[0],XI1[1]);\r
+               theta = Math.atan2(XI1[1],XI1[0]);\r
+               mag = Math.pow(r, -TUK-1.0);\r
+               ang = theta * (-TUK - 1.0);\r
+               realVar = mag * Math.cos(ang);\r
+               imagVar = mag * Math.sin(ang);\r
+               zmlt(XI1[0],XI1[1],XI1[0],XI1[1],cr,ci);\r
+               zmlt(realVar,imagVar,(cr[0]-1.0)*RR,ci[0]*RR,cr2,ci2);\r
+               FFG[0] = cr2[0];\r
+               FFG[1] = ci2[0];\r
+               zmlt(GG[I-1][0],GG[I-1][1],FFG[0],FFG[1],cr,ci);\r
+               REMND[0] = cr[0] + HH[I-1][0]*FFH;\r
+               REMND[1] = ci[0] + HH[I-1][1]*FFH;\r
+               TERM=zabs(REMND[0],REMND[1]);\r
+               MAXRM[0]=Math.max(MAXRM[0],TERM);\r
+           } // for (I=1; I <= NZZ; I++)\r
+       \r
+           if (MAXRM[0] < TOL[0]) {\r
+       \r
+               // ACCURACY IS ACHIEVED, BUT MAYBE TAU2 COULD BE INCREASED OR\r
+               // TAU1 DECREASED\r
+       \r
+               if (MAXRM[0] < HTOL) {\r
+       \r
+                   // TAU2 NEEDS INCREASING IF T1FXD (BUT THIS IS ONLY POSSIBLE IF \r
+                   // TAU2<TAU2I) OR TAU1 NEED DECREASING OTHERWISE (BUT THIS IS \r
+                   // ONLY POSSIBLE IF TAU1>TAU1I)\r
+       \r
+                   if (T1FXD && TAU2[0] < TAU2I) {\r
+                       LOWER=TAU2[0];\r
+                       TAU2[0]=0.5*(LOWER+UPPER);\r
+                       continue;\r
+                   }\r
+                   else if (!T1FXD && TAU1[0] > TAU1I) {\r
+                       UPPER=TAU1[0];\r
+                       TAU1[0]=0.5*(LOWER+UPPER);\r
+                       continue;\r
+                   }\r
+               } // if (MARXM[0] < HTOL)\r
+               break;\r
+           } // if (MAXRM[0] < TOL[0])\r
+           else {\r
+       \r
+               // ACCURACY NOT ACHIEVED AND TAU2 NEEDS DECREASING OR TAU1 NEEDS\r
+               // INCREASING.\r
+       \r
+               if (FIRST) {\r
+                   TOL[0]=HTOL;\r
+                   FIRST=false;\r
+               }\r
+               if (T1FXD) {\r
+                   UPPER=TAU2[0];\r
+                   TAU2[0]=0.5*(LOWER+UPPER);\r
+               }\r
+               else {\r
+                 LOWER=TAU1[0];\r
+                 TAU1[0]=0.5*(LOWER+UPPER);\r
+               }\r
+               continue;\r
+           } // else\r
+       } // while (true)\r
+       \r
+       // NORMAL EXIT\r
+       \r
+        IER[0]=0;\r
     } // private void DEPPL9\r
 \r
+    private double JACSUM(double X, int N, double A[], double B[], double H, double CO[]) {\r
+        // INTEGER N\r
+        // REAL X,A(*),B(*),H,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,CURR,NEXT;\r
+        int K;\r
+        double result;\r
+\r
+        if (N == 0) {\r
+            result=CO[0]/Math.sqrt(H);\r
+        }\r
+        else if (N > 0) {\r
+            PREV=CO[N];\r
+            CURR=CO[N-1]+(X-B[N-1])*PREV/A[N-1];\r
+            for (K=N-2; K >= 0; K--) {\r
+                NEXT=CO[K]+(X-B[K])*CURR/A[K]-A[K]*PREV/A[K+1];\r
+                PREV=CURR;\r
+                CURR=NEXT;\r
+            } // for (K=N-2; K >= 0; K--)\r
+            result=CURR/Math.sqrt(H);\r
+        }\r
+        else {\r
+            result=0.0;\r
+        }\r
+        return result;\r
+    } // private double JACSUM\r
+\r
+   private void DMPHYC(int NPTS, double PHYPT[][],double CANPT[][], boolean WANTM, int IER[]) {\r
+       \r
+       // INTEGER NPTS,IER\r
+       // INTEGER IGEOM(*),ISNPH(*),IQUPH(*)\r
+       // REAL RGEOM(*),RSNPH(*),RQUPH(*)\r
+       // COMPLEX CENTR\r
+       // COMPLEX PHYPT(*),CANPT(*),ZQUPH(*)\r
+       // LOGICAL INTER,WANTM\r
+       \r
+        // ......................................................................\r
+       \r
+       // 1.     DMPHYC\r
+       //           DOMAIN MAPPING FOR THE PHYSICAL --> CANONICAL MAP.\r
+       \r
+       // 2.     PURPOSE\r
+       //           GIVEN A VECTOR OF ARBITRARY POINTS IN THE PHYSICAL DOMAIN, \r
+       //           THIS ROUTINE COMPUTES THE VECTOR OF APPROXIMATE IMAGE POINTS\r
+       //           IN THE CANONICAL DOMAIN.\r
+                  \r
+       \r
+       // 3.     CALLING SEQUENCE\r
+       //           CALL DMPHYC(NPTS,PHYPT,CANPT,INTER,CENTR,IGEOM,RGEOM,ISNPH,\r
+       //                       RSNPH,IQUPH,RQUPH,ZQUPH,WANTM,IER)\r
+       \r
+       //        PARAMETERS\r
+       //         ON ENTRY\r
+       //            NPTS   - INTEGER\r
+       //                     THE NUMBER OF POINTS TO BE MAPPED.\r
+       \r
+       //            PHYPT  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST NPTS.  THIS IS\r
+       //                     THE VECTOR OF GIVEN POINTS IN THE PHYSICAL 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
+       //            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
+       //            IQUPH  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR IQUPH PREVIOUSLY SET UP BY\r
+       //                     GQPHYC.\r
+       \r
+       //            RQUPH  - REAL ARRAY\r
+       //                     THE REAL ARRAY PREVIOUSLY SET UP BY GQPHYC.\r
+       \r
+       //            ZQUPH  - COMPLEX ARRAY\r
+       //                     THE COMPLEX ARRAY PREVIOUSLY SET UP BY GQPHYC.\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
+       \r
+       //         ON EXIT\r
+       //            CANPT  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST NPTS.  CANPT(K)\r
+       //                     IS THE COMPUTED IMAGE IN THE CANONICAL DOMAIN OF\r
+       //                     THE GIVEN PHYSICAL POINT PHYPT(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 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
+       //             ROUTINES JAPHYC AND GQPHYC HAVE SUCCESSFULLY EXECUTED, \r
+       //             AND THAT MANY INPUT ARGUMENTS FOR DMPHYC ARE OUTPUT VALUES\r
+       //             FROM JAPHYC AND GQPHYC.\r
+       //           - THIS ROUTINE MAY BE USED FOR MAPPING POINTS ON THE BOUN-\r
+       //             DARY OF THE PHYSICAL DOMAIN, BUT THE ROUTINE BMPHYC WILL\r
+       //             BE SOMEWHAT MORE EFFICIENT FOR THIS CASE.\r
+       \r
+       // ......................................................................\r
+       //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+       //     LAST UPDATE: 6 JULY 1990\r
+       // ......................................................................\r
+            \r
+       //     LOCAL VARAIBLES\r
+       \r
+       int MNSUA,TNSUA;\r
+       // CHARACTER*6 IERTXT\r
+       \r
+       // EXTERNAL PHTCA1,IERTXT\r
+       \r
+       NARCS=ISNPH[0];\r
+       NQPTS=ISNPH[1];\r
+       TNSUA=ISNPH[2];\r
+       MNSUA=ISNPH[4];\r
+       MNEQN=ISNPH[5];\r
+       MQUPH=IQUPH[3];\r
+       \r
+       NJIND=NARCS+1;\r
+       TNGQP=NQPTS*NJIND;\r
+       \r
+       // SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC\r
+       \r
+       // SET UP THE POINTERS TO  ISNPH AND RSNPH, AS IN JAPHYC\r
+       \r
+       // SET UP POINTERS TO IQUPH AND RQUPH, AS IN GQPHYC\r
+       \r
+       // GET REQUIRED CANONICAL POINTS\r
+       \r
+       PHTCA1(NPTS,PHYPT,CANPT,NARCS,NQPTS,TNSUA,DGPOL,\r
+            JATYP,LOSUB,LQSBF,NPPQF,PARNT,\r
+            AICOF,ACOEF,BICOF,BCFSN,BCOEF,\r
+            H0VAL,HIVAL,HALEN,JACIN,RGEOM[1],\r
+            MIDPT,QUPTS,QUWTS,SOLUN,TPPQF,\r
+            TRRAD,VTARG,WPPQF,CENTR,FACTR,\r
+            ZPPQF,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 DMPHYC\r
+   \r
+   private void PHTCA1(int NPTS, double PHYPT[][], double CANPT[][], int NARCS, int NQPTS,\r
+       int TNSUA, int DGPOL[], int JATYP[], int LOSUB[], int LPQSB[], int NPPQF[],\r
+          int PARNT[], double A1COF[], double ACOEF[], double B1COF[], double BCFSN[],\r
+          double BCOEF[], double H0VAL[], double H1VAL[], double HALEN[], double JACIN[],\r
+          double LGTOL, double MIDPT[], double QUPTS[], double QUWTS[], double SOLUN[],\r
+          double TPPQF[], double TRRAD[], double VTARG[], double WPPQF[], double CENTR[],\r
+          double FACTR[], double ZPPQF[][], boolean INTER,int IER[]) {\r
+          // INTEGER IER,NPTS,NARCS,NQPTS,TNSUA\r
+       // INTEGER DGPOL(*),JATYP(*),LOSUB(*),LPQSB(*),NPPQF(*),PARNT(*)\r
+          // REAL LGTOL\r
+          // REAL A1COF(*),ACOEF(*),B1COF(*),BCFSN(*),BCOEF(*),H1VAL(*),\r
+          // +H0VAL(*),HALEN(*),JACIN(*),MIDPT(*),QUPTS(*),QUWTS(*),SOLUN(*),\r
+          // +TPPQF(*),TRRAD(*),VTARG(*),WPPQF(*)\r
+          //COMPLEX CENTR,FACTR\r
+          // COMPLEX CANPT(*),PHYPT(*),ZPPQF(*)\r
+          // LOGICAL INTER\r
+               \r
+               //     GIVEN THE ARRAY PHYPT OF NPTS POINTS IN THE PHYSICAL PLANE, THIS\r
+               //     ROUTINE COMPUTES THE ARRAY CANPT OF IMAGES IN THE CANONICAL\r
+               //     PLANE.\r
+               \r
+               //     IER=0     - NORMAL EXIT\r
+               //     IER=27    - LOCAL PARAMETER MXNQD NEEDS INCREASING\r
+               //     IER=28    - LOCAL PARAMETER MNCOF NEEDS INCREASING\r
+               \r
+               //     LOCAL VARIABLES\r
+               \r
+          final int MNCOF = 32;\r
+          final int MQIN1 = 11;\r
+          final int MXNQD = 80;\r
+          int AJT,DEG,I,IA,IP,K,J,J1,J2,JQ,JT,LIM,LOD,LOL,LOM,\r
+                   NQ,NQUAD,PT,QINTS;\r
+           /*double AISUM,ANGLE,ARGBR,ARGIN1,ARSUM,BETA,CURARG,DIST,HL,ISUM,\r
+                    +JACSUM,LIMIT,MD,MEAN,MINDS,NEWTL,PI,PTHTL,R1MACH,RR,RRB,RSUM,RT1,\r
+                    +RT2,SCO,SS,STARG,STRT1,STTH1,SUM1,THET1,THET2,TOLOU,TT,TXI,TUPI,WT\r
+                     COMPLEX BCF,CJACSU,CT,DPARFN,PARFUN,PSI,XI,DIFF1,DIFF2,\r
+                    +STDF1,ZXI,ZTOB1,ZZ\r
+                     LOGICAL FIRST\r
+                     EXTERNAL ARGIN1,CJACSU,DPARFN,JACSUM,PARFUN,PPSBI1,R1MACH,ZTOB1\r
+                     PARAMETER (MNCOF=32,MQIN1=11,MXNQD=80,PTHTL=1E-3)\r
+                     PARAMETER (LIMIT=2.3562E+0)\r
+                     REAL JACOF(MNCOF),TSPEC(MXNQD),WSPEC(MXNQD),XENPT(MQIN1)\r
+                     COMPLEX JCOFC(MNCOF),ZSPEC(MXNQD)\r
+               C\r
+                     NEWTL=SQRT(R1MACH(4))\r
+                     PI=4E+0*ATAN(1E+0)\r
+                     TUPI=2E+0*PI\r
+                     LOL=NARCS*NQPTS\r
+                     DO 300 IP=1,NPTS\r
+                       ZZ=PHYPT(IP)\r
+                       RSUM=0E+0\r
+                       ISUM=0E+0\r
+                       FIRST=.TRUE.\r
+                       DO 200 IA=1,TNSUA\r
+                         PT=PARNT(IA)\r
+                         JT=JATYP(IA)\r
+                         NQ=NPPQF(IA)\r
+                         K=LPQSB(IA)-1\r
+                         HL=HALEN(IA)\r
+                         MD=MIDPT(IA)\r
+                         ARSUM=0E+0\r
+                         AISUM=0E+0\r
+                         DO 100 JQ=1,NQ\r
+                           K=K+1\r
+                           DIFF2=ZZ-ZPPQF(K)\r
+                           RT2=MD+HL*TPPQF(K)\r
+                           DIST=ABS(DIFF2)\r
+                           IF (DIST .GE. TRRAD(K)) THEN\r
+                               WT=WPPQF(K)\r
+                               IF (WT .NE. 0E+0) THEN\r
+                                   ARSUM=ARSUM+WT*LOG(DIST)\r
+                                   IF (FIRST) THEN\r
+                                       CURARG=ATAN2(AIMAG(DIFF2),REAL(DIFF2))\r
+                                       THET2=CURARG\r
+                                       FIRST=.FALSE.\r
+                                       STARG=CURARG\r
+                                   ELSE\r
+               C                        CT=DIFF2/DIFF1\r
+               C                        CT=DIFF2*CONJG(DIFF1)\r
+               C                        ANGLE=ATAN2(AIMAG(CT),REAL(CT))\r
+                                       THET2=ATAN2(AIMAG(DIFF2),REAL(DIFF2))\r
+                                       ANGLE=THET2-THET1\r
+                                       IF (ANGLE .LE. -PI .OR. ANGLE .GT. PI) THEN\r
+                                         ANGLE=ANGLE-SIGN(TUPI,ANGLE)\r
+                                       ENDIF\r
+                                       IF (ABS(ANGLE) .GE. LIMIT) THEN\r
+                                           ANGLE=ARGIN1(RT1,RT2,PT,-DIFF1,-DIFF2,ZZ,\r
+                    +                                   LIMIT)\r
+                                       ENDIF\r
+                                       CURARG=CURARG+ANGLE\r
+                                   ENDIF\r
+                                   AISUM=CURARG*WT+AISUM\r
+                                   RT1=RT2\r
+                                   DIFF1=DIFF2\r
+                                   THET1=THET2\r
+                               ENDIF\r
+                           ELSE\r
+               C\r
+               C             ZZ IS TOO CLOSE TO ARC IA TO USE THE STANDARD RULE.\r
+               C             FIND THE QUADRATURE POINT NEAREST TO ZZ.\r
+               C\r
+                             J1=JQ\r
+                             MINDS=DIST\r
+                             TXI=TPPQF(K)\r
+                             ZXI=ZPPQF(K)\r
+               40            CONTINUE\r
+                             J1=J1+1\r
+                             IF (J1 .LE. NQ) THEN\r
+                               K=K+1\r
+                               DIFF2=ZZ-ZPPQF(K)\r
+                               DIST=ABS(DIFF2)\r
+                               IF (DIST .LT. MINDS) THEN\r
+                                 MINDS=DIST\r
+                                 TXI=TPPQF(K)\r
+                                 ZXI=ZPPQF(K)\r
+                                 GOTO 40\r
+                               ENDIF\r
+                             ENDIF\r
+               C\r
+               C             PRELIMINARIES\r
+               C\r
+                             IF (JT .GT. 0) THEN\r
+                               SS=1E+0\r
+                             ELSE\r
+                               SS=-1E+0\r
+                             ENDIF\r
+                             AJT=ABS(JT)\r
+                             BETA=JACIN(AJT)\r
+                             DEG=DGPOL(IA)\r
+                             IF (DEG+1 .GT. MNCOF) THEN\r
+                               IER=28\r
+                               RETURN\r
+                             ENDIF\r
+                             LOM=LOSUB(IA)\r
+                             LOD=(AJT-1)*NQPTS+1\r
+               C\r
+               C             NOW USE NEWTON'S METHOD TO ESTIMATE THE PARAMETRIC\r
+               C             PRE-IMAGE XI OF ZZ.\r
+               C\r
+                             XI=CMPLX(TXI)\r
+                             CT=MD+HL*XI\r
+                             DIFF2=(ZXI-ZZ)/(DPARFN(PT,CT)*HL)\r
+                             XI=XI-DIFF2\r
+               50            CONTINUE\r
+                             IF (ABS(DIFF2) .GT. NEWTL) THEN\r
+                               CT=MD+HL*XI\r
+                               DIFF2=(PARFUN(PT,CT)-ZZ)/(DPARFN(PT,CT)*HL)\r
+                               XI=XI-DIFF2\r
+                               GOTO 50\r
+                             ELSE\r
+               C\r
+               C               LAST ITERATION\r
+               C\r
+                               CT=MD+HL*XI\r
+                               DIFF2=(PARFUN(PT,CT)-ZZ)/(DPARFN(PT,CT)*HL)\r
+                               XI=XI-DIFF2\r
+                             ENDIF\r
+               C\r
+                             XI=SS*XI\r
+               C\r
+                             IF (ABS(AIMAG(XI)).LT.PTHTL .AND. ABS(REAL(XI)).LT.1E+\r
+                    +        0+PTHTL) THEN\r
+               C\r
+               C               ZZ IS PATHOLOGICALLY CLOSE TO ARC IA AND WE USE THE  \r
+               C               CONTINUATION OF THE BOUNDARY CORRESPONDENCE FUNCTION TO \r
+               C               ESTIMATE CANPT.\r
+               C\r
+                               PSI=CJACSU(XI,DEG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT),\r
+                    +                     BCFSN(LOM+1))\r
+                               PSI=ZTOB1(XI+1E+0,BETA+1E+0,JT,INTER)*\r
+                    +              (BCFSN(LOM)-(1E+0-XI)*PSI)\r
+                               IF (JT .GT. 0) THEN\r
+                                 BCF=VTARG(IA)\r
+                               ELSE\r
+                                 BCF=VTARG(IA+1)\r
+                               ENDIF\r
+                               BCF=BCF+SS*PSI\r
+                               CANPT(IP)=CEXP((0E+0,1E+0)*BCF)\r
+                               GOTO 300\r
+                             ELSE\r
+               C\r
+               C               SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS\r
+               C               PARTICULAR POINT ZZ.\r
+               C\r
+                               SCO=SS\r
+                               DO 55 J=1,DEG+1\r
+                                 J1=LOM+J-1\r
+                                 SCO=SCO*SS\r
+                                 JACOF(J)=SOLUN(J1)*SCO\r
+                                 JCOFC(J)=CMPLX(SOLUN(J1)*SCO)\r
+               55              CONTINUE\r
+                               CALL PPSBI1(XI,BETA,NQPTS,DEG,ACOEF(LOD),BCOEF(LOD),\r
+                    +                      H0VAL(AJT),JCOFC,LGTOL,TOLOU,XENPT,QINTS,\r
+                    +                      MQIN1,IER)\r
+                               IF (IER .GT. 0) THEN\r
+                                 RETURN\r
+                               ENDIF\r
+                               NQUAD=QINTS*NQPTS\r
+                               IF (NQUAD .GT. MXNQD) THEN\r
+                                 IER=27\r
+                                 RETURN\r
+                               ENDIF\r
+                               K=0\r
+                               SUM1=BETA+1E+0\r
+                               DO 70 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 60 J=1,NQPTS\r
+                                     J1=LOD+J-1\r
+                                     K=K+1\r
+                                     TT=(MEAN+RR*QUPTS(J1))\r
+                                     WSPEC(K)=RRB*QUWTS(J1)*JACSUM(TT,DEG,ACOEF(LOD),\r
+                    +                         BCOEF(LOD),H0VAL(AJT),JACOF)\r
+                                     TT=TT*SS\r
+                                     TSPEC(K)=MD+TT*HL\r
+                                     CT=CMPLX(TSPEC(K))\r
+                                     ZSPEC(K)=PARFUN(PT,CT)\r
+               60                  CONTINUE\r
+                                 ELSE\r
+                                   DO 65 J=1,NQPTS\r
+                                     J1=LOL+J\r
+                                     K=K+1\r
+                                     TT=(MEAN+RR*QUPTS(J1))\r
+                                     WSPEC(K)=RR*QUWTS(J1)*(1E+0+TT)**BETA*JACSUM(TT,\r
+                    +                         DEG,ACOEF(LOD),BCOEF(LOD),H0VAL(AJT),\r
+                    +                         JACOF)\r
+                                     TT=TT*SS\r
+                                     TSPEC(K)=MD+TT*HL\r
+                                     CT=CMPLX(TSPEC(K))\r
+                                     ZSPEC(K)=PARFUN(PT,CT)\r
+               65                  CONTINUE\r
+                                 ENDIF\r
+               70              CONTINUE\r
+                               IF (SS .LT. 0E+0) THEN\r
+                                 LIM=NQUAD\r
+                                 IF (MOD(LIM,2) .EQ. 0) THEN\r
+                                   LIM=LIM/2\r
+                                 ELSE\r
+                                   LIM=(LIM-1)/2\r
+                                 ENDIF\r
+                                 J1=0\r
+                                 J2=NQUAD+1\r
+                                 DO 75 J=1,LIM\r
+                                   J1=J1+1\r
+                                   J2=J2-1\r
+                                   TT=WSPEC(J1)\r
+                                   WSPEC(J1)=WSPEC(J2)\r
+                                   WSPEC(J2)=TT\r
+                                   TT=TSPEC(J1)\r
+                                   TSPEC(J1)=TSPEC(J2)\r
+                                   TSPEC(J2)=TT\r
+                                   CT=ZSPEC(J1)\r
+                                   ZSPEC(J1)=ZSPEC(J2)\r
+                                   ZSPEC(J2)=CT\r
+               75                CONTINUE\r
+                               ENDIF\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
+                               ARSUM=0E+0\r
+                               AISUM=0E+0\r
+                               IF (IA .EQ. 1) THEN\r
+                                 FIRST=.TRUE.\r
+                               ELSE\r
+                                 CURARG=STARG\r
+                                 RT1=STRT1\r
+                                 DIFF1=STDF1\r
+                                 THET1=STTH1\r
+                               ENDIF\r
+                               DO 80 K=1,NQUAD\r
+                                   WT=WSPEC(K)\r
+                                   DIFF2=ZZ-ZSPEC(K)\r
+                                   RT2=TSPEC(K)\r
+                                   DIST=ABS(DIFF2)\r
+                                   ARSUM=ARSUM+WT*LOG(DIST)\r
+                                   IF (FIRST) THEN\r
+                                     CURARG=ATAN2(AIMAG(DIFF2),REAL(DIFF2))\r
+                                     THET2=CURARG\r
+                                     FIRST=.FALSE.\r
+                                   ELSE\r
+               C                      CT=DIFF2/DIFF1\r
+               C                      CT=DIFF2*CONJG(DIFF1)\r
+               C                      ANGLE=ATAN2(AIMAG(CT),REAL(CT))\r
+                                     THET2=ATAN2(AIMAG(DIFF2),REAL(DIFF2))\r
+                                     ANGLE=THET2-THET1\r
+                                     IF (ANGLE .LE. -PI .OR. ANGLE .GT. PI) THEN\r
+                                         ANGLE=ANGLE-SIGN(TUPI,ANGLE)\r
+                                     ENDIF\r
+                                     IF (ABS(ANGLE) .GE. LIMIT) THEN\r
+                                         ANGLE=ARGIN1(RT1,RT2,PT,-DIFF1,-DIFF2,ZZ,\r
+                    +                                 LIMIT)\r
+                                     ENDIF\r
+                                     CURARG=CURARG+ANGLE\r
+                                   ENDIF\r
+                                   AISUM=CURARG*WT+AISUM\r
+                                   RT1=RT2\r
+                                   DIFF1=DIFF2\r
+                                   THET1=THET2\r
+               80              CONTINUE\r
+                               GOTO 150\r
+                             ENDIF              \r
+                           ENDIF\r
+               C\r
+               C         END OF QUADRATURE SUM LOOP \r
+               C\r
+               100       CONTINUE\r
+               C\r
+               150       CONTINUE\r
+                         RSUM=RSUM+ARSUM\r
+                         ISUM=ISUM+AISUM\r
+                         IF (JT .LT. 0) THEN\r
+               C\r
+               C           BRING THE ARGUMENT FORWARD TO THE CORNER POINT AND REPLACE\r
+               C           THE INCREMENTED CURARG VALUE BY AN INVERSE TANGENT\r
+               C           EVALUATION\r
+               C\r
+                           DIFF2=ZZ-PARFUN(PT,(1E+0,0E+0))\r
+                           RT2=1E+0\r
+                           THET2=ATAN2(AIMAG(DIFF2),REAL(DIFF2))\r
+                           ANGLE=THET2-THET1\r
+                           IF (ANGLE .LE. -PI .OR. ANGLE .GT. PI) THEN\r
+                               ANGLE=ANGLE-SIGN(TUPI,ANGLE)\r
+                           ENDIF\r
+                           IF (ABS(ANGLE) .GE. LIMIT) THEN\r
+                               ANGLE=ARGIN1(RT1,RT2,PT,-DIFF1,-DIFF2,ZZ,\r
+                    +                       LIMIT)\r
+                           ENDIF\r
+                           CURARG=CURARG+ANGLE\r
+                           ARGBR=ANINT((CURARG-THET2)/TUPI)\r
+                           CURARG=THET2+TUPI*ARGBR\r
+                           RT1=-1E+0\r
+                           DIFF1=DIFF2\r
+                           THET1=THET2\r
+                         ENDIF           \r
+                         STARG=CURARG\r
+                         STRT1=RT1\r
+                         STDF1=DIFF1\r
+                         STTH1=THET1\r
+               C\r
+               C       END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA\r
+               C\r
+               200     CONTINUE\r
+                       CT=CMPLX(RSUM,ISUM)\r
+                       CT=CEXP(CT)\r
+                       IF (INTER) THEN\r
+                         CANPT(IP)=(ZZ-CENTR)*FACTR/CT\r
+                       ELSE\r
+                         CANPT(IP)=CT*FACTR\r
+                       ENDIF\r
+               C\r
+               C     END OF MAP CALCULATION FOR FIELD POINT NUMBER IP\r
+               C\r
+               300   CONTINUE */\r
+               \r
+               IER[0]=0;\r
+        \r
+    } // private void PHTCA1\r
 \r
 \r
       /**\r