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

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

index 4c0113d..f554a8a 100644 (file)
@@ -9615,211 +9615,322 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           final int MNCOF = 32;\r
           final int MQIN1 = 11;\r
           final int MXNQD = 80;\r
+          int QINTS[] = new int[1];\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
+                   NQ,NQUAD,PT;\r
+          final double PTHTL = 1.0E-3;\r
+          final double LIMIT = 2.3562;\r
+          double TOLOU[] = new double[1];\r
+          double AISUM,ANGLE,ARGBR,ARSUM,BETA,CURARG,DIST,HL,ISUM,\r
+                    MD,MEAN,MINDS,NEWTL,PI,RR,RRB,RSUM,RT1,\r
+                    RT2,SCO,SS,STARG,STRT1,STTH1,SUM1,THET1,THET2,TT,TXI,TUPI,WT;\r
+          double BCF[] = new double[2];\r
+          double CT[] = new double[2];\r
+          double PSI[] = new double[2];\r
+          double XI[] = new double[2];\r
+          double DIFF1[] = new double[2];\r
+          double DIFF2[] = new double[2];\r
+          double STDF1[] = new double[2];\r
+          double ZXI[] = new double[2];\r
+          double ZZ[] = new double[2];\r
+       // COMPLEX BCF,CT,PSI,XI,DIFF1,DIFF2,STDF1,ZXI,ZZ\r
+          boolean FIRST;\r
+          // EXTERNAL ARGIN1,CJACSU,DPARFN,JACSUM,PARFUN,PPSBI1,ZTOB1\r
+          double JACOF[] = new double[MNCOF];\r
+          double TSPEC[] = new double[MXNQD];\r
+          double WSPEC[] = new double[MXNQD];\r
+          double XENPT[] = new double[MQIN1];\r
+          double JCOFC[][] = new double[MNCOF][2];\r
+          double ZSPEC[][] = new double[MXNQD][2];\r
+          // COMPLEX JCOFC(MNCOF),ZSPEC(MXNQD)\r
+          double mag;\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
+          double DIFF1IN[] = new double[2];\r
+          double DIFF2IN[] = new double[2];\r
+          double POUT[];\r
+          double DOUT[];\r
+          double A[];\r
+          double B[];\r
+          double CO[];\r
+          double X[] = new double[2];\r
+          double ZOUT[];\r
+          double AIN[];\r
+          double BIN[];\r
+          int I2;\r
+               \r
+          NEWTL=Math.sqrt(EPS);\r
+          PI=Math.PI;\r
+          TUPI=2.0*PI;\r
+          LOL=NARCS*NQPTS;\r
+          /*iploop: for (IP=1; IP <= NPTS; IP++) {\r
+              ZZ[0]=PHYPT[IP-1][0];\r
+              ZZ[1]=PHYPT[IP-1][1];\r
+                  RSUM=0.0;\r
+                  ISUM=0.0;\r
+                  FIRST= true;\r
+                  for (IA=1; IA <= TNSUA; IA++) {\r
+                      PT=PARNT[IA-1];\r
+                      JT=JATYP[IA-1];\r
+                      NQ=NPPQF[IA-1];\r
+                      K=LPQSB[IA-1]-1;\r
+                      HL=HALEN[IA-1];\r
+                      MD=MIDPT[IA-1];\r
+                      ARSUM=0.0;\r
+                      AISUM=0.0;\r
+                      for (JQ=1; JQ <= NQ; JQ++) {\r
+                          K=K+1;\r
+                          DIFF2[0]=ZZ[0]-ZPPQF[K-1][0];\r
+                          DIFF2[1]=ZZ[1]-ZPPQF[K-1][1];\r
+                          RT2=MD+HL*TPPQF[K-1];\r
+                          DIST=zabs(DIFF2[0],DIFF2[1]);\r
+                          if (DIST >= TRRAD[K-1]) {\r
+                              WT=WPPQF[K-1];\r
+                              if (WT != 0.0) {\r
+                                  ARSUM=ARSUM+WT*Math.log(DIST);\r
+                                  if (FIRST) {\r
+                                      CURARG=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                                      THET2=CURARG;\r
+                                      FIRST=false;\r
+                                      STARG=CURARG;\r
+                                  } // if (FIRST)\r
+                                  else { // !FIRST\r
+                                      // CT=DIFF2/DIFF1\r
+                                      // CT=DIFF2*CONJG(DIFF1)\r
+                                      // ANGLE=ATAN2(AIMAG(CT),REAL(CT))\r
+                                      THET2=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                                      ANGLE=THET2-THET1;\r
+                                      if (ANGLE <= -PI || ANGLE > PI) {\r
+                                         if (ANGLE > PI) {\r
+                                                 ANGLE = ANGLE - TUPI;\r
+                                         }\r
+                                         else {\r
+                                                 ANGLE = ANGLE + TUPI;\r
+                                         }\r
+                                      }\r
+                                      if (Math.abs(ANGLE) >= LIMIT) {\r
+                                          DIFF1IN[0] = -DIFF1[0];\r
+                                          DIFF1IN[1] = -DIFF1[1];\r
+                                          DIFF2IN[0] = -DIFF2[0];\r
+                                          DIFF2IN[1] = -DIFF2[1];\r
+                                          ANGLE=ARGIN1(RT1,RT2,PT,DIFF1IN,DIFF2IN,ZZ,LIMIT);\r
+                                      }\r
+                                      CURARG=CURARG+ANGLE;\r
+                                  } // else !FIRST\r
+                                  AISUM=CURARG*WT+AISUM;\r
+                                  RT1=RT2;\r
+                                  DIFF1[0]=DIFF2[0];\r
+                                  DIFF1[1]=DIFF2[1];\r
+                                  THET1=THET2;\r
+                          } // if (WT != 0.0)\r
+                      } // if (DIST >= TRADD[K-1])\r
+                      else { // DIST < TRADD[K-1]\r
+               \r
+                          // ZZ IS TOO CLOSE TO ARC IA TO USE THE STANDARD RULE.\r
+                          // FIND THE QUADRATURE POINT NEAREST TO ZZ.\r
+               \r
+                          J1=JQ;\r
+                          MINDS=DIST;\r
+                          TXI=TPPQF[K-1];\r
+                          ZXI[0]=ZPPQF[K-1][0];\r
+                          ZXI[1]=ZPPQF[K-1][1];\r
+                          while (true) {\r
+                              J1=J1+1;\r
+                              if (J1 <= NQ) {\r
+                                  K=K+1;\r
+                                  DIFF2[0]=ZZ[0]-ZPPQF[K-1][0];\r
+                                  DIFF2[1]=ZZ[1]-ZPPQF[K-1][1];\r
+                                  DIST=zabs(DIFF2[0],DIFF2[1]);\r
+                                  if (DIST < MINDS) {\r
+                                      MINDS=DIST;\r
+                                      TXI=TPPQF[K-1];\r
+                                      ZXI[0]=ZPPQF[K-1][0];\r
+                                      ZXI[1]=ZPPQF[K-1][1];\r
+                                      continue;\r
+                                  } // if (DIST < MINDS)\r
+                              } //  if (J1 <= NQ)\r
+                              break;\r
+                          } // while (true)\r
+               \r
+                          // PRELIMINARIES\r
+       \r
+                          if (JT > 0) {\r
+                              SS=1.0;\r
+                          }\r
+                          else {\r
+                              SS=-1.0;\r
+                          }\r
+                          AJT=Math.abs(JT);\r
+                          BETA=JACIN[AJT-1];\r
+                          DEG=DGPOL[IA-1];\r
+                          if (DEG+1 > MNCOF) {\r
+                              IER[0]=28;\r
+                              return;\r
+                          }\r
+                          LOM=LOSUB[IA-1];\r
+                          LOD=(AJT-1)*NQPTS+1;\r
+               \r
+                          // NOW USE NEWTON'S METHOD TO ESTIMATE THE PARAMETRIC\r
+                      // PRE-IMAGE XI OF ZZ.\r
+               \r
+                          XI[0]=TXI;\r
+                          XI[1] = 0.0;\r
+                          CT[0]=MD+HL*XI[0];\r
+                          CT[1]=HL*XI[1];\r
+                          DOUT = DPARFN(PT,CT);\r
+                          zdiv(ZXI[0]-ZZ[0],ZXI[1]-ZZ[1],DOUT[0]*HL,DOUT[1]*HL,cr,ci);\r
+                          DIFF2[0] = cr[0];\r
+                          DIFF2[1] = ci[0];\r
+                          XI[0]=XI[0]-DIFF2[0];\r
+                          XI[1]=XI[1]-DIFF2[1];\r
+                          while (true) {\r
+                              if (zabs(DIFF2[0],DIFF2[1]) > NEWTL) {\r
+                                  CT[0]=MD+HL*XI[0];\r
+                                  CT[1] = HL*XI[1];\r
+                                  POUT = PARFUN(PT,CT);\r
+                                  DOUT = DPARFN(PT,CT);\r
+                                  zdiv(POUT[0]-ZZ[0],POUT[1]-ZZ[1],DOUT[0]*HL,DOUT[1]*HL,cr,ci);\r
+                                  DIFF2[0] = cr[0];\r
+                                  DIFF2[1] = ci[0];\r
+                                  XI[0]=XI[0]-DIFF2[0];\r
+                                  XI[1]=XI[1]-DIFF2[1];\r
+                                  continue;\r
+                              } // if (zabs(DIFF2[0],DIFF2[1]) > NEWTL)\r
+                              else {\r
+               \r
+                                  // LAST ITERATION\r
+               \r
+                                  CT[0]=MD+HL*XI[0];\r
+                                  CT[1]=HL*XI[1];\r
+                                  POUT = PARFUN(PT,CT);\r
+                                  DOUT = DPARFN(PT,CT);\r
+                                  zdiv(POUT[0]-ZZ[0],POUT[1]-ZZ[1],DOUT[0]*HL,DOUT[1]*HL,cr,ci);\r
+                                  DIFF2[0] = cr[0];\r
+                                  DIFF2[1] = ci[0];\r
+                                  XI[0]=XI[0]-DIFF2[0];\r
+                                  XI[1]=XI[1]-DIFF2[1];\r
+                              } // else\r
+                             break;\r
+                          } // while (true)\r
+               \r
+                          XI[0]=SS*XI[0];\r
+                          XI[1]=SS*XI[1];\r
+               \r
+                          if (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL) {\r
+               \r
+                              // ZZ IS PATHOLOGICALLY CLOSE TO ARC IA AND WE USE THE  \r
+                              // CONTINUATION OF THE BOUNDARY CORRESPONDENCE FUNCTION TO \r
+                              // ESTIMATE CANPT.\r
+               \r
+                              A = new double[DEG-1];\r
+                              B = new double[DEG-1];\r
+                              CO = new double[DEG];\r
+                              for (I = 0; I < DEG-1; I++) {\r
+                                  A[I] = A1COF[LOD+I-1];\r
+                                  B[I] = B1COF[LOD+I-1];\r
+                              }\r
+                              for (I = 0; I < DEG; I++) {\r
+                                  CO[I] = BCFSN[LOM+I];\r
+                              }\r
+                                  PSI=CJACSU(XI,DEG-1,A,B,H1VAL[AJT-1],CO);\r
+                                  X[0] = XI[0] + 1.0;\r
+                                  X[1] = XI[1];\r
+                                  ZOUT = ZTOB1(X,BETA+1.0,JT,INTER);\r
+                                  zmlt(1.0-XI[0],-XI[1],PSI[0],PSI[1],cr,ci);\r
+                                  zmlt(ZOUT[0],ZOUT[1],BCFSN[LOM-1]-cr[0],-ci[0],cr2,ci2);\r
+                                  PSI[0] = cr2[0];\r
+                                  PSI[1] = ci2[0];\r
+                              if (JT > 0) {\r
+                                  BCF[0]=VTARG[IA-1];\r
+                                  BCF[1]=0.0;\r
+                              }\r
+                              else {\r
+                                  BCF[0]=VTARG[IA];\r
+                              }\r
+                              BCF[0]=BCF[0]+SS*PSI[0];\r
+                              BCF[1]=BCF[1]+SS*PSI[1];\r
+                              mag = Math.exp(-BCF[1]);\r
+                              CANPT[IP-1][0] = mag * Math.cos(BCF[0]);\r
+                              CANPT[IP-1][1] = mag * Math.sin(BCF[0]);\r
+                              continue iploop;\r
+                          } // if (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL)\r
+                          else { // ! (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL)\r
+               \r
+                              // SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS\r
+                              // PARTICULAR POINT ZZ.\r
+               \r
+                              SCO=SS;\r
+                              for (J=1; J <= DEG+1; J++) {\r
+                                  J1=LOM+J-1;\r
+                                  SCO=SCO*SS;\r
+                                  JACOF[J-1]=SOLUN[J1-1]*SCO;\r
+                                  JCOFC[J-1][0]=SOLUN[J1-1]*SCO;\r
+                                  JCOFC[J-1][1]= 0.0;\r
+                              } // for (J=1; J <= DEG+1; J++)\r
+                              AIN = new double[DEG];\r
+                              BIN = new double[DEG];\r
+                              for (I = 0; I < DEG; I++) {\r
+                                  AIN[I] = ACOEF[LOD+I-1];\r
+                                  BIN[I] = BCOEF[LOD+I-1];\r
+                              }\r
+                              PPSBI1(XI,BETA,NQPTS,DEG,AIN,BIN,\r
+                                         H0VAL[AJT-1],JCOFC,LGTOL,TOLOU,XENPT,QINTS,MQIN1,IER);\r
+                              if (IER[0] > 0) {\r
+                                  return;\r
+                              }\r
+                              NQUAD=QINTS[0]*NQPTS;\r
+                              if (NQUAD > MXNQD) {\r
+                                  IER[0]=27;\r
+                                  return;\r
+                              }\r
+                              K=0;\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
+                                           AIN = new double[DEG];\r
+                                           BIN = new double[DEG];\r
+                                           for (I2 = 0; I2 < DEG; I2++) {\r
+                                               AIN[I2] = ACOEF[LOD+I2-1];\r
+                                               BIN[I2] = BCOEF[LOD+I2-1];\r
+                                           }\r
+                                           WSPEC[K-1]=RRB*QUWTS[J1-1]*JACSUM(TT,DEG,AIN,BIN,H0VAL[AJT-1],JACOF);\r
+                                           TT=TT*SS;\r
+                                           TSPEC[K-1]=MD+TT*HL;\r
+                                           CT[0]=TSPEC[K-1];\r
+                                           CT[1] = 0.0;\r
+                                           ZSPEC[K-1]=PARFUN(PT,CT);\r
+                                       } // for (J=1; J <= NQPTS; J++)\r
+                                   } // if (I == 1)\r
+                                   else { // I != 1\r
+                                       for (J=1; J <= NQPTS; J++) {\r
+                                           J1=LOL+J;\r
+                                            K=K+1;\r
+                                            TT=(MEAN+RR*QUPTS[J1-1]);\r
+                                            AIN = new double[DEG];\r
+                                            BIN = new double[DEG];\r
+                                            for (I2 = 0; I2 < DEG; I2++) {\r
+                                                AIN[I2] = ACOEF[LOD+I2-1];\r
+                                                BIN[I2] = BCOEF[LOD+I2-1];\r
+                                            }\r
+                                            WSPEC[K-1]=RR*QUWTS[J1-1]*Math.pow((1.0+TT),BETA)*JACSUM(TT,\r
+                                             DEG,AIN,BIN,H0VAL[AJT-1],JACOF);\r
+                                            TT=TT*SS;\r
+                                            TSPEC[K-1]=MD+TT*HL;\r
+                                            CT[0]=TSPEC[K-1];\r
+                                            CT[1] = 0.0;\r
+                                            ZSPEC[K-1]=PARFUN(PT,CT);\r
+                                       } // for (J=1; J <= NQPTS; J++)\r
+                                   } // else I != 1\r
+                              } // for (I=1; I <= QINTS[0]; I++)\r
                                IF (SS .LT. 0E+0) THEN\r
                                  LIM=NQUAD\r
                                  IF (MOD(LIM,2) .EQ. 0) THEN\r
@@ -9888,12 +9999,12 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                                    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
+                          } // else ! (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL)              \r
+                        } //  else DIST < TRADD[K-1]\r
+               \r
+                        // END OF QUADRATURE SUM LOOP \r
+               \r
+                      } // for (JQ=1; JQ <= NQ; JQ++)\r
                C\r
                150       CONTINUE\r
                          RSUM=RSUM+ARSUM\r
@@ -9919,32 +10030,329 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                            ARGBR=ANINT((CURARG-THET2)/TUPI)\r
                            CURARG=THET2+TUPI*ARGBR\r
                            RT1=-1E+0\r
-                           DIFF1=DIFF2\r
+                           DIFF1[0]=DIFF2[0];\r
+                           DIFF1[1]=DIFF2[1];\r
                            THET1=THET2\r
                          ENDIF           \r
                          STARG=CURARG\r
                          STRT1=RT1\r
-                         STDF1=DIFF1\r
+                         STDF1[0]=DIFF1[0];\r
+                         STDF1[1]=DIFF1[1];\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
+                      // END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA\r
+               \r
+                  } // for (IA=1; IA <= TNSUA; IA++)\r
+                  CT[0] = RSUM;\r
+                  CT[1] = ISUM;\r
+                  mag = Math.exp(CT[0]);\r
+                  CT[0] = mag * Math.cos(CT[1]);\r
+                  CT[1] = mag * Math.sin(CT[1]);\r
+                  if (INTER) {\r
+                          zmlt(ZZ[0]-CENTR[0],ZZ[1]-CENTR[1],FACTR[0],FACTR[1],cr,ci);\r
+                          zdiv(cr[0],ci[0],CT[0],CT[1],cr2,ci2);\r
+                          CANPT[IP-1][0] = cr2[0];\r
+                          CANPT[IP-1][1] = ci2[0];\r
+                  }\r
+                  else {\r
+                          zmlt(CT[0],CT[1],FACTR[0],FACTR[1],cr,ci);\r
+                          CANPT[IP-1][0] = cr[0];\r
+                          CANPT[IP-1][1] = ci[0];\r
+                  }\r
+               \r
+               // END OF MAP CALCULATION FOR FIELD POINT NUMBER IP\r
+               \r
+           } // iploop: for (IP=1; IP <= NPTS; IP++)*/\r
                \r
                IER[0]=0;\r
         \r
     } // private void PHTCA1\r
+   \r
+   \r
+   private double ARGIN1(double RT1, double RT2, int PT, double DIFF1[], double DIFF2[],\r
+                  double ZZ[], double LIMIT) {\r
+       // INTEGER PT\r
+       // REAL RT1,RT2,LIMIT\r
+       // COMPLEX DIFF1,DIFF2,ZZ\r
+\r
+       // ZZ IS A GIVEN FIELD POINT AND DIFF1, DIFF2 ARE THE DIFFERENCES\r
+       // BETWEEN ZZ AND CONSECUTIVE POINTS ON THE BOUNDARY \r
+       // (DIFF1=PARFUN(PT,RT1)-ZZ, ZET2=PARFUN(PT,RT2)-ZZ).  THE\r
+       // PURPOSE OF THIS ROUTINE IS TO CALCULATE THE INCREASE IN THE \r
+       // ARGUMENT ARG(ZZ-Z) AS Z MOVES ALONG THE BOUNDARY FROM THE POINT\r
+       // WITH PARAMETER VALUE RT1 TO THE POINT WITH PARAMETER VALUE RT2.\r
+\r
+       // LOCAL VARIABLES\r
+\r
+       int NANGS,NINTS;\r
+       double ANGLE,T1,T2;\r
+       double D1[] = new double[2];\r
+       double D2[] = new double[2];\r
+       double V[] = new double[2];\r
+       // COMPLEX D1,D2,PARFUN,V\r
+       // EXTERNAL PARFUN\r
+       double PIN[] = new double[2];\r
+       double POUT[];\r
+       double result;\r
+       double cr[] = new double[1];\r
+       double ci[] = new double[1];\r
+       // LIMIT IS CURRENTLY SET TO 3*PI/4, APPROXIMATELY\r
+\r
+       T1=RT1;\r
+       T2=(RT1+RT2)*0.5;\r
+       D1[0]=DIFF1[0];\r
+       D1[1]=DIFF1[1];\r
+       PIN[0] = T2;\r
+       PIN[1] = 0.0;\r
+       POUT = PARFUN(PT,PIN);\r
+       D2[0] = POUT[0] - ZZ[0];\r
+       D2[1] = POUT[1] - ZZ[1];\r
+       NANGS=0;\r
+       NINTS=2;\r
+       result=0.0;\r
+\r
+       while (true) {\r
+          zmlt(D2[0],D2[1],D1[0],-D1[1],cr,ci);\r
+          V[0] = cr[0];\r
+          V[1] = ci[0];\r
+           ANGLE=Math.atan2(V[1],V[0]);\r
+           if (Math.abs(ANGLE) >= LIMIT) {\r
+               T2=(T1+T2)*0.5;\r
+               PIN[0] = T2;\r
+               PIN[1] = 0.0;\r
+               POUT = PARFUN(PT,PIN);\r
+               D2[0] = POUT[0] - ZZ[0];\r
+               D2[1] = POUT[1] - ZZ[1];\r
+               NINTS=NINTS+1;\r
+               continue;\r
+           } // if (Math.abs(ANGLE) >= LIMIT)\r
+           else {\r
+               result=result+ANGLE;\r
+               NANGS=NANGS+1;\r
+               if (NANGS != NINTS) {\r
+                   T1=T2;\r
+                   T2=RT2;\r
+                   D1[0]=D2[0];\r
+                   D1[1]=D2[1];\r
+                   D2[0]=DIFF2[0];\r
+                   D2[1]=DIFF2[1];\r
+                   continue;\r
+               } // if (NANGS != NINTS)\r
+               break;\r
+           } // else\r
+       } // while (true)\r
+       \r
+       return result;\r
+\r
+   } // private double ARGIN1\r
+\r
+   //COMPLEX FUNCTION CJACSU(X,N,A,B,H,CO)\r
+   private double[] CJACSU(double X[], int N, double A[], double B[], double H,\r
+                  double CO[]) {\r
+       // INTEGER N\r
+       // REAL A(*),B(*),H,CO(*)\r
+       // COMPLEX X\r
+\r
+       // ..TO CALCULATE SUMMATION{CO(K+1)*P(K,X)},K=0(1)N, WHERE P(K,X)\r
+       // ..DENOTES THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE K\r
+       // ..EVALUATED AT X, ARRAY CO STORES A GIVEN SET OF COEFFICIENTS,\r
+       // ..ARRAYS A,B STORE THE COEFFICIENTS IN THE THREE-TERM\r
+       // ..RECURRENCE FORMULA FOR THE JACOBI POLYNOMIALS (SEE ASONJ7)\r
+       // ..AND H IS THE SQUARED 2-NORM OF UNITY.\r
+\r
+       double PREV[] = new double[2];\r
+       double CURR[] = new double[2];\r
+       double NEXT[] = new double[2];\r
+          // COMPLEX PREV,CURR,NEXT\r
+       int K;\r
+       double result[] = new double[2];\r
+       double cr[] = new double[1];\r
+       double ci[] = new double[1];\r
+\r
+       if (N == 0) {\r
+           result[0]=CO[0]/Math.sqrt(H);\r
+           result[1]=0.0;\r
+       } // if (N == 0)\r
+       else if (N > 0) {\r
+           PREV[0]=CO[N];\r
+           PREV[1]=0.0;\r
+           zmlt(X[0]-B[N-1],X[1],PREV[0],PREV[1],cr,ci);\r
+           CURR[0] = CO[N-1] + cr[0]/A[N-1];\r
+           CURR[1] = ci[0]/A[N-1];\r
+           for (K=N-2; K >= 0; K--) {\r
+                  zmlt(X[0]-B[K],X[1],CURR[0],CURR[1],cr,ci);\r
+                  NEXT[0]=CO[K] + cr[0]/A[K] - A[K]*PREV[0]/A[K+1];\r
+                  NEXT[1]=ci[0]/A[K] - A[K]*PREV[1]/A[K+1];\r
+               PREV[0]=CURR[0];\r
+               PREV[1]=CURR[1];\r
+               CURR[0]=NEXT[0];\r
+               CURR[1]=NEXT[1];\r
+           } // for (K=N-2; K >= 0; K--)\r
+           result[0]=CURR[0]/Math.sqrt(H);\r
+           result[1]=CURR[1]/Math.sqrt(H);\r
+       } // else if (N > 0)\r
+       else {\r
+           result[0] = 0.0;\r
+           result[1] = 0.0;\r
+       }\r
+       return result;\r
+\r
+   } // private double[] CJACSU\r
+   \r
+   //COMPLEX FUNCTION ZTOB1(Z,B1,JT,IN)\r
+   private double[] ZTOB1(double Z[], double B1, int JT, boolean IN) {\r
+       // INTEGER JT\r
+       // REAL B1\r
+       // COMPLEX Z\r
+       // LOGICAL IN\r
+\r
+       // TO COMPUTE Z**B1 BUT CHOOSING THE BRANCH CUT TO POINT ALONG THE\r
+       // RAY WITH POLAR ANGLE DEFINED BY THE VARIABLE *CUT* BELOW.\r
+\r
+       // THIS IS A SPECIAL PURPOSE ROUTINE IN THAT IT IS ASSUMED THAT\r
+       // B1=1/ALHPA, WHERE ALPHA*PI IS THE INTERIOR ANGLE AT THE BRANCH\r
+       // POINT OF THE MAP : PHYSICAL --> CANONICAL.  HENCE THE USE OF THIS\r
+       // ROUTINE IS EFFECTIVELY RESTRCITED TO COMPUTING ONLY THE BOUNDARY\r
+       // CORRESPONDENCE FUNCTION FOR COMPLEX PARAMETERS.\r
+\r
+       double CUT;\r
+       double W[] = new double[2];\r
+       // COMPLEX W\r
+       final double PI = Math.PI;\r
+       final double TUPI = 2.0*PI;\r
+       double result[] = new double[2];\r
+       double mag;\r
+       double ang;\r
+\r
+       if (zabs(Z[0],Z[1]) == 0.0) {\r
+           result[0] = 0.0;\r
+           result[1] = 0.0;\r
+           return result;\r
+       }\r
+       W[0] = Math.log(zabs(Z[0],Z[1]));\r
+       W[1] = Math.atan2(Z[1],Z[0]);\r
+       CUT=PI*(0.5-B1)/B1;\r
+       if ((JT < 0.0 && IN) || (JT > 0.0 && !IN)) {\r
+           CUT=-CUT;\r
+           if (W[1] > CUT) {\r
+                  W[1] = W[1] - TUPI;\r
+           }\r
+       }\r
+       else if (W[1] < CUT) {\r
+          W[1] = W[1] + TUPI;\r
+       }\r
+       mag = Math.exp(W[0]*B1);\r
+       ang = W[1]*B1;\r
+       result[0] = mag * Math.cos(ang);\r
+       result[1] = mag * Math.sin(ang);\r
+       return result;\r
+\r
+   } // private double[] ZTOB1\r
+\r
+   private void PPSBI1(double ZZ[], double BETA, int NQUAD, int DGPOL, double ACOEF[], double BCOEF[],\r
+       double H0VAL, double SOLUN[][], double TOLIN, double TOLOU[], double XENPT[], int QINTS[],\r
+       int MQIN1, int IER[]) {\r
+          // INTEGER DGPOL,MQIN1,NQUAD,QINTS,IER\r
+          // REAL BETA,H0VAL,TOLIN,TOLOU\r
+          // REAL ACOEF(*),BCOEF(*),XENPT(*)\r
+          // COMPLEX SOLUN(*),ZZ\r
+               \r
+               //     CALCULATES THE NUMBER OF QUADRATURE INTERVALS (QINTS) REQUIRED\r
+               //     FOR THE COMPOSITE GAUSS-JACOBI/GAUSS-LEGENDRE ESTIMATION OF\r
+               \r
+               //          INTEGRAL  [(1+X)**BETA*FNPHI(X)*LOG(ZZ-X)*dX], \r
+               //         -1<=X<=1                                         \r
+               \r
+               //     WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY \r
+               //     CORRESPONDENCE DERIVATIVE - JACOBI WEIGHT QUOTIENT. \r
+                \r
+               //     ZZ IS ANY POINT IN THE PLANE.\r
+               \r
+               //     THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL,SOLUN ARE USED TO DEFINE\r
+               //     FNPHI AND ARE PASSED TO DEPPJ9 AND DEPPL9 FOR THIS PURPOSE.\r
+               \r
+               //     THE ENDPOINTS OF THE QUADRATURE INTERVALS ARE RETURNED IN VECTOR \r
+               //     XENPT, WITH XENPT(1)=-1<XENPT(2)<...<1=XENPT(QINTS+1).\r
+               \r
+               //     TOLOU RECORDS OUR ESTIMATE FOR THE MAXIMUM OF THE ABSOLUTE VALUES\r
+               //     OF THE REMAINDER ESTIMATES ON THE DELTA-CONTOUR.  WE REQUIRE THAT\r
+               //\r
+               //                  TOLOU <= TOLIN\r
+               \r
+               //     WITH THE TOLOU BEING REASONABLY CLOSE TO TOLIN.\r
+               \r
+               //     IER=0  - NORMAL EXIT\r
+               //     IER=29 - MQIN1 SHOULD BE INCREASED IN PHTCA1 (CHANGED TO ERROR NO.\r
+               //              40 IF CALLED FROM CINRAD OR ERROR NO. 49 IF CALLED FROM\r
+               //              CATPH4)\r
+                  \r
+               //     LOCAL VARIABLES\r
+               \r
+               double TAU[] = new double[1];\r
+               double TOL[] = new double[1];\r
+               double MAXRM[] = new double[1];\r
+           double RIGHT[] = new double[1];\r
+               boolean T1FXD,FIRST;\r
+               // EXTERNAL DEPPJ9,DEPPL9\r
+               double TAU1[] = new double[1];\r
+               double ZZIN[][] = new double[1][2];\r
+               \r
+               QINTS[0]=1;\r
+               XENPT[0]=-1.0;\r
+               TOL[0]=TOLIN;\r
+               TAU[0]=1.0;\r
+               ZZIN[0][0] = ZZ[0];\r
+               ZZIN[0][1] = ZZ[1];\r
+               DEPPJ9(ZZIN,1,BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL,SOLUN,TOL,MAXRM,IER);\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
+                   FIRST=true;\r
+                   DEPPL9(ZZIN,1,BETA,RIGHT,TAU,T1FXD,NQUAD,DGPOL,ACOEF,\r
+                   BCOEF,H0VAL,SOLUN,TOL,MAXRM,FIRST,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
+                         FIRST=true;\r
+                         TAU1[0] = XENPT[QINTS[0]-1];\r
+                         DEPPL9(ZZIN,1,BETA,TAU1,TAU,T1FXD,NQUAD,DGPOL,\r
+                                         ACOEF,BCOEF,H0VAL,SOLUN,TOL,MAXRM,FIRST,IER);\r
+                         XENPT[QINTS[0]-1] = TAU1[0];\r
+                         TOLOU[0]=TOLOU[0]+MAXRM[0];\r
+                         QINTS[0]=QINTS[0]+1;\r
+                         if (QINTS[0] >= MQIN1) {\r
+                           IER[0]=29;\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 PPSBI1\r
 \r
 \r
       /**\r