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

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

index 2db7415..1eba8bc 100644 (file)
@@ -10801,11 +10801,38 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
        final int MNDG = 20;\r
        final int MNQD = 128;\r
-          int AJT,AJTC,DG,DGC,I,I1,IC,JT,JTC,K,K1,LOD,LODC,LOL,LOM,\r
-                   LOS,NQUAD,PSA,PT,QINTS;\r
+       int PSA = 1;\r
+       int JT = 0;\r
+       int NQUAD = 0;\r
+       int LODC = 0;\r
+       int QINTS = 0;\r
+       int DG = 0;\r
+       int LOD = 0;\r
+       int DGC = 0;\r
+       int LOM = 0;\r
+       int JTC = 0;\r
+       int PT = 0;\r
+          int AJT,AJTC,I,I1,IC,K,K1,LOL,LOS;\r
           final double RFAC = 10.0;\r
-          double AA,BB,BC1,BETA,BETAC,CONST,H0,H0C,H1,HA,HB1,HL,LL,MD,MXPT,\r
-                    QHLEN,RRHS,RSLN,SJT,SJTC,TERM,TOLIW,TT,UU,XX;\r
+          double HB1 = 0.0;\r
+          double QHLEN = 0.0;\r
+          double BETAC = 0.0;\r
+          double AA = 0.0;\r
+          double BB = 0;\r
+          double CONST = 0.0;\r
+          double H0 = 0.0;\r
+          double H0C = 0.0;\r
+          double H1 = 0.0;\r
+          double HA = 0.0;\r
+          double HL = 0.0;\r
+          double BETA = 0.0;\r
+          double SJT = 0.0;\r
+          double SJTC = 0.0;\r
+          double MD = 0.0;\r
+          double TT = 0.0;\r
+          double BC1 = 0.0;\r
+          double MXPT = 0.0;\r
+          double LL,RRHS,RSLN,TERM,TOLIW,UU,XX;\r
           double JACOF[] = new double[MNDG];\r
           double QAB[] = new double[MNQD];\r
           double QWT[] = new double[MNQD];\r
@@ -10815,6 +10842,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           double NEW[][] = new double[MNDG][2];\r
           double OLD[][] = new double[MNDG][2];\r
           double RHOVL[][] = new double[MNQD][2];\r
+          int i;\r
           // COMPLEX NEW(MNDG),OLD(MNDG),RHOVL(MNQD)\r
           // EXTERNAL BISNEW,INVJCO,R1MACH\r
                \r
@@ -10822,10 +10850,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           LOL=(NJIND-1)*NQPTS+1;\r
           IC=1;\r
           LOS=1;\r
+          boolean do1 = true;\r
+          boolean do2 = true;\r
                \r
-          /*while (true) {\r
+          outer: while (true) {\r
                \r
-              if (IC > TNSUC) {\r
+            if (do1) {\r
+                  if (IC > TNSUC) {\r
                \r
                       // NORMAL EXIT\r
                \r
@@ -10853,245 +10884,611 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                   HL=HALEN[PSA-1];\r
                   MD=MIDPT[PSA-1];\r
                   PT=PARNT[PSA-1];\r
-                     DO 20 I=1,DG+1\r
-                       JACOF(I)=SOLUN(I+LOM-1)\r
-               20    CONTINUE\r
-                     DO 30 I=2,DG+1,2\r
-                       JACOF(I)=SJT*JACOF(I)\r
-               30    CONTINUE\r
-               C\r
-               C     INITIALISATION FOR ARC NUMBER IC ON CIRCLE\r
-               C\r
-               40    CONTINUE\r
-                     DO 50 I=1,NQPTS\r
-                       OLD(I)=0E+0\r
-               50    CONTINUE\r
-                     QINTS=1\r
-                     QHLEN=1E+0\r
-                     NQUAD=NQPTS\r
-                     DGC=NQPTS-1\r
-                     IF (DGC+1 .GT. MNDG) THEN\r
-                       IER=30\r
-                       RETURN\r
-                     ENDIF\r
-                     JTC=JTYPC(IC)\r
-                     AJTC=ABS(JTC)\r
-                     SJTC=REAL(SIGN(1,JTC))\r
-                     LODC=(AJTC-1)*NQPTS+1\r
-                     BETAC=JAINC(AJTC)\r
-                     H0C=H0VLC(AJTC)\r
-                     HA=(VARGC(IC+1)-VARGC(IC))*5E-1\r
-                     RSLN=HA/ERARC(PSA)\r
-                     MXPT=RSLN/RFAC\r
-                     BC1=BETAC+1E+0\r
-                     HB1=1E+0\r
-               C\r
-               C     SET UP RIGHT HAND SIDE *CONST* FOR THE BOUNDARY CORRESPONDENCE\r
-               C     EQUATION THAT WILL BE USED TO COMPUTE PHYSICAL PARAMETERS \r
-               C     CORRESPONDING TO GIVEN POINTS ON THE CIRCLE.\r
-               C\r
-                     IF (JT .LT. 0) THEN\r
-                       CONST=VTARG(PSA+1)-VARGC(IC+1)\r
-                     ELSE\r
-                       CONST=VARGC(IC)-VTARG(PSA)   \r
-                     ENDIF\r
-               C\r
-               C     SET UP AA,BB WHERE THE PHYSICAL ARC IS CORRESPONDS TO THE \r
-               C     PARAMETER INTERVAL [AA,BB].\r
-               C\r
-                     IF (PHPAS(IC+1) .LE. PHPAS(IC)) THEN\r
-                       BB=1E+0\r
-                     ELSE\r
-                       BB=PHPAS(IC+1)\r
-                     ENDIF\r
-                     AA=PHPAS(IC)\r
-               60    CONTINUE\r
-               C\r
-               C     SET UP THE (POSSIBLY) COMPOSITE QUADRATURE RULE BASED ON *QINTS*\r
-               C     SUBINTERVALS\r
-               C\r
-                     IF (NQUAD .GT. MNQD) THEN\r
-                       IER=31\r
-                       RETURN\r
-                     ENDIF\r
-                     DO 70 K1=1,NQPTS\r
-                       I1=LODC+K1-1\r
-                       QWT(K1)=HB1*QUWTC(I1)\r
-                       QAB(K1)=-1E+0+QHLEN*(1E+0+QUPTC(I1))\r
-               70    CONTINUE\r
-                     K1=NQPTS\r
-                     DO 90 K=2,QINTS\r
-                       DO 80 I=1,NQPTS\r
-                         K1=K1+1\r
-                         I1=LOL+I-1\r
-                         XX=2E+0*K-1E+0+QUPTC(I1)\r
-                         QWT(K1)=HB1*XX**BETAC*QUWTC(I1)\r
-                         QAB(K1)=-1E+0+QHLEN*XX\r
-               80      CONTINUE\r
-               90    CONTINUE\r
-               C\r
-               C     ESTIMATE THE JACOBI COEFFICIENTS FOR THE INVERSE DENSITY FOR\r
-               C     ARC NUMBER IC ON THE CIRCLE.\r
-               C\r
-                     CALL INVJCO(NEW,AICOF(LOD),AA,ACOEF(LOD),ACOFC(LODC),\r
-                    +     BICOF(LOD),BB,BCFSN(LOM),BCOEF(LOD),BCOFC(LODC),BETA,BETAC,\r
-                    +     CENTR,CONST,DGC,DG,H0,H0C,H1,HA,HL,IER,INTER,JACOF,JTC,MD,\r
-                    +     NEWTL,NQUAD,PT,QAB,QWT,RHOVL,SJT,SJTC,SVAL,TOLIW,TVAL,WORK)\r
-               C\r
-               C     CHECK THAT THE SIZE OF THE HIGHEST DEGREE COEFFICIENT IS SMALL\r
-               C     ENOUGH.\r
-               C\r
-                     TERM=ABS(NEW(NQPTS))*RIGLL(LODC+NQPTS-1)/LGTOL\r
-                     IF (TERM .GT. 1E+0) THEN\r
-               C\r
-               C       COEFFICIENT IS TOO LARGE - SUBDIVIDE CIRCULAR ARC AND RESTART.\r
-               C   \r
-               C       FIRST FIND THE LOCAL PHYSICAL PARAMETER *TT* CORRESPONDING TO\r
-               C       THE MIDPOINT OF THE CURRENT CIRCULAR ARC NUMBER IC.\r
-               C\r
-                       RRHS=CONST+HA\r
-                       LL=AA\r
-                       UU=BB\r
-                       CALL BISNEW(IER,LL,TT,UU,AICOF(LOD),ACOEF(LOD),BICOF(LOD),\r
-                    +              BCFSN(LOM),BCOEF(LOD),BETA,DG,H0,H1,JACOF,NEWTL,SJT,\r
-                    +              RRHS,TOLIW)\r
-                       IF (IER.GT.0) THEN\r
-                         RETURN\r
-                       ENDIF\r
-               C\r
-               C       NEXT UPDATE VARIOUS DATA ITEMS TO DESCRIBE THE NEW SUBDIVISION\r
-               C       OF THE CIRCLE.\r
-               C\r
-                       DO 100 I=TNSUC+1,IC+1,-1\r
-                         PHPAS(I+1)=PHPAS(I)\r
-                         VARGC(I+1)=VARGC(I)\r
-               100     CONTINUE\r
-                       PHPAS(IC+1)=TT\r
-                       VARGC(IC+1)=(VARGC(IC+1)+VARGC(IC))*5E-1\r
-                       DO 110 I=TNSUC,IC,-1\r
-                         PRNSA(I+1)=PRNSA(I)\r
-               110     CONTINUE\r
-                       DO 120 I=TNSUC,IC+1,-1\r
-                         JTYPC(I+1)=JTYPC(I)\r
-               120     CONTINUE\r
-                       IF (JTC .GT. 0) THEN\r
-                         JTYPC(IC+1)=NJIND\r
-                       ELSE\r
-                         JTYPC(IC+1)=JTC\r
-                         JTYPC(IC)=NJIND\r
-                       ENDIF\r
-                       TNSUC=TNSUC+1\r
-                       IF (TNSUC .GE. MNSUC) THEN\r
-                         IER=32\r
-                         RETURN\r
-                       ENDIF\r
-               C\r
-               C       START AGAIN WITH THE NEW REFINED ARC NUMBER IC\r
-               C\r
-                       GOTO 40\r
-                     ENDIF\r
-               C\r
-               C     ARC REFINEMENT DOES NOT SEEM TO BE REQUIRED.  EXAMINE THE \r
-               C     JACOBI COEFFICIENTS TO ESTIMATE THE DEGREE OF POLYNOMIAL\r
-               C     APPROXIMATION REQUIRED AND ALSO TEST FOR CONVERGENCE OF THE\r
-               C     SIGNIFICANT COEFFICIENTS.\r
-               C\r
-               130   CONTINUE\r
-                     DGC=DGC-1\r
-                     IF (DGC .LT. 0) THEN\r
-               C\r
-               C       ACCEPT THAT A POLYNOMIAL APPROXIMATION OF DEGREE ZERO WILL\r
-               C       DO FOR THIS ARC AND MOVE ON TO THE NEXT ARC.\r
-               C\r
-                       DGPOC(IC)=0\r
-                       LSUBC(IC)=LOS\r
-                       IF (LOS .GT. MNCOF) THEN\r
-                         IER=33\r
-                         RETURN\r
-                       ENDIF\r
-                       SOLNC(LOS)=NEW(1)\r
-                       LOS=LOS+1\r
-                       IC=IC+1\r
-                       GOTO 10\r
-                     ENDIF\r
-               C\r
-                     TERM=ABS(NEW(DGC+1))*RIGLL(LODC+DGC)/LGTOL\r
-                     IF (TERM .LE. 1E+0) THEN\r
-               C\r
-               C       THIS COEFFICIENT MAY BE IGNORED - CONSIDER THE COEFFICIENT FOR\r
-               C       NEXT LOWER DEGREE POLYNOMIAL.\r
-               C\r
-                       GOTO 130\r
-                     ENDIF\r
-               C\r
-               C     THE DEGREE IS POSSIBLY DGC+1; CHECK FOR CONVERGENCE OF THESE\r
-               C     COEFFICIENTS.\r
-               C\r
-                     I=DGC\r
-               140   CONTINUE\r
-                     TERM=ABS(NEW(I+1)-OLD(I+1))*RIGLL(LODC+I)/LGTOL\r
-                     IF (TERM .LE. 1E+0) THEN\r
-               C\r
-               C       CONVERGENCE FOR THIS TERM\r
-               C\r
-                       I=I-1\r
-                       IF (I .GT. 0) THEN\r
-               C\r
-               C         TAKE COEFFICIENT OF NEXT LOWER DEGREE.\r
-               C\r
-                         GOTO 140\r
-                       ELSE\r
-               C\r
-               C         ALL COEFFICIENTS HAVE CONVERGED.\r
-               C\r
-                         DGPOC(IC)=DGC+1\r
-                         LSUBC(IC)=LOS\r
-                         IF (LOS+DGC .GE. MNCOF) THEN\r
-                           IER=33\r
-                           RETURN\r
-                         ENDIF\r
-                         DO 150 I=1,DGC+2\r
-                           SOLNC(LOS+I-1)=NEW(I)\r
-               150       CONTINUE\r
-                         LOS=LOS+DGC+2\r
-                         IC=IC+1\r
-                         GOTO 10\r
-                       ENDIF\r
-                     ELSE\r
-               C\r
-               C       THIS TERM HASN'T CONVERGED - TRY AGAIN WITH REFINED QUADRATURE\r
-               C       RULE, IF RESOLUTION PERMITS.\r
-               C\r
-                       QINTS=QINTS*2\r
-                       NQUAD=QINTS*NQPTS\r
-                       IF (NQUAD .GE. MXPT) THEN\r
-               C\r
-               C         FURTHER REFINEMENT IS PRACTICALLY UNACCEPTABLE DUE TO LOCAL\r
-               C         CROWDING - ACCEPT CURRENT SOLUTION\r
-               C\r
-                         DGPOC(IC)=DGC+1\r
-                         LSUBC(IC)=LOS\r
-                         IF (LOS+DGC .GE. MNCOF) THEN\r
-                           IER=33\r
-                           RETURN\r
-                         ENDIF\r
-                         DO 160 I=1,DGC+2\r
-                           SOLNC(LOS+I-1)=NEW(I)\r
-               160       CONTINUE\r
-                         LOS=LOS+DGC+2\r
-                         IC=IC+1\r
-                         GOTO 10\r
-                       ENDIF\r
-                       QHLEN=QHLEN*5E-1\r
-                       HB1=QHLEN**BC1\r
-                       DGC=NQPTS-1\r
-                       DO 170 I=1,NQPTS\r
-                         OLD(I)=NEW(I)\r
-               170     CONTINUE\r
-                       GOTO 60\r
-                     ENDIF\r
-          } // while (true) */\r
-    } // private void JCFIM5 \r
+                  for (I=1; I <= DG+1; I++) {\r
+                      JACOF[I-1]=SOLUN[I+LOM-2];\r
+                  }\r
+                  for (I=2; I <= DG+1; I +=2) {\r
+                       JACOF[I-1]=SJT*JACOF[I-1];\r
+                  }\r
+            } // if (do1)\r
+            do1 = true;\r
+               \r
+                  // INITIALISATION FOR ARC NUMBER IC ON CIRCLE\r
+               \r
+                  while (true) {\r
+                          if (do2) {\r
+                      for (I=1; I <= NQPTS; I++) {\r
+                          OLD[I-1][0]=0.0;\r
+                          OLD[I-1][1] = 0.0;\r
+                      }\r
+                      QINTS=1;\r
+                      QHLEN=1.0;\r
+                      NQUAD=NQPTS;\r
+                      DGC=NQPTS-1;\r
+                      if (DGC+1 > MNDG) {\r
+                          IER[0]=30;\r
+                          return;\r
+                      }\r
+                      JTC=JTYPC[IC-1];\r
+                      AJTC=Math.abs(JTC);\r
+                      if (JTC >= 0.0) {\r
+                          SJTC = 1.0;\r
+                      }\r
+                      else {\r
+                          SJTC = -1.0;\r
+                      }\r
+                      LODC=(AJTC-1)*NQPTS+1;\r
+                      BETAC=JAINC[AJTC-1];\r
+                      H0C=H0VLC[AJTC-1];\r
+                      HA=(VARGC[IC]-VARGC[IC-1])*0.5;\r
+                      RSLN=HA/ERARC[PSA-1];\r
+                      MXPT=RSLN/RFAC;\r
+                      BC1=BETAC+1.0;\r
+                      HB1=1.0;\r
+               \r
+                      // SET UP RIGHT HAND SIDE *CONST* FOR THE BOUNDARY CORRESPONDENCE\r
+                      // EQUATION THAT WILL BE USED TO COMPUTE PHYSICAL PARAMETERS \r
+                      // CORRESPONDING TO GIVEN POINTS ON THE CIRCLE.\r
+               \r
+                      if (JT < 0) {\r
+                          CONST=VTARG[PSA]-VARGC[IC];\r
+                      }\r
+                      else {\r
+                          CONST=VARGC[IC-1]-VTARG[PSA-1];   \r
+                      }\r
+               \r
+                      // SET UP AA,BB WHERE THE PHYSICAL ARC IS CORRESPONDS TO THE \r
+                      // PARAMETER INTERVAL [AA,BB].\r
+               \r
+                      if (PHPAS[IC] <= PHPAS[IC-1]) {\r
+                          BB=1.0;\r
+                      }\r
+                      else {\r
+                          BB=PHPAS[IC];\r
+                      }\r
+                      AA=PHPAS[IC-1];\r
+                          } // if (do2)\r
+                          do2 = true;\r
+               \r
+                      // SET UP THE (POSSIBLY) COMPOSITE QUADRATURE RULE BASED ON *QINTS*\r
+                      // SUBINTERVALS\r
+               \r
+                      if (NQUAD > MNQD) {\r
+                          IER[0]=31;\r
+                          return;\r
+                      }\r
+                      for (K1=1; K1 <= NQPTS; K1++) {\r
+                          I1=LODC+K1-1;\r
+                          QWT[K1-1]=HB1*QUWTC[I1-1];\r
+                          QAB[K1-1]=-1.0+QHLEN*(1.0+QUPTC[I1-1]);\r
+                      } // for (K1=1; K1 <= NQPTS; K1++)\r
+                      K1=NQPTS;\r
+                      for (K=2; K <= QINTS; K++) {\r
+                          for (I=1; I <= NQPTS; I++) {\r
+                              K1=K1+1;\r
+                              I1=LOL+I-1;\r
+                              XX=2.0*K-1.0+QUPTC[I1-1];\r
+                              QWT[K1-1]=HB1*Math.pow(XX,BETAC)*QUWTC[I1-1];\r
+                              QAB[K1-1]=-1.0+QHLEN*XX;\r
+                          } // for (I=1; I <= NQPTS; I++)\r
+                      } // for (K=2; K <= QINTS; K++)\r
+               \r
+                      // ESTIMATE THE JACOBI COEFFICIENTS FOR THE INVERSE DENSITY FOR\r
+                      // ARC NUMBER IC ON THE CIRCLE.\r
+                      double A1COF[]= new double[DG-1];\r
+                      double B1COF[] = new double[DG-1];\r
+                      for (i = 0; i < DG-1; i++) {\r
+                          A1COF[i] = AICOF[LOD+i-1];\r
+                          B1COF[i] = BICOF[LOD+i-1];\r
+                      }\r
+                      double ACOEFIN[]= new double[DG];\r
+                      double BCOEFIN[] = new double[DG];\r
+                      for (i = 0; i < DG; i++) {\r
+                          ACOEFIN[i] = ACOEF[LOD+i-1];\r
+                          BCOEFIN[i] = BCOEF[LOD+i-1];\r
+                      }\r
+                      double ACOFCIN[] = new double[DGC];\r
+                      double BCOFCIN[] = new double[DGC];\r
+                      for (i = 0; i < DGC; i++) {\r
+                          ACOFCIN[i] = ACOFC[LODC+i-1];\r
+                          BCOFCIN[i] = BCOFC[LODC+i-1];\r
+                      }\r
+                      double BCFSNIN[] = new double[DG+1];\r
+                      for (i = 0; i < DG+1; i++) {\r
+                          BCFSNIN[i] = BCFSN[LOM+i-1];\r
+                      }\r
+                      INVJCO(NEW,A1COF,AA,ACOEFIN,ACOFCIN,\r
+                         B1COF,BB,BCFSNIN,BCOEFIN,BCOFCIN,BETA,BETAC,\r
+                         CENTR,CONST,DGC,DG,H0,H0C,H1,HA,HL,IER,INTER,JACOF,JTC,MD,\r
+                         NEWTL,NQUAD,PT,QAB,QWT,RHOVL,SJT,SJTC,SVAL,TOLIW,TVAL,WORK);\r
+               \r
+                      // CHECK THAT THE SIZE OF THE HIGHEST DEGREE COEFFICIENT IS SMALL\r
+                      // ENOUGH.\r
+               \r
+                      TERM=zabs(NEW[NQPTS-1][0],NEW[NQPTS-1][1])*RIGLL[LODC+NQPTS-2]/LGTOL;\r
+                      if (TERM > 1.0) {\r
+               \r
+                          // COEFFICIENT IS TOO LARGE - SUBDIVIDE CIRCULAR ARC AND RESTART.\r
+                  \r
+                          // FIRST FIND THE LOCAL PHYSICAL PARAMETER *TT* CORRESPONDING TO\r
+                          // THE MIDPOINT OF THE CURRENT CIRCULAR ARC NUMBER IC.\r
+               \r
+                          RRHS=CONST+HA;\r
+                          LL=AA;\r
+                          UU=BB;\r
+                          BISNEW(IER,LL,TT,UU,A1COF,ACOEFIN,B1COF,\r
+                                 BCFSNIN,BCOEFIN,BETA,DG,H0,H1,JACOF,NEWTL,SJT,RRHS,TOLIW);\r
+                          if (IER[0] > 0) {\r
+                              return;\r
+                          }\r
+               \r
+                          // NEXT UPDATE VARIOUS DATA ITEMS TO DESCRIBE THE NEW SUBDIVISION\r
+                          // OF THE CIRCLE.\r
+               \r
+                          for (I=TNSUC+1; I >= IC+1; I--) {\r
+                              PHPAS[I]=PHPAS[I-1];\r
+                              VARGC[I]=VARGC[I-1];\r
+                          } // for (I=TNSUC+1; I >= IC+1; I--)\r
+                          PHPAS[IC]=TT;\r
+                          VARGC[IC]=(VARGC[IC]+VARGC[IC-1])*0.5;\r
+                          for (I=TNSUC; I >= IC; I--) {\r
+                              PRNSA[I]=PRNSA[I-1];\r
+                          }\r
+                          for (I=TNSUC; I >= IC+1; I--) {\r
+                              JTYPC[I]=JTYPC[I-1];\r
+                          }\r
+                          if (JTC > 0) {\r
+                              JTYPC[IC]=NJIND;\r
+                          }\r
+                          else {\r
+                              JTYPC[IC]=JTC;\r
+                              JTYPC[IC-1]=NJIND;\r
+                          }\r
+                          TNSUC=TNSUC+1;\r
+                          if (TNSUC >= MNSUC) {\r
+                              IER[0]=32;\r
+                              return;\r
+                          }\r
+               \r
+                          // START AGAIN WITH THE NEW REFINED ARC NUMBER IC\r
+               \r
+                         continue;\r
+                     } // if (TERM > 1.0)\r
+                     else {\r
+                         break;\r
+                     }\r
+                  } // while (true)\r
+               \r
+                  // ARC REFINEMENT DOES NOT SEEM TO BE REQUIRED.  EXAMINE THE \r
+                  // JACOBI COEFFICIENTS TO ESTIMATE THE DEGREE OF POLYNOMIAL\r
+                  // APPROXIMATION REQUIRED AND ALSO TEST FOR CONVERGENCE OF THE\r
+                  // SIGNIFICANT COEFFICIENTS.\r
+               \r
+                  while (true) {\r
+                      DGC=DGC-1;\r
+                      if (DGC < 0) {\r
+               \r
+                          // ACCEPT THAT A POLYNOMIAL APPROXIMATION OF DEGREE ZERO WILL\r
+                          // DO FOR THIS ARC AND MOVE ON TO THE NEXT ARC.\r
+               \r
+                          DGPOC[IC-1]=0;\r
+                          LSUBC[IC-1]=LOS;\r
+                          if (LOS > MNCOF) {\r
+                              IER[0]=33;\r
+                              return;\r
+                          } // if (LOS > MNCOF)\r
+                          SOLNC[LOS-1][0]=NEW[0][0];\r
+                          SOLNC[LOS-1][1]=NEW[0][1];\r
+                          LOS=LOS+1;\r
+                          IC=IC+1;\r
+                          continue outer;\r
+                      } // if (DGC < 0)\r
+               \r
+                      TERM=zabs(NEW[DGC][0],NEW[DGC][1])*RIGLL[LODC+DGC-1]/LGTOL;\r
+                      if (TERM <= 1.0) {\r
+               \r
+                          // THIS COEFFICIENT MAY BE IGNORED - CONSIDER THE COEFFICIENT FOR\r
+                          // NEXT LOWER DEGREE POLYNOMIAL.\r
+               \r
+                          continue;\r
+                      }\r
+                      else {\r
+                          break;\r
+                      }\r
+                  } // while (true)\r
+               \r
+                  // THE DEGREE IS POSSIBLY DGC+1; CHECK FOR CONVERGENCE OF THESE\r
+                  // COEFFICIENTS.\r
+               \r
+                  I=DGC;\r
+                  while (true) {\r
+                      TERM=zabs(NEW[I][0]-OLD[I][0],NEW[I][1]-OLD[I][1])*RIGLL[LODC+I-1]/LGTOL;\r
+                     if (TERM <= 1.0) {\r
+               \r
+                         // CONVERGENCE FOR THIS TERM\r
+               \r
+                         I=I-1;\r
+                         if (I > 0) {\r
+               \r
+                             // TAKE COEFFICIENT OF NEXT LOWER DEGREE.\r
+               \r
+                             continue;\r
+                         } // if (I > 0)\r
+                         else { // I <= 0\r
+               \r
+                             // ALL COEFFICIENTS HAVE CONVERGED.\r
+               \r
+                             DGPOC[IC-1]=DGC+1;\r
+                             LSUBC[IC-1]=LOS;\r
+                             if (LOS+DGC >= MNCOF) {\r
+                                 IER[0]=33;\r
+                                 return;\r
+                             }\r
+                             for (I=1; I <= DGC+2; I++) {\r
+                                 SOLNC[LOS+I-2][0]=NEW[I-1][0];\r
+                                 SOLNC[LOS+I-2][1]=NEW[I-1][1];\r
+                             } // for (I=1; I <= DGC+2; I++)\r
+                             LOS=LOS+DGC+2;\r
+                             IC=IC+1;\r
+                             continue outer;\r
+                         } // else I <= 0\r
+                     } // if (TERM <= 1.0)\r
+                     else { // TERM > 1.0\r
+               \r
+                         // THIS TERM HASN'T CONVERGED - TRY AGAIN WITH REFINED QUADRATURE\r
+                         // RULE, IF RESOLUTION PERMITS.\r
+               \r
+                         QINTS=QINTS*2;\r
+                         NQUAD=QINTS*NQPTS;\r
+                         if (NQUAD >= MXPT) {\r
+               \r
+                             // FURTHER REFINEMENT IS PRACTICALLY UNACCEPTABLE DUE TO LOCAL\r
+                             // CROWDING - ACCEPT CURRENT SOLUTION\r
+               \r
+                             DGPOC[IC-1]=DGC+1;\r
+                             LSUBC[IC-1]=LOS;\r
+                             if (LOS+DGC >= MNCOF) {\r
+                                 IER[0]=33;\r
+                                 return;\r
+                             }\r
+                             for (I=1; I <= DGC+2; I++) {\r
+                                 SOLNC[LOS+I-2][0]=NEW[I-1][0];\r
+                                 SOLNC[LOS+I-2][1]=NEW[I-1][1];\r
+                             }\r
+                             LOS=LOS+DGC+2;\r
+                             IC=IC+1;\r
+                             continue outer;\r
+                         } // if (NQUAD >= MXPT)\r
+                         QHLEN=QHLEN*0.5;\r
+                         HB1=Math.pow(QHLEN,BC1);\r
+                         DGC=NQPTS-1;\r
+                         for (I=1; I <= NQPTS; I++) {\r
+                             OLD[I-1][0]=NEW[I-1][0];\r
+                             OLD[I-1][1]=NEW[I-1][1];\r
+                         } // for (I=1; I <= NQPTS; I++\r
+                       do1 = false;\r
+                       do2 = false;\r
+                       continue outer;\r
+                       } // while (true)\r
+              } // else TERM > 1.0\r
+          } // outer: while (true)\r
+    } // private void JCFIM5\r
+   \r
+   private void INVJCO(double SOLNC[][], double A1COF[], double AA, double ACOEF[], double ACOFC[],\r
+       double B1COF[], double BB, double BCFSN[], double BCOEF[], double BCOFC[], double BETA, double BETAC,\r
+       double CENTR[], double CONST, int DGPOC, int DGPOL, double H0VAL, double H0VLC, double H1VAL,\r
+       double HAANG, double HALEN, int IER[], boolean INTER, double JACOF[], int JTYPC, double MIDPT,\r
+       double NEWTL, int NQUAD, int PARNT, double QUPTC[], double QUWTC[], double RHOVL[][], double SJT,\r
+       double SJTC, double SVAL[], double TOLIW, double TVAL[], double WORK[]) {\r
+          // INTEGER DGPOC,DGPOL,IER,JTYPC,NQUAD,PARNT\r
+          // REAL AA,BB,BETA,BETAC,CONST,H0VAL,H0VLC,H1VAL,HAANG,HALEN,MIDPT,\r
+          // +NEWTL,SJT,SJTC,TOLIW\r
+          // REAL A1COF(*),ACOEF(*),ACOFC(*),B1COF(*),BCFSN(*),BCOEF(*),\r
+          // +BCOFC(*),JACOF(*),QUPTC(*),QUWTC(*),SVAL(*),TVAL(*),WORK(*)\r
+          // COMPLEX CENTR\r
+          // COMPLEX SOLNC(*),RHOVL(*)\r
+          // LOGICAL INTER\r
+               \r
+          // COMPUTES THE JACOBI COEFFICIENT VECTOR  *SOLNC* FOR THE INVERSE  \r
+          // DENSITY FUNCTION FOR THE PARTICULAR ARC SPECIFIED BY THE OTHER\r
+          // PARAMETERS AND USING THE *NQUAD* POINT RULE STORED IN *QUPTC* AND\r
+          // *QUWTS*.\r
+               \r
+          // IER=0 - NORMAL EXIT\r
+               \r
+          // LOCAL VARIABLES\r
+               \r
+          int I,K;\r
+          double LL,RRHS,UU;\r
+          // EXTERNAL BISNEW,RHOFN\r
+               \r
+          for (I=1; I <= NQUAD; I++) {\r
+              SVAL[I-1]=SJTC*QUPTC[I-1];\r
+          }\r
+               \r
+          // GET LOCAL PHYSICAL PARAMETER VALUES *TVAL* CORRESPONDING TO \r
+          // QUADRATURE PARAMETERS *SVAL* ON CIRCULAR ARC\r
+               \r
+          for (I=1; I <= NQUAD; I++) {\r
+              RRHS=CONST+HAANG*(1.0+SJT*SVAL[I-1]);\r
+                  LL=AA;\r
+                  UU=BB;\r
+                  BISNEW(IER,LL,TVAL[I-1],UU,A1COF,ACOEF,B1COF,BCFSN,BCOEF,\r
+                   BETA,DGPOL,H0VAL,H1VAL,JACOF,NEWTL,SJT,RRHS,TOLIW);\r
+                  if (IER[0] > 0) {\r
+                      return;\r
+                  }\r
+          } // for (I=1; I <= NQUAD; I++)\r
+               \r
+          // GET VALUES OF DENSITY *RHOVL* CORRESPONDING TO *SVAL* AND\r
+          // *TVAL* \r
+               \r
+          RHOFN(IER,RHOVL,ACOEF,BCOEF,BETA,BETAC,CENTR,DGPOL,H0VAL,\r
+                    HAANG,HALEN,INTER,MIDPT,NQUAD,PARNT,SJT,JACOF,SVAL,TVAL);  \r
+               \r
+               if (IER[0] > 0) {\r
+                   return;\r
+               }\r
+               \r
+               for (I=1; I <= 1+DGPOC; I++) {\r
+                   SOLNC[I-1][0]=0.0;\r
+                   SOLNC[I-1][1] = 0.0;\r
+               } // for (I=1; I <= 1+DGPOC; I++)\r
+               \r
+               WORK[0]=1.0/Math.sqrt(H0VLC);\r
+               for (I=1; I <= NQUAD; I++) {\r
+                   JAPAR7(WORK,QUPTC[I-1],ACOFC,BCOFC,DGPOC);\r
+                   for (K=1; K <= 1+DGPOC; K++) {\r
+                       SOLNC[K-1][0]=SOLNC[K-1][0]+QUWTC[I-1]*RHOVL[I-1][0]*WORK[K-1];\r
+                       SOLNC[K-1][1]=SOLNC[K-1][1]+QUWTC[I-1]*RHOVL[I-1][1]*WORK[K-1];\r
+                   } // for (K=1; K <= 1+DGPOC; K++)\r
+               } // for (I=1; I <= NQUAD; I++)\r
+               \r
+               if (JTYPC < 0) {\r
+                   for (K=2; K <= 1+DGPOC; K +=2) {\r
+                         SOLNC[K-1][0]=-SOLNC[K-1][0];\r
+                         SOLNC[K-1][1]=-SOLNC[K-1][1];\r
+                   } // for (K=2; K <= 1+DGPOC; K +=2)\r
+               } // if (JTYPC < 0)\r
+               \r
+               // NORMAL EXIT\r
+               \r
+        IER[0]=0;\r
+               \r
+    } // private void INVJCO\r
+   \r
+   private void BISNEW(int IER[], double LL, double TT, double UU, double A1COF[], double ACOEF[],\r
+       double B1COF[], double BCFSN[], double BCOEF[], double BETA, int DEG, double H0VAL, double H1VAL,\r
+       double JACOF[], double NEWTL, double SJT, double RRHS, double TOLIW) {\r
+          // INTEGER DEG,IER\r
+          //  REAL A1COF(*),ACOEF(*),B1COF(*),BCFSN(*),BCOEF(*),BETA,H0VAL,\r
+          // +H1VAL,JACOF,LL,NEWTL,TT,UU,SJT,RRHS,TOLIW\r
+          // Should be JACOF(*)\r
+                    \r
+               //     A MIXTURE OF BISECTION AND NEWTON'S METHOD TO SOLVE THE NON-LINEAR\r
+               //     BOUNDARY CORRESPONDENCE EQUATION\r
+               \r
+               //                 THETA(TT) = CONST\r
+               \r
+               //     FOR REAL PARAMETER TT GIVEN REAL CONST; SEE RB#50 P134.  THE\r
+               //     INTERVAL (LL,UU) SHOULD BRACKET TT. \r
+               \r
+               //     IER=0   -  NORMAL EXIT\r
+               //     IER=34  -  FUNCTION HAS SAME SIGN AT LL AND UU\r
+               //     IER=35  -  ZERO FUNCTION DERIVATIVE DETECTED\r
+               \r
+               //     LOCAL VARIABLES\r
+          final int MNITS = 10;\r
+          final int NBSCT = 3;\r
+          int NITS,STEPS;\r
+          double DFT,EPS,FL,FT,FU,JACSUM,RDIFF,TUPI;\r
+          double CO[];\r
+          // EXTERNAL JACSUM\r
+          int i;\r
+               \r
+          TUPI=2.0*Math.PI;\r
+          CO = new double[DEG];\r
+          for (i = 0; i < DEG; i++) {\r
+                  CO[i] = BCFSN[1+i];\r
+          }\r
+          FL=JACSUM(SJT*LL,DEG-1,A1COF,B1COF,H1VAL,CO);\r
+          FL=BCFSN[0]-(1.0-SJT*LL)*FL;\r
+          FL=Math.pow((1.0+SJT*LL),(1.0+BETA))*FL-RRHS;\r
+               \r
+          for (i = 0; i < DEG; i++) {\r
+                  CO[i] = BCFSN[1+i];\r
+          }\r
+          FU=JACSUM(SJT*UU,DEG-1,A1COF,B1COF,H1VAL,CO);\r
+          FU=BCFSN[0]-(1.0-SJT*UU)*FU;\r
+          FU=Math.pow((1.0+SJT*UU),(1.0+BETA))*FU-RRHS;\r
+               \r
+          if (FL*FU > 0.0) {\r
+              IER[0]=34;\r
+                  return;\r
+          }\r
+               \r
+          //  ENTER NEWTON ITERATION MODE\r
+               \r
+          outer: while (true) {\r
+              TT=(UU+LL)*0.5;\r
+                  NITS=0;\r
+                  middle: while (true) {\r
+                          for (i = 0; i < DEG; i++) {\r
+                                  CO[i] = BCFSN[1+i];\r
+                          }\r
+                      FT=JACSUM(SJT*TT,DEG-1,A1COF,B1COF,H1VAL,CO);\r
+                      FT=BCFSN[0]-(1.0-SJT*TT)*FT;\r
+                      FT=Math.pow((1.0+SJT*TT),(1.0+BETA))*FT-RRHS;\r
+                      DFT=JACSUM(SJT*TT,DEG,ACOEF,BCOEF,H0VAL,JACOF);\r
+                      DFT=TUPI*Math.pow((1.0+SJT*TT),BETA)*DFT*SJT;\r
+                      if (DFT == 0.0) {\r
+                          IER[0]=35;\r
+                          return;\r
+                      }\r
+                      RDIFF=FT/DFT;\r
+                      TT=TT-RDIFF;\r
+                      NITS=NITS+1;\r
+                      if (Math.abs(RDIFF) < NEWTL) {\r
+               \r
+                          // NEWTON ITERATIONS HAVE CONVERGED FOR TT\r
+               \r
+                          IER[0]=0;\r
+                          return;\r
+                      }\r
+                      else if (TT <= LL || TT >= UU || NITS == MNITS) {\r
+               \r
+                          // PERFORM NBSCT BISECTION STEPS\r
+               \r
+                          STEPS=0;\r
+                          inner: while (true) {\r
+                              EPS=(UU-LL)*0.5;\r
+                              if (EPS < TOLIW) {\r
+               \r
+                                  // BISECTION ITERATIONS HAVE CONVERGED FOR TT\r
+               \r
+                                  IER[0]=0;\r
+                                  return;\r
+                              } // if (EPS < TOLIW)\r
+                              TT=(UU+LL)*0.5;\r
+                              for (i = 0; i < DEG; i++) {\r
+                                                  CO[i] = BCFSN[1+i];\r
+                                          }\r
+                              FT=JACSUM(SJT*TT,DEG-1,A1COF,B1COF,H1VAL,CO);\r
+                              FT=BCFSN[0]-(1.0-SJT*TT)*FT;\r
+                              FT=Math.pow((1.0+SJT*TT),(1.0+BETA))*FT-RRHS;\r
+               \r
+                              if (FT*FL < 0.0) {\r
+                                  UU=TT;\r
+                                  FU=FT;\r
+                              }\r
+                              else {\r
+                                  LL=TT;\r
+                                  FL=FT;\r
+                              }\r
+                              STEPS=STEPS+1;\r
+                              if (STEPS == NBSCT) {\r
+               \r
+                                  // RE-START NEWTON MODE\r
+               \r
+                                  continue outer;\r
+                              }\r
+                              else {\r
+               \r
+                                  // CONTINUE WITH BISECTION\r
+               \r
+                                  continue inner;\r
+                              }\r
+                          } // inner: while (true)\r
+                      } // else if (TT <= LL || TT >= UU || NITS == MNITS)\r
+                      else {\r
+               \r
+                          // CONTINUE WITH NEWTON MODE\r
+               \r
+                          continue middle;\r
+                      }\r
+                  } // middle: while (true);\r
+          } // outer: while (true)\r
+    } // private void BISNEW\r
+   \r
+   private void RHOFN(int IER[], double RHOVL[][], double ACOEF[], double BCOEF[], double BETA,\r
+       double BETAC, double CENTR[], int DGPOL, double H0VAL, double HAANG, double HALEN, boolean INTER,\r
+       double MIDPT, int NVALS, int PARNT, double SJT, double JACOF[], double SVAL[], double TVAL[]) {\r
+          // INTEGER IER,DGPOL,NVALS,PARNT\r
+          // REAL ACOEF(*),BCOEF(*),BETA,BETAC,H0VAL,HAANG,HALEN,\r
+          // +MIDPT,SJT,JACOF(*),SVAL(*),TVAL(*)\r
+          // COMPLEX CENTR,RHOVL(*)\r
+          //LOGICAL INTER\r
+               \r
+          //     GIVEN THE ARRAY *SVAL* OF PARAMETER VALUES OF POINTS ON A \r
+          //     CIRCULAR ARC AND THE ARRAY *TVAL* OF LOCAL PARAMETERS OF THE\r
+          //     CORRESPONDING POINTS ON THE PHYSICAL SUBARC, TO DETERMINE THE\r
+          //     ARRAY *RHOVL* OF VALUES OF THE FUNCTION *RHO* (SEE #50, p115)\r
+          //     AT THESE PARAMETER VALUES.\r
+               \r
+          //     THE FIRST ELEMENT IN VECTOR *JACOF* MUST BE THE FIRST COMPUTED\r
+          //     JACOBI COEFFICIENT FOR THE RELEVANT PHYSICAL ARC WITH SIGN CHANGES\r
+          //     APPROPRIATE TO THE JACOBI TYPE OF THE ARC.\r
+               \r
+          //     THE FIRST ELEMENTS IN VECTORS *ACOEF* AND *BCOEF* MUST BE THE \r
+          //     FIRST THREE-TERM RECURRENCE COEFFICIENTS FOR THE RELEVANT PHYSICAL\r
+          //     ARC.\r
+               \r
+          //     IER=0  - NORMAL EXIT\r
+          //     IER=36 - AN ELEMENT OF ARRAY *SVAL* IS EITHER +1 OR -1, A\r
+          //              POSSIBILITY NOT ALLOWED BY THE CURRENT CODE.\r
+               \r
+          //     LOCAL VARIABLES\r
+               \r
+          int I;\r
+          double PHI,TT,TUPI;\r
+          double C1[] = new double[2];\r
+          double CT[] = new double[2];\r
+          // COMPLEX C1,CT,DPARFN,PARFUN\r
+          // EXTERNAL DPARFN,JACSUM,PARFUN\r
+          double cr[] = new double[1];\r
+          double ci[] = new double[1];\r
+          double denom;\r
+          double TTIN[] = new double[2];\r
+          double DOUT[];\r
+          double POUT[];\r
+               \r
+          TUPI=2.0*Math.PI;\r
+               \r
+          for (I=1; I <= NVALS; I++) {\r
+              TT=SJT*TVAL[I-1];\r
+                  PHI=JACSUM(TT,DGPOL,ACOEF,BCOEF,H0VAL,JACOF);\r
+                  RHOVL[I-1][0]=TUPI*Math.pow((1.0+TT),BETA)*PHI;\r
+                  RHOVL[I-1][1] = 0.0;\r
+                  // AT THIS POINT RHOVL STORES THE BOUNDARY CORRESPONDENCE\r
+                  // DERIVATIVE.\r
+          } // for (I=1; I <= NVALS; I++)\r
+               \r
+          for (I=1; I <= NVALS; I++) {\r
+              if (1.0+SJT*SVAL[I-1] == 0.0) {\r
+                      IER[0]=36;\r
+                      return;\r
+              }\r
+              else {\r
+                  zdiv(HAANG,0.0,RHOVL[I-1][0],RHOVL[I-1][1],cr,ci);\r
+                  denom = Math.pow((1.0+SJT*SVAL[I-1]), BETAC);\r
+                  RHOVL[I-1][0] = cr[0]/denom;\r
+                  RHOVL[I-1][1] = ci[0]/denom;\r
+              }\r
+          } // for (I=1; I <= NVALS; I++) \r
+               \r
+          C1[0] = 0.0;\r
+          if (INTER) {\r
+                  C1[1] = 1.0/TUPI;\r
+          }\r
+          else {\r
+                  C1[1] = -1.0/TUPI;\r
+          }\r
+               \r
+          for (I=1; I <= NVALS; I++) {\r
+           TT=MIDPT+HALEN*TVAL[I-1];\r
+           TTIN[0] = TT;\r
+           TTIN[1] = 0.0;\r
+           DOUT = DPARFN(PARNT,TTIN);\r
+           zmlt(DOUT[0],DOUT[1],RHOVL[I-1][0]*HALEN,RHOVL[I-1][1]*HALEN,cr,ci);\r
+           CT[0] = cr[0];\r
+           CT[1] = ci[0];\r
+                  POUT = PARFUN(PARNT,TTIN);\r
+                  zdiv(CT[0],CT[1],POUT[0]-CENTR[0],POUT[1]-CENTR[1],cr,ci);\r
+                  CT[0] = cr[0];\r
+                  CT[1] = ci[0];\r
+                  zmlt(CT[0],CT[1],C1[0],C1[1],cr,ci);\r
+                  RHOVL[I-1][0] = cr[0];\r
+                  RHOVL[I-1][1] = ci[0];\r
+          } // for (I=1; I <= NVALS; I++)\r
+       \r
+          // NORMAL EXIT\r
+               \r
+           IER[0]=0;\r
+               \r
+    } // private void RHOFN\r
+\r
+\r
+\r
 \r
       /**\r
        * zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.\r