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

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

index f554a8a..2db7415 100644 (file)
@@ -211,6 +211,47 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 //                     BOUNDARY.\r
        private double FACTR[] = new double[2]; // At first location of ZQUPH\r
        private double ZPPQF[][] = new double[MQUPH][2];\r
+       // From JACANP:\r
+       //            ISNCA  - INTEGER ARRAY\r
+       //                     AN INTEGER VECTOR OF SIZE AT LEAST\r
+       //                        4*IBNDS(1) + 6 ;\r
+       //                     ISNCA MAINLY STORES POINTERS TO ACCESS RSNCA AND\r
+       //                     ZSNCA.\r
+       private int ISNCA[] = new int[6];\r
+       private int DGPOC[] = new int[IBNDS[0]];\r
+    private int JTYPC[] = new int[IBNDS[0]];\r
+    private int LSUBC[] = new int[IBNDS[0]];\r
+    private int PRNSA[] = new int[IBNDS[0]];\r
+    //            RSNCA  - REAL ARRAY\r
+       //                     A REAL VECTOR OF SIZE AT LEAST\r
+       //                         2*IBNDS(1) + (4 + 6*NQPTS)*(NARCS + 1) + 2, \r
+       //                     WHERE NARCS, NQPTS ARE INPUT ARGUMENTS TO JAPHYC.\r
+       //                     (NOTE: NARCS=IGEOM(1), NQPTS=IGEOM(2))\r
+       //                     STORES DATA RELATING TO THREE-TERM RECURRENCE\r
+       //                     SCHEMES, ELEMENTARY GAUSS-JACOBI QUADRATURE RULES,\r
+       //                     AND THE ARGUMENTS OF SUB-ARC ENDPOINTS ON THE UNIT\r
+       //                     DISC.\r
+    private double RSNCA[] = new double[1];\r
+    private double ACOFC[] = new double[TNGQP];\r
+    private double BCOFC[] = new double[TNGQP];\r
+    private double AICOC[] = new double[TNGQP];\r
+    private double BICOC[] = new double[TNGQP];\r
+    private double QUPTC[] = new double[TNGQP];\r
+    private double QUWTC[] = new double[TNGQP];\r
+    private double H0VLC[] = new double[NJIND];\r
+    private double HIVLC[] = new double[NJIND];\r
+    private double JAINC[] = new double[NJIND];\r
+    private double COARG[] = new double[NJIND];\r
+    private double PHPAS[] = new double[IBNDS[0]];\r
+    private double VARGC[] = new double[IBNDS[0]];\r
+       //            CSNCA  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST 2*IBNDS(2) + 1;\r
+       //                     STORES THE JACOBI COEFFICIENTS FOR THE COMPLEX\r
+       //                     (INVERSE) BOUNDARY CORRESPONDENCE FUNCTION AND\r
+       //                     ITS DERIVATIVE.\r
+    private double CSNCA[] = new double[2]; // At first location of CSNCA\r
+    private double BFSNC[][] = new double[IBNDS[1]][2];\r
+    private double SOLNC[][] = new double[IBNDS[1]][2];\r
        public SymmsIntegralMapping() {\r
                \r
        }\r
@@ -9621,9 +9662,15 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           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 CURARG = 0.0;\r
+          double RT1 = 0.0;\r
+          double STARG = 0.0;\r
+          double STRT1 = 0.0;\r
+          double STTH1 = 0.0;\r
+          double THET1 = 0.0;\r
+          double AISUM,ANGLE,ARGBR,ARSUM,BETA,DIST,HL,ISUM,\r
+                    MD,MEAN,MINDS,NEWTL,PI,RR,RRB,RSUM,\r
+                    RT2,SCO,SS,SUM1,THET2,TT,TXI,TUPI,WT;\r
           double BCF[] = new double[2];\r
           double CT[] = new double[2];\r
           double PSI[] = new double[2];\r
@@ -9660,12 +9707,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           double AIN[];\r
           double BIN[];\r
           int I2;\r
+          double PIN[] = new double[2];\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
+          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
@@ -9680,7 +9728,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                       MD=MIDPT[IA-1];\r
                       ARSUM=0.0;\r
                       AISUM=0.0;\r
-                      for (JQ=1; JQ <= NQ; JQ++) {\r
+                      jqloop: 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
@@ -9931,114 +9979,142 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                                        } // 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
-                                   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
+                              if (SS < 0.0) {\r
+                                  LIM=NQUAD;\r
+                                  if ((LIM%2) == 0) {\r
+                                      LIM=LIM/2;\r
+                                  }\r
+                                  else {\r
+                                      LIM=(LIM-1)/2;\r
+                                  }\r
+                                  J1=0;\r
+                                  J2=NQUAD+1;\r
+                                  for (J=1; J <= LIM; J++) {\r
+                                      J1=J1+1;\r
+                                      J2=J2-1;\r
+                                      TT=WSPEC[J1-1];\r
+                                      WSPEC[J1-1]=WSPEC[J2=1];\r
+                                      WSPEC[J2-1]=TT;\r
+                                      TT=TSPEC[J1-1];\r
+                                      TSPEC[J1-1]=TSPEC[J2-1];\r
+                                      TSPEC[J2-1]=TT;\r
+                                      CT[0]=ZSPEC[J1-1][0];\r
+                                      CT[1]=ZSPEC[J1-1][1];\r
+                                      ZSPEC[J1-1][0]=ZSPEC[J2-1][0];\r
+                                      ZSPEC[J1-1][1]=ZSPEC[J2-1][1];\r
+                                      ZSPEC[J2-1][0]=CT[0];\r
+                                      ZSPEC[J2-1][1]=CT[1];\r
+                                  } // for (J=1; J <= LIM; J++) \r
+                              } // if (SS < 0.0)\r
+               \r
+                              // THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS\r
+                              // AND POINTS WSPEC AND ZSPEC.  NOW ESTIMATE THE INTEGRAL.\r
+               \r
+                              ARSUM=0.0;\r
+                              AISUM=0.0;\r
+                              if (IA == 1) {\r
+                                  FIRST=true;\r
+                              }\r
+                              else {\r
+                                  CURARG=STARG;\r
+                                  RT1=STRT1;\r
+                                  DIFF1[0]=STDF1[0];\r
+                                  DIFF1[1]=STDF1[1];\r
+                                  THET1=STTH1;\r
+                              }\r
+                              for (K=1; K <= NQUAD; K++) {\r
+                                  WT=WSPEC[K-1];\r
+                                  DIFF2[0]=ZZ[0]-ZSPEC[K-1][0];\r
+                                  DIFF2[1]=ZZ[1]-ZSPEC[K-1][1];\r
+                                  RT2=TSPEC[K-1];\r
+                                  DIST=zabs(DIFF2[0],DIFF2[1]);\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
+                                  } // if (FIRST)\r
+                                  else {\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
+                                      } // if (ANGLE <= -PI || ANGLE > PI) \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
+                              } // for (K=1; K <= NQUAD; K++)\r
+                              break jqloop;\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
-                         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[0]=DIFF2[0];\r
-                           DIFF1[1]=DIFF2[1];\r
-                           THET1=THET2\r
-                         ENDIF           \r
-                         STARG=CURARG\r
-                         STRT1=RT1\r
-                         STDF1[0]=DIFF1[0];\r
-                         STDF1[1]=DIFF1[1];\r
-                         STTH1=THET1\r
+                      } // jqloop: for (JQ=1; JQ <= NQ; JQ++)\r
+               \r
+                      RSUM=RSUM+ARSUM;\r
+                      ISUM=ISUM+AISUM;\r
+                      if (JT < 0) {\r
+               \r
+                          // BRING THE ARGUMENT FORWARD TO THE CORNER POINT AND REPLACE\r
+                          // THE INCREMENTED CURARG VALUE BY AN INVERSE TANGENT\r
+                          // EVALUATION\r
+               \r
+                          PIN[0] = 1.0;\r
+                          PIN[1] = 0.0;\r
+                          POUT = PARFUN(PT,PIN);\r
+                          DIFF2[0] = ZZ[0] - POUT[0];\r
+                          DIFF2[1] = ZZ[1] - POUT[1];\r
+                          RT2=1.0;\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
+                          } // if (ANGLE <= -PI || ANGLE > PI)\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
+                          ARGBR=(int)Math.round((CURARG-THET2)/TUPI);\r
+                          CURARG=THET2+TUPI*ARGBR;\r
+                          RT1=-1.0;\r
+                          DIFF1[0]=DIFF2[0];\r
+                          DIFF1[1]=DIFF2[1];\r
+                          THET1=THET2;\r
+                      } // if (JT < 0)          \r
+                      STARG=CURARG;\r
+                      STRT1=RT1;\r
+                      STDF1[0]=DIFF1[0];\r
+                      STDF1[1]=DIFF1[1];\r
+                      STTH1=THET1;\r
                \r
                       // END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA\r
                \r
@@ -10353,7 +10429,669 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         IER[0]=0;\r
                \r
    } // private void PPSBI1\r
+   \r
+   private void JACANP(double RWORK[], int IER[]) {\r
+               \r
+       // INTEGER CHNL,IER\r
+          // INTEGER IBNDS(*),ISNCA(*),ISNPH(*),IGEOM(*),IQUPH(*)\r
+          // REAL RSNCA(*),RSNPH(*),RGEOM(*),RQUPH(*),RWORK\r
+          // COMPLEX CENTR\r
+          // COMPLEX ZSNCA(*),ZQUPH(*)\r
+          // LOGICAL INTER\r
+               \r
+               // ......................................................................\r
+               \r
+               // 1.     JACANP\r
+               //           COMPUTATION OF PIECEWISE ORTHOGONAL JACOBI POLYNOMIAL \r
+               //           APPROXIMATIONS TO THE COMPLEX BOUNDARY CORRESPONDENCE DERIV-\r
+               //           ATIVE FOR THE MAP: CANONICAL--> PHYSICAL.\r
+               \r
+               // 2.     PURPOSE\r
+               //           THE MAIN PURPOSE IS TO CALCULATE THE COEFFICIENTS IN THE\r
+               //           PIECEWISE ORTHOGONAL JACOBI POLYNOMIAL APPROXIMATIONS TO THE\r
+               //           COMPLEX BOUNDARY CORRESPONDENCE DERIVATIVE FOR THE CONFORMAL\r
+               //           MAP OF THE CANONICAL DOMAIN ONTO THE SIMPLY-CONNECTED PHYS-\r
+               //           ICAL DOMAIN.  \r
+               //           THE METHOD ATTEMPTS TO COMPUTE TRUNCATED FOURIER-JACOBI APP-\r
+               //           ROXIMATIONS BY DIRECT QUADRATURE ESTIMATION OF THE FOURIER-\r
+               //           JACOBI COEFFICIENTS.  IF DECAY OF THESE COEFFICIENTS ISN'T\r
+               //           SUFFICIENTLY RAPID ON A GIVEN SUB-ARC OF THE DISC, THEN\r
+               //           THE SUB-ARC IS DIVIDED.\r
+               //           A NUMBER OF DATA ARRAYS ASSOCIATED WITH THE POLYNOMIAL\r
+               //           APPROXIMATIONS ARE ALSO COMPUTED AND MAY BE USED FOR SUBSE-\r
+               //           QUENT PROCESSING.  IN ADDITION TO BEING RETURNED AS \r
+               //           PARAMETERS OF THE SUBROUTINE THESE ARE ALSO AUTOMATICALLY\r
+               //           OUTPUT TO DATA FILES.\r
 \r
+               \r
+               // 3.     CALLING SEQUENCE\r
+               //           CALL JACANP(IBNDS,INTER,CENTR,IGEOM,RGEOM,ISNPH,RSNPH,IQUPH,\r
+               //                       RQUPH,ZQUPH,RWORK,CHNL,ISNCA,RSNCA,ZSNCA,IER)\r
+               \r
+               //        PARAMETERS\r
+               //         ON ENTRY\r
+               //            IBNDS  - INTEGER ARRAY\r
+               //                     INTEGER VECTOR OF SIZE AT LEAST 2.\r
+               //                     IBNDS(K), K=1,2 DEFINE VARIOUS UPPER LIMITS\r
+               //                     THAT HAVE BEEN SET IN THE CALLING PROGRAM AND\r
+               //                     WHICH CONTROL THE SIZES OF THE ARRAYS ISNCA,RSNCA,\r
+               //                     ZSNCA.\r
+               //                     IBNDS(1) - THE MAXIMUM NUMBER OF SUB-ARCS ALLOWED\r
+               //                                ON THE UNIT DISC.\r
+               //                     IBNDS(2) - THE MAXIMUM TOTAL NUMBER OF JACOBI CO-\r
+               //                                EFFICIENTS ALLOWED.\r
+               //                                (IBNDS(2) <= IBNDS(1)*NQPTS WHERE\r
+               //                                 NQPTS = IGEOM(2))\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
+               //            RWORK  - REAL ARRAY\r
+               //                     REAL WORKING VECTOR OF SIZE AT LEAST \r
+               //                              (NARCS + 1)*NQPTS \r
+               //                     WHERE NARCS, NQPTS ARE INPUT ARGUMENTS TO JAPHYC.\r
+               //                     (NOTE: NARCS=IGEOM(1), NQPTS=IGEOM(2))\r
+               \r
+               //             CHNL  - INTEGER\r
+               //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+               //                     WRITING THE FILES <JBNM>ca, <JBNM>cl.\r
+               \r
+               //         ON EXIT\r
+               //            ISNCA  - INTEGER ARRAY\r
+               //                     AN INTEGER VECTOR OF SIZE AT LEAST\r
+               //                        4*IBNDS(1) + 6 ;\r
+               //                     ISNCA MAINLY STORES POINTERS TO ACCESS RSNCA AND\r
+               //                     ZSNCA.\r
+               \r
+               //            RSNCA  - REAL ARRAY\r
+               //                     A REAL VECTOR OF SIZE AT LEAST\r
+               //                         2*IBNDS(1) + (4 + 6*NQPTS)*(NARCS + 1) + 2, \r
+               //                     WHERE NARCS, NQPTS ARE INPUT ARGUMENTS TO JAPHYC.\r
+               //                     (NOTE: NARCS=IGEOM(1), NQPTS=IGEOM(2))\r
+               //                     STORES DATA RELATING TO THREE-TERM RECURRENCE\r
+               //                     SCHEMES, ELEMENTARY GAUSS-JACOBI QUADRATURE RULES,\r
+               //                     AND THE ARGUMENTS OF SUB-ARC ENDPOINTS ON THE UNIT\r
+               //                     DISC.\r
+               \r
+               //            CSNCA  - COMPLEX ARRAY\r
+               //                     A COMPLEX VECTOR OF SIZE AT LEAST 2*IBNDS(2) + 1;\r
+               //                     STORES THE JACOBI COEFFICIENTS FOR THE COMPLEX\r
+               //                     (INVERSE) BOUNDARY CORRESPONDENCE FUNCTION AND\r
+               //                     ITS DERIVATIVE.\r
+               \r
+               //            IER    - INTEGER\r
+               //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+               //                     A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
+               //                     WRITTEN ON THE STANDARD OUTPUT CHANNEL AND THE\r
+               //                     LISTING FILE <JBNM>cl.\r
+               //                     IER=0 - NORMAL EXIT.\r
+               //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+               //                             BE SELF EXPLANATORY.\r
+                 \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 SOME INPUT ARGUMENTS FOR JACANP ARE OUTPUT VALUES\r
+               //             FROM JAPHYC AND GQPHYC.\r
+               //           - THE DATA WHICH MAY BE REQUIRED FOR LATER PROCESSING BY\r
+               //             OTHER CONFPACK ROUTINES IS WRITTEN ON THE FIEL <JBNM>ca,\r
+               //             WHERE <JBNM> IS COLLECTED FROM THE FILE jbnm PREVIOUSLY \r
+               //             CREATED BY JAPHYC.\r
+               //           - A SUMMARY LISTING OF ACTIONS TAKEN IS AUTOMATICALLY\r
+               //             WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+               \r
+               // ......................................................................\r
+               //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+               //     LAST UPDATE: 3 JULY 1990\r
+               // ......................................................................\r
+                    \r
+               //     LOCAL VARAIBLES\r
+               //\r
+               \r
+               //**** POINTERS USED TO PROCESS ARRAYS\r
+               \r
+               // INTEGER ACOEF,ACOFC,AICOC,AICOF,BCFSN,BCOEF,BCOFC,BFSNC,BICOC,\r
+               // BICOF,COARG,DGPOC,DGPOL,ERARC,H0VAL,H0VLC,HALEN,HIVAL,HIVLC,JACIN,\r
+               // JAINC,JATYP,JTYPC,LOSUB,LQSBF,LSUBC,MIDPT,NPPQF,PARNT,PHPAS,PRNSA,\r
+               // QUPTC,QUPTS,QUWTC,QUWTS,SOLNC,SOLUN,TPPQF,TRRAD,VARGC,VTARG,WPPQF,\r
+               // ZPPQF\r
+               \r
+               // OTHER SCALAR VARIABLES\r
+               \r
+               int I,MNCOF,MNEQN,MNSUA,MNSUC,NEQNS,NJCOG,PT,TNSUA,TNSUC;\r
+               double COLSC[] = new double[]{-1.0};\r
+               double INNRAD,LGTOL,PI,SUPER,THET0,NEWTL;\r
+               double DCAP0[] = new double[2];\r
+               double FACTR[] = new double[2];\r
+               //COMPLEX DCAP0,FACTR,CINRAD\r
+               \r
+               // EXTERNAL BCFSNG,CINRAD,IGNLVL,JCFIM5,OPQUD1,R1MACH,OUPTCA,OUPTCL,\r
+               // WRHEAD,WRTAIL\r
+               \r
+               // OUTPUT CONFPACK HEADING\r
+               \r
+               WRHEAD(3,0,null);\r
+               \r
+               NEWTL=Math.sqrt(EPS);\r
+               PI=Math.PI;\r
+               \r
+               NARCS=IGEOM[0];\r
+               NQPTS=IGEOM[1];\r
+               MNSUA=IGEOM[3];\r
+               TNSUA=ISNPH[2];\r
+               NEQNS=ISNPH[3];\r
+               MNEQN=ISNPH[5];\r
+               SUPER=RGEOM[0];\r
+               LGTOL=RGEOM[1];\r
+               MNCOF=IBNDS[1];\r
+               MNSUC=IBNDS[0];\r
+               MQUPH=IQUPH[3];\r
+               NJIND=NARCS+1;\r
+               TNGQP=NQPTS*NJIND;\r
+               \r
+               // ASSIGN FIXED DATA TO ISNCA, RSNCA\r
+               \r
+               ISNCA[0]=NARCS;\r
+               ISNCA[1]=NQPTS;\r
+               ISNCA[4]=MNSUC;\r
+               ISNCA[5]=MNCOF;\r
+               RSNCA[0]=LGTOL;\r
+               \r
+               // SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC\r
+               \r
+               // SET UP POINTERS TO ELEMENTS IN ISNPH AND RSNPH,AS IN JAPHYC\r
+               \r
+               // SET UP POINTERS TO ELEMENTS IN IQUPH AND RQUPH, AS IN GQPHYC\r
+               \r
+               // SET UP POINTERS TO ELEMENTS IN ISNCA, RSNCA AND ZSNCA\r
+               \r
+           // INITIALISE JACOBI INDECES *JAINC* FOR THE INVERSE MAP\r
+               \r
+               for (I=1; I <= NARCS; I++) {\r
+                   JAINC[I-1]=-JACIN[I-1]/(1.0+JACIN[I-1]);\r
+               }\r
+               JAINC[NJIND-1]=0.0;\r
+               \r
+               //     SET UP GAUSS-JACOBI QUADRATURE DATA FOR INVERSE MAP AND STORE IN\r
+               //     ARRAYS *QUPTC* AND *QUWTC*.  SET THE CORRESPONDING THREE TERM\r
+               //     RECURRENCE COEFFICIENTS AND STORE IN *ACOFC*, *BCOFC*.  DETERMINE\r
+               //     THE ZEROTH MOMENTS OF THE JACOBI DISTRIBUTIONS AND STORE IN \r
+               //     *H0VLC*. ALSO SET UP THE DATA *AICOC*,*BICOC* AND *HICOC* FOR THE\r
+               //     INTEGRATED POLYNOMIALS NEEDED FOR PROCESSING AFTER THIS MODULE.\r
+               \r
+               OPQUD1(NJIND,NQPTS,JAINC,ACOFC,BCOFC,\r
+                    H0VLC,AICOC,BICOC,HIVLC,QUPTC,\r
+                    QUWTC,RWORK,IER);\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(3,0,IER[0],null);\r
+                       return;\r
+               }\r
+               System.out.println("BASIC GAUSS QUADRATURE DATA DONE:");\r
+               \r
+               // SET UP THE ARRAY *RWORK* OF REFERENCE IGNORE LEVELS.\r
+               \r
+               IGNLVL(RWORK,COLSC,ACOFC,BCOFC,H0VLC,\r
+                    JAINC,NJIND,NQPTS,IER);\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(3,0,IER[0],null);\r
+                       return;\r
+               }\r
+               \r
+               //     INITIALISE THE DATA RELATING TO THE DESCRIPTION OF THE DISCRET-\r
+               //     ISATION AND SOLUTION HOUSEKEEPING ON THE UNIT CIRCLE\r
+               \r
+               TNSUC=TNSUA;\r
+                    \r
+               for (I=1; I <= TNSUA; I++) {\r
+                   DGPOC[I-1]=NQPTS-1;\r
+                   JTYPC[I-1]=JATYP[I-1];\r
+                   PHPAS[I-1]=-1.0;\r
+                   PRNSA[I-1]=I;\r
+                   VARGC[I-1]=VTARG[I-1];\r
+               } // for (I=1; I <= TNSUA; I++)\r
+               VARGC[TNSUA]=VTARG[TNSUA];\r
+               PHPAS[TNSUA]=-1.0;\r
+               COARG[0]=VARGC[0];\r
+               PT=1;\r
+               for (I=2; I <= TNSUA; I++) {\r
+                   if (PARNT[I-1] != PT) {\r
+                       PT=PARNT[I-1];\r
+                       COARG[PT-1]=VARGC[I-1];\r
+                   }\r
+               } // for (I=2; I <= TNSUA; I++)\r
+               COARG[NARCS]=COARG[0]+2.0*PI;\r
+               \r
+               // SET UP THE JACOBI COEFFICIENTS OF THE COMPLEX DENSITY FOR\r
+               // THE INVERSE MAP\r
+               \r
+               JCFIM5(DGPOC,IER,JTYPC,LSUBC,\r
+                    PHPAS,PRNSA,SOLNC,TNSUC,VARGC,\r
+                    AICOF,ACOEF,ACOFC,BICOF,BCFSN,\r
+                    BCOEF,BCOFC,CENTR,DGPOL,ERARC,\r
+                    H0VAL,H0VLC,HIVAL,HALEN,INTER,\r
+                    JACIN,JAINC,JATYP,LGTOL,LOSUB,\r
+                    MIDPT,MNCOF,MNSUC,NJIND,NQPTS,PARNT,QUPTC,\r
+                    QUWTC,RWORK,SOLUN,NEWTL,VTARG);\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(3,0,IER[0],null);\r
+                       return;\r
+               }\r
+               NJCOG=LSUBC[TNSUC-1]+DGPOC[TNSUC-1];\r
+               ISNCA[2]=TNSUC;\r
+               ISNCA[3]=NJCOG;\r
+               System.out.println("JACOBI COEFFICIENTS DONE:");\r
+               \r
+               // COMPUTE THE COMPLEX INNER RADIUS FOR INTERIOR DOMAINS OR THE\r
+               // COMPLEX CAPACITY FOR EXTERIOR DOMAINS,  STORING RELEVANT VALUE\r
+               // IN *DCAP0*\r
+               \r
+               THET0=VTARG[0];\r
+               FACTR[0] = Math.cos(THET0);\r
+               FACTR[1] = Math.sin(THET0);\r
+               /*      DCAP0=CINRAD(NARCS,NQPTS,TNSUA,ISNPH(DGPOL),ISNPH(JATYP),\r
+                    +ISNPH(LOSUB),IQUPH(LQSBF),IQUPH(NPPQF),IGEOM(PARNT),RSNPH(ACOEF),\r
+                    +RSNPH(BCOEF),RSNPH(H0VAL),RGEOM(HALEN),RSNPH(JACIN),LGTOL,\r
+                    +RGEOM(MIDPT),RSNPH(QUPTS),RSNPH(QUWTS),RSNPH(SOLUN),RQUPH(TPPQF),\r
+                    +RQUPH(TRRAD),RQUPH(WPPQF),CENTR,FACTR,ZQUPH(ZPPQF),IER)\r
+                     IF (IER .GT. 0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+               C\r
+                     IF (INTER) THEN\r
+                       INNRAD=ABS(DCAP0)\r
+                       WRITE(*,2) 'INNER RADIUS:',INNRAD\r
+                     ELSE\r
+                       DCAP0=EXP(-RSNPH(SOLUN+NEQNS-1))*DCAP0/ABS(DCAP0)\r
+                     ENDIF\r
+                     ZSNCA(1)=DCAP0\r
+               C\r
+               C     GET THE COEFFICIENTS BFSNC FOR THE COMPLEX BOUNDARY CORRESPONDENCE\r
+               C     FUNCTION FOR THE INVERSE MAP.\r
+               C\r
+                     CALL BCFSNG(TNSUC,ISNCA(DGPOC),ISNCA(JTYPC),ISNCA(LSUBC),\r
+                    +RSNCA(H0VLC),RSNCA(JAINC),ZSNCA(BFSNC),ZSNCA(SOLNC))\r
+               C\r
+                     CALL OUPTCL(ISNCA(DGPOC),ISNCA(JTYPC),LGTOL,ISNCA(LSUBC),NQPTS,\r
+                    +CHNL,IGEOM(PARNT),ISNCA(PRNSA),RWORK,ZSNCA(SOLNC),TNSUC,INTER,\r
+                    +INNRAD,IER)\r
+               C\r
+               C     OUTPUT ALL RESULTS REQUIRED FOR SUBSEQUENT PROCESSING\r
+               C\r
+                     CALL OUPTCA(ISNCA,RSNCA,ZSNCA,CHNL)\r
+               C\r
+               999   CONTINUE\r
+               C\r
+                     CALL WRTAIL(3,0,IER)\r
+               C */\r
+    } // private void JACANP\r
+\r
+\r
+   private void JCFIM5(int DGPOC[], int IER[], int JTYPC[], int LSUBC[], double PHPAS[],\r
+       int PRNSA[], double SOLNC[][], int TNSUC, double VARGC[], double AICOF[], double  ACOEF[],\r
+       double ACOFC[], double BICOF[], double BCFSN[], double BCOEF[], double BCOFC[],\r
+       double CENTR[], int DGPOL[], double ERARC[], double H0VAL[], double H0VLC[], double HIVAL[],\r
+       double HALEN[], boolean INTER, double JACIN[], double JAINC[], int JATYP[], double LGTOL,\r
+       int LOSUB[], double MIDPT[], int MNCOF, int MNSUC, int NJIND,int NQPTS, int PARNT[],\r
+       double QUPTC[], double QUWTC[], double RIGLL[], double SOLUN[], double NEWTL, \r
+       double VTARG[]) {\r
+               \r
+       // INTEGER DGPOC(*),DGPOL(*),IER,JATYP(*),JTYPC(*),LOSUB(*),LSUBC(*),\r
+          // +MNCOF,MNSUC,NJIND,NQPTS,PARNT(*),PRNSA(*),TNSUC\r
+          // REAL AICOF(*),ACOEF(*),ACOFC(*),BICOF(*),BCFSN(*),BCOEF(*),\r
+          // +BCOFC(*),ERARC(*),H0VAL(*),H0VLC(*),HIVAL(*),HALEN(*),JACIN(*),\r
+          // +JAINC(*),LGTOL,MIDPT(*),QUPTC(*),QUWTC(*),PHPAS(*),RIGLL(*),\r
+          // +SOLUN(*),NEWTL,VARGC(*),VTARG(*)\r
+          // COMPLEX CENTR,SOLNC(*)\r
+          // LOGICAL INTER\r
+               \r
+               //     TO CARRY OUT A DYNAMIC ESTIMATION OF THE JACOBI COEFFICIENTS OF\r
+               //     THE INVERSE COMPLEX DENSITY FUNCTIONS *RHO*; SEE #50 p115 et seq.\r
+               \r
+               //     IER=0  - NORMAL EXIT\r
+               //     IER=30 - LOCAL PARAMETER *MNDG* BELOW NEEDS INCREASING \r
+               //     IER=31 - LOCAL PARAMETER *MNQD* BELOW NEEDS INCREASING \r
+               //     IER=32 - THE SUBROUTINE PARAMETER *TNSUC* HAS REACHED ITS MAXIMUM\r
+               //              PERMITTED VALUE *MNSUC*; *MNSUC* SHOULD BE INCREASED IN \r
+               //              CALLING PROGRAM\r
+               //     IER=33 - THE REQUIRED TOTAL NUMBER OF JACOBI COEFFICIENTS FOR THE\r
+               //              INVERSE BOUNDARY DENSITY EXCEEDS THE LIMIT *MNCOF*; \r
+               //              *MNCOF* SHOULD BE INCREASED IN THE CALLING PROGRAM.\r
+               \r
+               //     LOCAL VARIABLES\r
+               \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
+          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 JACOF[] = new double[MNDG];\r
+          double QAB[] = new double[MNQD];\r
+          double QWT[] = new double[MNQD];\r
+          double SVAL[] = new double[MNQD];\r
+          double TVAL[] = new double[MNQD];\r
+          double WORK[] = new double[MNDG];\r
+          double NEW[][] = new double[MNDG][2];\r
+          double OLD[][] = new double[MNDG][2];\r
+          double RHOVL[][] = new double[MNQD][2];\r
+          // COMPLEX NEW(MNDG),OLD(MNDG),RHOVL(MNQD)\r
+          // EXTERNAL BISNEW,INVJCO,R1MACH\r
+               \r
+          TOLIW=1.0*EPS;\r
+          LOL=(NJIND-1)*NQPTS+1;\r
+          IC=1;\r
+          LOS=1;\r
+               \r
+          /*while (true) {\r
+               \r
+              if (IC > TNSUC) {\r
+               \r
+                      // NORMAL EXIT\r
+               \r
+                      IER[0]=0;\r
+                      return;\r
+              } // if (IC > TNSUC)\r
+               \r
+                  // INITIALISATION FOR PARENT PHYSICAL SUBARC\r
+               \r
+                  PSA=PRNSA[IC-1];\r
+                  DG=DGPOL[PSA-1];\r
+                  JT=JATYP[PSA-1];\r
+                  AJT=Math.abs(JT);\r
+                  if (JT >= 0) {\r
+                          SJT = 1.0;\r
+                  }\r
+                  else {\r
+                          SJT = -1.0;\r
+                  }\r
+                  LOD=(AJT-1)*NQPTS+1;\r
+                  BETA=JACIN[AJT-1];\r
+                  H0=H0VAL[AJT-1];\r
+                  H1=HIVAL[AJT-1];\r
+                  LOM=LOSUB[PSA-1];\r
+                  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
 \r
       /**\r
        * zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.\r