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

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

index 69d4b93..bcb2924 100644 (file)
@@ -181,6 +181,36 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        private double P0VAL;\r
        private double SCALE;\r
        private int TYPE;\r
+       \r
+       // From GQPHYC:\r
+       // IQUPH  - INTEGER ARRAY\r
+       //                     AN INTEGER VECTOR OF SIZE AT LEAST 2*IBNDS(1) + 4,\r
+       //                     WHERE IBNDS(1) (=IGEOM(4)) IS THE VALUE PREVIOUSLY\r
+       //                     SUPPLIED TO JAPHYC; IQUPH STORES POINTERS TO \r
+       //                     ACCESS RQUPH AND ZQUPH.\r
+       private int IQUPH[] = new int[4];\r
+       private int LQSBF[] = new int[IBNDS[0]];\r
+       private int NPPQF[] = new int[IBNDS[0]];\r
+//  MQUPH  - INTEGER\r
+//                     THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
+//                     IN THE FINAL GLOBAL RULE.  (THE VALUE OF THIS\r
+//                     ARGUMENT IS LINKED TO THOSE OF ARGUMENTS NQPTS\r
+//                     AND IBNDS(1) PREVIOUSLY SUPPLIED TO JAPHYC VIA\r
+//                     MQUPH <= (MQIN1-1)*NQPTS*IBNDS(1))\r
+       private int MQUPH;\r
+       // RQUPH  - REAL ARRAY\r
+       //                     A REAL VECTOR OF SIZE AT LEAST 3*MQUPH + 1; STORES\r
+       //                     THE REAL QUADRATURE DATA.\r
+       private double RQUPH[] = new double[1];\r
+       private double TPPQF[] = new double[MQUPH];\r
+       private double TRRAD[] = new double[MQUPH];\r
+       private double WPPQF[] = new double[MQUPH];\r
+//  ZQUPH  - COMPLEX ARRAY\r
+//                     A COMPLEX VECTOR OF SIZE AT LEAST MQUPH + 1; \r
+//                     STORES THE QUADRATURE POINTS ON THE PHYSICAL \r
+//                     BOUNDARY.\r
+       private double FACTR[] = new double[2]; // At first location of ZQUPH\r
+       private double ZPPQF[][] = new double[MQUPH][2];\r
        public SymmsIntegralMapping() {\r
                \r
        }\r
@@ -2346,11 +2376,10 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        return result;\r
    }\r
 \r
-   private void JAPHYC(String JBNM, String HEAD, double MAXER,\r
+   private void JAPHYC(double MAXER,\r
                      int ISYGP, boolean INCST,\r
                      int RFARC, double RFARG[], int TSTNG,int OULVL,\r
                      double MATRX[][][],\r
-                     int OCH,\r
                      int IER[]) {\r
                \r
                      //INTEGER IBNDS(*),IGEOM(*),ISNPH(*),IWORK(*)\r
@@ -2660,9 +2689,10 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                //     LOCAL VARAIBLES\r
                \r
           int TNSUA[] = new int[1];\r
-          int I,IMXER,INDEG,J,L,MDGPO,MNJXS,MNQUA,MNSUA,MQIN1,NCOLL,\r
-                    NEFF,NEQNS,NROWS,NTEST,\r
-                    ORDSG;\r
+          int NTEST[] = new int[1];\r
+          int IMXER[] = new int[1];\r
+          int I,INDEG,J,MDGPO,MNJXS,MNQUA,MNSUA,MQIN1,NCOLL,\r
+                    NEFF,NEQNS,NROWS,ORDSG;\r
                int SOLCO = 0;\r
                int QIERC[] = new int[7];\r
                int QIERR[] = new int[7];\r
@@ -2675,8 +2705,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                double MQERR[] = new double[1];\r
                double RCOND[] = new double[1];\r
                double ESTOL[] = new double[1];\r
+               double MXERM[] = new double[1];\r
                double AQTOL,CONST,GAQTL,GLGTL,GRQTL,GSUPE,GTGTE,LGTOL,\r
-                    MCHEP,MXERM,PI,RQTOL,SSUPE,\r
+                    MCHEP,PI,RQTOL,SSUPE,\r
                     TGTER,TOLNR;\r
                double SOLUNB[][] = new double[MNEQN][1];\r
                double A[][] = new double[MNEQN][MNEQN];\r
@@ -2688,31 +2719,12 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                boolean ACCPT[] = new boolean[1];\r
                boolean NUQTL[] = new boolean[1];\r
                boolean GACPT,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
@@ -2823,7 +2835,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                //**** LIST THE INPUT ARGUMENTS AND ASSOCIATED QUANTITIES\r
                \r
-               RSLT80(JBNM,HEAD,GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
+               RSLT80(GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
                     INCST,INDEG,RFARC,RFARG[0],CENTR,JACIN,LNSEG,\r
                     TSTNG,OULVL,IBNDS,MNEQN);\r
                PI=Math.PI;\r
@@ -3277,11 +3289,11 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                System.out.println();\r
                System.out.println("ERRORS IN MODULUS STARTED:");\r
                \r
-               /*TSJAC3(LOTES,HITES,COLPR,ZCOLL,\r
+               TSJAC3(LOTES,HITES,COLPR,ZCOLL,\r
                  NQPTS,NTEST,ORDSG,TNSUA[0],TSTNG,DGPOL,JATYP,\r
                  PARNT,AICOF,BICOF,DIAG,HALEN,\r
                  JACIN,MIDPT,SDIAG,IER);\r
-                 if (IER[0] > 0) {\r
+           if (IER[0] > 0) {\r
                          if (SOLCO >= 1) {\r
                                        \r
                                    // ****   COMPUTE THE BOUNDARY CORRESPONDENCE COEFFICIENTS BCFSN AND THE\r
@@ -3309,7 +3321,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                          }\r
                  }\r
                  TESMD9(WORK2,A,SOLUN,MNEQN,NCOLL,\r
-                    NTEST,NQPTS,TNSUA[0],JATYP,PARNT,DGPOL,\r
+                    NTEST[0],NQPTS,TNSUA[0],JATYP,PARNT,DGPOL,\r
                     LOSUB,HISUB,LOTES,HITES,NQUAD,\r
                     LOQSB,TOLNR,MIDPT,HALEN,H0VAL,\r
                     COLSC,ACOEF,BCOEF,COLPR,QCOMX,\r
@@ -3317,57 +3329,48 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                     QIERR,MQERR,JACIN,A1COF,B1COF,AQTOL,RQTOL,\r
                     AQCOF,BQCOF,CQCOF,MXERM,IMXER,ZMXER,\r
                     ERARC,ORDSG,REFLN);\r
+                 for (I = 0; I < MNEQN; I++) {\r
+                         for (J = 0; J < MNEQN; J++) {\r
+                                 MATRX[I][J][1] = A[I][J];\r
+                         }\r
+                 }\r
                \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
+                 System.out.println("ERRORS IN MODULUS DONE:");\r
+                 System.out.println("MAXIMUM ERROR AT TEST POINTS: " + MXERM[0]);\r
+                 for (I=0; I <= 6; I++) {\r
+                     QIERC[I]=QIERC[I]+QIERR[I];\r
+                 }\r
+                 if ((OULVL%2) == 1) {\r
+                     RSLT84(WORK2,TNSUA[0],MXERM[0],ZMXER,IMXER[0],LOTES,\r
+                     HITES,QIERC,PARNT,ORDSG);\r
+                 }\r
+                 else {\r
+                     RSLT83(ERARC,TNSUA[0],MXERM[0],ZMXER,IMXER[0],QIERC,\r
+                     PARNT,ORDSG);\r
+                 }\r
+               \r
+                 // RESCALE SOLUTIONS TO OBTAIN STANDARD JACOBI COEFFICIENTS\r
+                 RESCAL(NQPTS,TNSUA[0],LOSUB,HISUB,JATYP,SOLUN,COLSC);\r
                  \r
-                     CALL WRTAIL(1,0,IER)*/\r
+                 if (SOLCO >= 1) {\r
+                               \r
+                           // ****   COMPUTE THE BOUNDARY CORRESPONDENCE COEFFICIENTS BCFSN AND THE\r
+                           // ****   ARGUMENTS OF ALL SUBARC END POINTS ON THE UNIT DISC,\r
+                           // ****   AS REQUIRED BY SUBSEQUENT PROCESSING ROUTINES.\r
+                       \r
+                           BCFVTF(BCFSN,VTARG,DGPOL,JATYP,\r
+                               LOSUB,PARNT,RFARC,TNSUA[0],H0VAL,JACIN,\r
+                               RFARG[0],SOLUN);\r
+                       \r
+                           // ****   OUTPUT DATA REQUIRED FOR POST-PROCESSING.\r
+                       \r
+                           IGEOM[2]=TNSUA[0];\r
+                           ISNPH[2]=TNSUA[0];\r
+                           ISNPH[3]=NEQNS;\r
+               } // if (SOLCO >= 1)\r
+                         \r
+               WRTAIL(1,0,IER[0],null);\r
+                   return;     \r
                \r
     } // private void JAPHYC\r
 \r
@@ -3472,7 +3475,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 \r
     } // private void ANGLE7\r
     \r
-    private void RSLT80(String JBNM, String HEAD, double SUPER, double MAXER, double AQTOL, boolean INTER, int NARCS, int ORDSG,\r
+    private void RSLT80(double SUPER, double MAXER, double AQTOL, boolean INTER, int NARCS, int ORDSG,\r
         int NQPTS, boolean INCST, int INDEG, int RFARC, double RFARG,double CENTR[], double BETA[], boolean LINEAR[],\r
         int TSTNG, int OULVL, int IBNDS[], int MNEQN) {\r
        \r
@@ -3489,10 +3492,6 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
        int I;\r
        final String TXT1 = "REQUESTED ACCURACY UNREALISTIC";\r
-\r
-\r
-       Preferences.debug("HEAD\n",Preferences.DEBUG_ALGORITHM);\r
-        Preferences.debug("JOB NAME : " + JBNM + "\n",Preferences.DEBUG_ALGORITHM);\r
        \r
        if (INTER) {\r
            Preferences.debug("INTERIOR DOMAIN WITH " + NARCS + " ARCS\n",Preferences.DEBUG_ALGORITHM);\r
@@ -7329,7 +7328,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 \r
 \r
     private void TSJAC3(int LOTES[], int HITES[], double TESPR[], double ZTEST[][], int NQPTS,\r
-        int NTEST, int ORDSG, int TNSUA, int TSTNG, int DGPOL[], int JATYP[], int PARNT[],\r
+        int NTEST[], int ORDSG, int TNSUA, int TSTNG, int DGPOL[], int JATYP[], int PARNT[],\r
         double AICOF[], double BICOF[], double DIAG[], double HALEN[], double JACIN[],\r
         double MIDPT[], double SDIAG[], int IER[]) {\r
        \r
@@ -7429,18 +7428,18 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                } // if (D > 0)\r
            } // for (I=1; I <= TSFBS; I++) \r
         } // else TSTNG == 1\r
-       NTEST=HITES[TSFBS-1];\r
+       NTEST[0]=HITES[TSFBS-1];\r
        if (ORDSG > 1) {\r
-           NTEST=NTEST+1;\r
-           LOTES[TSFBS]=NTEST;\r
-           HITES[TSFBS]=NTEST;\r
-           TESPR[NTEST-1]=-1.0;\r
+           NTEST[0]=NTEST[0]+1;\r
+           LOTES[TSFBS]=NTEST[0];\r
+           HITES[TSFBS]=NTEST[0];\r
+           TESPR[NTEST[0]-1]=-1.0;\r
            P=PARNT[TSFBS];\r
            PIN[0] = -1.0;\r
            PIN[1] = 0.0;\r
            POUT = PARFUN(P,PIN);\r
-           ZTEST[NTEST-1][0] = POUT[0];\r
-           ZTEST[NTEST-1][1] = POUT[1];\r
+           ZTEST[NTEST[0]-1][0] = POUT[0];\r
+           ZTEST[NTEST[0]-1][1] = POUT[1];\r
        }\r
        \r
        IER[0]=0;\r
@@ -7454,8 +7453,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         double BCOEF[], double TESPR[], double QCOMX[], double QCOMW[], double CENTR[], \r
         double ZTEST[][], boolean INTER, boolean LNSEG[], double WORK[], int QIERR[],\r
         double MQERR[], double JACIN[], double A1COF[], double B1COF[], double AQTOL,\r
-        double RQTOL, double AQCOF[], double BQCOF[], double CQCOF[], double MXERM,\r
-        int IMXER, double ZMXER[], double ERARC[], int ORDSG, boolean REFLN) {\r
+        double RQTOL, double AQCOF[], double BQCOF[], double CQCOF[], double MXERM[],\r
+        int IMXER[], double ZMXER[], double ERARC[], int ORDSG, boolean REFLN) {\r
        \r
        //INTEGER MNEQN,NCOLL,NTEST,NQPTS,TNSUA,JATYP(*),PARNT(*),DGPOL(*),\r
        //+LOSUB(*),HISUB(*),LOTES(*),HITES(*),NQUAD(*),LOQSB(*),QIERR(0:6),\r
@@ -7486,6 +7485,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        // COMPLEX GTQ,ZQ,ZT,PARFUN,DPARFN\r
        // COMMON /FNDEF/BETA,A1,B1,P0VAL,SCALE,TYPE\r
        // EXTERNAL DPARFN,FNVAL,JAPAR7,PARFUN,QAWS,R1MACH\r
+       double PIN[] = new double[2];\r
+       double POUT[];\r
        qaws qmod;\r
        int routine = Integration2.DQAWSE;\r
        /**\r
@@ -7509,7 +7510,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
        // NOW SELECT THE INTEGRATION ARC\r
        \r
-       /*for (IA=1; IA <= TNSUA; IA++) {\r
+       for (IA=1; IA <= TNSUA; IA++) {\r
 \r
        \r
            // INITIALISE DATA FOR THIS ARC\r
@@ -7657,247 +7658,1064 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                            IER[0] = qmod.getErrorStatus();\r
                            QIERR[IER[0]]=QIERR[IER[0]]+1;\r
                    }\r
-                     WORK(2)=U0+U1\r
-                     MQERR=MAX(MQERR,ABER0+ABER1)\r
-       C\r
-       C             NOW USE THE (WEAKLY) STABLE FORWARD RECURRENCE SCHEME FOR \r
-       C             WORK(I),I=3,NQPTS (WITH SCALE FACTOR FOR WORK(2))\r
+                   WORK[1]=U0+U1;\r
+                   MQERR[0]=Math.max(MQERR[0],ABER0+ABER1);\r
+       \r
+                   // NOW USE THE (WEAKLY) STABLE FORWARD RECURRENCE SCHEME FOR \r
+                   // WORK(I),I=3,NQPTS (WITH SCALE FACTOR FOR WORK(2))\r
+       \r
+                   CURR=WORK[1];\r
+                   PREV=SCALE;\r
+                   for (I=1; I <= DG-1; I++) {\r
+                       J=LOD+I-1;\r
+                       NEXT=(AQCOF[J-1]*TT-BQCOF[J-1])*CURR-CQCOF[J-1]*PREV;\r
+                       WORK[I+1]=NEXT;\r
+                       PREV=CURR;\r
+                       CURR=NEXT;\r
+                   } // for (I=1; I <= DG-1; I++)\r
+       \r
+                   // ASSIGN CORRECT SCALE FACTORS.\r
+       \r
+                   for (I=3; I <= DG+1; I++) {\r
+                       J=LOD+I-1;\r
+                       WORK[I-1]=WORK[I-1]*COLSC[J-1]/SCALE;\r
+                   } // for (I=3; I <= DG+1; I++)\r
+               } // if (DG > 0)\r
+       \r
+               SG=1.0;\r
+               for (J=LOCLM; J <= HICLM; J++) {\r
+                   MATRX[ROW-1][J-1]=MATRX[ROW-1][J-1]+SG*WORK[J-LOCLM];\r
+                   SG=SG*SJI;\r
+               } // for (J=LOCLM; J <= HICLM; J++)\r
+       \r
+           } // for (IT=FIRST; IT <= LAST; IT++)\r
+       \r
+           // INITIALISE SOME DATA FOR THE NON-SINGULAR INTEGRALS\r
+       \r
+           WORK[0]=1.0/RH;\r
+           NQ=NQUAD[AJI-1];\r
+           LOQ=LOQSB[AJI-1];\r
+           if (IA == TNSUA) {\r
+               I1=2;\r
+           }\r
+           else {\r
+               I1=1;\r
+           }\r
+       \r
+           for (IQ=1; IQ <= NQ; IQ++) {\r
+               I=LOQ+IQ-1;\r
+               TQ=QCOMX[I-1];\r
+               WQ=QCOMW[I-1];\r
+               double AA[] = new double[DG];\r
+               double BB[] = new double[DG];\r
+               for (I = 0; I < DG; I++) {\r
+                       AA[I] = ACOEF[LOD+I-1];\r
+                       BB[I] = BCOEF[LOD+I-1];\r
+               }\r
+               JAPAR7(WORK,TQ,AA,BB,DG);\r
+               if (JI < 0) {\r
+                   TQ=-TQ;\r
+               }\r
+               GTQ[0] = MD+HL*TQ;\r
+               GTQ[1] = 0.0;\r
+               ZQ=PARFUN(PT,GTQ);\r
+       \r
+               // ACCUMULATE THE ELEMENTS ABOVE THE DIAGONAL BLOCK\r
+       \r
+       \r
+               for (IT=I1; IT <= FIRST-1; IT++) {\r
+                   ZT[0]=ZTEST[IT-1][0];\r
+                   ZT[1]=ZTEST[IT-1][1];\r
+                   DD=zabs(ZT[0]-ZQ[0],ZT[1]-ZQ[1]);\r
+                   LDD=Math.log(DD)*WQ;\r
+                   SG=1.0;\r
+                   for (J1=LOCLM; J1 <= HICLM; J1++) {\r
+                       J=J1-LOCLM+1;\r
+                       I=J1-LOCLM+LOD;\r
+                       MATRX[IT-1][J1-1]=MATRX[IT-1][J1-1]+SG*WORK[J-1]*LDD*COLSC[I-1];\r
+                       SG=SG*SJI;\r
+                   } // for (J1=LOCLM; J1 <= HICLM; J1++)\r
+               } // for (IT=I1; IT <= FIRST-1; IT++)\r
+       \r
+               // ACCUMULATE THE ELEMENTS BELOW THE DIAGONAL BLOCK\r
+       \r
+               for (IT=LAST+1; IT <= NTEST; IT++) {\r
+                   ZT[0]=ZTEST[IT-1][0];\r
+                   ZT[1]=ZTEST[IT-1][1];\r
+                   DD=zabs(ZT[0]-ZQ[0],ZT[1]-ZQ[1]);\r
+                   LDD=Math.log(DD)*WQ;\r
+                   SG=1.0;\r
+                   for (J1=LOCLM; J1 <= HICLM; J1++) {\r
+                       J=J1-LOCLM+1;\r
+                       I=J1-LOCLM+LOD;\r
+                       MATRX[IT-1][J1-1]=MATRX[IT-1][J1-1]+SG*WORK[J-1]*LDD*COLSC[I-1];\r
+                       SG=SG*SJI;\r
+                   } // for (J1=LOCLM; J1 <= HICLM; J1++)\r
+               } // for (IT=LAST+1; IT <= NTEST; IT++)\r
+       \r
+                // ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
+               // THE DIAGONAL BLOCK FOR THE NON-LINE-SEGMENT CASE.\r
+       \r
+               if (!LNSEG[IA-1]) {\r
+                   for (IT=FIRST; IT <= LAST; IT++) {\r
+                       if (ORDSG > 1 && IA == (TSFBS+1) && IT == LAST) { \r
+                           continue;\r
+                       }\r
+                       if (IT == LAST) {\r
+                           TT=1.0;\r
+                       }\r
+                       else {\r
+                           TT=TESPR[IT-1];\r
+                       }\r
+       \r
+                       if (IT > NTEST) {\r
+                           ROW=1;\r
+                       }\r
+                       else {\r
+                           ROW=IT;\r
+                       }\r
+       \r
+                       DD=Math.abs(TT-TQ);\r
+                       if (DD <= TOLNR) {\r
+                               POUT = DPARFN(PT,GTQ);\r
+                           DD=zabs(POUT[0],POUT[1])*HL;\r
+                       }\r
+                       else {\r
+                           if (IT > NTEST) {\r
+                               ZT[0]=ZTEST[0][0];\r
+                               ZT[1]=ZTEST[0][1];\r
+                           }\r
+                           else {\r
+                               ZT[0]=ZTEST[IT-1][0];\r
+                               ZT[1]=ZTEST[IT-1][1];\r
+                           }\r
+                           DD=zabs(ZT[0]-ZQ[0],ZT[1]-ZQ[1])/DD;\r
+                           if (DD < TOLNR) {\r
+                               POUT = DPARFN(PT,GTQ);\r
+                               DD=zabs(POUT[0],POUT[1])*HL;\r
+                           }\r
+                       }\r
+                       LDD=Math.log(DD)*WQ;\r
+                       SG=1.0;\r
+                       for (J1=LOCLM; J1 <= HICLM; J1++) {\r
+                           J=J1-LOCLM+1;\r
+                           I=J1-LOCLM+LOD;\r
+                           MATRX[ROW-1][J1-1]=MATRX[ROW-1][J1-1]+SG*WORK[J-1]*LDD*COLSC[I-1];\r
+                           SG=SG*SJI;\r
+                       } // for (J1=LOCLM; J1 <= HICLM; J1++)\r
+                   } // for (IT=FIRST; IT <= LAST; IT++)\r
+               } // if (!LNSEG[IA-1])\r
+           } // for (IQ=1; IQ <= NQ; IQ++)\r
+       \r
+           // ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
+           // THE DIAGONAL BLOCK FOR THE LINE-SEGMENT CASE.\r
+       \r
+           if (LNSEG[IA-1]) {\r
+               PIN[0] = MD;\r
+               PIN[1] = 0.0;\r
+               POUT = DPARFN(PT,PIN);\r
+               ZT[0] = POUT[0]*HL;\r
+               ZT[1] = POUT[1]*HL;\r
+               C0=zabs(ZT[0],ZT[1]);\r
+               C0=RH*Math.log(C0)*COLSC[LOD-1];\r
+               for (IT=FIRST; IT <= LAST; IT++) {\r
+                   if (ORDSG > 1 && IA == (TSFBS+1) && IT == LAST) {\r
+                       continue;\r
+                   }\r
+                   if (IT > NTEST) {\r
+                       ROW=1;\r
+                   }\r
+                   else {\r
+                       ROW=IT;\r
+                   }\r
+                   MATRX[ROW-1][LOCLM-1]=MATRX[ROW-1][LOCLM-1]+C0;\r
+               } // for (IT=FIRST; IT <= LAST; IT++)\r
+           } // if (LNSEG[IA-1])\r
+       } // for (IA=1; IA <= TNSUA; IA++)\r
+       \r
+       // SET UP THE LAST COLUMN\r
+       \r
+       for (I=1; I <= NTEST; I++) {\r
+           MATRX[I-1][NEQNS-1]=1.0;\r
+       }\r
+       \r
+       // COMPUTE MATRIX-VECTOR PRODUCT\r
+       \r
+       for (I=1; I <= NTEST; I++) {\r
+           SUM=0.0;\r
+           for (J=1; J <= NEQNS; J++) {\r
+                 SUM=SUM+MATRX[I-1][J-1]*SOLUN[J-1];\r
+           }\r
+           ERMOD[I-1]=SUM;\r
+       } // for (I=1; I <= NTEST; I++)\r
+       \r
+       // FORM THE ERROR IN MODULUS\r
+       \r
+       if (INTER) {\r
+           for (I=1; I <= NTEST; I++) {\r
+               SUM=Math.exp(ERMOD[I-1]);\r
+               ERMOD[I-1]=Math.abs(1.0-zabs(ZTEST[I-1][0]-CENTR[0],ZTEST[I-1][1]-CENTR[1])/SUM);\r
+           }\r
+       }\r
+       else {\r
+           for (I=1; I <= NTEST; I++) {\r
+               SUM=Math.exp(ERMOD[I-1]);\r
+               ERMOD[I-1]=Math.abs(1.0-SUM);\r
+           }\r
+       }\r
+       \r
+       // FIND MAXIMUM ERROR IN MODULUS AND THE POINT AND THE ARC AT WHICH \r
+       // IT OCCURS\r
+       \r
+       MXERM[0]=0.0;\r
+           for (IA=1; IA <= TSFBS; IA++) {\r
+               FIRST=LOTES[IA-1];\r
+               LAST=HITES[IA-1];\r
+               MD=0.0;\r
+               for (IT=FIRST; IT <= LAST; IT++) {\r
+                   if (ERMOD[IT-1] > MD) {\r
+                       MD=ERMOD[IT-1];\r
+                   }\r
+                   if (MD > MXERM[0]) {\r
+                       MXERM[0]=MD;\r
+                       IMXER[0]=IA;\r
+                       ZMXER[0]=ZTEST[IT-1][0];\r
+                       ZMXER[1]=ZTEST[IT-1][1];\r
+                   }\r
+               } // for (IT=FIRST; IT <= LAST; IT++)\r
+               if (IA == TSFBS && ORDSG == 1) {\r
+                   IT=1;\r
+               }\r
+               else {\r
+                   IT=LAST+1;\r
+               }\r
+               if (ERMOD[IT-1] > MD) {\r
+                   MD=ERMOD[IT-1];\r
+               }\r
+               if (MD > MXERM[0]) {\r
+                   MXERM[0]=MD;\r
+                   IMXER[0]=IA;\r
+                   ZMXER[0]=ZTEST[IT-1][0];\r
+                   ZMXER[1]=ZTEST[IT-1][1];\r
+               } // if (MD > MXERM[0])\r
+               ERARC[IA-1]=MD;\r
+           } // for (IA=1; IA <= TSFBS; IA++)\r
+       \r
+           // IF REGION IS SYMMETRIC, FILL UP THE WHOLE ERARC VECTOR USING\r
+           // SYMMETRY\r
+       \r
+           if (ORDSG > 1) {\r
+               if (REFLN) {\r
+                   for (IA=1; IA <= TSFBS; IA++) {\r
+                       IB=2*TSFBS+1-IA;\r
+                       ERARC[IB-1]=ERARC[IA-1];\r
+                   } // for (IA=1; IA <= TSFBS; IA++)\r
+                   ORDRG=ORDSG/2;\r
+                   for (I=2; I <= ORDRG; I++) {\r
+                       I1=(I-1)*TSFBS*2;\r
+                       for (IA=1; IA <= TSFBS*2; IA++) {\r
+                           IB=I1+IA;\r
+                           ERARC[IB-1]=ERARC[IA-1];\r
+                       } // for (IA=1; IA <= TSFBS*2; IA++)\r
+                   } //  for (I=2; I <= ORDRG; I++)\r
+               } // if (REFLN)\r
+               else {\r
+                   for (I=2; I <= ORDSG; I++) {\r
+                       I1=(I-1)*TSFBS;\r
+                       for (IA=1; IA <= TSFBS; IA++) {\r
+                           IB=I1+IA;\r
+                           ERARC[IB-1]=ERARC[IA-1];\r
+                       } // for (IA=1; IA <= TSFBS; IA++)\r
+                   } // for (I=2; I <= ORDSG; I++)\r
+               } // else\r
+           } // if (ORDSG > 1)\r
+\r
+    } // private void TESMD9\r
+    \r
+    private void RSLT83(double ERARC[], int TNSUA, double MXERM, double ZMXER[], int IMXER,\r
+        int QIERC[], int PARNT[], int ORDSG) {\r
+       // INTEGER TNSUA,IMXER,ORDSG,OC\r
+       // INTEGER PARNT(*),QIERC(0:6)\r
+       // REAL MXERM,ERARC(*)\r
+       // COMPLEX ZMXER\r
+       \r
+       // LOCAL VARIABLES\r
+       \r
+       int I,TSFBS;\r
+       String QTEXT[] = new String[7];\r
+        //CHARACTER QTEXT(0:6)*22,LINE*72\r
+       final String LINE="_________________________________________________________________";\r
+       \r
+       QTEXT[0]="...........NORMAL EXIT";\r
+       QTEXT[1]=".....MAX. SUBDIVISIONS";\r
+       QTEXT[2]="....ROUNDOFF DETECTION";\r
+       QTEXT[3]=".........BAD INTEGRAND";\r
+       QTEXT[6]=".........INVALID INPUT";\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug(LINE+"\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("RESULTS FROM ERROR IN MODULUS TESTS",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("QAWS TERMINATIONS WITH......\n",Preferences.DEBUG_ALGORITHM); \r
+       for (I=0; I <= 6; I++) {\r
+           if (QIERC[I] > 0) {\r
+               Preferences.debug(QTEXT[I]+ " " + QIERC[I] + "\n",Preferences.DEBUG_ALGORITHM); \r
+           }\r
+       } // for (I=0; I <= 6; I++)\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("SUB ARC  PARENT ARC     MAX ERROR IN MODULUS\n",Preferences.DEBUG_ALGORITHM); \r
+       TSFBS=TNSUA/ORDSG;\r
+       for (I=1; I <= TSFBS; I++) {\r
+               Preferences.debug(" " + I + "       " + PARNT[I-1] + "          " + ERARC[I-1] + "\n",Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       \r
+       Preferences.debug("MAXIMUM ERROR IN MODULUS IS " + MXERM + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("THIS OCCURS AT ("+ZMXER[0]+","+ZMXER[1]+") ON ARC " + IMXER + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+    } // private void RSLT83\r
+\r
+\r
+    private void RSLT84(double ERMOD[], int TNSUA, double MXERM, double ZMXER[], int IMXER,\r
+        int LOTES[], int HITES[], int QIERC[], int PARNT[], int ORDSG) {\r
+       // INTEGER TNSUA,IMXER,ORDSG,OC\r
+       // INTEGER LOTES(*),HITES(*),PARNT(*),QIERC(0:6)\r
+       // REAL MXERM,ERMOD(*)\r
+       // COMPLEX ZMXER\r
+       \r
+       // LOCAL VARIABLES\r
+        \r
+       int I,J,L,H,TSFBS;\r
+       double S;\r
+        String QTEXT[] = new String[7];\r
+        //CHARACTER QTEXT(0:6)*22,LINE*72\r
+       final String LINE="_________________________________________________________________";\r
+       \r
+       QTEXT[0]="...........NORMAL EXIT";\r
+       QTEXT[1]=".....MAX. SUBDIVISIONS";\r
+       QTEXT[2]="....ROUNDOFF DETECTION";\r
+       QTEXT[3]=".........BAD INTEGRAND";\r
+       QTEXT[6]=".........INVALID INPUT";\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug(LINE+"\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("RESULTS FROM ERROR IN MODULUS TESTS",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("QAWS TERMINATIONS WITH......\n",Preferences.DEBUG_ALGORITHM); \r
+       for (I=0; I <= 6; I++) {\r
+           if (QIERC[I] > 0) {\r
+               Preferences.debug(QTEXT[I]+ " " + QIERC[I] + "\n",Preferences.DEBUG_ALGORITHM); \r
+           }\r
+       } // for (I=0; I <= 6; I++)\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("ERROR IN MODULI AT TEST POINTS\n",Preferences.DEBUG_ALGORITHM);\r
+       TSFBS=TNSUA/ORDSG;\r
+       for (I=1; I <= TSFBS; I++) {\r
+           Preferences.debug("SUB ARC = " + I + " ON PARENT ARC " + PARNT[I-1] + "\n",Preferences.DEBUG_ALGORITHM);\r
+           Preferences.debug("N         ERROR IN MOD\n",Preferences.DEBUG_ALGORITHM);\r
+           L=LOTES[I-1];\r
+           H=HITES[I-1];\r
+           for (J=L; J <= H; J++) {\r
+               S=ERMOD[J-1];\r
+               Preferences.debug(J + "   " + S + "\n",Preferences.DEBUG_ALGORITHM);\r
+           } // for (J=L; J <= H; J++)\r
+           if (I == TSFBS && ORDSG == 1) {\r
+               J=1;\r
+           }\r
+           else {\r
+               J=H+1;\r
+           }\r
+           S=ERMOD[J-1];\r
+           Preferences.debug(J + "   " + S + "\n",Preferences.DEBUG_ALGORITHM);\r
+       } // for (I=1; I <= TSFBS; I++)\r
+       \r
+       Preferences.debug("MAXIMUM ERROR IN MODULUS IS " + MXERM + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("THIS OCCURS AT ("+ZMXER[0]+","+ZMXER[1]+") ON ARC " + IMXER + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+    } // private void RSLT84\r
+    \r
+    private void RESCAL(int NQPTS,int TNSUA, int LOSUB[], int HISUB[], int JATYP[],\r
+               double SOLUN[], double COLSC[]) {\r
+        // INTEGER NQPTS,TNSUA\r
+        // INTEGER LOSUB(*),HISUB(*),JATYP(*)\r
+        // REAL SOLUN(*),COLSC(*)\r
+\r
+        // TO RESCALE THE SOLUTION VECTOR SOLUN SO AS TO OBTAIN THE \r
+        // STANDARD JACOBI COEFFICIENTS.\r
+\r
+        // LOCAL VARIABLES\r
+\r
+        int H,I,J,J1,JT,L,LOD;\r
+\r
+        for (I=1; I <= TNSUA; I++) {\r
+            L=LOSUB[I-1];\r
+            H=HISUB[I-1];\r
+            JT=Math.abs(JATYP[I-1]);\r
+            LOD=(JT-1)*NQPTS+1;\r
+            for (J=L; J <= H; J++) {\r
+                J1=LOD+J-L;\r
+                SOLUN[J-1]=SOLUN[J-1]*COLSC[J1-1];\r
+            } // for (J=L; J <= H; J++)\r
+        } // for (I=1; I <= TNSUA; I++)\r
+\r
+    } // private void RESCAL\r
+\r
+    private void GQPHYC(int MQIN1, int IER[]) {\r
+       \r
+       // INTEGER MQIN1,MQUPH,CHNL,IER\r
+       // INTEGER IQUPH(*),IGEOM(*),ISNPH(*)\r
+       // REAL RQUPH(*),RGEOM(*),RSNPH(*),RWORK(*)\r
+       // COMPLEX CENTR\r
+       // COMPLEX ZQUPH(*)\r
+       // LOGICAL INTER\r
+       \r
+       // ......................................................................\r
+       //\r
+       // 1.     GQPHYC\r
+       //           COMPUTES A GLOBAL QUADRATURE RULE FOR APPROXIMATING THE\r
+       //           BOUNDARY INTEGRAL REPRESENTATION OF THE MAP : PHYSICAL -->\r
+       //           CANONICAL.\r
+       \r
+       \r
+       // 2.     PURPOSE\r
+       //           THE ROUTINE SETS UP THE BOUNDARY QUADRATURE POINTS AND\r
+       //           CORRESPONDING WEIGHTS FOR A COMPOSITE GAUSS-JACOBI/GAUSS-\r
+       //           LEGENDRE RULE FOR ESTIMATING THE BOUNDARY INTEGRAL THAT\r
+       //           APPEARS IN THE RESPRESENTATION FOR THE CONFORMAL MAP OF THE\r
+       //           PHYSICAL DOMAIN ONTO THE CANONICAL DOMAIN.  THIS QUADRATURE\r
+       //           RULE IS USED IN THE STANDARD NON-SINGULAR CASE WHEN THE\r
+       //           FIELD POINT IN THE PHYSICAL DOMAIN DOES NOT LIE CLOSE TO THE\r
+       //           BOUNDARY. \r
+                  \r
+       \r
+       // 3.     CALLING SEQUENCE\r
+       //           CALL GQPHYC(MQIN1,MQUPH,INTER,CENTR,IGEOM,RGEOM,ISNPH,RSNPH,\r
+       //                       RWORK,CHNL,IQUPH,RQUPH,ZQUPH,IER)\r
+       \r
+       //        PARAMETERS\r
+       //         ON ENTRY\r
+       //            MQIN1  - INTEGER\r
+       //                     DEFINES THE NUMBER OF PANELS ALLOWED IN A \r
+       //                     COMPOSITE RULE.  SPECIFICALLY, MQIN1 = 1 + (THE\r
+       //                     MAXIMUM NUMBER OF PANELS IN A COMPOSITE RULE FOR\r
+       //                     A SINGLE SUB-ARC ON THE BOUNDARY)\r
+       \r
+       //            MQUPH  - INTEGER\r
+       //                     THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
+       //                     IN THE FINAL GLOBAL RULE.  (THE VALUE OF THIS\r
+       //                     ARGUMENT IS LINKED TO THOSE OF ARGUMENTS NQPTS\r
+       //                     AND IBNDS(1) PREVIOUSLY SUPPLIED TO JAPHYC VIA\r
+       //                     MQUPH <= (MQIN1-1)*NQPTS*IBNDS(1))\r
+       \r
+       //            INTER  - LOGICAL\r
+       //                     TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE \r
+       //                     OTHERWISE.\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
+       \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
+       //            RWORK  - REAL ARRAY\r
+       //                     A WORKING VECTOR OF SIZE AT LEAST MQIN1.\r
+       \r
+       //             CHNL  - INTEGER\r
+       //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+       //                     WRITING THE FILE <JBNM>pq.\r
+       \r
+       //         ON EXIT\r
+       //            IQUPH  - INTEGER ARRAY\r
+       //                     AN INTEGER VECTOR OF SIZE AT LEAST 2*IBNDS(1) + 4,\r
+       //                     WHERE IBNDS(1) (=IGEOM(4)) IS THE VALUE PREVIOUSLY\r
+       //                     SUPPLIED TO JAPHYC; IQUPH STORES POINTERS TO \r
+       //                     ACCESS RQUPH AND ZQUPH.\r
+       \r
+       //            RQUPH  - REAL ARRAY\r
+       //                     A REAL VECTOR OF SIZE AT LEAST 3*MQUPH + 1; STORES\r
+       //                     THE REAL QUADRATURE DATA.\r
+       \r
+       //            ZQUPH  - COMPLEX ARRAY\r
+       //                     A COMPLEX VECTOR OF SIZE AT LEAST MQUPH + 1; \r
+       //                     STORES THE QUADRATURE POINTS ON THE PHYSICAL \r
+       //                     BOUNDARY.\r
+                            \r
+       //            IER    - INTEGER\r
+       //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+       //                     A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
+       //                     WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+       //                     IER=0 - NORMAL EXIT.\r
+       //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+       //                             BE SELF EXPLANATORY.\r
+       \r
+                           \r
+       // 4.     SUBROUTINES OR FUNCTIONS NEEDED\r
+       //              - THE CONFPACK LIBRARY.\r
+       //              - THE REAL FUNCTION R1MACH.\r
+       //              - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN.\r
+       \r
+       \r
+       // 5.     FURTHER COMMENTS\r
+       //           - NOTE THAT THIS ROUTINE CAN ONLY BE USED  A F T E R  THE  \r
+       //             ROUTINE JAPHYC HAS SUCCESSFULLY EXECUTED, AND THAT SOME  \r
+       //             INPUT ARGUMENTS FOR GQPHYC ARE OUTPUT VALUES FROM JAPHYC.\r
+       //           - THE GLOBAL QUADRATURE DATA ARE AUTOMATICALLY OUTPUT ON THE\r
+       //             FILE <JBNM>pq, WHERE <JBNM> IS COLLECTED FROM THE FILE \r
+       //             jbnm PREVIOUSLY 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
+       final int NINTS = 5;\r
+       int ACOEF,BCOEF,DGPOL,FACTR,H0VAL,HALEN,JACIN,JATYP,LOSUB,\r
+            LQSBF,MIDPT,NPPQF,PARNT,QUPTS,QUWTS,SOLUN,TPPQF,TRRAD,VTARG,WPPQF,\r
+            ZPPQF;\r
+       int MNEQN,MNSUA,NARCS,NEQNS,NJIND,NQPTS,TNGQP,TNPQP,\r
+            TNSUA;\r
+       final double DLETA = 0.2;\r
+       double LGTOL,SUPER,THET0,TOLOU;\r
+       double CT[] = new double[2];\r
+       //COMPLEX CT\r
+       \r
+        // EXTERNAL POPQF1,OUPTPQ,WRHEAD,WRTAIL\r
+       \r
+       // WRITE HEADING TO STANDARD OUTPUT CHANNEL\r
+       \r
+       WRHEAD(2,0,null);\r
+       \r
+       NARCS=ISNPH[0];\r
+       NQPTS=ISNPH[1];\r
+       TNSUA=ISNPH[2];\r
+       NEQNS=ISNPH[3];\r
+       MNSUA=ISNPH[4];\r
+       MNEQN=ISNPH[5];\r
+       NJIND=NARCS+1;\r
+       TNGQP=NJIND*NQPTS;\r
+       SUPER=RGEOM[0];\r
+       LGTOL=RGEOM[1];\r
+       \r
+       IQUPH[1]=TNSUA;\r
+       IQUPH[2]=MNSUA;\r
+       IQUPH[3]=MQUPH;\r
+       \r
+       // COPY POINTERS FROM JAPHYC\r
+       \r
+       // SET UP POINTERS FOR QUADRATURE DATA.\r
+       \r
+       \r
+       /*RQUPH[0]=SOLUN[NEQNS-1];\r
+       \r
+       System.out.println("QUADRATURE RULES STARTED:");\r
+       POPQF1(NPPQF,LQSBF,TNPQP,TOLOU,TPPQF,\r
+            TRRAD,WPPQF,ZPPQF,MQUPH,MQIN1,NARCS,NINTS,\r
+            NQPTS,TNSUA,DGPOL,JATYP,LOSUB,PARNT,\r
+            DELTA,LGTOL,ACOEF,BCOEF,H0VAL,HALEN,\r
+            JACIN,MIDPT,QUPTS,QUWTS,SOLUN,\r
+            RWORK,IER);\r
+             WRITE(*,10) 'QUADRATURE RULES DONE:'\r
        C\r
-                     CURR=WORK(2)\r
-                     PREV=SCALE\r
-                     DO 10 I=1,DG-1\r
-                         J=LOD+I-1\r
-                         NEXT=(AQCOF(J)*TT-BQCOF(J))*CURR-CQCOF(J)*PREV\r
-                         WORK(I+2)=NEXT\r
-                         PREV=CURR\r
-                         CURR=NEXT\r
-       10            CONTINUE\r
+             IF (IER .GT. 0) THEN\r
+               GOTO 999\r
+             ENDIF\r
        C\r
-       C             ASSIGN CORRECT SCALE FACTORS.\r
+             IQUPH(1)=TNPQP\r
        C\r
-                     DO 20 I=3,DG+1\r
-                         J=LOD+I-1\r
-                         WORK(I)=WORK(I)*COLSC(J)/SCALE\r
-       20            CONTINUE\r
-               } // if (DG > 0)\r
+       C**** SET UP THE CONSTANT FACTOR FOR THE MAPPING FORMULA\r
        C\r
-                 SG=1E+0\r
-                 DO 30 J=LOCLM,HICLM\r
-                   MATRX(ROW,J)=MATRX(ROW,J)+SG*WORK(J-LOCLM+1)\r
-                   SG=SG*SJI\r
-       30        CONTINUE\r
+             THET0=RGEOM(VTARG)\r
+             IF (INTER) THEN\r
+               ZQUPH(FACTR)=CMPLX(COS(THET0),SIN(THET0))\r
+             ELSE\r
+               ZQUPH(FACTR)=(1E+0,0E+0)\r
+               CALL DMPHYC(1,CENTR,CT,INTER,CENTR,IGEOM,RGEOM,ISNPH,RSNPH,\r
+            +              IQUPH,RQUPH,ZQUPH,.TRUE.,IER)\r
+               IF (IER .GT. 0) THEN\r
+                 GOTO 999\r
+               ENDIF\r
+               CT=CMPLX(RQUPH(1),THET0-AIMAG(LOG(CT)))\r
+               ZQUPH(FACTR)=CEXP(CT)\r
+             ENDIF\r
        C\r
-           } // for (IT=FIRST; IT <= LAST; IT++)\r
+             CALL OUPTPQ(IQUPH,RQUPH,ZQUPH,CHNL)\r
+       999   CONTINUE\r
        C\r
-       C       INITIALISE SOME DATA FOR THE NON-SINGULAR INTEGRALS\r
+       C**** WRITE CLOSING MESSAGE TO STANDARD OUTPUT CHANNEL\r
        C\r
-               WORK(1)=1E+0/RH\r
-               NQ=NQUAD(AJI)\r
-               LOQ=LOQSB(AJI)\r
-               IF (IA.EQ.TNSUA) THEN\r
-                   I1=2\r
-               ELSE\r
-                   I1=1\r
+             CALL WRTAIL(2,0,IER)\r
+       C*/\r
+    } // private void GQPHYC\r
+\r
+    private void POPQF1(int NPPQF[], int LPQSB[], int TNPQP[], double TOLOU, double TPPQF[],\r
+        double TRRAD[], double WPPQF[], double ZPPQF[][], int MNQUA, int MQIN1, int NARCS,\r
+        int NINTS, int NQPTS,int TNSUA, int DGPOL[], int JATYP[], int LOSUB[], int PARNT[],\r
+       double DELTA, double LGTOL, double ACOEF[], double BCOEF[], double H0VAL[], double HALEN[],\r
+       double JACIN[], double MIDPT[], double QUPTS[], double QUWTS[], double SOLUN[],\r
+        double XENPT[], int IER[]) {\r
+       // INTEGER IER,MQIN1,MNQUA,NARCS,NINTS,NQPTS,TNPQP,TNSUA\r
+       // INTEGER DGPOL(*),JATYP(*),LOSUB(*),LPQSB(*),NPPQF(*),\r
+       // +PARNT(*)\r
+       // REAL DELTA,LGTOL,TOLOU\r
+       // REAL ACOEF(*),BCOEF(*),HALEN(*),H0VAL(*),JACIN(*),MIDPT(*),\r
+       // +QUPTS(*),QUWTS(*),SOLUN(*),TPPQF(*),TRRAD(*),WPPQF(*),XENPT(*)\r
+       // COMPLEX ZPPQF(*)\r
+       \r
+       //     THE MAIN PURPOSE OF THIS ROUTINE IS TO SET UP THE REAL ARRAY\r
+       //     WPPQF OF QUADRATURE WEIGHTS, THE COMPLEX ARRAY ZPPQF OF\r
+       //     QUADRATURE POINTS ON THE PHYSICAL BOUNDARY AND THE REAL ARRAY\r
+       //     TPPQF OF LOCAL PARAMETER VALUES CORRESPONDING TO ZPPQF; ALL THESE\r
+       //     DATA DEFINE THE POST-PROCESSING COMPOSITE GAUSSIAN\r
+       //     QUADRATURE RULES FOR THE ESTIMATE OF THE MAP F (PHYSICAL ONTO\r
+       //     CANONICAL) WHEN Z IS NOT TOO CLOSE TO THE BOUNDARY.\r
+       \r
+       //     THIS ROUTINE ALSO RETURNS THE ARRAY TRRAD OF SO-CALLED TRANSITION\r
+       //     RADII;  THE IDEA IS THAT IF A PHYSICAL FIELD POINT Z SATISFIES\r
+       \r
+       //                 ABS(Z-ZPPQF(I)) < TRRAD(I)\r
+       \r
+       //     THEN THE QUADRATURE RULE PRODUCED BY THIS ROUTINE IS PROBABLY\r
+       //     NOT SUFFICIENTLY ACCURATE ON THE PARTICULAR ARC ON WHICH\r
+       //     ZPPQF(I) LIES.  THIS ROUTINE ALSO ASSIGNS FICTICIOUS QUADRATURE \r
+       //     POINTS WITH WEIGHTS SET TO ZERO.  THESE POINTS ARE INSERTED TO \r
+       //     ENSURE THAT IF A FIELD POINT Z IS CLOSE TO THE PHYSICAL BOUNDARY\r
+       //     THEN THE ABOVE INEQUALITY SHOULD BE SATISFIED FOR AT LEAST ONE\r
+       //     "GENERALISED" QUADRATURE POINT.\r
+       \r
+       //     THE ARRAY ELEMENT NPPQF(I) RECORDS THE NUMBER OF QUADRATURE\r
+       //     POINTS FOR THE SUBARC NUMBER I, I=1,...,TNSUA.  THE WEIGHTS,\r
+       //     QUADRATURE POINTS AND LOCAL QUADRATURE POINT PARAMETERS FOR ARC \r
+       //     NUMBER I ARE STORED IN WPPQF, ZPPQF AND TPPQF STARTING AT THE \r
+       //     ELEMENT WITH INDEX LPQSB(I).\r
+       \r
+       //     IER=0 - NORMAL EXIT \r
+       //     IER=22- THE TOTAL NUMBER OF QUADRATURE POINTS REQUESTED FOR \r
+       //             THE WHOLE BOUNDARY EXCEEDS THE LIMIT DEFINED BY THE\r
+       //             INPUT PARAMETER MNQUA; MNQUA MUST BE INCREASED IN A\r
+       //             HIGHER LEVEL ROUTINE.\r
+       //     IER=23- THE LOCAL PARAMETER MNCOF SHOULD BE INCREASED\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MNCOF = 32;\r
+       int AJT,DEG,HI,HI1,I,I1,J,J1,J2,JT,K,LIM,LOD,LOL,LOM,PT,QINTS;\r
+       final double RHO = 0.13;\r
+       double BETA,DIST,HL,JACSUM,MD,MEAN,MXDIS,RR,RRB,SS,SUM1,TT,SCO;\r
+       double CTT[] = new double[2];\r
+       // COMPLEX CTT,DPARFN,PARFUN\r
+       double JACOF[] = new double[MNCOF];\r
+       double JCOFC[][] = new double[MNCOF][2];\r
+       // COMPLEX JCOFC(MNCOF)\r
+       // EXTERNAL DPARFN,JACSUM,PARFUN,PPSBI7\r
+       \r
+       MXDIS=DELTA;\r
+       HI=0;\r
+       LOL=NARCS*NQPTS;\r
+       /*for (I1=1; I1 <= TNSUA; I1++) {\r
+           JT=JATYP[I1-1];\r
+           if (JT > 0) {\r
+               SS=1.0;\r
+           }\r
+           else {\r
+               SS=-1.0;\r
+           }\r
+           AJT=Math.abs(JT);\r
+           BETA=JACIN[AJT-1];\r
+           DEG=DGPOL[I1-1];\r
+           if (DEG+1 > MNCOF) {\r
+               IER[0]=23;\r
+               return;\r
+           }\r
+           LOM=LOSUB[I1-1];\r
+           LOD=(AJT-1)*NQPTS+1;\r
+           PT=PARNT[I1-1];\r
+           HL=HALEN[I1-1];\r
+           MD=MIDPT[I1-1];\r
+       \r
+           SCO=SS;\r
+           for (J=1; J <= DEG+1; J++) {\r
+               J1=LOM+J-1;\r
+               SCO=SCO*SS;\r
+               JACOF[J-1]=SOLUN[J1-1]*SCO;\r
+               JCOFC[J-1][0]=SOLUN[J1-1]*SCO;\r
+               JCOFC[J-1][1] = 0.0;\r
+           } // for (J=1; J <= DEG+1; J++)\r
+       \r
+           PPSBI7(DELTA,NINTS,BETA,NQPTS,DEG,ACOEF(LOD),\r
+            BCOEF(LOD),H0VAL(AJT),JCOFC,LGTOL,TOLOU,XENPT,\r
+            QINTS,MQIN1,IER);\r
+               IF (IER .GT. 0) THEN\r
+                 RETURN\r
                ENDIF\r
-       C\r
-               DO 100 IQ=1,NQ\r
-                 I=LOQ+IQ-1\r
-                 TQ=QCOMX(I)\r
-                 WQ=QCOMW(I)\r
-                 CALL JAPAR7(WORK,TQ,ACOEF(LOD),BCOEF(LOD),DG)\r
-                 IF (JI .LT. 0) THEN\r
-                   TQ=-TQ\r
+               NPPQF(I1)=QINTS*NQPTS\r
+               LPQSB(I1)=HI+1\r
+               HI1=HI+NPPQF(I1)\r
+               IF (HI1 .GT. MNQUA) THEN\r
+                 IER=22\r
+                 RETURN\r
+               ENDIF\r
+               K=HI\r
+               SUM1=BETA+1E+0\r
+               DO 30 I=1,QINTS\r
+                 RR=(XENPT(I+1)-XENPT(I))*5E-1\r
+                 MEAN=(XENPT(I+1)+XENPT(I))*5E-1\r
+                 IF (I .EQ. 1) THEN\r
+                   RRB=RR**SUM1\r
+                   DO 10 J=1,NQPTS\r
+                     J1=LOD+J-1\r
+                     K=K+1\r
+                     TT=(MEAN+RR*QUPTS(J1))\r
+                     WPPQF(K)=RRB*QUWTS(J1)*JACSUM(TT,DEG,ACOEF(LOD),\r
+            +                 BCOEF(LOD),H0VAL(AJT),JACOF)\r
+                     TT=TT*SS\r
+                     TPPQF(K)=TT\r
+                     CTT=CMPLX(MD+TT*HL)\r
+                     ZPPQF(K)=PARFUN(PT,CTT)\r
+                     TRRAD(K)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
+       10          CONTINUE\r
+                 ELSE\r
+                   DO 20 J=1,NQPTS\r
+                     J1=LOL+J\r
+                     K=K+1\r
+                     TT=(MEAN+RR*QUPTS(J1))\r
+                     WPPQF(K)=RR*QUWTS(J1)*(1E+0+TT)**BETA*JACSUM(TT,DEG,\r
+            +                 ACOEF(LOD),BCOEF(LOD),H0VAL(AJT),JACOF)\r
+                     TT=TT*SS\r
+                     TPPQF(K)=TT\r
+                     CTT=CMPLX(MD+TT*HL)\r
+                     ZPPQF(K)=PARFUN(PT,CTT)\r
+                     TRRAD(K)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
+       20          CONTINUE\r
                  ENDIF\r
-                 GTQ=CMPLX(MD+HL*TQ)\r
-                 ZQ=PARFUN(PT,GTQ)\r
-       C\r
-       C         ACCUMULATE THE ELEMENTS ABOVE THE DIAGONAL BLOCK\r
+       30      CONTINUE\r
+               IF (SS .LT. 0E+0) THEN\r
+                 LIM=NPPQF(I1)\r
+                 IF (MOD(LIM,2) .EQ. 0) THEN\r
+                   LIM=LIM/2\r
+                 ELSE\r
+                   LIM=(LIM-1)/2\r
+                 ENDIF\r
+                 J1=LPQSB(I1)-1\r
+                 J2=HI1+1\r
+                 DO 40 J=1,LIM\r
+                   J1=J1+1\r
+                   J2=J2-1\r
+                   TT=WPPQF(J1)\r
+                   WPPQF(J1)=WPPQF(J2)\r
+                   WPPQF(J2)=TT\r
+                   TT=TRRAD(J1)\r
+                   TRRAD(J1)=TRRAD(J2)\r
+                   TRRAD(J2)=TT\r
+                   TT=TPPQF(J1)\r
+                   TPPQF(J1)=TPPQF(J2)\r
+                   TPPQF(J2)=TT\r
+                   CTT=ZPPQF(J1)\r
+                   ZPPQF(J1)=ZPPQF(J2)\r
+                   ZPPQF(J2)=CTT\r
+       40        CONTINUE\r
+               ENDIF\r
        C\r
+       C       NEXT WE INSERT ANY NECESSARY FICTICIOUS QUADRATURE POINTS\r
        C\r
-                 DO 50 IT=I1,FIRST-1\r
-                   ZT=ZTEST(IT)\r
-                   DD=ABS(ZT-ZQ)\r
-                   LDD=LOG(DD)*WQ\r
-                   SG=1E+0\r
-                   DO 40 J1=LOCLM,HICLM\r
-                     J=J1-LOCLM+1\r
-                     I=J1-LOCLM+LOD\r
-                     MATRX(IT,J1)=MATRX(IT,J1)+SG*WORK(J)*LDD*COLSC(I)\r
-                     SG=SG*SJI\r
-       40          CONTINUE\r
+               J1=LPQSB(I1)\r
+               IF (TPPQF(J1)+1E+0 .GT. RHO*MXDIS) THEN\r
+                 J2=HI1\r
+                 DO 50 I=J2,J1,-1\r
+                   WPPQF(I+1)=WPPQF(I)\r
+                   TRRAD(I+1)=TRRAD(I)\r
+                   TPPQF(I+1)=TPPQF(I)\r
+                   ZPPQF(I+1)=ZPPQF(I)\r
        50        CONTINUE\r
+                 HI1=J2+1\r
+                 NPPQF(I1)=NPPQF(I1)+1\r
+                 WPPQF(J1)=0E+0\r
+                 TPPQF(J1)=-1E+0\r
+                 CTT=CMPLX(MD-HL)\r
+                 ZPPQF(J1)=PARFUN(PT,CTT)\r
+                 TRRAD(J1)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
+               ENDIF\r
        C\r
-       C         ACCUMULATE THE ELEMENTS BELOW THE DIAGONAL BLOCK\r
-       C\r
-                 DO 70 IT=LAST+1,NTEST\r
-                   ZT=ZTEST(IT)\r
-                   DD=ABS(ZT-ZQ)\r
-                   LDD=LOG(DD)*WQ\r
-                   SG=1E+0\r
-                   DO 60 J1=LOCLM,HICLM\r
-                     J=J1-LOCLM+1\r
-                     I=J1-LOCLM+LOD\r
-                     MATRX(IT,J1)=MATRX(IT,J1)+SG*WORK(J)*LDD*COLSC(I)\r
-                     SG=SG*SJI\r
-       60          CONTINUE\r
-       70        CONTINUE\r
-       C\r
-       C         ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
-       C         THE DIAGONAL BLOCK FOR THE NON-LINE-SEGMENT CASE.\r
-       C\r
-                 IF (.NOT.LNSEG(IA)) THEN\r
-                   DO 90 IT=FIRST,LAST\r
-                     IF (ORDSG.GT.1 .AND. IA.EQ.(TSFBS+1) .AND. IT.EQ.LAST) \r
-            +        GOTO 90\r
-                     IF (IT .EQ. LAST) THEN\r
-                         TT=1E+0\r
-                     ELSE\r
-                         TT=TESPR(IT)\r
-                     ENDIF\r
-       C\r
-                     IF (IT .GT. NTEST) THEN\r
-                         ROW=1\r
-                     ELSE\r
-                         ROW=IT\r
-                     ENDIF\r
-       C\r
-                     DD=ABS(TT-TQ)\r
-                     IF (DD .LE. TOLNR) THEN\r
-                         DD=ABS(DPARFN(PT,GTQ))*HL\r
-                     ELSE\r
-                         IF (IT .GT. NTEST) THEN\r
-                             ZT=ZTEST(1)\r
-                         ELSE\r
-                             ZT=ZTEST(IT)\r
-                         ENDIF\r
-                         DD=ABS(ZT-ZQ)/DD\r
-                         IF (DD.LT.TOLNR) DD=ABS(DPARFN(PT,GTQ))*HL\r
-                     ENDIF\r
-                     LDD=LOG(DD)*WQ\r
-                     SG=1E+0\r
-                     DO 80 J1=LOCLM,HICLM\r
-                       J=J1-LOCLM+1\r
-                       I=J1-LOCLM+LOD\r
-                       MATRX(ROW,J1)=MATRX(ROW,J1)+SG*WORK(J)*LDD*COLSC(I)\r
-                       SG=SG*SJI\r
-       80            CONTINUE\r
-       90          CONTINUE\r
+       60      CONTINUE\r
+               J1=J1+1\r
+               IF (J1 .LT. HI1) THEN\r
+                 DIST=TPPQF(J1+1)-TPPQF(J1)\r
+                 IF (DIST .GT. MXDIS) THEN\r
+                   J2=HI1\r
+                   DO 70 I=J2,J1+1,-1\r
+                     WPPQF(I+1)=WPPQF(I)\r
+                     TRRAD(I+1)=TRRAD(I)\r
+                     TPPQF(I+1)=TPPQF(I)\r
+                     ZPPQF(I+1)=ZPPQF(I)\r
+       70          CONTINUE\r
+                   HI1=J2+1\r
+                   NPPQF(I1)=NPPQF(I1)+1\r
+                   WPPQF(J1+1)=0E+0\r
+                   TPPQF(J1+1)=(TPPQF(J1)+TPPQF(J1+2))*5E-1\r
+                   CTT=CMPLX(MD+HL*TPPQF(J1+1))\r
+                   ZPPQF(J1+1)=PARFUN(PT,CTT)\r
+                   TRRAD(J1+1)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
                  ENDIF\r
-       100     CONTINUE\r
-       C\r
-       C       ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
-       C       THE DIAGONAL BLOCK FOR THE LINE-SEGMENT CASE.\r
-       C\r
-               IF (LNSEG(IA)) THEN\r
-                 ZT=DPARFN(PT,CMPLX(MD))*HL\r
-                 C0=ABS(ZT)\r
-                 C0=RH*LOG(C0)*COLSC(LOD)\r
-                 DO 110 IT=FIRST,LAST\r
-                   IF (ORDSG.GT.1 .AND. IA.EQ.(TSFBS+1) .AND. IT.EQ.LAST) \r
-            +      GOTO 110\r
-                   IF (IT .GT. NTEST) THEN\r
-                       ROW=1\r
-                   ELSE\r
-                       ROW=IT\r
-                   ENDIF\r
-                   MATRX(ROW,LOCLM)=MATRX(ROW,LOCLM)+C0\r
-       110       CONTINUE\r
+                 GOTO 60\r
+               ELSE IF (1E+0-TPPQF(J1) .GT. RHO*MXDIS) THEN\r
+                 J1=J1+1\r
+                 HI1=J1\r
+                 NPPQF(I1)=NPPQF(I1)+1\r
+                 WPPQF(J1)=0E+0\r
+                 TPPQF(J1)=1E+0\r
+                 CTT=CMPLX(MD+HL)\r
+                 ZPPQF(J1)=PARFUN(PT,CTT)\r
+                 TRRAD(J1)=HL*DELTA*ABS(DPARFN(PT,CTT))\r
                ENDIF\r
-       } // for (IA=1; IA <= TNSUA; IA++)\r
-       C\r
-       C     SET UP THE LAST COLUMN\r
-       C\r
-             DO 130 I=1,NTEST\r
-               MATRX(I,NEQNS)=1E+0\r
-       130   CONTINUE\r
-       C\r
-       C     COMPUTE MATRIX-VECTOR PRODUCT\r
+               HI=HI1\r
+       } // for (I1=1; I1 <= TNSUA; I1++)*/\r
+       \r
+       TNPQP[0]=HI;\r
+       \r
+       IER[0]=0;\r
+       \r
+    } // private void POPQF1\r
+    \r
+    private void PPSBI7(double DELTA, int NINTS, double BETA, int NQUAD, int DGPOL, double ACOEF[],\r
+        double BCOEF[], double H0VAL, double SOLUN[][], double TOLIN, double TOLOU, double XENPT[],\r
+        int QINTS[], int MQIN1, int IER[]) {\r
+       // INTEGER DGPOL,MQIN1,NQUAD,QINTS,IER,NINTS\r
+       // REAL BETA,H0VAL,TOLIN,TOLOU,DELTA\r
+       // REAL ACOEF(*),BCOEF(*),XENPT(*)\r
+       // COMPLEX SOLUN(*)\r
+       \r
+       //     CALCULATES THE NUMBER OF QUADRATURE INTERVALS (QINTS) REQUIRED\r
+       //     FOR THE COMPOSITE GAUSS-JACOBI/GAUSS-LEGENDRE ESTIMATION OF\r
+       \r
+       //          INTEGRAL  [(1+X)**BETA*FNPHI(X)*LOG(ZZ(J)-X)*dX], \r
+       //         -1<=X<=1                                         \r
+       \r
+       //     WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY \r
+       //     CORRESPONDENCE DERIVATIVE / JACOBI WEIGHT QUOTIENT. \r
+        \r
+       //     ZZ IS ANY POINT ON A "DELTA-CONTOUR" IN THE UPPER HALF PLANE,\r
+       //     THIS CONTOUR BEING DEFINED BY THE PARAMETER DELTA.  TEST VALUES \r
+       //     FOR ZZ ARE ASSIGNED IN THE SUBROUTINES DEPPJ8 AND DEPPL8.\r
+       \r
+       //     THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO \r
+       //     DEFINE FNPHI AND ARE PASSED TO DEPPJ8 AND DEPPL8 FOR THIS PURPOSE.\r
+       \r
+       //     THE ENDPOINTS OF THE QUADRATURE INTERVALS ARE RETURNED IN VECTOR \r
+       //     XENPT, WITH XENPT(1)=-1<XENPT(2)<...<1=XENPT(QINTS+1).\r
+       \r
+       //     TOLOU RECORDS OUR ESTIMATE FOR THE MAXIMUM OF THE ABSOLUTE VALUES\r
+       //     OF THE REMAINDER ESTIMATES ON THE DELTA-CONTOUR.  WE REQUIRE THAT\r
+       \r
+       //                  TOLOU <= TOLIN\r
+       \r
+       //     WITH THE TOLOU BEING REASONABLY CLOSE TO TOLIN.\r
+       \r
+       //     IER=0 - NORMAL EXIT\r
+       //     IER=24- THE REQUIRED NUMBER OF QUADRATURE INTERVALS EXCEEDS THAT\r
+       //             SPECIFIED BY THE GLOBAL PARAMETER MQIN1; MQIN1 MUST BE\r
+       //             INCREASED (ERROR NUMBER 43 IF CALLED FROM POPQG1)\r
+          \r
+       //     LOCAL VARIABLES\r
+       \r
+       int INTS;\r
+       double TAU[] = new double[1];\r
+       double MAXRM[] = new double[1];\r
+       double TOL,RIGHT;\r
+       boolean T1FXD;\r
+       // EXTERNAL DEPPJ8,DEPPL8\r
+       \r
+       QINTS[0]=1;\r
+       XENPT[0]=-1.0;\r
+       TOL=TOLIN;\r
+       INTS=NINTS;\r
+       DEPPJ8(BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL,\r
+            SOLUN,TOL,MAXRM,INTS,DELTA,IER);\r
+       /*      IF (IER .GT. 0) THEN\r
+               RETURN\r
+             ENDIF\r
+             TOLOU=MAXRM\r
+             XENPT(2)=TAU\r
        C\r
-             DO 150 I=1,NTEST\r
-               SUM=0E+0\r
-               DO 140 J=1,NEQNS\r
-                 SUM=SUM+MATRX(I,J)*SOLUN(J)\r
-       140     CONTINUE\r
-               ERMOD(I)=SUM\r
-       150   CONTINUE\r
+             IF (XENPT(2) .LT. 1E+0) THEN\r
+               QINTS=2\r
+               T1FXD=.FALSE.\r
+               TAU=1E+0\r
+               RIGHT=-1E+0\r
+               INTS=NINTS\r
+               CALL DEPPL8(BETA,RIGHT,TAU,T1FXD,NQUAD,DGPOL,ACOEF,\r
+            +              BCOEF,H0VAL,SOLUN,TOL,MAXRM,INTS,DELTA,IER)\r
+               IF (IER .GT. 0) THEN\r
+                 RETURN\r
+               ENDIF\r
+               TOLOU=TOLOU+MAXRM\r
+               T1FXD=.TRUE.\r
        C\r
-       C     FORM THE ERROR IN MODULUS\r
+       100     CONTINUE\r
        C\r
-             IF (INTER) THEN\r
-                 DO 160 I=1,NTEST\r
-                   SUM=EXP(ERMOD(I))\r
-                   ERMOD(I)=ABS(1E+0-ABS(ZTEST(I)-CENTR)/SUM)\r
-       160       CONTINUE\r
-             ELSE\r
-                 DO 170 I=1,NTEST\r
-                   SUM=EXP(ERMOD(I))\r
-                   ERMOD(I)=ABS(1E+0-SUM)\r
-       170       CONTINUE\r
+               IF (XENPT(QINTS) .GT. RIGHT) THEN\r
+                 XENPT(QINTS)=5E-1*(XENPT(QINTS)+RIGHT)\r
+                 XENPT(QINTS+1)=1E+0\r
+               ELSE\r
+                 TAU=1E+0\r
+                 INTS=NINTS\r
+                 CALL DEPPL8(BETA,XENPT(QINTS),TAU,T1FXD,NQUAD,DGPOL,\r
+            +                ACOEF,BCOEF,H0VAL,SOLUN,TOL,MAXRM,INTS,DELTA,IER) \r
+                 IF (IER .GT. 0) THEN\r
+                   RETURN\r
+                 ENDIF\r
+                 TOLOU=TOLOU+MAXRM\r
+                 QINTS=QINTS+1\r
+                 IF (QINTS .GE. MQIN1) THEN\r
+                   IER=24\r
+                   RETURN\r
+                 ENDIF\r
+                 XENPT(QINTS)=TAU\r
+                 GOTO 100\r
+               ENDIF\r
              ENDIF\r
        C\r
-       C     FIND MAXIMUM ERROR IN MODULUS AND THE POINT AND THE ARC AT WHICH \r
-       C     IT OCCURS\r
+             IER=0\r
+       C */\r
+    } // private void PPSBI7\r
+    \r
+    private void DEPPJ8(double BETA, double TAU[], int NQUAD, int DGPOL, double ACOEF[],\r
+        double BCOEF[], double H0VAL, double SOLUN[][], double TOL, double MAXRM[], int NINTS,\r
+        double DELTA,int IER[]) {\r
+        //INTEGER NQUAD,IER,DGPOL,NINTS\r
+       //REAL BETA,TAU,TOL,MAXRM,ACOEF(*),BCOEF(*),H0VAL,DELTA\r
+       //COMPLEX SOLUN(*)\r
+       \r
+       // WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN\r
+       //     USING AN NQUAD-POINT GAUSS-JACOBI RULE TO ESTIMATE THE INTEGRALS\r
+       \r
+       //            INTEGRAL  [(1+X)**BETA*FNPHI(X)*LOG(ZZ-X)*dX]\r
+       //           -1<=X<=TAU                                       \r
+       \r
+       //     WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY \r
+       //     CORRESPONDENCE DERIVATIVE / JACOBI WEIGHT QUOTIENT. \r
+        \r
+       //     ZZ IS ANY POINT ON A "DELTA-CONTOUR" IN THE UPPER HALF PLANE,\r
+       //     THIS CONTOUR BEING DEFINED BY THE PARAMETE DELTA.  TEST VALUES \r
+       //     FOR ZZ ARE ASSIGNED IN THE BODY OF THE ROUTINE AND THE LOCAL ARRAY\r
+       //     XIVAL IS USED FOR STORING THESE TEST VALUES.\r
+       \r
+       //     THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO \r
+       //     DEFINE FNPHI AND ARE PASSED TO DEPPJ9 FOR THIS PURPOSE.\r
+       \r
+       //     MAXRM RECORDS THE MAXIMUM OF THE ABSOLUTE VALUES OF THE REMAINDER\r
+       //     ESTIMATES ASSOCIATED WITH THE DELTA-CONTOUR.\r
+       \r
+       //     ON ENTRY NINTS IS THE NUMBER OF INITIAL TEST INTERVALS TO BE USED\r
+       //     ON THE UPPER HALF DELTA-CONTOUR; ON EXIT NINTS GIVES THE FINAL\r
+       //     NUMBEROF INTERVALS REQUIRED FOR CONVERGENCE TO TAU.\r
+       \r
+       //     THE PURPOSE OF THIS ROUTINE IS TO DETERMINE A VALUE FOR TAU\r
+       //     SUCH THAT\r
+       \r
+       //                         MAXRM < TOL \r
+       \r
+       //     AND THAT, IF POSSIBLE, \r
+       \r
+       //                   0.5*TOL <= MAXRM < TOL.\r
+       \r
+       //     IER=0 - NORMAL EXIT\r
+       //     IER=25- THE LOCAL ARRAY BOUND PARAMETER MNXI NEEDS INCREASING.\r
+                    \r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MNXI = 100;\r
+       int NXI;\r
+       double SS,PI,SS1,SS2,SS3,DP,LEN,TAUI,SINC,STAR,RXI,DT,HDP;\r
+       double XI[] = new double[2];\r
+       //COMPLEX XI\r
+       double XIVAL[][] = new double[MNXI][2];\r
+       // COMPLEX XIVAL(MNXI)\r
+       // EXTERNAL DEPPJ9\r
+       \r
+       //     INITIALISATION\r
+       \r
+       PI=Math.PI;\r
+             /*TAUI=1E+0\r
+             DP=DELTA*PI\r
+             DT=DELTA*2E+0\r
+             HDP=5E-1*DP\r
+             LEN=2E+0+DP\r
+             SINC=LEN/NINTS\r
+             STAR=SINC\r
        C\r
-             MXERM=0E+0\r
-             DO 190 IA=1,TSFBS\r
-               FIRST=LOTES(IA)\r
-               LAST=HITES(IA)\r
-               MD=0E+0\r
-               DO 180 IT=FIRST,LAST\r
-                 IF (ERMOD(IT) .GT. MD) THEN\r
-                     MD=ERMOD(IT)\r
-                 ENDIF\r
-                 IF (MD .GT. MXERM) THEN\r
-                     MXERM=MD\r
-                     IMXER=IA\r
-                     ZMXER=ZTEST(IT)\r
-                 ENDIF\r
-       180     CONTINUE\r
-               IF (IA.EQ.TSFBS .AND. ORDSG.EQ.1) THEN\r
-                 IT=1\r
+       C     START OF LOOP FOR INTERVAL HALVING\r
+       C\r
+       10    CONTINUE\r
+             NXI=0\r
+             SS1=-STAR+SINC\r
+             SS2=LEN\r
+             DO 30 SS=SS1,SS2,SINC\r
+               IF (SS .LT. HDP) THEN\r
+                   SS3=SS/DELTA\r
+                   XI=1E+0+DELTA*CMPLX(COS(SS3),SIN(SS3))\r
+               ELSE IF (SS .GT. 2E+0+HDP) THEN\r
+                   SS3=(SS-2E+0-HDP)/DELTA\r
+                   XI=-1E+0+DELTA*(0E+0,1E+0)*CMPLX(COS(SS3),SIN(SS3))\r
                ELSE\r
-                 IT=LAST+1\r
-               ENDIF\r
-               IF (ERMOD(IT) .GT. MD) THEN\r
-                   MD=ERMOD(IT)\r
+                   SS3=SS-HDP\r
+                   XI=CMPLX(1E+0-SS3,DELTA)\r
                ENDIF\r
-               IF (MD .GT. MXERM) THEN\r
-                   MXERM=MD\r
-                   IMXER=IA\r
-                   ZMXER=ZTEST(IT)\r
+               RXI=REAL(XI)\r
+               IF (RXI .LT. (TAUI+DT)) THEN\r
+                 NXI=NXI+1\r
+                 IF (NXI .GT. MNXI) THEN\r
+                   IER=25\r
+                   RETURN\r
+                 ENDIF\r
+                 XIVAL(NXI)=XI\r
                ENDIF\r
-               ERARC(IA)=MD\r
-       190   CONTINUE\r
-       C\r
-       C     IF REGION IS SYMMETRIC, FILL UP THE WHOLE ERARC VECTOR USING\r
-       C     SYMMETRY\r
+       30    CONTINUE\r
        C\r
-             IF (ORDSG.GT.1) THEN\r
-             IF (REFLN) THEN\r
-               DO 200 IA=1,TSFBS\r
-                 IB=2*TSFBS+1-IA\r
-                 ERARC(IB)=ERARC(IA)\r
-       200     CONTINUE\r
-               ORDRG=ORDSG/2\r
-               DO 220 I=2,ORDRG\r
-                 I1=(I-1)*TSFBS*2\r
-                 DO 210 IA=1,TSFBS*2\r
-                   IB=I1+IA\r
-                   ERARC(IB)=ERARC(IA)\r
-       210       CONTINUE\r
-       220     CONTINUE\r
-             ELSE\r
-               DO 240 I=2,ORDSG\r
-                 I1=(I-1)*TSFBS\r
-                 DO 230 IA=1,TSFBS\r
-                   IB=I1+IA\r
-                   ERARC(IB)=ERARC(IA)\r
-       230       CONTINUE\r
-       240     CONTINUE\r
+             IF (NXI .EQ. 0) THEN\r
+               SINC=STAR\r
+               STAR=5E-1*STAR\r
+               NINTS=NINTS*2\r
+               GOTO 10\r
              ENDIF\r
+       C      \r
+             TAU=TAUI\r
+             CALL DEPPJ9(XIVAL,NXI,BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL,\r
+            +            SOLUN(1),TOL,MAXRM,IER) \r
+             IF (IER .GT. 0) THEN\r
+               RETURN\r
              ENDIF\r
-       C*/\r
-    } // private void TESMD9\r
+       C\r
+             IF (TAU .NE. TAUI) THEN\r
+               TAUI=TAU\r
+               SINC=STAR\r
+               STAR=5E-1*STAR\r
+               NINTS=NINTS*2\r
+               GOTO 10\r
+             ENDIF\r
+       C\r
+             IER=0\r
+       C\r
+             END*/\r
+    } // private void DEPPJ8\r
 \r
 \r
       /**\r