Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 26 Dec 2017 22:21:23 +0000 (22:21 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 26 Dec 2017 22:21:23 +0000 (22:21 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15324 ba61647d-9d00-f842-95cd-605cb4296b96

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

index b46bbe1..c1b089c 100644 (file)
@@ -1888,9 +1888,11 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     // ......................................................................C\r
     //     LOCAL VARIABLES\r
 \r
-    int I,IA,L;\r
+    int I,IA;\r
     int IMX = 0;\r
-    double A1,DIFF,ERR,HH,MINC,R1MACH,RMAX,RMEAN,RMIN,T,TINC,TOL1,TMX,TSD;\r
+    double TINC = 0.0;\r
+    double TMX = 0.0;\r
+    double A1,DIFF,ERR,HH,MINC,RMAX,RMEAN,RMIN,T,TOL1,TSD;\r
     double TT[] = new double[2];\r
     //REAL TT(2)\r
     double C1[] = new double[2];\r
@@ -1902,7 +1904,6 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     //COMPLEX C1,C2,CENTR,ZZ0,DZZ,NDZZ;\r
     double ZZ[][] = new double[2][2];\r
     //COMPLEX ZZ(2)\r
-    String OFL;\r
     //CHARACTER OFL*6\r
     final int MNARC = 200;\r
     final double DTOL = 1.0E-1;\r
@@ -1912,6 +1913,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     double PIN[] = new double[2];\r
     double PAROUT[];\r
     int zzindex;\r
+    double cr[] = new double[1];\r
+    double ci[] = new double[1];\r
     //EXTERNAL DPARFN,LINSEG,PARFUN,R1MACH,WRHEAD,WRTAIL,ZDPARF\r
     //COMPLEX PARFUN,DPARFN,ZDPARF\r
 \r
@@ -2020,7 +2023,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 \r
    MXDIF=0.0;\r
    zzindex = 0;\r
-   /*for (IA=1; IA <= NARCS; IA++) {\r
+   for (IA=1; IA <= NARCS; IA++) {\r
        TT[0]=-1.0;\r
        PIN[0] = TT[0];\r
        PIN[1] = 0.0;\r
@@ -2033,99 +2036,119 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        }\r
        FIRST = true;\r
        WARND= false;\r
-20      CONTINUE\r
-C\r
-C****   TEST THE COMPATIBILTY OF PARFUN AND DPARFN BY ESTIMATING DPARFN\r
-C****   NUMERICALLY AT BOTH REAL AND COMPLEX PARAMETER VALUES.\r
-C\r
-     DO 30 I=1,2\r
-       IF (I.EQ.1) THEN\r
-         C1=CMPLX(TT(1))\r
-       ELSE\r
-         C1=CMPLX(TT(1),MINC)\r
-       ENDIF\r
-       DZZ=DPARFN(IA,C1)\r
-       NDZZ=ZDPARF(IA,C1)\r
-       A1=ABS(DZZ)\r
-C\r
-       IF (A1.EQ.0E+0) THEN\r
-         IER=60\r
-         WRITE(*,*)\r
-         WRITE(*,*) '              ***DPARFN=(0.,0.)***'\r
-         WRITE(*,*) '                             ARC:',IA \r
-         WRITE(*,*) '    STANDARDISED PARAMETER VALUE:',TT(1) \r
-         GOTO 999\r
-       ENDIF\r
-C\r
-       IF (A1.LE.TOL1 .AND. .NOT.WARND) THEN\r
-         WRITE(*,4) '*** W A R N I N G  ***'\r
-         WRITE(*,2) 'PATHOLOGICALLY SMALL DERIVATIVE ON ARC',IA\r
-         WARND=.TRUE.\r
-       ENDIF\r
-C\r
-       IF (FIRST) THEN\r
-         TINC=RMEAN/A1\r
-         TINC=MAX(TINC,MINC)\r
-         FIRST=.FALSE.\r
-       ENDIF\r
-C\r
-       ERR=ABS(1E+0-NDZZ/DZZ)\r
-       IF (ERR.GT.MXDIF) THEN\r
-         MXDIF=ERR\r
-         IMX=IA\r
-         TMX=TT(1)\r
-       ENDIF\r
-30      CONTINUE\r
-C\r
-     IF (.NOT.LNSEG(IA)) THEN\r
-C\r
-C****     DETERMINE THE NEXT BOUNDARY POINT TO BE PLOTTED\r
-C\r
-40        CONTINUE\r
-       TT(2)=TT(1)+TINC\r
-       IF (TT(2) .GE. 1E+0) THEN\r
-         TT(2)=1E+0\r
-         ATEND=.TRUE.\r
-       ELSE\r
-         ATEND=.FALSE.\r
-       ENDIF\r
-C\r
-       ZZ(2)=PARFUN(IA,CMPLX(TT(2)))\r
-       DIFF=ABS(ZZ(2)-ZZ(1))\r
-       IF (DIFF.EQ.0E+0 .AND. .NOT.ATEND) THEN\r
-         TINC=MAX(MINC,2*TINC)\r
-         GOTO 40\r
-       ENDIF\r
-C\r
-       IF (DIFF.GT.RMAX .OR. (DIFF.LT.RMIN .AND. .NOT.ATEND)) THEN\r
-         TINC=RMEAN*TINC/DIFF\r
-         TINC=MAX(TINC,MINC)\r
-         GOTO 40\r
-       ENDIF\r
-C\r
-       WRITE(CHNL,'(2E16.7)') ZZ(2)\r
-       IF (.NOT. ATEND) THEN\r
-         ZZ(1)=ZZ(2)\r
-         TT(1)=TT(2)\r
-         GOTO 20\r
-       ENDIF\r
-     ENDIF\r
+       while (true) {\r
+          \r
+           //****   TEST THE COMPATIBILTY OF PARFUN AND DPARFN BY ESTIMATING DPARFN\r
+           //****   NUMERICALLY AT BOTH REAL AND COMPLEX PARAMETER VALUES.\r
+\r
+                for (I=1; I <= 2; I++) {\r
+                    if (I == 1) {\r
+                        C1[0]=TT[0];\r
+                        C1[1] = 0.0;\r
+                    }\r
+                    else {\r
+                        C1[0]=TT[0];\r
+                     C1[1] = MINC;\r
+                    }\r
+                    DZZ=DPARFN(IA,C1);\r
+                    NDZZ=ZDPARF(IA,C1);\r
+                    A1=zabs(DZZ[0],DZZ[1]);\r
+       \r
+                    if (A1 == 0.0) {\r
+                        IER[0]=60;\r
+                        System.out.println();\r
+                        System.out.println("              ***DPARFN=(0.,0.)***");\r
+                        System.out.println("                             ARC: " + IA); \r
+                        System.out.println(" STANDARDISED PARAMETER VALUE: " + TT[0]);\r
+                        WRTAIL(7,0,IER[0],null);\r
+                        return;\r
+                    } // if (A1 == 0.0)\r
+       \r
+                    if (A1 <=TOL1 && !WARND) {\r
+                        System.out.println("*** W A R N I N G  ***");\r
+                        System.out.println("PATHOLOGICALLY SMALL DERIVATIVE ON ARC" + IA);\r
+                        WARND=true;\r
+                    } // if (A1 <=TOL1 && !WARND)\r
+       \r
+                    if (FIRST) {\r
+                        TINC=RMEAN/A1;\r
+                        TINC=Math.max(TINC,MINC);\r
+                        FIRST=false;\r
+                    } // if (FIRST)\r
+       \r
+                    zdiv(1.0-NDZZ[0],-NDZZ[1],DZZ[0],DZZ[1],cr,ci);\r
+                    ERR = zabs(cr[0],ci[0]);\r
+                    if (ERR > MXDIF) {\r
+                        MXDIF=ERR;\r
+                        IMX=IA;\r
+                        TMX=TT[0];\r
+                    } // if (ERR > MXDIF)\r
+                } // for (I=1; I <= 2; I++)\r
+       \r
+            if (!LNSEG[IA-1]) {\r
+       \r
+                //DETERMINE THE NEXT BOUNDARY POINT TO BE PLOTTED\r
+                while (true) {\r
+                    TT[1]=TT[0]+TINC;\r
+                    if (TT[1] >= 1.0) {\r
+                        TT[1]=1.0;\r
+                        ATEND=true;\r
+                    }\r
+                    else {\r
+                        ATEND=false;\r
+                    }\r
+              \r
+                    PIN[0] = TT[1];\r
+                    PIN[1] = 0.0;\r
+                    ZZ[1] = PARFUN(IA,PIN);\r
+                    DIFF=zabs(ZZ[1][0]-ZZ[0][0],ZZ[1][1]-ZZ[0][1]);\r
+                    if (DIFF == 0.0 && !ATEND) {\r
+                        TINC=Math.max(MINC,2*TINC);\r
+                        continue;\r
+                    } // if (DIFF == 0.0 && !ATEND)\r
+       \r
+                    if (DIFF > RMAX || (DIFF < RMIN && !ATEND)) {\r
+                        TINC=RMEAN*TINC/DIFF;\r
+                        TINC=Math.max(TINC,MINC);\r
+                        continue;\r
+                    } // if (DIFF > RMAX || (DIFF < RMIN && !ATEND))\r
+                    break;\r
+                } // while (true)\r
+       \r
+                zzset[zzindex][0] = ZZ[1][0];\r
+                zzset[zzindex++][1] = ZZ[1][1];\r
+                if (!ATEND) {\r
+                    ZZ[0][0]=ZZ[1][0];\r
+                    ZZ[0][1] = ZZ[1][1];\r
+                    TT[0]=TT[1];\r
+              }\r
+              else {\r
+                  break;\r
+              }\r
+            } // if (!LNSEG[IA-1])\r
+            else {\r
+                break;\r
+            }\r
+       } // while (true)\r
 \r
    } // for (IA=1; IA <= NARCS; IA++)\r
-   IF (LNSEG(NARCS)) WRITE(CHNL,'(2E16.7)') ZZ0\r
-   CLOSE(CHNL)   \r
-C\r
-   IF (MXDIF .GT. DTOL) THEN\r
-     WRITE(*,*)\r
-     WRITE(*,2) 'POSSIBLE PARFUN/DPARFN INCONSISTECY ON ARC:',IMX\r
-     WRITE(*,3) 'OCCURS AT STANDARDISED PARAMETER VALUE:',TMX\r
-     WRITE(*,3) 'RELATIVE FINITE DIFF ERROR:',MXDIF \r
-   ELSE\r
-     WRITE(*,*)\r
-     WRITE(*,1) 'PARFUN AND DPARFN ARE CONSISTENT:'\r
-   ENDIF\r
-C\r
-999   CALL WRTAIL(7,0,IER)*/\r
+   if (LNSEG[NARCS-1]) {\r
+          zzset[zzindex][0] = ZZ0[0];\r
+       zzset[zzindex++][1] = ZZ0[1];\r
+   }   \r
+\r
+   if (MXDIF > DTOL) {\r
+       System.out.println();\r
+       System.out.println("POSSIBLE PARFUN/DPARFN INCONSISTECY ON ARC: " + IMX);\r
+       System.out.println("OCCURS AT STANDARDISED PARAMETER VALUE: " + TMX);\r
+       System.out.println("RELATIVE FINITE DIFF ERROR: " + MXDIF);\r
+   }\r
+   else {\r
+       System.out.println();\r
+       System.out.println("PARFUN AND DPARFN ARE CONSISTENT:");\r
+   }\r
+\r
+   WRTAIL(7,0,IER[0],null);\r
 \r
    } // private void TSTPLOT\r
    \r
@@ -2138,7 +2161,7 @@ C
 \r
        final int NPTS = 9;\r
           int IA,J,NINTS;\r
-       double DIFF,HH,MXDIF,TOL,R1MACH;\r
+       double DIFF,HH,MXDIF,TOL;\r
        double SUM[] = new double[2];\r
        double TT[] = new double[2];\r
        // COMPLEX SUM,TT\r
        } // for (IA=1; IA <= NARCS; IA++)\r
 \r
    } // private void LINSEG\r
+   \r
 \r
+   //COMPLEX FUNCTION ZDPARF(I,T)\r
+   private double[] ZDPARF(int I, double T[]) {\r
+       //COMPLEX T\r
 \r
-      \r
+       //**** NUMERICAL ESTIMATION OF THE DERIVATIVE OF THE PARAMETRIC FUNCTION\r
+       //**** USING 2- OR 4-POINT TRAPEZOIDAL RULE ESTIMATES IN CAUCHY'S\r
+       //**** FORMULA.  THE 2-POINT ESTIMATE IS THE STANDARD CENTRAL DIFFERENCE\r
+       //**** IN THE REAL AXIS DIRECTION.\r
+\r
+       double EPSZ;\r
+       //final double IM[] = new double[]{0.0,1.0};\r
+       final boolean FOUR = false;\r
+       double SUM[] = new double[2];\r
+       //COMPLEX IM,SUM\r
+       double POUT1[];\r
+       double POUT2[];\r
+       double result[] = new double[2];\r
+       double PIN[] = new double[2];\r
+\r
+       //EXTERNAL PARFUN\r
+       // COMPLEX PARFUN\r
+\r
+       EPSZ = Math.pow(EPS, 0.3333);\r
+       PIN[0] = T[0] + EPSZ;\r
+       PIN[1] = T[1];\r
+       POUT1 = PARFUN(I,PIN);\r
+       PIN[0] = T[0] - EPSZ;\r
+       PIN[1] = T[1];\r
+       POUT2 = PARFUN(I,PIN);\r
+       SUM[0] = (POUT1[0] - POUT2[0])/2.0/EPSZ;\r
+       SUM[1] = (POUT1[1] - POUT2[1])/2.0/EPSZ;\r
+\r
+       if (FOUR) {\r
+           PIN[0] = T[0];\r
+           PIN[1] = T[1]+EPSZ;\r
+           POUT1 = PARFUN(I,PIN);\r
+           PIN[0] = T[0];\r
+           PIN[1] = T[1]-EPSZ;\r
+           POUT2 = PARFUN(I,PIN);\r
+            result[0] = SUM[0]/2.0 + (POUT1[1] - POUT2[1])/4.0/EPSZ;\r
+            result[1] = SUM[1]/2.0 - (POUT1[0] - POUT2[0])/4.0/EPSZ;\r
+       }\r
+       else {\r
+            result[0]=SUM[0];\r
+            result[1] = SUM[1];\r
+       }\r
+       return result;\r
+   }\r
+\r
+   private void JAPHYC(String JBNM, String HEAD, double MAXER,boolean INTER, int NARCS,\r
+                     int ISYGP, int NQPTS, boolean INCST,\r
+                     int RFARC, double RFARG[], double CENTR[], int TSTNG,int OULVL, int IBNDS[],\r
+                     int MNEQN, double MATRX[][][], int IWORK[], double RWORK[],\r
+                     double ZWORK[][], boolean LWORK[], int OCH, int IGEOM[], double RGEOM[],\r
+                     int ISNPH[], double RSNPH[], int IER[]) {\r
+               \r
+                     //INTEGER IBNDS(*),IGEOM(*),ISNPH(*),IWORK(*)\r
+                     //REAL RGEOM(*),MATRX(MNEQN,MNEQN,2),RSNPH(*),RWORK(*)\r
+                     //COMPLEX CENTR\r
+                     //COMPLEX ZWORK(*)\r
+                     //LOGICAL LWORK(*)\r
+                     //CHARACTER JBNM*4,HEAD*72\r
+               \r
+               // ......................................................................\r
+               \r
+               // 1.     JAPHYC\r
+               //           COMPUTATION OF PIECEWISE ORTHOGONAL JACOBI POLYNOMIAL \r
+               //           APPROXIMATIONS TO THE BOUNDARY CORRESPONDENCE DERIVATIVE FOR\r
+               //           THE MAP:PHYSICAL --> CANONICAL.\r
+               \r
+               // 2.     PURPOSE\r
+               //           THE MAIN PURPOSE IS TO CALCULATE THE COEFFICIENTS IN THE\r
+               //           PIECEWISE ORTHOGONAL JACOBI POLYNOMIAL APPROXIMATIONS TO THE\r
+               //           BOUNDARY CORRESPONDENCE DERIVATIVE FOR THE CONFORMAL MAP OF\r
+               //           A GIVEN SIMPLY CONNECTED PHYSICAL DOMAIN (WITH PIECEWISE\r
+               //           ANALYTIC BOUNDARY) ONTO A CANONICAL DOMAIN (WITH UNIT CIRCLE\r
+               //           AS BOUNDARY).  AN INTERIOR PHYSICAL DOMAIN IS MAPPED TO THE \r
+               //           UNIT DISC, AN EXTERIOR PHYSICAL DOMAIN TO THE COMPLEMENT OF\r
+               //           THE CLOSED UNIT DISC.\r
+               //           THE METHOD USED IS AN ADAPTIVE COLLOCATION SOLUTION OF\r
+               //           SYMM'S INTEGRAL EQUATION.  \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
+               // 3.     CALLING SEQUENCE\r
+               //           CALL JAPHYC(JBNM,HEAD,MAXER,INTER,NARCS,ISYGP,NQPTS,INCST,\r
+               //                       RFARC,RFARG,CENTR,TSTNG,OULVL,IBNDS,MATRX,IWORK,\r
+               //                       RWORK,ZWORK,LWORK,OCH,IGEOM,RGEOM,ISNPH,RSNPH,\r
+               //                       IER)\r
+               \r
+               //        PARAMETERS\r
+               //         ON ENTRY\r
+               //            JBNM   - CHARACTER*4\r
+               //                     THE JOB NAME.  THIS IS USED TO CREATE THREE OUT-\r
+               //                     PUT FILES WITH FILENAMES\r
+               \r
+               //                         <JBNM>pl, <JBNM>gm, <JBNM>ph,\r
+               \r
+               //                     WHERE <JBNM> DENOTES THE VALUE OF VARIABLE JBNM\r
+               //                     WITH ANY TRAILING SPACES DELETED.  THE FIRST OF\r
+               //                     THESE IS A LISTING FILE RECORDING THE PROGRESS\r
+               //                     AND RESULTS OF THE CALCULATION FOR LATER READING\r
+               //                     BY THE USER.  THE TWO FILES <JBNM>gm AND <JBNM>ph\r
+               //                     ARE DATA FILES, NOT REALLY INTENDED TO BE READ\r
+               //                     BY THE USER.\r
+               //                     THE VALUE OF JBNM IS ALSO THE ONLY ITEM IN A\r
+               //                     FOURTH OUTPUT FILE NAMED (LITERALLY) jbnm.\r
+               \r
+               //            HEAD   - CHARACTER*72\r
+               //                     A HEADING FOR THE PROBLEM, TO APPEAR ON THE\r
+               //                     LISTING FILE <JBNM>pl.\r
+               \r
+               //            MAXER  - REAL\r
+               //                     RELATIVE ACCURACY REQUESTED FOR THE CONFORMAL MAP;\r
+               //                     THIS IS THE SAME AS THE ABSOLUTE ACCURACY ON THE\r
+               //                     BOUNDARY OF THE PHYSICAL DOMAIN.\r
+               \r
+               //            INTER  - LOGICAL\r
+               //                     TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE \r
+               //                     OTHERWISE.\r
+               \r
+               //            NARCS  - INTEGER\r
+               //                     THE NUMBER OF ANALYTIC ARCS THAT MAKE UP THE \r
+               //                     W H O L E BOUNDARY OF THE PHYSICAL DOMAIN.\r
+               \r
+               //            ISYGP  - INTEGER\r
+               //                     THE MAGNITUDE OF ISYGP IS THE ORDER OF THE\r
+               //                     SYMMETRY GROUP OF THE PHYSICAL DOMAIN.\r
+               //                     ISYGP.EQ.1  -THE SYMMETRY GROUP HAS ONLY ONE ELE-\r
+               //                                  MENT,THE IDENTITY TRANSFORMATION;  IN\r
+               //                                  OTHER WORDS, THE DOMAIN HAS 'NO \r
+               //                                  SYMMETRY'.\r
+               //                     ISYGP.GT.1  -THE SYMMETRY GROUP CONTAINS ONLY\r
+               //                                  PROPER (IN-PLANE) ROTATIONS; IN OTHER\r
+               //                                  WORDS, THE DOMAIN HAS ONLY ROTATIONAL\r
+               //                                  SYMMETRIES.\r
+               //                     ISYGP.LT.-1 -THE SYMMETRY GROUP CONTAINS IMPROPER\r
+               //                                  (OUT-OF-PLANE) ROATIONS; IN OTHER\r
+               //                                  WORDS, THE DOMAIN HAS REFLECTIONAL\r
+               //                                  SYMMETRY AND MAY ALSO HAVE ROTATIONAL\r
+               //                                  SYMMETRIES.\r
+               //                     AN INPUT VALUE OF -1 OR 0 IS TREATED AS IF IT WERE\r
+               //                     1.\r
+               \r
+               //            NQPTS  - INTEGER\r
+               //                     PLAYS TWO ROLES.\r
+               //                     1. THE NUMBER OF QUADRATURE POINTS TO BE USED IN \r
+               //                        AN ELEMENTARY GAUSS-JACOBI RULE;  COMPOSITE \r
+               //                        RULES ARE CONSTRUCTED FROM PANELS OF NQPTS-\r
+               //                        POINT RULES.\r
+               //                     2. THE MAXIMUM DEGREE OF POLYNOMIAL APPROXIMATION\r
+               //                        IS FIXED AT NQPTS-1.\r
+               //                     NQPTS SHOULD BE REASONABLY LARGE; A PRACTICAL RULE\r
+               //                     OF THUMB IS THAT IF MACHINE PRECISION IS X*1E-N,\r
+               //                    1<X<10, THEN NQPTS=N+1.\r
+               \r
+               //            INCST  - LOGICAL\r
+               //                     IF INCST IS TRUE THEN AN INCREMENTAL STRATEGY IS\r
+               //                     USED TO TRY TO ACHIEVE THE ACCURACY SPECIFIED BY\r
+               //                     MAXER; VERY ROUGHLY SPEAKING, THIS MEANS THAT THE \r
+               //                     METHOD SUCCESSIVELY ACHIEVES THE TARGET ACCURACIES\r
+               //                     1E-1,1E-2,...UNTIL MAXER HAS BEEN ACHIEVED.  IF \r
+               //                     THE PROBLEM IS THOUGHT TO BE EITHER PARTICULARLY\r
+               //                     DIFFICULT OR PARTICULARLY SIMPLE, THEN INCST \r
+               //                     SHOULD BE SET TO .TRUE.  FOR PROBLEMS OF 'AVERAGE'\r
+               //                     DIFFICULTY, SETTING INCST TO .FALSE. IS USUALLY \r
+               //                     MORE EFFICIENT.\r
+               \r
+               //            RFARC  - INTEGER\r
+               //                     THE REFERENCE ARC USED TO DEFINE THE ORIENTATION \r
+               //                     THAT IS GIVEN TO THE MAP.  THE CONVENTION IS THAT \r
+               //                     THE POINT AT THE START OF ANALYTIC ARC NUMBER \r
+               //                     RFARC IS MAPPED TO THE POINT WITH ARGUMENT \r
+               //                     RFARG*PI ON THE UNIT DISC.\r
+               \r
+               //            RFARG  - REAL\r
+               //                     THE REFERENCE ARGUMENT/PI USED TO DEFINE THE \r
+               //                     ORIENTATION THAT IS GIVEN TO THE MAP.  SEE RFARC \r
+               //                     ABOVE.\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.\r
+               //                     IN CASE ABS(ISYGP).GT.1 THEN CENTR MUST ALSO BE\r
+               //                     A CENTRE OF SYMMETRY FOR THE PHYSICAL DOMAIN.\r
+               \r
+               //            TSTNG  - INTEGER\r
+               //                     EITHER 0 OR 1.\r
+               //                     ON SUCCESSFUL COMPLETION OF THE NUMERICAL SOLUTION\r
+               //                     OF SYMM'S EQUATION, A MODULE IS PROVIDED FOR\r
+               //                     TESTING THE ERROR IN THE MODULUS OF THE COMPUTED\r
+               //                     MAP ON THE BOUNDARY OF THE DOMAIN.\r
+               //                     TSTNG=0 - TEST ONLY AT SUB-ARC END POINTS\r
+               //                     TSTNG=1 - IN ADDITION TO TESTING AT SUB-ARC END\r
+               //                               POINTS TEST ALSO AT INTERIOR POINTS\r
+               //                               ON EACH SUB-ARC.\r
+               \r
+               //            OULVL  - INTEGER\r
+               //                     EITHER 0,1,2,3,4 OR 5.\r
+               //                     CONTROLS THE AMOUNT OF OUTPUT IN THE LISTING FILE\r
+               //                     <JBNM>pl.\r
+               //                     OULVL=0 - OUTPUT A SOLUTION SUMMARY AT EACH STAGE\r
+               //                               IN THE ADAPTIVE PROCESS AND A SHORT\r
+               //                               SUMMARY OF THE ERRORS IN MODULUS.\r
+               //                     OULVL=1 - AS 0, BUT ALSO OUTPUT A DETAILED LIST OF\r
+               //                               THE ERRORS IN MODULUS.\r
+               //                     OULVL=2 - AS 0, BUT ALSO OUTPUT FULL DETAILS OF \r
+               //                               THE FINAL COMPUTED JACOBI COEFFICIENTS \r
+               //                               ON SUCCESSFUL COMPLETION.\r
+               //                     OULVL=3 - AS 2, BUT ALSO OUTPUT A DETAILED LIST OF\r
+               //                               THE ERRORS IN MODULUS.\r
+               //                     OULVL=4 - OUTPUT FULL DETAILS OF THE OF THE COMPU-\r
+               //                               TED JACOBI COEFFICIENTS AT EVERY STAGE\r
+               //                               IN THE ADAPTIVE PROCESS AND A SHORT\r
+               //                               SUMMARY OF THE ERRORS IN MODULUS.\r
+               //                     OULVL=5 - AS 4, BUT ALSO OUTPUT A DETAILED LIST OF\r
+               //                               THE ERRORS IN MODULUS.\r
+               \r
+               //            IBNDS  - INTEGER ARRAY\r
+               //                     INTEGER VECTOR OF SIZE AT LEAST 5.\r
+               //                     IBNDS(K), K=1,2,3,4,5, DEFINE VARIOUS UPPER LIMITS\r
+               //                     THAT HAVE BEEN SET IN THE CALLING PROGRAM AND \r
+               //                     WHICH CONTROL THE SIZES OF THE ARRAYS IGEOM,RGEOM,\r
+               //                     MATRX,ISNPH,RSNPH,IWORK,RWORK,ZWORK,LWORK. \r
+               //                     THEIR MEANINGS ARE AS FOLLOWS:\r
+               //                     IBNDS(1) - THE MAXIMUM NUMBER OF SUB-ARCS ALLOWED.\r
+               //                     IBNDS(2) - THE MAXIMUM NUMBER OF JACOBI INDECES\r
+               //                                ALLOWED (WHICH IS ALSO THE 1 + THE\r
+               //                                MAXIMUM NUMBER OF CORNERS ALLOWED ON \r
+               //                                PHYSICAL BOUNDARY).\r
+               //                     IBNDS(3) - 1 + THE MAXIMUM NUMBER OF PANELS\r
+               //                                ALLOWED IN A SINGLE COMPOSITE GAUSSIAN\r
+               //                                RULE.\r
+               //                     IBNDS(4) - THE MAXIMUM TOTAL NUMBER OF QUADRATURE \r
+               //                                POINTS ALLOWED OVER ALL COMPOSITE \r
+               //                                GAUSSIAN RULES.\r
+               //                                (IBNDS(4)<=(IBNDS(3)-1)*IBNDS(2)*NQPTS)\r
+               \r
+               //            MNEQN  - INTEGER\r
+               //                     THE MAXIMUM NUMBER OF EQUATIONS ALLOWED IN THE \r
+               //                     LINEAR ALGEBRAIC SYSTEM RESULTING FROM THE \r
+               //                     COLLOCATION METHOD. (MNEQN <= 1+IBNDS(1)*NQPTS)\r
+               \r
+               //            MATRX  - REAL ARRAY\r
+               //                     A 3-DIMENSIONAL MATRIX OF SIZE \r
+               //                          MNEQN X MNEQN X 2 .\r
+               //                     (IN THE ADAPTIVE PROCESS, MATRX(*,*,2) WILL STORE \r
+               //                     THE COEFFICIENT MATRIX OF THE CURRENT COLLOCATION \r
+               //                     SYSTEM AND MATRX(*,*,1) WILL STORE THE COEFFICIENT\r
+               //                     MATRIX OF THE PREVIOUS SYSTEM)\r
+               \r
+               //            IWORK  - INTEGER ARRAY\r
+               //                     A WORKING VECTOR OF SIZE AT LEAST\r
+               //                        8*IBNDS(1)+MNEQN+2*IBNDS(2) .\r
+               \r
+               //            RWORK  - REAL ARRAY\r
+               //                     A WORKING VECTOR OF SIZE AT LEAST\r
+               //                       (4 + 3*NQPTS + 5*IBNDS(2))*NQPTS + 2*IBNDS(1) +\r
+               //                       2*MNEQN + IBNDS(3) + 5*IBNDS(2) + 2*IBNDS(4)\r
+               \r
+               //            ZWORK  - COMPLEX ARRAY\r
+               //                     A WORKING VECTOR OF SIZE AT LEAST\r
+               //                         MNEQN + 2*IBNDS(2)\r
+               \r
+               //            LWORK  - LOGICAL ARRAY\r
+               //                     A WORKING VECTOR OF SIZE AT LEAST\r
+               //                         3*IBNDS(1) + IBNDS(2)\r
+               \r
+               //            OCH    - INTEGER\r
+               //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+               //                     WRITING THE FILES <JBNM>pl, <JBNM>gm, <JBNM>ph.\r
+               \r
+               //         ON EXIT\r
+               //            RFARG  - REAL\r
+               //                     EXIT VALUE IS PI*(ENTRY VALUE)\r
+               \r
+               //            IGEOM  - INTEGER ARRAY\r
+               //                     A VECTOR OF SIZE AT LEAST \r
+               //                          IBNDS(1) + 4;\r
+               //                     STORES DATA RELATING TO THE ARC SUBDIVISIONS THAT\r
+               //                     HAVE TAKEN PLACE.\r
+               \r
+               //            RGEOM  - REAL ARRAY\r
+               //                     A VECTOR OF SIZE AT LEAST \r
+               //                          3*IBNDS(1)+2;\r
+               //                     STORES DATA RELATING TO THE ARC SUBDIVISIONS THAT\r
+               //                     HAVE TAKEN PLACE AND THE ARGUMENTS OF SUB-ARC END\r
+               //                     POINTS ON THE UNIT DISC.\r
+               \r
+               //            ISNPH  - INTEGER ARRAY\r
+               //                     A SOLUTION VECTOR OF SIZE AT LEAST \r
+               //                          3*IBNDS(1)+6;\r
+               //                     STORES DATA DEFINING THE FINAL POLYNOMIAL DEGREES\r
+               //                     ON THE SUB-ARCS, THE JACOBI INDEX FOR EACH SUB-ARC\r
+               //                     AND POINTERS TO THE SOLUTIONS STORED IN RSNPH.\r
+               \r
+               //            RSNPH  - REAL ARRAY\r
+               //                     A SOLUTION VECTOR OF SIZE AT LEAST\r
+               //                         IBNDS(1)+2*MNEQN+3*IBNDS(2)*(1+2*NQPTS);\r
+               //                     STORES DATA DEFINING THREE-TERM RECURRENCE\r
+               //                     SCHEMES, ELEMENTARY GAUSS-JACOBI QUADRATURE RULES,\r
+               //                     THE JACOBI COEFFICIENTS FOR THE BOUNDARY\r
+               //                     CORRESPONDENCE FUNCTION AND ITS DERIVATIVE AND\r
+               //                     THE ERRORS IN MODULUS ON EACH SUB-ARC.\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>pl.\r
+               //                     IER=0 - NORMAL EXIT.\r
+               //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+               //                             BE SELF EXPLANATORY.\r
+                                                 \r
+               \r
+               // 4.     SUBROUTINES OR FUNCTIONS NEEDED\r
+               //              - THE CONFPACK LIBRARY.\r
+               //              - THE REAL FUNCTION R1MACH, WHICH IS A MACHINE CONSTANTS \r
+               //                ROUTINE OBTAINED FROM THE PORT LIBRARY.   \r
+               //                IT MUST BE ADJUSTED TO SUIT EACH PARTICULAR MACHINE.  \r
+               //                IF YOUR MACHINE ISN'T LISTED IN R1MACH THEN YOU'LL  \r
+               //                HAVE TO WRITE YOUR OWN VERSION, BUT NOTE THAT CONFPACK \r
+               //                ONLY USES R1MACH(1), R1MACH(2) AND R1MACH(4).\r
+               //              - THE FOLLOWING LINPACK ROUTINES:\r
+               //                    ISAMAX   SASUM    SAXPY   SDOT    SGECO\r
+               //                    SGEFA    SGEDI    SGESL   SSCAL   SSWAP\r
+               //              - THE FOLLOWING QUADPACK ROUTINES:\r
+               //                    QAWS     QAWSE    QC25S   QCHEB   QK15W\r
+               //                    QMAC     QMOMO    QSORT   QWGTS\r
+               //              - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN\r
+               //                WHICH DEFINE THE PARAMETRIC FUNCTION FOR THE PHYSICAL\r
+               //                BOUNDARY AND THE DERIVATIVE OF THE PARAMETRIC FUNCTION.\r
+               //                THE PARAMETRIC FUNCTION DEFINING THE K'TH ANALYTIC ARC\r
+               //                SHOULD HAVE THE SUBROUTINE HEADING\r
+               \r
+               //                    COMPLEX FUNCTION PARFUN(K,T)\r
+               //                    INTEGER K\r
+               //                    COMPLEX T\r
+               \r
+               //                WITH THE REAL PARAMETER INTERVAL -1 < REAL(T) < +1\r
+               //                BEING MAPPED TO THE PHYSICAL ARC.  A SIMILAR HEADING\r
+               //                SHOULD BE GIVEN FOR THE DERIVATIVE DPARFN.  THE PRE-\r
+               //                PROCESSING PROGRAM PARGEN IS AVAILABLE TO HELP WITH\r
+               //                THE CREATION OF PARFUN AND DPARFN.\r
+               \r
+               \r
+               // 5.     FURTHER COMMENTS\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: 15 JULY 1990\r
+               // ......................................................................\r
+                    \r
+               //     LOCAL VARAIBLES\r
+               \r
+               //**** POINTERS FOR IGEOM,RGEOM,ISNPH,RSNPH\r
+               \r
+               int ACOEF,AICOF,BCFSN,BCOEF,BICOF,DGPOL,ERARC,H0VAL,HALEN,\r
+                    HIVAL,JACIN,JATYP,LOSUB,MIDPT,PARNT,QUPTS,QUWTS,SOLUN,VTARG;\r
+               \r
+               //**** POINTERS FOR IWORK,RWORK,ZWORK,LWORK\r
+               \r
+               int A1COF,AQCOF,AXION,B1COF,BQCOF,COLPR,COLSC,CQCOF,DIAG,\r
+                    HIOLD,HISUB,HITES,ICOPY,IPIVT,LCOPY,LNSEG,LOOLD,LOQSB,LOTES,NEWDG,\r
+                    NEWHL,NEWQU,NQUAD,PNEWQ,QCOMW,QCOMX,RCOPY,RIGLL,SDIAG,TOLOU,WORK2,\r
+                    WORKQ,WORK,WORKT,XENPT,XIDST,XIVAL,ZCOLL;\r
+               \r
+               //**** OTHER LOCAL VARIABLES\r
+               \r
+               int I,IMXER,INDEG,J,L,MDGPO,MNJXS,MNQUA,MNSUA,MQIN1,NCOLL,\r
+                    NEFF,NEQNS,NJIND,NROWS,NTEST,\r
+                    TNSUA,TNGQP,ORDSG;\r
+               int SOLCO = 0;\r
+               int QIERC[] = new int[7];\r
+               int QIERR[] = new int[7];\r
+               //int QIERC(0:6),QIERR(0:6)\r
+               //DATA QIERC/7*0E+0/\r
+               \r
+               final double SFACT = 0.8;\r
+               final double QFACT = 0.1;\r
+               double AQTOL,CONST,GAQTL,GLGTL,GRQTL,GSUPE,GTGTE,LGTOL,ESTOL,\r
+                    MCHEP,MCQER,MQERR,MXERM,PI,R1MACH,RCOND,RQTOL,SSUPE,\r
+                    TGTER,TOLNR;\r
+               \r
+               double ZMXER[] = new double[2];\r
+               //COMPLEX ZMXER\r
+               \r
+               final boolean INIBT = true;\r
+               boolean ACCPT,GACPT,NUQTL,REFLN;\r
+               \r
+               String OFL;\r
+               //CHARACTER OFL*6\r
+       \r
+               //EXTERNAL AXION1,ANGLE7,ASQUC7,BCFVTF,CPJAC3,CSCAL3,ICOQR1,IGNLVL,\r
+               //         LINSEG,LNSY11,OPQUD1,OUPTGM,OUPTPH,R1MACH,RECON,RESCAL,RSLT80,\r
+               //         RSLT71,RSLT72,RSLT83,RSLT84,SETIGL,SGECO,SGEDI,SGESL,TESMD9,\r
+               //         TSJAC3,UPCOQ1,UPJAC1,WRHEAD,WRTAIL\r
+               \r
+               //C**** DEFINE SOME OUTPUT FORMATS\r
+               \r
+               //1     FORMAT(A45)\r
+               //3     FORMAT(A45,2X,E9.2)\r
+               \r
+               //**** NAME AND OPEN THE MAIN LISTING FILE AND OUTPUT THE JOBNAME TO FILE\r
+               //**** jbnm.\r
+               \r
+               //OPEN(OCH,FILE='jbnm')\r
+               //WRITE(OCH,'(A4)') JBNM\r
+               //CLOSE(OCH)\r
+               //L=INDEX(JBNM,' ')-1\r
+               //IF (L.EQ.-1) L=4\r
+               //OFL=JBNM(1:L)//'pl'\r
+               //OPEN(OCH,FILE=OFL)\r
+               \r
+               //**** OUTPUT CONFPACK HEADING\r
+               \r
+               WRHEAD(1,0,null);\r
+               //WRHEAD(1,OCH)\r
+               \r
+               if (NQPTS < 1) {\r
+                   IER[0]=3;\r
+                   WRTAIL(1,0,IER[0],null);\r
+                   return;\r
+               } \r
+               \r
+               //**** INITIALISE SOME VARIABLES\r
+               \r
+               if (ISYGP == 0 || ISYGP == -1) {\r
+                   ORDSG=1;\r
+                   REFLN=false;\r
+               }\r
+               else {\r
+                   ORDSG=Math.abs(ISYGP);\r
+                   REFLN=(ISYGP < -1);\r
+               }\r
+               \r
+               if ((NARCS % ORDSG) != 0) {\r
+                   IER[0]=55;\r
+                   WRTAIL(1,0,IER[0],null);\r
+                   return;\r
+               }\r
+               \r
+               SOLCO=0;\r
+               NEFF=0;\r
+               MCHEP=EPS;\r
+               TOLNR=Math.sqrt(MCHEP);\r
+               NJIND=NARCS+1;\r
+               TNGQP=NQPTS*NJIND;\r
+               MDGPO=NQPTS-1;\r
+               MNSUA=IBNDS[0];\r
+               MNJXS=IBNDS[1];\r
+               MQIN1=IBNDS[2];\r
+               MNQUA=IBNDS[3];\r
+               if (2*NARCS > MNSUA) {\r
+                   IER[0]=1;\r
+                   WRTAIL(1,0,IER[0],null);\r
+                   return;\r
+               }\r
+               if (NARCS+1 > MNJXS) {\r
+                   IER[0]=2;\r
+                   WRTAIL(1,0,IER[0],null);\r
+                   return;            \r
+               }\r
+               if (TSTNG != 1) {\r
+                   TSTNG=0;\r
+               }\r
+               GSUPE=MAXER;\r
+               GTGTE=GSUPE*SFACT;\r
+               GAQTL=QFACT*GTGTE;\r
+               if (GAQTL < 5*MCHEP) {\r
+                   GAQTL=5*MCHEP;\r
+                   GTGTE=GAQTL/QFACT;\r
+                   GSUPE=GTGTE/SFACT;\r
+               }\r
+               GLGTL=Math.log(1+GTGTE);\r
+               GRQTL=GAQTL;\r
+               IGEOM[0]=NARCS;\r
+               IGEOM[1]=NQPTS;\r
+               IGEOM[3]=MNSUA;\r
+               ISNPH[0]=NARCS;\r
+               ISNPH[1]=NQPTS;\r
+               ISNPH[4]=MNSUA;\r
+               ISNPH[5]=MNEQN;\r
+               RGEOM[0]=GSUPE;\r
+               RGEOM[1]=GLGTL;\r
+               \r
+               //**** SET UP THE POINTERS TO ELEMENTS IN ARRAYS IGEOM AND RGEOM\r
+               \r
+               PARNT=5;\r
+               HALEN=3;\r
+               MIDPT=MNSUA+3;\r
+               VTARG=2*MNSUA+3; \r
+               \r
+               //**** SET UP THE POINTERS TO ELEMENTS IN ARRAYS ISNPH AND RSNPH\r
+               \r
+               DGPOL=7;\r
+               JATYP=MNSUA+7;\r
+               LOSUB=2*MNSUA+7;\r
+               ACOEF=1;\r
+               BCOEF=TNGQP+1;\r
+               AICOF=2*TNGQP+1;\r
+               BICOF=3*TNGQP+1;\r
+               QUPTS=4*TNGQP+1;\r
+               QUWTS=5*TNGQP+1;\r
+               H0VAL=6*TNGQP+1;\r
+               HIVAL=NJIND+6*TNGQP+1;\r
+               JACIN=2*NJIND+6*TNGQP+1;\r
+               ERARC=3*NJIND+6*TNGQP+1;\r
+               BCFSN=MNSUA+3*NJIND+6*TNGQP+1;\r
+               SOLUN=MNEQN+MNSUA+3*NJIND+6*TNGQP+1;\r
+               \r
+               //**** SET UP THE POINTERS TO ELEMENTS IN ARRAYS IWORK,RWORK,ZWORK AND \r
+               //**** LWORK\r
+               \r
+               IPIVT=1;\r
+               LOQSB=MNEQN+1;\r
+               NQUAD=MNJXS+MNEQN+1;\r
+               HISUB=2*MNJXS+MNEQN+1;\r
+               LOTES=MNSUA+2*MNJXS+MNEQN+1;\r
+               HITES=2*MNSUA+2*MNJXS+MNEQN+1;\r
+               AXION=3*MNSUA+2*MNJXS+MNEQN+1;\r
+               NEWDG=4*MNSUA+2*MNJXS+MNEQN+1;\r
+               ICOPY=5*MNSUA+2*MNJXS+MNEQN+1;\r
+               LOOLD=6*MNSUA+2*MNJXS+MNEQN+1;\r
+               HIOLD=7*MNSUA+2*MNJXS+MNEQN+1;\r
+               WORK2=1;\r
+               COLPR=MNEQN+1;\r
+               A1COF=2*MNEQN+1;\r
+               B1COF=MNJXS+2*MNEQN+1;\r
+               TOLOU=2*MNJXS+2*MNEQN+1;\r
+               XIDST=3*MNJXS+2*MNEQN+1;\r
+               XENPT=5*MNJXS+2*MNEQN+1;\r
+               QCOMX=MQIN1+5*MNJXS+2*MNEQN+1;\r
+               QCOMW=MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               RCOPY=2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               NEWHL=MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               AQCOF=2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               BQCOF=TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               CQCOF=2*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               COLSC=3*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               RIGLL=4*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               WORK=5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               DIAG=2*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               SDIAG=3*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               WORKT=4*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+2*MNEQN+1;\r
+               WORKQ=2*NQPTS*NQPTS+4*NQPTS+5*TNGQP+2*MNSUA+2*MNQUA+MQIN1+5*MNJXS+\r
+                     2*MNEQN+1;\r
+               ZCOLL=1;\r
+               XIVAL=MNEQN+1;\r
+               NEWQU=1;\r
+               LCOPY=MNJXS+1;\r
+               PNEWQ=MNSUA+MNJXS+1;\r
+               LNSEG=2*MNSUA+MNJXS+1;\r
+               \r
+               //**** ASSIGN THE JACOBI INDECES FOR EACH ARC.\r
+               \r
+               //ANGLE7(RSNPH(JACIN),NARCS,INTER);\r
+               double BE[] = new double[NARCS];\r
+               ANGLE7(BE, NARCS, INTER);\r
+               for (I = 0; I < NARCS; I++) {\r
+                       RSNPH[JACIN + I - 1] = BE[I];\r
+               }\r
+               RSNPH[JACIN+NJIND-2]=0.0;\r
+               \r
+               //**** SET SUB-TOLERANCES AND INDEG\r
+               \r
+               if (INCST && GSUPE <= 3.16E-2) {\r
+               \r
+                   //****   FOLLOW INCREMENTAL STRATEGY\r
+               \r
+                   SSUPE=0.1;\r
+                   TGTER=SSUPE*SFACT;\r
+                   AQTOL=TGTER*QFACT;\r
+                   LGTOL=Math.log(1.0+TGTER);\r
+                   RQTOL=AQTOL;\r
+                   INDEG=Math.min(3,NQPTS-1);\r
+               }\r
+               else {\r
+               \r
+                   //****   SUB-TOLERANCES SAME AS GLOBAL TOLERANCES, INDEG DETERMINED\r
+                   // ****   ACCORDING TO ACCURACY REQUESTED\r
+               \r
+                   SSUPE=GSUPE;\r
+                   TGTER=GTGTE;\r
+                   AQTOL=GAQTL;\r
+                   LGTOL=GLGTL;\r
+                   RQTOL=GRQTL;\r
+                   INDEG=(int)Math.round((-Math.log10(TGTER)))+2;\r
+                   INDEG=Math.min(INDEG,NQPTS-1);\r
+               }\r
+               \r
+               //**** ASSIGN THE LOGICAL LINE SEGMENT TYPE TO EACH ARC.\r
+               \r
+               /*      CALL LINSEG(LWORK(LNSEG),NARCS)\r
+               C\r
+               C**** LIST THE INPUT ARGUMENTS AND ASSOCIATED QUANTITIES\r
+               C\r
+                     CALL RSLT80(JBNM,HEAD,GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
+                    +INCST,INDEG,RFARC,RFARG,CENTR,RSNPH(JACIN),LWORK(LNSEG),\r
+                    +TSTNG,OULVL,IBNDS,MNEQN,OCH)\r
+                     PI=4E+0*ATAN(1E+0)\r
+                     RFARG=RFARG*PI\r
+               C\r
+               C**** SET UP THE GAUSS-JACOBI AND GAUSS-LEGENDRE QUADRATURE DATA AND \r
+               C**** STORE IN ARRAYS QUPTS AND QUWTS.  SET UP THREE TERM RECURRENCE\r
+               C**** COEFFICIENTS AND STORE IN ACOEF, BCOEF.  DETERMINE ZEROTH\r
+               C**** MOMENTS OF JACOBI DISTRIBUTIONS AND STORE IN H0VAL. \r
+               C**** ALSO SET UP THREE TERM RECURRENCE COEFFICIENTS AND ZEROTH MOMENTS\r
+               C**** FOR THE INTEGRATED POLYNOMIALS, STORING RESULTS IN AICOF,BICOF\r
+               C**** AND HIVAL.\r
+               C      \r
+                     CALL OPQUD1(NJIND,NQPTS,RSNPH(JACIN),RSNPH(ACOEF),RSNPH(BCOEF),\r
+                    +RSNPH(H0VAL),RSNPH(AICOF),RSNPH(BICOF),RSNPH(HIVAL),RSNPH(QUPTS),\r
+                    +RSNPH(QUWTS),RWORK(WORK),IER)\r
+                     IF (IER .GT. 0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+                     J=1-NQPTS\r
+                     DO 10 I=1,NJIND\r
+                       J=J+NQPTS\r
+                       RWORK(A1COF+I-1)=RSNPH(ACOEF+J-1)\r
+                       RWORK(B1COF+I-1)=RSNPH(BCOEF+J-1)\r
+               10    CONTINUE  \r
+                     WRITE(*,1) 'BASIC GAUSS QUADRATURE DATA DONE:'\r
+               C\r
+               C**** SET UP THE COEFFICIENTS IN THE THREE TERM RECURRENCE FORMULAE\r
+               C**** FOR THE PRINCIPAL SINGULAR INTEGRALS ASSOCIATED WITH THE VARIOUS\r
+               C**** JACOBI WEIGHT FUNCTIONS AND THEIR ORTHONORMAL POLYNOMIALS; STORE\r
+               C**** THESE COEFFICIENTS IN AQCOF, BQCOF AND CQCOF\r
+               C\r
+                     CALL ASQUC7(RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),RSNPH(JACIN),\r
+                    +NJIND,NQPTS)\r
+                     WRITE(*,1) 'DATA FOR SINGULAR INTEGRALS DONE:'\r
+               C\r
+               C**** SET UP THE A PRIORI COLUMN SCALE FACTORS, STORED IN COLSC.\r
+               C\r
+                     CALL CSCAL3(RWORK(COLSC),NQPTS,NJIND,RSNPH(ACOEF),RSNPH(BCOEF),\r
+                    +RSNPH(H0VAL),RSNPH(QUPTS),RSNPH(QUWTS),RSNPH(JACIN),RWORK(WORK),\r
+                    +RWORK(WORKT),RWORK(WORKQ))\r
+               C\r
+               C**** SET UP THE ARRAY RIGLL OF REFERENCE IGNORE LEVELS.\r
+               C\r
+                     CALL IGNLVL(RWORK(RIGLL),RWORK(COLSC),RSNPH(ACOEF),RSNPH(BCOEF),\r
+                    +RSNPH(H0VAL),RSNPH(JACIN),NJIND,NQPTS,IER)\r
+                     IF (IER .GT. 0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+               C\r
+               C**** SET UP THE ARRAY OF COLLOCATION POINTS PARAMETER VALUES, COLPR,\r
+               C**** THE ARRAY OF COLLOCATION POINTS ZCOLL AND THE ARRAYS LOSUB AND \r
+               C**** HISUB NEEDED TO ACCESS COLPR AND ZCOLL CORRECTLY.  INITIALISE\r
+               C**** DGPOL AND UPDATE LNSEG FOR ARC HALVING. \r
+               C\r
+                     CALL CPJAC3(NARCS,NQPTS,INDEG,ISNPH(DGPOL),RSNPH(JACIN),\r
+                    +RSNPH(ACOEF),RSNPH(BCOEF),RWORK(DIAG),RWORK(SDIAG),TNSUA,\r
+                    +ISNPH(LOSUB),IWORK(HISUB),ISNPH(JATYP),IGEOM(PARNT),RGEOM(MIDPT),\r
+                    +RGEOM(HALEN),RWORK(COLPR),ZWORK(ZCOLL),LWORK(LNSEG),IWORK(LOOLD),\r
+                    +IWORK(HIOLD),EPS,IER,INIBT)\r
+                     IF (IER .GT. 0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+                     NCOLL=IWORK(HISUB+TNSUA-1)\r
+                     NEQNS=NCOLL+1\r
+                     NROWS=NCOLL/ORDSG+1\r
+                     IF (NEQNS .GT. MNEQN) THEN\r
+                       IER=8\r
+                       GOTO 999\r
+                     ENDIF\r
+                     WRITE(*,1) 'COLLOCATION POINT CHOICE DONE:'\r
+               C\r
+               C**** SET UP THE COMPOSITE GAUSSIAN QUADRATURE RULES, STORING ABSCISSAE\r
+               C**** AND WEIGHTS IN QCOMX AND QCOMW.  SET UP ARRAYS NQUAD,LOQSB\r
+               C**** NEEDED TO ACCESS THESE DATA.  RECORD MAXIMUM QUADRATURE ERRORS\r
+               C**** FOR COLUMN SCALED INTEGRALS IN ARRAY TOLOU.\r
+               C\r
+                     CALL ICOQR1(NARCS,NJIND,NQPTS,MDGPO,MQIN1,AQTOL,RSNPH(QUPTS),\r
+                    +RSNPH(QUWTS),RSNPH(JACIN),RGEOM(MIDPT),RGEOM(HALEN),RSNPH(ACOEF),\r
+                    +RSNPH(BCOEF),RSNPH(H0VAL),RWORK(COLSC),IWORK(NQUAD),IWORK(LOQSB),\r
+                    +RWORK(QCOMX),RWORK(QCOMW),MNQUA,RWORK(TOLOU),MCQER,RWORK(XENPT),\r
+                    +ZWORK(XIVAL),RWORK(XIDST),IER)\r
+                     NUQTL=.FALSE.\r
+                     IF (IER .GT. 0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+                     WRITE(*,1) 'COMPOSITE GAUSSIAN RULES DONE:'\r
+               C\r
+               C**** SET UP LINEAR ALGEBRAIC SYSTEM.\r
+               C\r
+               23    CONTINUE\r
+                     SOLCO=SOLCO+1\r
+                     WRITE(*,24) '********SOLUTION',SOLCO,'********',NROWS,'EQUATIONS'\r
+               24    FORMAT(/,T18,A,1X,I2,1X,A,/,T25,I3,1X,A)\r
+               C\r
+                     CALL LNSY11(MATRX,RSNPH(SOLUN),MNEQN,NCOLL,ORDSG,REFLN,NQPTS,\r
+                    +TNSUA,ISNPH(JATYP),IGEOM(PARNT),ISNPH(DGPOL),ISNPH(LOSUB),\r
+                    +IWORK(HISUB),IWORK(NQUAD),IWORK(LOQSB),TOLNR,RGEOM(MIDPT),\r
+                    +RGEOM(HALEN),RSNPH(H0VAL),RWORK(COLSC),RSNPH(ACOEF),RSNPH(BCOEF),\r
+                    +RWORK(COLPR),RWORK(QCOMX),RWORK(QCOMW),CENTR,ZWORK(ZCOLL),INTER,\r
+                    +LWORK(LNSEG),RWORK(WORK),QIERR,MQERR,RSNPH(JACIN),RWORK(A1COF),\r
+                    +RWORK(B1COF),AQTOL,RQTOL,RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),\r
+                    +IWORK(LOOLD),IWORK(HIOLD))\r
+               C\r
+                     DO 25 I=0,6\r
+                       QIERC(I)=QIERC(I)+QIERR(I)\r
+               25    CONTINUE\r
+                     WRITE(*,1) 'LINEAR SYSTEM SET UP DONE:'\r
+               C\r
+               C**** SOLVE LINEAR SYSTEM BY GAUSSIAN ELIMINATION USING LINPACK\r
+               C\r
+                     CALL SGECO(MATRX(1,1,2),MNEQN,NROWS,IWORK(IPIVT),RCOND,\r
+                    +RWORK(WORK2))\r
+                     IF (RCOND .EQ. 0E+0) THEN\r
+                       IER=15\r
+                       SOLCO=SOLCO-1\r
+                       GOTO 999\r
+                     ENDIF    \r
+                     CALL SGESL(MATRX(1,1,2),MNEQN,NROWS,IWORK(IPIVT),RSNPH(SOLUN),0)\r
+                     NEFF=NEFF+NROWS**3\r
+                     WRITE(*,1) 'LINEAR SYSTEM SOLUTION DONE:'\r
+               C\r
+               C**** RECONSTITUTE FULL SOLUTION VECTOR\r
+               C\r
+                     IF (ORDSG.GT.1) THEN\r
+                       CALL RECON(ORDSG,REFLN,NCOLL,TNSUA,ISNPH(LOSUB),IWORK(HISUB),\r
+                    +  RSNPH(SOLUN))\r
+                     ENDIF\r
+                     CONST=RSNPH(SOLUN+NEQNS-1)\r
+               C\r
+               C**** SET UP THE ARRAY WORK2 OF ACTUAL COEFFICIENT IGNORE LEVELS\r
+               C\r
+                     CALL SETIGL(RWORK(WORK2),IWORK(HISUB),ISNPH(JATYP),ISNPH(LOSUB),\r
+                    +NQPTS,RWORK(RIGLL),TNSUA)\r
+               C\r
+               C**** DETERMINE THE ACTIONS THAT HAVE TO BE TAKEN ON EACH ARC\r
+               C\r
+                     CALL AXION1(IWORK(AXION),IWORK(NEWDG),RSNPH(SOLUN),MDGPO,TNSUA,\r
+                    +ISNPH(DGPOL),ISNPH(LOSUB),IWORK(HISUB),RWORK(RIGLL),LGTOL,ACCPT,\r
+                    +RSNPH(JACIN),ISNPH(JATYP),NJIND,RWORK(NEWHL),ESTOL,IER)\r
+                     ESTOL=ESTOL/SFACT\r
+                     IF (IER.GT.0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+                     WRITE(*,1) 'DECISIONS FOR EACH ARC DONE:'\r
+                     WRITE(*,3) 'EFFECTIVE STOPPING TOLERANCE:',ESTOL\r
+                     IF (ACCPT .AND. ESTOL.LE.GSUPE) THEN\r
+                       GACPT=.TRUE.\r
+                     ELSE\r
+                       GACPT=.FALSE.\r
+                     ENDIF\r
+               C\r
+                     IF (GACPT) THEN\r
+               C\r
+               C****   OUTPUT RESULTS\r
+               C\r
+                       IF (OULVL .LT. 4) THEN\r
+                         CALL RSLT72(QIERC,RCOND,CONST,NROWS,ISNPH(DGPOL),ISNPH(JATYP),\r
+                    +                IGEOM(PARNT),TNSUA,INTER,MQERR,MCQER,IWORK(AXION),\r
+                    +                IWORK(NEWDG),NJIND,IWORK(NQUAD),RWORK(TOLOU),\r
+                    +                LGTOL,SOLCO,OCH)\r
+                       ELSE \r
+                         CALL RSLT71(QIERC,RCOND,RSNPH(SOLUN),NEQNS,ISNPH(LOSUB),\r
+                    +                IWORK(HISUB),RWORK(COLSC),NQPTS,ISNPH(JATYP),\r
+                    +                IGEOM(PARNT),TNSUA,INTER,MQERR,MCQER,RWORK(WORK2),\r
+                    +                IWORK(AXION),IWORK(NEWDG),NJIND,RSNPH(JACIN),\r
+                    +                IWORK(NQUAD),RWORK(TOLOU),LGTOL,SOLCO,OCH)\r
+                       ENDIF\r
+                       WRITE(OCH,*) 'EFFECTIVE STOPPING TOLERANCE :',ESTOL\r
+                       NEFF=NINT(REAL(NEFF)**3.3333333E-1)\r
+                       WRITE(*,54) '****THE SOLUTION IS ACCEPTED****'\r
+                       WRITE(*,55) 'EFFECTIVE SIZE OF ALL SYSTEMS:',NEFF\r
+                       IF (INTER) THEN\r
+                         WRITE(*,3) 'ZERO:',CONST\r
+                       ELSE\r
+                         WRITE(*,56) 'CAPACITY:',EXP(-CONST)\r
+                       ENDIF\r
+               54      FORMAT(/,T17,A)\r
+               55      FORMAT(/,A45,I4)\r
+               56      FORMAT(A45,E16.8)\r
+\r
+                       WRITE(OCH,*)\r
+                       WRITE(OCH,*) '****THE SOLUTION IS ACCEPTED****'\r
+                       WRITE(OCH,*) 'EFFECTIVE SIZE OF ALL SYSTEMS : ',NEFF\r
+                       WRITE(OCH,*) \r
+                     ELSE\r
+                       IF (ACCPT .OR. ESTOL.LE.SSUPE) THEN\r
+               C\r
+               C****     SOLUTION AT INTERMEDIATE ACCURACY IS ACCEPTED; SET TOLERANCES\r
+               C****     FOR GREATER ACCURACY AND RE-ASSESS UPDATING ACTIONS BEFORE \r
+               C****     CONTINUING\r
+               C\r
+                         SSUPE=1E-1*MIN(SSUPE,ESTOL)\r
+                         TGTER=SFACT*SSUPE\r
+                         IF (TGTER .LE. 2E+0*GTGTE) THEN\r
+                           TGTER=GTGTE\r
+                         ENDIF\r
+                         AQTOL=TGTER*QFACT\r
+                         NUQTL=.TRUE.\r
+                         LGTOL=LOG(1E+0+TGTER)\r
+                         RQTOL=AQTOL\r
+                         I=NINT(-LOG10(TGTER))+2\r
+                         INDEG=MIN(I,MDGPO)\r
+               C\r
+               C****     DETERMINE THE ACTIONS THAT HAVE TO BE TAKEN ON EACH ARC\r
+               C\r
+                         CALL AXION1(IWORK(AXION),IWORK(NEWDG),RSNPH(SOLUN),MDGPO,\r
+                    +                TNSUA,ISNPH(DGPOL),ISNPH(LOSUB),IWORK(HISUB),\r
+                    +                RWORK(RIGLL),LGTOL,ACCPT,RSNPH(JACIN),\r
+                    +                ISNPH(JATYP),NJIND,RWORK(NEWHL),ESTOL,IER)\r
+                         ESTOL=ESTOL/SFACT\r
+                         IF (IER.GT.0) THEN\r
+                           GOTO 999\r
+                         ENDIF\r
+                         WRITE(*,1) 'DECISIONS FOR EACH ARC RE-DONE:'\r
+                       ENDIF\r
+               C\r
+               C****   OUTPUT RESULTS\r
+               C\r
+                       IF (OULVL .LT. 4) THEN\r
+                         CALL RSLT72(QIERC,RCOND,CONST,NROWS,ISNPH(DGPOL),ISNPH(JATYP),\r
+                    +                IGEOM(PARNT),TNSUA,INTER,MQERR,MCQER,IWORK(AXION),\r
+                    +                IWORK(NEWDG),NJIND,IWORK(NQUAD),RWORK(TOLOU),\r
+                    +                LGTOL,SOLCO,OCH)\r
+                       ELSE \r
+                         CALL RSLT71(QIERC,RCOND,RSNPH(SOLUN),NEQNS,ISNPH(LOSUB),\r
+                    +                IWORK(HISUB),RWORK(COLSC),NQPTS,ISNPH(JATYP),\r
+                    +                IGEOM(PARNT),TNSUA,INTER,MQERR,MCQER,RWORK(WORK2),\r
+                    +                IWORK(AXION),IWORK(NEWDG),NJIND,RSNPH(JACIN),\r
+                    +                IWORK(NQUAD),RWORK(TOLOU),LGTOL,SOLCO,OCH)\r
+                       ENDIF\r
+                       WRITE(OCH,*) 'EFFECTIVE STOPPING TOLERANCE :',ESTOL\r
+                       IF (RCOND .LT. 5E+0*MCHEP) THEN\r
+                         IER=16\r
+                         GOTO 999\r
+                       ELSE IF (RCOND .LT. AQTOL) THEN\r
+                         NUQTL=.TRUE.\r
+                         AQTOL=1E-1*RCOND\r
+                         IF (AQTOL .LT. 5E+0*MCHEP) AQTOL=5E+0*MCHEP\r
+                       ENDIF\r
+               C\r
+               C****   IMPLEMENT UPDATING PROCEDURES.\r
+               C****   FIRST UPDATE THE COLLOCATION PARAMETERS AND OTHER DATA\r
+               C****   RELATING TO SUB-ARC DEFINITIONS.\r
+               C\r
+                       CALL UPJAC1(NQPTS,NJIND,INDEG,IWORK(AXION),ISNPH(DGPOL),\r
+                    +  IWORK(NEWDG),RSNPH(ACOEF),RSNPH(BCOEF),RWORK(DIAG),\r
+                    +  RWORK(SDIAG),TNSUA,MNSUA,ISNPH(LOSUB),IWORK(HISUB),\r
+                    +  ISNPH(JATYP),IGEOM(PARNT),RGEOM(MIDPT),RGEOM(HALEN),\r
+                    +  RWORK(COLPR),ZWORK(ZCOLL),LWORK(LNSEG),LWORK(PNEWQ),EPS,IER,\r
+                    +  RWORK(WORK),RWORK(NEWHL),RWORK(RCOPY),IWORK(ICOPY),\r
+                    +  LWORK(LCOPY),IWORK(LOOLD),IWORK(HIOLD))\r
+                       IF (IER .GT. 0) THEN\r
+                         GOTO 999\r
+                       ENDIF\r
+                       WRITE(*,1) 'ARC REFINEMENTS DONE:'\r
+                       NCOLL=IWORK(HISUB+TNSUA-1)\r
+                       NEQNS=NCOLL+1\r
+                       NROWS=NCOLL/ORDSG+1\r
+                       IF (NEQNS .GT. MNEQN) THEN\r
+                         IER=18\r
+                         GOTO 999\r
+                       ENDIF\r
+               C\r
+               C****   NEXT UPDATE THE COMPOSITE QUADRATURE RULES\r
+               C\r
+                       CALL UPCOQ1(NARCS,NJIND,NQPTS,MDGPO,MQIN1,AQTOL,RSNPH(QUPTS),\r
+                    +  RSNPH(QUWTS),RSNPH(JACIN),RGEOM(MIDPT),RGEOM(HALEN),\r
+                    +  RSNPH(ACOEF),RSNPH(BCOEF),RSNPH(H0VAL),RWORK(COLSC),\r
+                    +  IWORK(NQUAD),IWORK(LOQSB),RWORK(QCOMX),RWORK(QCOMW),\r
+                    +  MNQUA,RWORK(TOLOU),MCQER,RWORK(XENPT),ZWORK(XIVAL),\r
+                    +  RWORK(XIDST),TNSUA,LWORK(PNEWQ),LWORK(NEWQU),ISNPH(JATYP),\r
+                    +  IGEOM(PARNT),NUQTL,IER)\r
+                       IF (IER .GT. 0) THEN\r
+                         GOTO 999\r
+                       ENDIF\r
+                       WRITE(*,1) 'QUADRATURE UPDATES DONE:'\r
+                       GOTO 23\r
+                     ENDIF \r
+               C\r
+                     IF (OULVL.EQ.2 .OR. OULVL.EQ.3) THEN\r
+                       CALL RSLT71(QIERC,RCOND,RSNPH(SOLUN),NEQNS,ISNPH(LOSUB),\r
+                    +  IWORK(HISUB),RWORK(COLSC),NQPTS,ISNPH(JATYP),IGEOM(PARNT),TNSUA,\r
+                    +  INTER,MQERR,MCQER,RWORK(WORK2),IWORK(AXION),IWORK(NEWDG),NJIND,\r
+                    +  RSNPH(JACIN),IWORK(NQUAD),RWORK(TOLOU),LGTOL,SOLCO,OCH)\r
+                     ENDIF\r
+               C\r
+               C**** ESTIMATE MAXIMUM ERROR IN MODULUS.\r
+               C\r
+                       WRITE(*,*)\r
+                       WRITE(*,1) 'ERRORS IN MODULUS STARTED:'\r
+               C\r
+                     CALL TSJAC3(IWORK(LOTES),IWORK(HITES),RWORK(COLPR),ZWORK(ZCOLL),\r
+                    +NQPTS,NTEST,ORDSG,TNSUA,TSTNG,ISNPH(DGPOL),ISNPH(JATYP),\r
+                    +IGEOM(PARNT),RSNPH(AICOF),RSNPH(BICOF),RWORK(DIAG),RGEOM(HALEN),\r
+                    +RSNPH(JACIN),RGEOM(MIDPT),RWORK(SDIAG),IER)\r
+                     IF (IER .GT. 0) THEN\r
+                       GOTO 999\r
+                     ENDIF\r
+               C\r
+                     CALL TESMD9(RWORK(WORK2),MATRX(1,1,2),RSNPH(SOLUN),MNEQN,NCOLL,\r
+                    +NTEST,NQPTS,TNSUA,ISNPH(JATYP),IGEOM(PARNT),ISNPH(DGPOL),\r
+                    +ISNPH(LOSUB),IWORK(HISUB),IWORK(LOTES),IWORK(HITES),IWORK(NQUAD),\r
+                    +IWORK(LOQSB),TOLNR,RGEOM(MIDPT),RGEOM(HALEN),RSNPH(H0VAL),\r
+                    +RWORK(COLSC),RSNPH(ACOEF),RSNPH(BCOEF),RWORK(COLPR),RWORK(QCOMX),\r
+                    +RWORK(QCOMW),CENTR,ZWORK(ZCOLL),INTER,LWORK(LNSEG),RWORK(WORK),\r
+                    +QIERR,MQERR,RSNPH(JACIN),RWORK(A1COF),RWORK(B1COF),AQTOL,RQTOL,\r
+                    +RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),MXERM,IMXER,ZMXER,\r
+                    +RSNPH(ERARC),ORDSG,REFLN)\r
+               C\r
+                     WRITE(*,1) 'ERRORS IN MODULUS DONE:'\r
+                     WRITE(*,3) 'MAXIMUM ERROR AT TEST POINTS:',MXERM\r
+                     DO 60 I=0,6\r
+                       QIERC(I)=QIERC(I)+QIERR(I)\r
+               60    CONTINUE\r
+                     IF (MOD(OULVL,2) .EQ. 1) THEN\r
+                       CALL RSLT84(RWORK(WORK2),TNSUA,MXERM,ZMXER,IMXER,IWORK(LOTES),\r
+                    +  IWORK(HITES),QIERC,IGEOM(PARNT),ORDSG,OCH)\r
+                     ELSE\r
+                       CALL RSLT83(RSNPH(ERARC),TNSUA,MXERM,ZMXER,IMXER,QIERC,\r
+                    +  IGEOM(PARNT),ORDSG,OCH)\r
+                     ENDIF\r
+               C\r
+               C**** RESCALE SOLUTIONS TO OBTAIN STANDARD JACOBI COEFFICIENTS\r
+                     CALL RESCAL(NQPTS,TNSUA,ISNPH(LOSUB),IWORK(HISUB),ISNPH(JATYP),\r
+                    +RSNPH(SOLUN),RWORK(COLSC))\r
+\r
+               999   CONTINUE\r
+               C\r
+                     CALL WRTAIL(1,OCH,IER)\r
+                     CLOSE(OCH)\r
+               C\r
+                     IF (SOLCO .GE. 1) THEN\r
+               C\r
+               C****   COMPUTE THE BOUNDARY CORRESPONDENCE COEFFICIENTS BCFSN AND THE\r
+               C****   ARGUMENTS OF ALL SUBARC END POINTS ON THE UNIT DISC,\r
+               C****   AS REQUIRED BY SUBSEQUENT PROCESSING ROUTINES.\r
+               C\r
+                       CALL BCFVTF(RSNPH(BCFSN),RGEOM(VTARG),ISNPH(DGPOL),ISNPH(JATYP),\r
+                    +  ISNPH(LOSUB),IGEOM(PARNT),RFARC,TNSUA,RSNPH(H0VAL),RSNPH(JACIN),\r
+                    +  RFARG,RSNPH(SOLUN))\r
+               C\r
+               C****   OUTPUT DATA REQUIRED FOR POST-PROCESSING.\r
+               C\r
+                       IGEOM(3)=TNSUA\r
+                       ISNPH(3)=TNSUA\r
+                       ISNPH(4)=NEQNS\r
+               C\r
+                       OFL=JBNM(1:L)//'gm'\r
+                       OPEN(OCH,FILE=OFL)\r
+                       CALL OUPTGM(IGEOM,RGEOM,CENTR,INTER,OCH)\r
+                       CLOSE(OCH)\r
+               C\r
+                       OFL=JBNM(1:L)//'ph'\r
+                       OPEN(OCH,FILE=OFL)\r
+                       CALL OUPTPH(ISNPH,RSNPH,OCH)\r
+                       CLOSE(OCH)\r
+                     ENDIF\r
+               C  \r
+                     CALL WRTAIL(1,0,IER)\r
+               C */\r
+    } // private void JAPHYC\r
+\r
+    private void ANGLE7(double BE[], int NA, boolean IN) {\r
+        // REAL BE(*)\r
+\r
+        // **** BE=JACIN,NA=NARCS,IN=INTER\r
+\r
+        // **** TO COMPUTE THE ARRAY OF JACOBI INDECES CORRESPONDING TO THE\r
+        // **** CORNER ANGLES ON THE BOUNDARY\r
+\r
+        // **** LOCAL VARIABLES\r
+        int K,B0,B1,B2;\r
+        double X,Y,ANG,PI,R1MACH,EPSA,XI,APP;\r
+        double U[] = new double[2];\r
+        double V[] = new double[2];\r
+        double DIN[] = new double[2];\r
+        double absu;\r
+        double absv;\r
+        double cr[] = new double[1];\r
+        double ci[] = new double[1];\r
+        //COMPLEX U,V,DPARFN\r
+        //EXTERNAL DPARFN\r
+\r
+        PI=Math.PI;\r
+        EPSA=Math.sqrt(EPS);\r
+        for (K=1; K <= NA; K++) {\r
+               DIN[0] = 1.0;\r
+               DIN[1] = 0.0;\r
+            U=DPARFN(K,DIN);\r
+            absu = zabs(U[0],U[1]);\r
+            U[0] = -U[0]/absu;\r
+            U[1] = -U[1]/absu;\r
+            DIN[0] = -1.0;\r
+               DIN[1] = 0.0;\r
+            if (K == NA) {\r
+                V=DPARFN(1,DIN);\r
+            }\r
+            else {\r
+                V=DPARFN(K+1,DIN);\r
+            }\r
+            absv = zabs(V[0],V[1]);\r
+            V[0] = V[0]/absv;\r
+            V[1] = V[1]/absv;\r
+            zdiv(U[0],U[1],V[0],V[1],cr,ci);\r
+            V[0]= cr[0];\r
+            V[1] = ci[0];\r
+            X = V[0];\r
+            Y = V[1];\r
+            ANG=Math.atan2(Y,X);\r
+            if (ANG < 0) {\r
+                ANG=ANG+2.0*PI;\r
+            }\r
+            ANG=ANG/PI;\r
+            if (!IN) {\r
+                ANG=2.0-ANG;\r
+            }\r
+            ANG=-1.0+1.0/ANG;\r
+\r
+            // ****   TRY TO DETECT SIMPLE RATIONAL INDECES AND FORCE BEST REAL\r
+            // ****   APPROXIMATIONS\r
+\r
+            if (Math.abs(ANG) < EPSA) {\r
+                ANG=0.0;\r
+            }\r
+            else {\r
+                XI=Math.abs(ANG);\r
+                B0=(int)(XI);\r
+                XI=XI-(double)(B0);\r
+                if (Math.abs(XI) < EPSA) {\r
+                    APP=(double)(B0);\r
+                }\r
+                else {\r
+                    XI=1.0/XI;\r
+                    B1=(int)(XI);\r
+                    XI=XI-(double)(B1);\r
+                    if (Math.abs(XI) < EPSA) {\r
+                        APP=(double)(1+B0*B1)/(double)(B1);\r
+                    }\r
+                    else {\r
+                        XI=1.0/XI;\r
+                        B2=(int)(XI);\r
+                        APP=(double)(B0*(1+B1*B2)+B2)/(double)(1+B1*B2);   \r
+                    } // else\r
+                } // else\r
+                if (ANG < 0.0) {\r
+                       APP = -APP;\r
+                }\r
+                if (Math.abs(ANG-APP) < EPSA) {\r
+                       ANG=APP;\r
+                }\r
+            } // else\r
+\r
+            if (K == NA) {\r
+                BE[0]=ANG;\r
+            }\r
+            else {\r
+                BE[K]=ANG;\r
+            }\r
+  \r
+        } // for (K=1; K <= NA; K++)\r
+\r
+    } // private void ANGLE7\r
\r
       /**\r
        * zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.\r
        * \r
@@ -2234,5 +3325,30 @@ C
 \r
           return;\r
       }\r
+      \r
+      /**\r
+       * complex divide c = a/b.\r
+       * \r
+       * @param ar double\r
+       * @param ai double\r
+       * @param br double\r
+       * @param bi double\r
+       * @param cr double[]\r
+       * @param ci double[]\r
+       */\r
+      public void zdiv(final double ar, final double ai, final double br, final double bi, final double[] cr,\r
+              final double[] ci) {\r
+          double bm, cc, cd, ca, cb;\r
+\r
+          bm = 1.0 / zabs(br, bi);\r
+          cc = br * bm;\r
+          cd = bi * bm;\r
+          ca = ( (ar * cc) + (ai * cd)) * bm;\r
+          cb = ( (ai * cc) - (ar * cd)) * bm;\r
+          cr[0] = ca;\r
+          ci[0] = cb;\r
+\r
+          return;\r
+      }\r
 \r
 }
\ No newline at end of file