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

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

index 9dd318c..bc95cf6 100644 (file)
@@ -2672,18 +2672,21 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                final double SFACT = 0.8;\r
                final double QFACT = 0.1;\r
                double MCQER[] = new double[1];\r
-               double RCOOND[] = new double[1];\r
+               double MQERR[] = new double[1];\r
+               double RCOND[] = new double[1];\r
                double ESTOL[] = new double[1];\r
                double AQTOL,CONST,GAQTL,GLGTL,GRQTL,GSUPE,GTGTE,LGTOL,\r
-                    MCHEP,MQERR,MXERM,PI,RQTOL,SSUPE,\r
+                    MCHEP,MXERM,PI,RQTOL,SSUPE,\r
                     TGTER,TOLNR;\r
+               double SOLUNB[][] = new double[MNEQN][1];\r
                \r
                double ZMXER[] = new double[2];\r
                //COMPLEX ZMXER\r
                \r
                final boolean INIBT = true;\r
                boolean ACCPT[] = new boolean[1];\r
-               boolean GACPT,NUQTL,REFLN;\r
+               boolean NUQTL[] = new boolean[1];\r
+               boolean GACPT,REFLN;\r
                \r
                String OFL;\r
                //CHARACTER OFL*6\r
@@ -2821,7 +2824,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                RSLT80(JBNM,HEAD,GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
                     INCST,INDEG,RFARC,RFARG[0],CENTR,JACIN,LNSEG,\r
-                    TSTNG,OULVL,IBNDS,MNEQN,OCH);\r
+                    TSTNG,OULVL,IBNDS,MNEQN);\r
                PI=Math.PI;\r
                RFARG[0]=RFARG[0]*PI;\r
                \r
@@ -2904,7 +2907,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                     BCOEF,H0VAL,COLSC,NQUAD,LOQSB,\r
                     QCOMX,QCOMW,MNQUA,TOLOU,MCQER,XENPT,\r
                     XIVAL,XIDST,IER);\r
-               NUQTL=false;\r
+               NUQTL[0]=false;\r
                if (IER[0] > 0) {\r
                        WRTAIL(1,0,IER[0],null);\r
                    return;    \r
@@ -2913,9 +2916,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                //**** SET UP LINEAR ALGEBRAIC SYSTEM.\r
                \r
-/*         while (true) {\r
+           while (true) {\r
                    SOLCO=SOLCO+1;\r
-                   System.out.prinln("********SOLUTION "+ SOLCO+ " ******** " + NROWS + " EQUATIONS");\r
+                   System.out.println("********SOLUTION "+ SOLCO+ " ******** " + NROWS + " EQUATIONS");\r
                \r
                    LNSY11(MATRX,SOLUN,MNEQN,NCOLL,ORDSG,REFLN,NQPTS,\r
                     TNSUA[0],JATYP,PARNT,DGPOL,LOSUB,\r
@@ -2964,7 +2967,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                                \r
                                    BCFVTF(BCFSN,VTARG,DGPOL,JATYP,\r
                                        LOSUB,PARNT,RFARC,TNSUA[0],H0VAL,JACIN,\r
-                                       RFARG,SOLUN);\r
+                                       RFARG[0],SOLUN);\r
                                \r
                                    // ****   OUTPUT DATA REQUIRED FOR POST-PROCESSING.\r
                                \r
@@ -2979,7 +2982,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    // LINPACK SGESL is equivalent to LAPACK DGETRS\r
                    // Solves a general system of linear equations, after factorization by SGECO or SGEFA\r
                    // SOLUN has the solution vector on exit.\r
-                   le2.dgetrs('N', NROWS, 1, A, MNEQN, IPIVT, SOLUN, MNEQN, info);\r
+                   for (I = 0; I < MNEQN; I++) {\r
+                       SOLUNB[I][0] = SOLUN[I];\r
+                   }\r
+                   le2.dgetrs('N', NROWS, 1, A, MNEQN, IPIVT, SOLUNB, MNEQN, info);\r
+                   for (I = 0; I < MNEQN; I++) {\r
+                       SOLUN[I] = SOLUNB[I][0];\r
+                   }\r
                    //SGESL(MATRX(1,1,2),MNEQN,NROWS,IPIVT,SOLUN,0)\r
                    for (I = 0; I < MNEQN; I++) {\r
                        for (J = 0; J < MNEQN; J++) {\r
@@ -3016,7 +3025,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                                        \r
                                            BCFVTF(BCFSN,VTARG,DGPOL,JATYP,\r
                                                LOSUB,PARNT,RFARC,TNSUA[0],H0VAL,JACIN,\r
-                                               RFARG,SOLUN);\r
+                                               RFARG[0],SOLUN);\r
                                        \r
                                            // ****   OUTPUT DATA REQUIRED FOR POST-PROCESSING.\r
                                        \r
@@ -3030,7 +3039,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                     } // if (IER[0] > 1)\r
                     System.out.println("DECISIONS FOR EACH ARC DONE:");\r
                     System.out.println("EFFECTIVE STOPPING TOLERANCE: " + ESTOL[0]);\r
-                    if (ACCPT[0] && ESTOL <= GSUPE) {\r
+                    if (ACCPT[0] && ESTOL[0] <= GSUPE) {\r
                         GACPT=true;\r
                     }\r
                     else {\r
@@ -3042,133 +3051,221 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                          // OUTPUT RESULTS\r
                \r
                          if (OULVL < 4) {\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
+                             RSLT72(QIERC,RCOND[0],CONST,NROWS,DGPOL,JATYP,\r
+                                    PARNT,TNSUA[0],INTER,MQERR[0],MCQER[0],AXION,\r
+                                    NEWDG,NJIND,NQUAD,TOLOU,LGTOL,SOLCO);\r
                          }\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
+                             RSLT71(QIERC,RCOND[0],SOLUN,NEQNS,LOSUB,\r
+                                    HISUB,COLSC,NQPTS,JATYP,\r
+                                    PARNT,TNSUA[0],INTER,MQERR[0],MCQER[0],WORK2,\r
+                                    AXION,NEWDG,NJIND,JACIN,\r
+                                    NQUAD,TOLOU,LGTOL,SOLCO);\r
                          }\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
-                       break;\r
+                         NEFF=(int)Math.round(Math.pow((double)(NEFF),3.3333333E-1));\r
+                         System.out.println("\n****THE SOLUTION IS ACCEPTED****");\r
+                         System.out.println("\nEFFECTIVE SIZE OF ALL SYSTEMS: " + NEFF);\r
+                         if (INTER) {\r
+                             System.out.println("ZERO " + CONST);\r
+                         }\r
+                         else {\r
+                             System.out.println("CAPACITY: " + (Math.exp(-CONST)));\r
+                         }\r
+               \r
+                  Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+                  Preferences.debug("****THE SOLUTION IS ACCEPTED****\n",Preferences.DEBUG_ALGORITHM);\r
+                         Preferences.debug("EFFECTIVE SIZE OF ALL SYSTEMS: " + NEFF + "\n\n");\r
+                         break;\r
                      } // if (GACPT)\r
                      else { // !GACPT\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
-                       continue;\r
+                         if (ACCPT[0] || ESTOL[0] <= SSUPE) {\r
+               \r
+                             // SOLUTION AT INTERMEDIATE ACCURACY IS ACCEPTED; SET TOLERANCES\r
+                             // FOR GREATER ACCURACY AND RE-ASSESS UPDATING ACTIONS BEFORE \r
+                             // CONTINUING\r
+               \r
+                             SSUPE=0.1*Math.min(SSUPE,ESTOL[0]);\r
+                             TGTER=SFACT*SSUPE;\r
+                             if (TGTER <= 2.0*GTGTE) {\r
+                                 TGTER=GTGTE;\r
+                             }\r
+                             AQTOL=TGTER*QFACT;\r
+                             NUQTL[0]=true;\r
+                             LGTOL=Math.log(1.0+TGTER);\r
+                             RQTOL=AQTOL;\r
+                             I=(int)Math.round(-Math.log10(TGTER))+2;\r
+                             INDEG=Math.min(I,MDGPO);\r
+               \r
+                             // DETERMINE THE ACTIONS THAT HAVE TO BE TAKEN ON EACH ARC\r
+               \r
+                             AXION1(AXION,NEWDG,SOLUN,MDGPO,\r
+                                    TNSUA[0],DGPOL,LOSUB,HISUB,\r
+                                    RIGLL,LGTOL,ACCPT,JACIN,\r
+                                    JATYP,NJIND,NEWHL,ESTOL,IER);\r
+                             ESTOL[0]=ESTOL[0]/SFACT;\r
+                             if (IER[0] > 0) {\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
+                                    } // if (IER[0] > 1)\r
+                            System.out.println("DECISIONS FOR EACH ARC RE-DONE:");\r
+                         } //  if (ACCPT[0] || ESTOL[0] <= SSUPE)\r
+               \r
+                         // OUTPUT RESULTS\r
+               \r
+                         if (OULVL < 4) {\r
+                             RSLT72(QIERC,RCOND[0],CONST,NROWS,DGPOL,JATYP,\r
+                                    PARNT,TNSUA[0],INTER,MQERR[0],MCQER[0],AXION,\r
+                                    NEWDG,NJIND,NQUAD,TOLOU,\r
+                                    LGTOL,SOLCO);\r
+                         }\r
+                         else {\r
+                             RSLT71(QIERC,RCOND[0],SOLUN,NEQNS,LOSUB,\r
+                                    HISUB,COLSC,NQPTS,JATYP,\r
+                                    PARNT,TNSUA[0],INTER,MQERR[0],MCQER[0],WORK2,\r
+                                    AXION,NEWDG,NJIND,JACIN,\r
+                                    NQUAD,TOLOU,LGTOL,SOLCO);\r
+                         }\r
+                         System.out.println("EFFECTIVE STOPPING TOLERANCE : " + ESTOL[0]);\r
+                         if (RCOND[0] < 5.0*MCHEP) {\r
+                             IER[0]=16;\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
+                         } // if (RCOND[0] < 5.0*MCHEP)\r
+                         else if (RCOND[0] < AQTOL) {\r
+                             NUQTL[0]=true;\r
+                             AQTOL=0.1*RCOND[0];\r
+                             if (AQTOL < 5.0*MCHEP) AQTOL=5.0*MCHEP;\r
+                         }\r
+               \r
+                         // IMPLEMENT UPDATING PROCEDURES.\r
+                         // FIRST UPDATE THE COLLOCATION PARAMETERS AND OTHER DATA\r
+                         // RELATING TO SUB-ARC DEFINITIONS.\r
+               \r
+                         UPJAC1(NQPTS,NJIND,INDEG,AXION,DGPOL,\r
+                         NEWDG,ACOEF,BCOEF,DIAG,\r
+                         SDIAG,TNSUA[0],MNSUA,LOSUB,HISUB,\r
+                         JATYP,PARNT,MIDPT,HALEN,\r
+                         COLPR,ZCOLL,LNSEG,PNEWQ,EPS,IER,\r
+                         WORK,NEWHL,RCOPY,ICOPY,\r
+                         LCOPY,LOOLD,HIOLD);\r
+                         if (IER[0] > 0) {\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
+                         } // if (IER[0] > 0)\r
+                         System.out.println("ARC REFINEMENTS DONE:");\r
+                         NCOLL=HISUB[TNSUA[0]-1];\r
+                         NEQNS=NCOLL+1;\r
+                         NROWS=NCOLL/ORDSG+1;\r
+                         if (NEQNS > MNEQN) {\r
+                             IER[0]=18;\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
+                         } // if (NEQNS > MNEQN)\r
+               \r
+                         // NEXT UPDATE THE COMPOSITE QUADRATURE RULES\r
+               \r
+                         UPCOQ1(NARCS,NJIND,NQPTS,MDGPO,MQIN1,AQTOL,QUPTS,\r
+                         QUWTS,JACIN,MIDPT,HALEN,\r
+                         ACOEF,BCOEF,H0VAL,COLSC,\r
+                         NQUAD,LOQSB,QCOMX,QCOMW,\r
+                         MNQUA,TOLOU,MCQER,XENPT,XIVAL,\r
+                         XIDST,TNSUA[0],PNEWQ,NEWQU,JATYP,\r
+                         PARNT,NUQTL,IER);\r
+                         if (IER[0] > 0) {\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
+                         } // if (IER[0] > 0)\r
+                         System.out.println("QUADRATURE UPDATES DONE:");\r
+                         continue;\r
                      } // else !GACPT \r
            } // while (true)\r
-               C\r
-                     IF (OULVL.EQ.2 .OR. OULVL.EQ.3) THEN\r
+               \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
@@ -3354,7 +3451,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     \r
     private void RSLT80(String JBNM, String HEAD, 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,int OCH) {\r
+        int TSTNG, int OULVL, int IBNDS[], int MNEQN) {\r
        \r
        //INTEGER IBNDS(*)\r
        //REAL BETA(*)\r
@@ -5504,103 +5601,96 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     } // private void DELEG7\r
 \r
 \r
-            /* SUBROUTINE RSLT71(QIERC,RCOND,SOLUN,NEQNS,LOSUB,HISUB,COLSC,\r
-            +NQPTS,JATYP,PARNT,TNSUA,INTER,MQERR,MCQER,CINFN,ACTIN,\r
-            +NEWDG,NJIND,JACIN,NQUAD,TOLOU,LGTOL,SOLCO,OUCH1)\r
-             INTEGER NEQNS,TNSUA,OUCH1,NQPTS,NJIND,NEWDG(*),NQUAD(*),LOSUB(*),\r
-            +HISUB(*),QIERC(0:6),JATYP(*),PARNT(*),ACTIN(*),SOLCO\r
-             REAL SOLUN(*),RCOND,COLSC(*),MQERR,MCQER,LGTOL,\r
-            +CINFN(*),JACIN(*),TOLOU(*)\r
-             LOGICAL INTER\r
-             CHARACTER QTEXT(0:6)*22,LINE*72\r
-             PARAMETER (LINE='_________________________________________________\r
-            +________________')\r
-       C \r
-       C     LOCAL VARIABLES\r
-       C\r
-             INTEGER I,J,JI,K,L,LOD,N,H\r
-             REAL S,CAP\r
-       C\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
-       C\r
-             WRITE(OUCH1,*) LINE\r
-             WRITE(OUCH1,*) '             SOLUTION NUMBER =',SOLCO\r
-             WRITE(OUCH1,*) '                       NEQNS =',NEQNS \r
-             WRITE(OUCH1,*) 'RECIPROCAL COND NO. ESTIMATE =',RCOND\r
-             WRITE(OUCH1,*) '   CONDITION NO. LOWER BOUND =',1E+0/RCOND\r
-       C\r
-             WRITE(OUCH1,*) \r
-             WRITE(OUCH1,997) 'JACOBI INDEX','POINTS','TOLERANCE ACHIEVED'\r
-             DO 10 I=1,NJIND\r
-               WRITE(OUCH1,998) I,NQUAD(I),TOLOU(I)\r
-       10    CONTINUE\r
-       C\r
-             WRITE(OUCH1,*) \r
-             WRITE(OUCH1,*) 'QAWS TERMINATIONS WITH......' \r
-             DO 20 I=0,6\r
-               IF (QIERC(I) .GT. 0) THEN\r
-                 WRITE(OUCH1,1000) QTEXT(I),QIERC(I)\r
-               ENDIF\r
-       20    CONTINUE\r
-       C\r
-             WRITE(OUCH1,*) \r
-             WRITE(OUCH1,999) '              MAXIMUM QAWS ERROR =',MQERR\r
-             WRITE(OUCH1,999) 'MAXIMUM COMPOSITE GAUSSIAN ERROR =',MCQER\r
-             WRITE(OUCH1,*) \r
-             DO 40 I=1,TNSUA\r
-                 WRITE(OUCH1,*)\r
-                 WRITE(OUCH1,*) 'SUB ARC =',I,' ON PARENT ARC',PARNT(I)\r
-                 WRITE(OUCH1,990) 'N','SCALED SOLUN','UNSCALED SOLUN','IGNORE L\r
-            +EVEL'\r
-                 L=LOSUB(I)\r
-                 H=HISUB(I)\r
-                 JI=ABS(JATYP(I))\r
-                 LOD=(JI-1)*NQPTS+1\r
-                 DO 30 J=L,H\r
-                     N=J-L\r
-                     K=LOD+N\r
-                     S=SOLUN(J)\r
-                     WRITE(OUCH1,991) N,S,S*COLSC(K),LGTOL/CINFN(J)\r
-       30        CONTINUE\r
-                 IF (ACTIN(I) .EQ. -1) THEN\r
-                     WRITE(OUCH1,*)'ACTION: REDUCE DEGREE TO ',NEWDG(I),' ***'\r
-                 ELSE IF (ACTIN(I) .EQ. 0) THEN\r
-                     WRITE(OUCH1,*)'ACTION: NONE            ***'\r
-                 ELSE IF (ACTIN(I) .EQ. 1) THEN\r
-                     WRITE(OUCH1,*)'ACTION: INCREASE DEGREE TO ',NEWDG(I)\r
-                 ELSE\r
-                     WRITE(OUCH1,*)'ACTION: SUBDIVIDE THIS ARC'\r
-                 ENDIF\r
-       40    CONTINUE\r
-       C\r
-             WRITE(OUCH1,*) 'KAPPA =',SOLUN(NEQNS)\r
-             IF (.NOT.INTER) THEN\r
-                 CAP=EXP(-SOLUN(NEQNS))\r
-                 WRITE(OUCH1,*) 'CAPACITY = ',CAP\r
-             ENDIF\r
-       C\r
-       990   FORMAT(A,T7,A,T26,A,T44,A)\r
-       991   FORMAT(I3,T6,E15.8,T25,E15.8,T44,E10.3)\r
-       992   FORMAT(E15.8)\r
-       993   FORMAT(I3,T8,E15.8,'  (',E14.7,',',E14.7,')')\r
-       994   FORMAT(A,T8,A,T34,A)\r
-       995   FORMAT(A,T6,A,T23,A,T36,A)\r
-       996   FORMAT(I2,T6,E14.7,T23,F10.5,T36,E14.7)\r
-       997   FORMAT(A,T24,A,T40,A)\r
-       998   FORMAT(T5,I3,T26,I3,T45,E9.2)\r
-       999   FORMAT(A,E10.2)\r
-       1000  FORMAT(A,1X,I5)\r
-       C\r
-    }*/\r
+    private void RSLT71(int QIERC[], double RCOND,double SOLUN[], int NEQNS, int LOSUB[],\r
+        int HISUB[], double COLSC[], int NQPTS, int JATYP[], int PARNT[], int TNSUA, boolean INTER,\r
+        double MQERR, double MCQER, double CINFN[], int ACTIN[], int NEWDG[], int NJIND,\r
+        double JACIN[], int NQUAD[], double TOLOU[], double LGTOL, int SOLCO) {\r
+       //INTEGER NEQNS,TNSUA,OUCH1,NQPTS,NJIND,NEWDG(*),NQUAD(*),LOSUB(*),\r
+       //+HISUB(*),QIERC(0:6),JATYP(*),PARNT(*),ACTIN(*),SOLCO\r
+       //REAL SOLUN(*),RCOND,COLSC(*),MQERR,MCQER,LGTOL,\r
+       //+CINFN(*),JACIN(*),TOLOU(*)\r
+       //LOGICAL INTER\r
+       String QTEXT[] = new String[7];\r
+       final String LINE = "_________________________________________________________________";\r
+       //CHARACTER QTEXT(0:6)*22,LINE*72\r
+             \r
+       // LOCAL VARIABLES\r
+       \r
+       int I,J,JI,K,L,LOD,N,H;\r
+       double S,CAP;\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(LINE+"\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("             SOLUTION NUMBER = " + SOLCO + "\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("                       NEQNS = " + NEQNS + "\n", Preferences.DEBUG_ALGORITHM); \r
+       Preferences.debug("RECIPROCAL COND NO. ESTIMATE = " + RCOND + "\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("   CONDITION NO. LOWER BOUND = " + (1.0/RCOND) + "\n", Preferences.DEBUG_ALGORITHM);\r
+       \r
+       Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("JACOBI INDEX           POINTS          TOLERANCE ACHIEVED\n", Preferences.DEBUG_ALGORITHM);\r
+       for (I=1; I <= NJIND; I++) {\r
+           Preferences.debug("     "+I+"                 "+NQUAD[I-1]+"           "+TOLOU[I-1]+"\n", Preferences.DEBUG_ALGORITHM);\r
+       } // for (I=1; I <= NJIND; I++)+\r
+       \r
+       \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
+       }      \r
+       \r
+       Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("              MAXIMUM QAWS ERROR = " + MQERR + "\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM COMPOSITE GAUSSIAN ERROR = " + MCQER + "\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+             \r
+       for (I=1; I <= TNSUA; I++) {\r
+               Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+           Preferences.debug("SUB ARC = " + I + " ON PARENT ARC " + PARNT[I-1] + "\n", Preferences.DEBUG_ALGORITHM);\r
+           Preferences.debug("N     SCALED SOLUN       UNSCALED SOLUN    IGNORE LEVEL\n", Preferences.DEBUG_ALGORITHM);\r
+           L=LOSUB[I-1];\r
+           H=HISUB[I-1];\r
+           JI=Math.abs(JATYP[I-1]);\r
+           LOD=(JI-1)*NQPTS+1;\r
+           for (J=L; J <= H; J++) {\r
+               N=J-L;\r
+               K=LOD+N;\r
+               S=SOLUN[J-1];\r
+               Preferences.debug(N+"  "+S+"    "+(S*COLSC[K-1])+"    "+(LGTOL/CINFN[J-1])+ "\n", Preferences.DEBUG_ALGORITHM);\r
+           } // for (J=L; J <= H; J++)\r
+           if (ACTIN[I-1] == -1) {\r
+               Preferences.debug("ACTION: REDUCE DEGREE TO " + NEWDG[I-1] + " ***\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+           else if (ACTIN[I-1] == 0) {\r
+               Preferences.debug("ACTION: NONE            ***\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+           else if (ACTIN[I-1] == 1) {\r
+               Preferences.debug("ACTION: INCREASE DEGREE TO " + NEWDG[I-1]+ "\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+           else {\r
+               Preferences.debug("ACTION: SUBDIVIDE THIS ARC\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+       Preferences.debug("KAPPA = " + SOLUN[NEQNS-1] + "\n", Preferences.DEBUG_ALGORITHM);\r
+       if (!INTER) {\r
+           CAP=Math.exp(-SOLUN[NEQNS-1]);\r
+           Preferences.debug("CAPACITY = " + CAP + "\n", Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       \r
+    } // private void RSLT71\r
     \r
     private void RSLT72(int QIERC[], double RCOND, double GAMMA, int NEQNS, int DGPOL[],\r
         int JATYP[], int PARNT[], int TNSUA, boolean INTER, double MQERR, double MCQER,\r
         int ACTIN[], int NEWDG[], int NJIND, int NQUAD[], double TOLOU[], double LGTOL,\r
-       int SOLCO, int OUCH1) {\r
+       int SOLCO) {\r
        // INTEGER NEQNS,TNSUA,OUCH1,NJIND,NEWDG(*),NQUAD(*),QIERC(0:6),\r
        // PARNT(*),ACTIN(*),DGPOL(*),JATYP(*),SOLCO\r
        // REAL GAMMA,RCOND,MQERR,MCQER,LGTOL,TOLOU(*)\r
@@ -5626,56 +5716,50 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        Preferences.debug("RECIPROCAL COND NO. ESTIMATE = " + RCOND + "\n", Preferences.DEBUG_ALGORITHM);\r
        Preferences.debug("   CONDITION NO. LOWER BOUND = " + (1.0/RCOND) + "\n", Preferences.DEBUG_ALGORITHM);\r
        \r
-       /*      WRITE(OUCH1,*) \r
-             WRITE(OUCH1,997) 'JACOBI INDEX','POINTS','TOLERANCE ACHIEVED'\r
-             DO 10 I=1,NJIND\r
-               WRITE(OUCH1,998) I,NQUAD(I),TOLOU(I)\r
-       10    CONTINUE\r
-       C\r
-             WRITE(OUCH1,*) \r
-             WRITE(OUCH1,*) 'QAWS TERMINATIONS WITH......' \r
-             DO 20 I=0,6\r
-               IF (QIERC(I) .GT. 0) THEN\r
-                 WRITE(OUCH1,1000) QTEXT(I),QIERC(I)\r
-               ENDIF\r
-       20    CONTINUE\r
-       C\r
-             WRITE(OUCH1,*) \r
-             WRITE(OUCH1,999) '              MAXIMUM QAWS ERROR =',MQERR\r
-             WRITE(OUCH1,999) 'MAXIMUM COMPOSITE GAUSSIAN ERROR =',MCQER\r
-             WRITE(OUCH1,*) \r
-             WRITE(OUCH1,992) 'SUB ARC','PARENT ARC','TYPE','CURRENT DEGREE','\r
-            +ACTION'\r
-             DO 40 I=1,TNSUA\r
-               IF (ACTIN(I) .EQ. -1) THEN\r
-                 WRITE(OUCH1,993) I,PARNT(I),JATYP(I),DGPOL(I),'REDUCE TO ',NEW\r
-            +                     DG(I)\r
-               ELSE IF (ACTIN(I) .EQ. 0) THEN\r
-                 WRITE(OUCH1,994) I,PARNT(I),JATYP(I),DGPOL(I),'NONE'\r
-               ELSE IF (ACTIN(I) .EQ. 1) THEN\r
-                 WRITE(OUCH1,993) I,PARNT(I),JATYP(I),DGPOL(I),'INCREASE TO ',N\r
-            +                     EWDG(I)\r
-               ELSE\r
-                 WRITE(OUCH1,994) I,PARNT(I),JATYP(I),DGPOL(I),'SUBDIVIDE'\r
-               ENDIF\r
-       40    CONTINUE\r
-       C\r
-             WRITE(OUCH1,*) 'KAPPA =',GAMMA\r
-             IF (.NOT.INTER) THEN\r
-                 CAP=EXP(-GAMMA)\r
-                 WRITE(OUCH1,*) 'CAPACITY = ',CAP\r
-             ENDIF\r
-       C\r
-       990   FORMAT(A,T7,A,T26,A,T44,A)\r
-       991   FORMAT(I3,T6,E15.8,T25,E15.8,T44,E10.3)\r
-       992   FORMAT(A,T10,A,T24,A,T34,A,T53,A)\r
-       993   FORMAT(T2,I2,T14,I3,T25,I3,T40,I2,T53,A,1X,I2)\r
-       994   FORMAT(T2,I2,T14,I3,T25,I3,T40,I2,T53,A)\r
-       997   FORMAT(A,T24,A,T40,A)\r
-       998   FORMAT(T5,I3,T26,I3,T45,E9.2)\r
-       999   FORMAT(A,E10.2)\r
-       1000  FORMAT(A,1X,I5)\r
-       */\r
+       Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("JACOBI INDEX           POINTS          TOLERANCE ACHIEVED\n", Preferences.DEBUG_ALGORITHM);\r
+       for (I=1; I <= NJIND; I++) {\r
+           Preferences.debug("     "+I+"                 "+NQUAD[I-1]+"           "+TOLOU[I-1]+"\n", Preferences.DEBUG_ALGORITHM);\r
+       } // for (I=1; I <= NJIND; I++)+\r
+       \r
+       \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
+       }\r
+       \r
+       Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("              MAXIMUM QAWS ERROR = " + MQERR + "\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM COMPOSITE GAUSSIAN ERROR = " + MCQER + "\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("SUB ARC   PARENT ARC   TYPE     CURRENT DEGREE     ACTION\n", Preferences.DEBUG_ALGORITHM);\r
+       for (I=1; I <= TNSUA; I++) {\r
+           if (ACTIN[I-1] == -1) {\r
+                 Preferences.debug("  "+I+"         "+PARNT[I-1]+"        "+JATYP[I-1]+"            "+\r
+                                   DGPOL[I-1]+" REDUCE TO " + NEWDG[I-1] + "\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+           else if (ACTIN[I-1] == 0) {\r
+               Preferences.debug("  "+I+"         "+PARNT[I-1]+"        "+JATYP[I-1]+"            "+\r
+                        DGPOL[I-1]+" NONE " + "\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+           else if (ACTIN[I-1] == 1) {\r
+               Preferences.debug("  "+I+"         "+PARNT[I-1]+"        "+JATYP[I-1]+"            "+\r
+                        DGPOL[I-1]+" INCREASE TO " + NEWDG[I-1] + "\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+           else {\r
+               Preferences.debug("  "+I+"         "+PARNT[I-1]+"        "+JATYP[I-1]+"            "+\r
+                        DGPOL[I-1]+" SUBDIVIDE " + "\n", Preferences.DEBUG_ALGORITHM);\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+       Preferences.debug("KAPPA = " + GAMMA + "\n", Preferences.DEBUG_ALGORITHM);\r
+       if (!INTER) {\r
+           CAP=Math.exp(-GAMMA);\r
+           Preferences.debug("CAPACITY = " + CAP + "\n", Preferences.DEBUG_ALGORITHM);\r
+       }\r
     } // private void RSLT72\r
 \r
     \r
@@ -6593,6 +6677,632 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         IER[0]=0;\r
 \r
     } // private void STATS1\r
+    \r
+    private void UPJAC1(int NQPTS,int NJIND,int INDEG,int AXION[],int DGPOL[],\r
+        int NEWDG[],double ACOEF[],double BCOEF[],double DIAG[],double SDIAG[],\r
+        int TNSUA,int MNSUA,int LOSUB[],int HISUB[],int JATYP[],int PARNT[],\r
+        double MIDPT[],double HALEN[],double COLPR[], double ZCOLL[][],\r
+        boolean LNSEG[],boolean PNEWQ[],double EPS,int IER[],double WORK[],\r
+        double NEWHL[],double RCOPY[],int ICOPY[],boolean LCOPY[], int LOOLD[],\r
+       int HIOLD[]) {\r
+       //INTEGER NQPTS,INDEG,TNSUA,MNSUA,IER,NJIND\r
+       //INTEGER DGPOL(*),LOSUB(*),HISUB(*),JATYP(*),PARNT(*),ICOPY(*),\r
+       //+AXION(*),NEWDG(*),LOOLD(*),HIOLD(*)\r
+       //REAL EPS,ACOEF(*),BCOEF(*),DIAG(*),SDIAG(*),WORK(*),MIDPT(*),\r
+       //+HALEN(*),COLPR(*),RCOPY(*),NEWHL(*)\r
+       //COMPLEX ZCOLL(*)\r
+       //LOGICAL LNSEG(*),LCOPY(*),PNEWQ(*)\r
+       \r
+       // TO UPDATE THE COLLOCATION PARAMETERS (STORED IN COLPR), THE \r
+       // COLLOCATION POINTS ON THE PHYSICAL BOUNDARY (STORED IN ZCOLL) \r
+       // AND THE ARRAYS LOSUB AND HISUB NEEDED TO ACCESS THIS DATA \r
+       // CORRECTLY.  \r
+       // ALSO TO UPDATE/DETERMINE THE ARRAYS\r
+       //   JATYP - THE JACOBI INDEX TYPE OF EACH SUBARC\r
+       //   PARNT - THE PARENT ARC OF EACH SUBARC\r
+       //   MIDPT - THE GLOBAL PARAMETRIC MIDPOINT OF EACH SUBARC\r
+       //   HALEN - THE GLOBAL PARAMETRIC HALF-LENGTH OF EACH SUBARC\r
+       //   DGPOL - THE POLYNOMIAL DEGREE ON EACH SUBARC\r
+       //   LNSEG - THE LINE SEGMENT BOOLEAN FOR EACH SUBARC\r
+       //   PNEWQ - BOOLEAN INDICATING POSSIBLE NEW QUADRATURE FOR SUBARC\r
+       //   IER=0 - NORMAL TERMINATION\r
+       //   IER=7 - FAILURE IN IMTQLH\r
+       //   IER=17- NUMBER OF SUBARCS REQUIRED EXCEEDS MNSUA\r
+       \r
+       // LOCAL VARIABLES\r
+       \r
+       int IFAIL[] = new int[1];\r
+       int D,D1,FIRST,I,J,K,K1,K2,P,PREV,JT,NTNSA,J1,J2;\r
+       double S,TC,MD,HH,F1,F2;\r
+       double PIN[] = new double[2];\r
+       double POUT[];\r
+       //COMPLEX PARFUN\r
+       final boolean USEIN = true;\r
+       boolean LS;\r
+       // EXTERNAL PARFUN,IMTQLH\r
+       \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           MD=MIDPT[I-1];\r
+           if (AXION[I-1] < 2) {\r
+               J=J+1;\r
+               if (J > MNSUA) {\r
+                   IER[0]=17;\r
+                   return;\r
+               }\r
+               RCOPY[J-1]=MD;\r
+           } // if (AXION[I-1] < 2)\r
+           else { // AXION[I-1] >= 2\r
+               JT=JATYP[I-1];\r
+               if (JT > 0) {\r
+                   F1=1.0-NEWHL[I-1];\r
+               }\r
+               else {\r
+                   F1=NEWHL[I-1];\r
+               }\r
+               F2=1.0-F1;\r
+               HH=HALEN[I-1];\r
+               J=J+1;\r
+               if (J > MNSUA) {\r
+                   IER[0]=17;\r
+                   return;\r
+               }\r
+               RCOPY[J-1]=MD-F1*HH;\r
+               J=J+1;\r
+               if (J > MNSUA) {\r
+                   IER[0]=17;\r
+                   return;\r
+               }\r
+               RCOPY[J-1]=MD+F2*HH;\r
+           } // else AXION[I-1] >= 2\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       NTNSA=J;\r
+       for (I=1; I <= NTNSA; I++) {\r
+           MIDPT[I-1]=RCOPY[I-1];\r
+       }\r
+            \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           HH=HALEN[I-1];      \r
+           if (AXION[I-1] < 2) {\r
+               J=J+1;\r
+               RCOPY[J-1]=HH;\r
+           }\r
+           else {\r
+               JT=JATYP[I-1];\r
+               if (JT > 0) {\r
+                   F1=NEWHL[I-1];\r
+               }\r
+               else {\r
+                   F1=1.0-NEWHL[I-1];\r
+               }\r
+               F2=1.0-F1;\r
+               J=J+1;\r
+               RCOPY[J-1]=F1*HH;\r
+               J=J+1;\r
+               RCOPY[J-1]=F2*HH;\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       for (I=1; I <= NTNSA; I++) {\r
+           HALEN[I-1]=RCOPY[I-1];\r
+       } // for (I=1; I <= NTNSA; I++)\r
+             \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           JT=JATYP[I-1];\r
+           if (AXION[I-1] < 2) {\r
+               J=J+1;\r
+               ICOPY[J-1]=JT;\r
+           }\r
+           else { \r
+               if (JT < NJIND) {\r
+                   if (JT > 0) {\r
+                       J1=JT;\r
+                       J2=NJIND;\r
+                   }\r
+                   else {\r
+                       J1=NJIND;\r
+                       J2=JT;\r
+                   }\r
+               }\r
+               else {\r
+                   J1=NJIND;\r
+                   J2=J1;\r
+               }\r
+               J=J+1;\r
+               ICOPY[J-1]=J1;\r
+               J=J+1;\r
+               ICOPY[J-1]=J2;\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       for (I=1; I <= NTNSA; I++) {\r
+           JATYP[I-1]=ICOPY[I-1];\r
+       }\r
+       \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           if (AXION[I-1] < 2) {\r
+               J=J+1;\r
+               ICOPY[J-1]=PARNT[I-1];\r
+           }\r
+           else {\r
+               J=J+1;\r
+               ICOPY[J-1]=PARNT[I-1];\r
+               J=J+1;\r
+               ICOPY[J-1]=PARNT[I-1];\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       for (I=1; I <= NTNSA; I++) {\r
+           PARNT[I-1]=ICOPY[I-1];\r
+       }\r
+             \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           LS=LNSEG[I-1];       \r
+           if (AXION[I-1] < 2) {\r
+               J=J+1;\r
+               LCOPY[J-1]=LS;\r
+           }\r
+           else {\r
+               J=J+1;\r
+               LCOPY[J-1]=LS;\r
+               J=J+1;\r
+               LCOPY[J-1]=LS;\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       for (I=1; I <= NTNSA; I++) {\r
+           LNSEG[I-1]=LCOPY[I-1];\r
+       }\r
+         \r
+       if (USEIN) { \r
+       \r
+           // USE INDEG ON SUBDIVIDED ARCS\r
+          \r
+           J=0;\r
+           for (I=1; I <= TNSUA; I++) {\r
+               if (AXION[I-1] < 2) {\r
+                   J=J+1;\r
+                   DGPOL[J-1]=NEWDG[I-1];\r
+               }\r
+               else {\r
+                   J=J+1;\r
+                   DGPOL[J-1]=INDEG;\r
+                   J=J+1;\r
+                   DGPOL[J-1]=INDEG;\r
+               }\r
+           } // for (I=1; I <= TNSUA; I++)\r
+       } // if (USEIN)\r
+       else {\r
+       \r
+           // USE CURRENT DEGREE ON SUBDIVIDED ARCS\r
+       \r
+           J=0;\r
+           for (I=1; I <= TNSUA; I++) {\r
+               if (AXION[I-1] < 2) {\r
+                   J=J+1;\r
+                   ICOPY[J-1]=NEWDG[I-1];\r
+               }\r
+               else {\r
+                   J=J+1;\r
+                   ICOPY[J-1]=DGPOL[I-1];\r
+                   J=J+1;\r
+                   ICOPY[J-1]=DGPOL[I-1];\r
+               }\r
+           } // for (I=1; I <= TNSUA; I++) \r
+           for (I=1; I <= NTNSA; I++) {\r
+               DGPOL[I-1]=ICOPY[I-1];\r
+           }\r
+       } // else !USEIN\r
+       \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           if (AXION[I-1] == 2) {\r
+               J=J+1;\r
+               LOOLD[J-1]=0;\r
+               HIOLD[J-1]=-1;     \r
+               J=J+1;\r
+               LOOLD[J-1]=0;\r
+               HIOLD[J-1]=-1;\r
+           }\r
+           else if (AXION[I-1] == 0) {\r
+               J=J+1;\r
+               LOOLD[J-1]=LOSUB[I-1];\r
+               HIOLD[J-1]=HISUB[I-1];\r
+           }\r
+           else {     \r
+               J=J+1;\r
+               LOOLD[J-1]=0;\r
+               HIOLD[J-1]=-1;\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+       J=0;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           if (AXION[I-1] < 2) {\r
+               J=J+1;\r
+               LCOPY[J-1]=true;\r
+           }\r
+           else {\r
+               J=J+1;\r
+               LCOPY[J-1]=false;\r
+               J=J+1;\r
+               LCOPY[J-1]=false;\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+       if (LCOPY[0] && LCOPY[NTNSA-1]) {\r
+           if (LCOPY[1]) {\r
+               PNEWQ[0]=false;\r
+           }\r
+           else {\r
+               PNEWQ[0]=true;\r
+           }\r
+           if (LCOPY[NTNSA-2]) {\r
+               PNEWQ[NTNSA-1]=false;\r
+           }\r
+           else {\r
+               PNEWQ[NTNSA-1]=true;\r
+           }\r
+       }\r
+       else {\r
+           PNEWQ[0]=true;\r
+           PNEWQ[NTNSA-1]=true;\r
+       }\r
+       \r
+       J=NTNSA-1;\r
+       for (I=2; I <= J; I++) {\r
+           if (LCOPY[I-2] && LCOPY[I-1] && LCOPY[I]) {\r
+               PNEWQ[I-1]=false;\r
+           }\r
+           else {\r
+               PNEWQ[I-1]=true;\r
+           }\r
+       } // for (I=2; I <= J; I++)\r
+       \r
+       TNSUA=NTNSA;\r
+       LOSUB[0]=1;\r
+       HISUB[0]=1+DGPOL[0];\r
+       for (I=2; I <= TNSUA; I++) {\r
+           LOSUB[I-1]=HISUB[I-2]+1;\r
+           HISUB[I-1]=LOSUB[I-1]+DGPOL[I-1];\r
+       } // for (I=2; I <= TNSUA; I++)\r
+       \r
+       for (I=1; I <= TNSUA; I++) {\r
+           J=JATYP[I-1];\r
+           P=PARNT[I-1];\r
+           D=DGPOL[I-1];\r
+           D1=D+1;\r
+           if (J > 0) {\r
+               S=1.0;\r
+           }\r
+           else {\r
+               S=-1.0;\r
+               J=-J;\r
+           }\r
+           PREV=(J-1)*NQPTS;\r
+           FIRST=LOSUB[I-1];\r
+           for (K=1; K <= D1; K++) {\r
+               WORK[K-1]=0.0;\r
+               K1=PREV+K;\r
+               DIAG[K-1]=BCOEF[K1-1];\r
+               if (K == 1) {\r
+                   SDIAG[K-1]=0.0;\r
+               }\r
+               else {\r
+                   SDIAG[K-1]=ACOEF[K1-2];\r
+               }\r
+           } // for (K=1; K <= D1; K++) \r
+           WORK[0]=1.0;\r
+           IMTQLH(D1,DIAG,SDIAG,IFAIL);\r
+           if (IFAIL[0] > 0) {\r
+               IER[0]=7;\r
+               return;\r
+           }\r
+           for (K=1; K <= D1; K++) {\r
+               TC=S*DIAG[K-1];\r
+               K2=FIRST+K-1;\r
+               COLPR[K2-1]=TC;\r
+               TC=MIDPT[I-1]+HALEN[I-1]*TC;\r
+               PIN[0] = TC;\r
+               PIN[1] = 0.0;\r
+               POUT = PARFUN(P,PIN);\r
+               ZCOLL[K2-1][0] = POUT[0];\r
+               ZCOLL[K2-1][1] = POUT[1];\r
+           } // for (K=1; K <= D1; K++)\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+       // NORMAL EXIT\r
+       \r
+        IER[0]=0;\r
+       \r
+    } // private void UPJAC1\r
+    \r
+    private void UPCOQ1(int NARCS,int NJIND,int NQPTS,int MDGPO,int MQIN1,double AQTOL,\r
+        double QUPTS[], double QUWTS[],double JACIN[],double MIDPT[],double HALEN[],\r
+        double ACOEF[],double BCOEF[],double H0VAL[],double COLSC[],int NQUAD[],\r
+        int LOQSB[],double QCOMX[],double QCOMW[],int MNQUA,double TOLOU[], double MCQER[],\r
+        double XENPT[], double XIVAL[][], double XIDST[],int TNSUA, boolean PNEWQ[],\r
+        boolean NEWQU[],int JATYP[],int PARNT[],boolean NUQTL[],int IER[]) {\r
+       //INTEGER NARCS,NQPTS,MDGPO,MQIN1,TNSUA,IER,NQUAD(*),LOQSB(*),\r
+       //+NJIND,JATYP(*),PARNT(*),MNQUA\r
+       //REAL AQTOL,QUPTS(*),QUWTS(*),JACIN(*),MIDPT(*),HALEN(*),ACOEF(*),\r
+       //+BCOEF(*),H0VAL(*),COLSC(*),QCOMX(*),QCOMW(*),TOLOU(*),XENPT(*),\r
+       //+XIDST(*),MCQER\r
+       //LOGICAL NUQTL\r
+       //LOGICAL PNEWQ(*),NEWQU(*)\r
+       //COMPLEX XIVAL(*)\r
+       \r
+       //     THE PURPOSE OF THIS ROUTINE IS TO UPDATE THE ABSCISSAE \r
+       //     (QCOMX) AND WEIGHTS (QCOMW) FOR THE COMPOSITE GAUSSIAN RULES \r
+       //     FOR THE ESTIMATION OF\r
+       \r
+       //       INTEGRAL  [(1+X)**BETA*P(X,I)*LOG|ZZ-X|*dX], I=0,1,...,MDGPO.\r
+       //      -1<=X<=1                                      J=1,NZZ\r
+       \r
+       //     HERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I\r
+       //     ASSOCIATED WITH THE WEIGHT (1+X)**BETA AND ZZ IS ANY COLLOCATION\r
+       //     POINT PREIMAGE NOT ON [-1,1].  BETA TAKES ON THE VARIOUS VALUES\r
+       //     DEFINED BY ARRAY JACIN.  THE ROUTINE ALSO COMPUTES\r
+       \r
+       //     NQUAD - NQUAD(I) IS THE NUMBER OF QUADRATURE POINTS IN THE\r
+       //             COMPOSITE RULE FOR BETA=JACIN(I).\r
+       //     LOQSB - THE ABSCISSAE AND WEIGHTS OF THE COMPOSITE RULE FOR\r
+       //             BETA=JACIN(I) ARE STORED IN ARRAYS QCOMX AND QCOMW IN\r
+       //             THE POSITIONS LOQSB(I) TO LOQSB(I)+NQUAD(I)-1 INCLUSIVE.\r
+       //     XIDST,\r
+       //     XIVAL - XIVAL(2*I-1) STORES THE COLLOCATION PREIMAGE THOUGHT\r
+       //             TO BE NEAREST TO -1 AND XIDST(2*I-1) STORES ITS DISTANCE\r
+       //             FROM -1; SIMILARLY, XIVAL(2*I) STORES THE PREIMAGE\r
+       //             THOUGHT TO BE NEAREST TO +1 AND XIDST(2*I) ITS DISTANCE\r
+       //             FROM +1. THE PREIMAGES ARE WITH RESPECT TO\r
+       //             THE PARAMETRIC FUNCTIONS DEFINING THE SUBARCS WHICH\r
+       //             MEET AT THE PHYSICAL CORNER WHERE BETA=JACIN(I).\r
+       //     TOLOU - TOLOU(I) IS THE ESTIMATED MAXIMUM ERROR OVER ALL\r
+       //             COLLOCATION POINTS IN USING THE COMPOSITE RULE\r
+       //             FOR BETA=JACIN(I).\r
+       //     IER  - IER=0 FOR NORMAL TERMINATION.\r
+       //            IER=19 THE REQUIRED TOTAL NUMBER OF COMPOSITE QUADRATURE\r
+       //                   POINTS EXCEEDS THE LIMIT MNQUA.\r
+       //            IER=20 THE PARAMETER MQIN1 NEEDS INCREASING; MQIN1-1 IS \r
+       //                   THE MAXIMUM ALLOWED NUMBER OF SUBINTERVALS IN OUR\r
+       //                   COMPOSITE GAUSSIAN RULE. (SAME AS IER=11, BUT 20\r
+       //                   IDICATES DURING REFINEMENT PROCESS)\r
+       \r
+       //     ALL THE ABOVE (APART FROM IER) SHOULD HAVE EXISTING VALUES ON \r
+       //     INPUT WHICH ARE UPDATED BY THE NEW VERSIONS ON OUTPUT.\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       int I,I0,I1,I2,J,K,JI,JI0,JI1,JI2,P0,P1,P2,HI,LO,QINTS,NQ,DIFF,TNCQP;\r
+       double DST[] = new double[2];\r
+       double BETA,H1,M1,T0,T2,SUM1,RR,RRB,MEAN,RXI,IXI;\r
+       final double ONE[] = new double[]{1.0,0.0};\r
+       double ZZ[] = new double[2];\r
+       double Z0[] = new double[2];\r
+       double Z2[] = new double[2];\r
+       double XI[][] = new double[2][2];\r
+       //COMPLEX ONE,ZZ,Z0,Z2,XI(2),PARFUN,DPARFN\r
+       //PARAMETER (ONE=(1E+0,0E+0))\r
+        //EXTERNAL PARFUN,DPARFN,SUBIN7\r
+       double PIN[] = new double[2];\r
+       double POUT[];\r
+       double DOUT[];\r
+       double cr[] = new double[1];\r
+       double ci[] = new double[1];\r
+       \r
+       //**** NEWQU(J) IS TRUE IF THE QUADRATURE RULE FOR THE J'TH JACOBI INDEX\r
+       //**** NEEDS UPDATING.  FIRST SET NEWQU FOR THE CASE WHERE A NEW\r
+       //**** QUADRATURE TOLERANCE HAS BEEN FIXED; THIS IS INDEPENDENT OF ANY\r
+       //**** POSSIBLE ARC SUBDIVISIONS.  IF THE PURE LEGENDRE RULE DOESN'T\r
+       //**** ALREADY EXIST THEN IT DOESN'T HAVE TO BE UPDATED.\r
+       \r
+       for (J=1; J <= NARCS; J++) {\r
+           NEWQU[J-1]=NUQTL[0];\r
+       }\r
+       NEWQU[NJIND-1]=(NUQTL[0] && (NQUAD[NJIND-1] > 0));\r
+       \r
+       //**** NEXT OVERWRITE NEWQU TO PICK UP THOSE CASES WHERE A BOUNDARY\r
+       //**** SUBDIVISION HAS OCCURRED AND UPDATE THE NEAR POINT VECTOR XIVAL.\r
+       \r
+       for (I1=1; I1 <= TNSUA; I1++) {\r
+           if (PNEWQ[I1-1]) {\r
+               if (I1 == 1) {\r
+                   I0=TNSUA;\r
+               }\r
+               else { \r
+                   I0=I1-1;\r
+               }\r
+       \r
+               if (I1 == TNSUA) {\r
+                   I2=1;\r
+               }\r
+               else {\r
+                   I2=I1+1;\r
+               }\r
+       \r
+               JI0=JATYP[I0-1];\r
+               JI1=JATYP[I1-1];\r
+               JI2=JATYP[I2-1];\r
+       \r
+               if (JI0 > 0) {\r
+                   T0=QUPTS[JI0*NQPTS-1];\r
+               }\r
+               else {\r
+                   JI0=-JI0;\r
+                   T0=-QUPTS[(JI0-1)*NQPTS];\r
+               }\r
+       \r
+               if (JI2 > 0) {\r
+                   T2=QUPTS[(JI2-1)*NQPTS];\r
+               }\r
+               else {\r
+                   JI2=-JI2;\r
+                   T2=-QUPTS[JI2*NQPTS-1];\r
+               }\r
+       \r
+               T0=MIDPT[I0-1]+T0*HALEN[I0-1];\r
+               T2=MIDPT[I2-1]+T2*HALEN[I2-1];\r
+               P0=PARNT[I0-1];\r
+               P1=PARNT[I1-1];\r
+               P2=PARNT[I2-1];\r
+               PIN[0] = T0;\r
+               PIN[1] = 0.0;\r
+               Z0=PARFUN(P0,PIN);\r
+               PIN[0] = T2;\r
+               PIN[1] = 0.0;\r
+               Z2=PARFUN(P2,PIN);\r
+               H1=HALEN[I1-1];\r
+               M1=MIDPT[I1-1];\r
+               ZZ[0]=M1-H1;\r
+               ZZ[1] = 0.0;\r
+               POUT = PARFUN(P1,ZZ);\r
+               DOUT = DPARFN(P1,ZZ);\r
+               zdiv(POUT[0] - Z0[0], POUT[1] - Z0[1],DOUT[0],DOUT[1],cr,ci);\r
+               XI[0][0] = -1.0 - cr[0]/H1;\r
+               XI[0][1] = -ci[0]/H1;\r
+               ZZ[0] = M1+H1;\r
+               ZZ[1] = 0.0;\r
+               POUT = PARFUN(P1,ZZ);\r
+               DOUT = DPARFN(P1,ZZ);\r
+               zdiv(POUT[0]-Z2[0], POUT[1]-Z2[1],DOUT[0],DOUT[1],cr,ci);\r
+               XI[1][0] = 1.0 - cr[0]/H1;\r
+               XI[1][1] = -ci[0]/H1;\r
+       \r
+               if (JI1 < 0) {\r
+                   Z0[0]=XI[0][0];\r
+                   Z0[1]=XI[0][1];\r
+                   XI[0][0]=-XI[1][0];\r
+                   XI[0][1]=-XI[1][1];\r
+                   XI[1][0]=-Z0[0];\r
+                   XI[1][1]=-Z0[1];\r
+                   JI1=-JI1;\r
+               } // if (JI1 < 0)\r
+       \r
+               for (J=1; J <= 2; J++) {\r
+                   RXI=XI[J-1][0];\r
+                   IXI=XI[J-1][1];\r
+                   if (-1.0 <= RXI && RXI <= 1.0) {\r
+                       DST[J-1]=Math.abs(IXI);\r
+                   }\r
+                   else if (RXI < -1.0) {\r
+                       DST[J-1]=zabs(XI[J-1][0]+1.0,XI[J-1][1]);\r
+                   }\r
+                   else {\r
+                       DST[J-1]=zabs(XI[J-1][0]-1.0,XI[J-1][1]);\r
+                   }\r
+               } // for (J=1; J <= 2; J++)\r
+       \r
+               J=2*JI1-2;\r
+               for (I=1; I <= 2; I++) {\r
+                   J=J+1;\r
+                   if (DST[I-1] < XIDST[J-1]) {\r
+                       NEWQU[JI1-1]=true;\r
+                       XIVAL[J-1][0]=XI[I-1][0];\r
+                       XIVAL[J-1][1]=XI[I-1][1];\r
+                       XIDST[J-1]=DST[I-1];\r
+                   }\r
+               } // for (I=1; I <= 2; I++) \r
+           } // if (PNEWQ[I1-1])\r
+       } // for (I1=1; I1 <= TNSUA; I1++)\r
+       \r
+       //     FOR THOSE INDECES FOR WHICH NEWQU IS TRUE WE NOW SET UP THE NEW\r
+       //     COMPOSITE GAUSSIAN QUADRATURE DATA.\r
+       \r
+       TNCQP=LOQSB[NJIND-1]+NQUAD[NJIND-1]-1;\r
+       HI=0;\r
+       /*for (JI=1; JI <= NJIND; JI++) {\r
+           NQ=NQUAD[JI-1];\r
+           if (NEWQU[JI-1]) {\r
+               LO=(JI-1)*NQPTS+1;\r
+               BETA=JACIN[JI-1];\r
+               I2=2*JI-1;\r
+               SUBIN7(XIVAL(I2),2,BETA,MDGPO,NQPTS,ACOEF(LO),BCOEF(LO),\r
+                     H0VAL(JI),COLSC(LO),AQTOL,TOLOU(JI),XENPT,QINTS,MQIN1,\r
+                     IER);\r
+                 IF (IER .GT. 0) THEN\r
+                   IF (IER .EQ. 11) THEN\r
+                     IER=20\r
+                   ENDIF\r
+                   RETURN\r
+                 ENDIF\r
+       C\r
+                 DIFF=QINTS*NQPTS-NQ\r
+                 IF (TNCQP+DIFF .GT. MNQUA) THEN\r
+                   IER=19\r
+                   RETURN\r
+                 ENDIF\r
+                 I1=HI+NQ+1\r
+       C \r
+       C         IF DIFF IS NON-ZERO WE MUST MAKE SPACE IN ARRAYS QCOMX AND\r
+       C         QCOMW TO RECEIVE THE NEW DATA.\r
+       C\r
+                 IF (DIFF .GT. 0) THEN\r
+                   DO 40 I=TNCQP,I1,-1\r
+                     J=I+DIFF\r
+                     QCOMX(J)=QCOMX(I)\r
+                     QCOMW(J)=QCOMW(I)\r
+       40          CONTINUE\r
+                 ELSE IF (DIFF .LT. 0) THEN\r
+                   DO 50 I=I1,TNCQP\r
+                     J=I+DIFF\r
+                     QCOMX(J)=QCOMX(I)\r
+                     QCOMW(J)=QCOMW(I)\r
+       50          CONTINUE\r
+                 ENDIF\r
+       C\r
+       C         NOW SET UP THE NEW RULE AND STORE DATA IN QCOMX, QCOMW\r
+       C\r
+                 TNCQP=TNCQP+DIFF\r
+                 NQUAD(JI)=NQ+DIFF\r
+                 LOQSB(JI)=HI+1\r
+                 SUM1=BETA+1E+0\r
+                 K=HI\r
+                 DO 80 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
+                     LO=LO-1\r
+                     DO 60 J=1,NQPTS\r
+                       K=K+1\r
+                       QCOMX(K)=MEAN+RR*QUPTS(LO+J)\r
+                       QCOMW(K)=RRB*QUWTS(LO+J)\r
+       60            CONTINUE\r
+                   ELSE\r
+                     LO=NARCS*NQPTS\r
+                     DO 70 J=1,NQPTS\r
+                       K=K+1\r
+                       QCOMX(K)=MEAN+RR*QUPTS(LO+J)\r
+                       QCOMW(K)=RR*QUWTS(LO+J)*(1E+0+QCOMX(K))**BETA\r
+       70            CONTINUE\r
+                   ENDIF\r
+       80        CONTINUE\r
+                 HI=HI+NQUAD(JI)\r
+           } // if (NEWQU[JI-1])\r
+           else {\r
+       C\r
+       C         HERE WE DO NOTHING OTHER THAN UPDATE SOME SUBSCRIPTS.\r
+       C\r
+                 LOQSB(JI)=HI+1\r
+                 HI=HI+NQ\r
+           }\r
+       C\r
+       } // for (JI=1; JI <= NJIND; JI++)\r
+       C\r
+             MCQER=0E+0\r
+             DO 100 I=1,NJIND\r
+               MCQER=MAX(MCQER,TOLOU(I))\r
+       100   CONTINUE\r
+       C\r
+             NUQTL=.FALSE.\r
+       C\r
+       C     NORMAL TERMINATION\r
+       C\r
+             IER=0\r
+       C */\r
+    } // private void UPOCOQ1\r
+\r
+\r
 \r
       /**\r
        * zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.\r