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

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

index ba88cb1..9dd318c 100644 (file)
@@ -6,6 +6,8 @@ import java.io.RandomAccessFile;
 import java.util.Scanner;\r
 \r
 import gov.nih.mipav.model.structures.ModelImage;\r
+import gov.nih.mipav.model.structures.jama.GeneralizedEigenvalue;\r
+import gov.nih.mipav.model.structures.jama.LinearEquations2;\r
 import gov.nih.mipav.view.MipavUtil;\r
 import gov.nih.mipav.view.Preferences;\r
 \r
@@ -89,33 +91,33 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     Scanner input = new Scanner(System.in);\r
     private double zzset[][] = new double[400][2];\r
     private int IBNDS[] = new int[5];\r
-    private int IGEOM[] = new int[4];\r
-    private int PARNT[] = new int[IBNDS[0]];\r
-    private double RGEOM[] = new double[2];\r
-    private double HALEN[]= new double[IBNDS[0]];\r
-       private double MIDPT[] = new double[IBNDS[0]];\r
-       private double VTARG[] = new double[IBNDS[0]]; \r
-    private int ISNPH[] = new int[6];\r
-    private int DGPOL[] = new int[IBNDS[0]];\r
-       private int JATYP[] = new int[IBNDS[0]];\r
-       private int LOSUB[] = new int[IBNDS[0]];\r
+    private int IGEOM[] = new int[4]; // Written in original JAPHYC\r
+    private int PARNT[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
+    private double RGEOM[] = new double[2]; // Written in original JAPHYC\r
+    private double HALEN[]= new double[IBNDS[0]]; // Written in original JAPHYC\r
+       private double MIDPT[] = new double[IBNDS[0]]; // Written in original JAPHYC\r
+       private double VTARG[] = new double[IBNDS[0]]; // Written in original JPAHYC\r
+    private int ISNPH[] = new int[6]; // Written in original JAPHYC\r
+    private int DGPOL[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
+       private int JATYP[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
+       private int LOSUB[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
        // Parts of RSNPH\r
        private int NQPTS;\r
        private int NJIND = NARCS + 1;\r
        private int TNGQP = NQPTS * NJIND;\r
        private int MNEQN;\r
-       private double ACOEF[] = new double[TNGQP];\r
-       private double BCOEF[] = new double[TNGQP];\r
-       private double AICOF[] = new double[TNGQP];\r
-       private double BICOF[] = new double[TNGQP];\r
-       private double QUPTS[] = new double[TNGQP];\r
-       private double QUWTS[] = new double[TNGQP];\r
-       private double H0VAL[] = new double[NJIND];\r
-       private double HIVAL[] = new double[NJIND];\r
-       private double JACIN[] = new double[NJIND];\r
-       private double ERARC[] = new double[IBNDS[0]];\r
-       private double BCFSN[] = new double[MNEQN];\r
-       private double SOLUN[] = new double[MNEQN];\r
+       private double ACOEF[] = new double[TNGQP]; // Written in original JAPHYC\r
+       private double BCOEF[] = new double[TNGQP]; // Written in original JAPHYC\r
+       private double AICOF[] = new double[TNGQP]; // Written in original JAPHYC\r
+       private double BICOF[] = new double[TNGQP]; // Written in original JAPHYC\r
+       private double QUPTS[] = new double[TNGQP]; // Written in original JAPHYC\r
+       private double QUWTS[] = new double[TNGQP]; // Written in original JAPHYC\r
+       private double H0VAL[] = new double[NJIND]; // Written in original JAPHYC\r
+       private double HIVAL[] = new double[NJIND]; // Written in original JAPHYC\r
+       private double JACIN[] = new double[NJIND]; // Written in original JAPHYC\r
+       private double ERARC[] = new double[IBNDS[0]]; // Written in original JAPHYC\r
+       private double BCFSN[] = new double[MNEQN]; // Written in original JAPHYC\r
+       private double SOLUN[] = new double[MNEQN]; // Written in original JAPHYC\r
        // Parts of IWORK\r
        private int IPIVT[] = new int[MNEQN];\r
        private int LOQSB[] = new int[IBNDS[1]];\r
@@ -159,6 +161,19 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        private boolean PNEWQ[] = new boolean[IBNDS[0]];\r
        private boolean LNSEG[] = new boolean[IBNDS[0]];\r
        \r
+//  CENTR  - COMPLEX\r
+//                     THE POINT IN THE PHYSICAL PLANE THAT IS TO BE\r
+//                     MAPPED TO THE CENTRE OF THE UNIT DISC.  FOR\r
+//                     EXTERIOR DOMAINS CENTR MUST BE SOME POINT IN THE\r
+//                     COMPLEMENTARY INTERIOR  PHYSICAL DOMAIN.\r
+//                     IN CASE ABS(ISYGP).GT.1 THEN CENTR MUST ALSO BE\r
+//                     A CENTRE OF SYMMETRY FOR THE PHYSICAL DOMAIN.\r
+       private double CENTR[] = new double[2]; // Written in original JAPHYC\r
+//  INTER  - LOGICAL\r
+//                     TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE \r
+//                     OTHERWISE.\r
+       private boolean INTER; // Written in in original JAPHYC\r
+       \r
        //COMMON /FNDEF/\r
        private double BETA;\r
        private double A1;\r
@@ -2331,9 +2346,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        return result;\r
    }\r
 \r
-   private void JAPHYC(String JBNM, String HEAD, double MAXER,boolean INTER,\r
+   private void JAPHYC(String JBNM, String HEAD, double MAXER,\r
                      int ISYGP, boolean INCST,\r
-                     int RFARC, double RFARG[], double CENTR[], int TSTNG,int OULVL,\r
+                     int RFARC, double RFARG[], int TSTNG,int OULVL,\r
                      double MATRX[][][],\r
                      int OCH,\r
                      int IER[]) {\r
@@ -2657,15 +2672,18 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                final double SFACT = 0.8;\r
                final double QFACT = 0.1;\r
                double MCQER[] = new double[1];\r
-               double AQTOL,CONST,GAQTL,GLGTL,GRQTL,GSUPE,GTGTE,LGTOL,ESTOL,\r
-                    MCHEP,MQERR,MXERM,PI,R1MACH,RCOND,RQTOL,SSUPE,\r
+               double RCOOND[] = 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
                     TGTER,TOLNR;\r
                \r
                double ZMXER[] = new double[2];\r
                //COMPLEX ZMXER\r
                \r
                final boolean INIBT = true;\r
-               boolean ACCPT,GACPT,NUQTL,REFLN;\r
+               boolean ACCPT[] = new boolean[1];\r
+               boolean GACPT,NUQTL,REFLN;\r
                \r
                String OFL;\r
                //CHARACTER OFL*6\r
@@ -2895,7 +2913,7 @@ 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
                \r
@@ -2916,63 +2934,126 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    //**** SOLVE LINEAR SYSTEM BY GAUSSIAN ELIMINATION USING LINPACK\r
                    // LAPACK equivalent of LINPACK SGECO are DLANGE, DGETRF, and DGECON\r
                    // LU factorization and condition estimation of a general matrix.\r
-                   SGECO(MATRX(1,1,2),MNEQN,NROWS,IWORK(IPIVT),RCOND,RWORK(WORK2));\r
-                     IF (RCOND .EQ. 0E+0) THEN\r
-                       IER=15\r
-                       SOLCO=SOLCO-1\r
-                       GOTO 999\r
-                     ENDIF    \r
-                     CALL SGESL(MATRX(1,1,2),MNEQN,NROWS,IWORK(IPIVT),RSNPH(SOLUN),0)\r
-                     NEFF=NEFF+NROWS**3\r
-                     WRITE(*,1) 'LINEAR SYSTEM SOLUTION DONE:'\r
-               C\r
-               C**** RECONSTITUTE FULL SOLUTION VECTOR\r
-               C\r
-                     IF (ORDSG.GT.1) THEN\r
-                       CALL RECON(ORDSG,REFLN,NCOLL,TNSUA,ISNPH(LOSUB),IWORK(HISUB),\r
-                    +  RSNPH(SOLUN))\r
-                     ENDIF\r
-                     CONST=RSNPH(SOLUN+NEQNS-1)\r
-               C\r
-               C**** SET UP THE ARRAY WORK2 OF ACTUAL COEFFICIENT IGNORE LEVELS\r
-               C\r
-                     CALL SETIGL(RWORK(WORK2),IWORK(HISUB),ISNPH(JATYP),ISNPH(LOSUB),\r
-                    +NQPTS,RWORK(RIGLL),TNSUA)\r
-               C\r
-               C**** DETERMINE THE ACTIONS THAT HAVE TO BE TAKEN ON EACH ARC\r
-               C\r
-                     CALL AXION1(IWORK(AXION),IWORK(NEWDG),RSNPH(SOLUN),MDGPO,TNSUA,\r
-                    +ISNPH(DGPOL),ISNPH(LOSUB),IWORK(HISUB),RWORK(RIGLL),LGTOL,ACCPT,\r
-                    +RSNPH(JACIN),ISNPH(JATYP),NJIND,RWORK(NEWHL),ESTOL,IER)\r
-                     ESTOL=ESTOL/SFACT\r
-                     IF (IER.GT.0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-                     WRITE(*,1) 'DECISIONS FOR EACH ARC DONE:'\r
-                     WRITE(*,3) 'EFFECTIVE STOPPING TOLERANCE:',ESTOL\r
-                     IF (ACCPT .AND. ESTOL.LE.GSUPE) THEN\r
-                       GACPT=.TRUE.\r
-                     ELSE\r
-                       GACPT=.FALSE.\r
-                     ENDIF\r
-               C\r
+                   double A[][] = new double[MNEQN][MNEQN];\r
+                   for (I = 0; I < MNEQN; I++) {\r
+                       for (J = 0; J < MNEQN; J++) {\r
+                               A[I][J] = MATRX[I][J][1];\r
+                       }\r
+                   }\r
+                   LinearEquations2 le2 = new LinearEquations2();\r
+           GeneralizedEigenvalue ge = new GeneralizedEigenvalue();\r
+           double anorm;\r
+           int iwork[] = new int[NROWS];\r
+           int info[] = new int[1];\r
+           anorm = ge.dlange('1', MNEQN, MNEQN, A, MNEQN, WORK2);\r
+           le2.dgetrf(MNEQN, MNEQN, A, MNEQN, IPIVT, info);\r
+           if (info[0] > 0) {\r
+               MipavUtil.displayError("In dgetrf factor U is exactly singular");\r
+           }\r
+           le2.dgecon('1', NROWS, A, MNEQN, anorm, RCOND, WORK2, iwork, info);\r
+           // Reciprocal condition number of A in 1-norm.\r
+                   //SGECO(MATRX(1,1,2),MNEQN,NROWS,IPIVT,RCOND,WORK2);\r
+                   if (RCOND[0] == 0.0) {\r
+                       IER[0]=15;\r
+                       SOLCO=SOLCO-1;\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,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] == 0.0) \r
+                   // 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
+                   //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
+                               MATRX[I][J][1] = A[I][J];\r
+                       }\r
+                   }\r
+                   NEFF=NEFF+NROWS*NROWS*NROWS;\r
+                   System.out.println("LINEAR SYSTEM SOLUTION DONE:");\r
+               \r
+                   // **** RECONSTITUTE FULL SOLUTION VECTOR\r
+               \r
+                   if (ORDSG > 1) {\r
+                       RECON(ORDSG,REFLN,NCOLL,TNSUA[0],LOSUB,HISUB,SOLUN);\r
+                   } // if (ORDSG > 1)\r
+                   CONST=SOLUN[NEQNS-1];\r
+               \r
+                   // **** SET UP THE ARRAY WORK2 OF ACTUAL COEFFICIENT IGNORE LEVELS\r
+               \r
+                   SETIGL(WORK2,HISUB,JATYP,LOSUB,\r
+                    NQPTS,RIGLL,TNSUA[0]);\r
+               \r
+                    // **** DETERMINE THE ACTIONS THAT HAVE TO BE TAKEN ON EACH ARC\r
+               \r
+                    AXION1(AXION,NEWDG,SOLUN,MDGPO,TNSUA[0],\r
+                    DGPOL,LOSUB,HISUB,RIGLL,LGTOL,ACCPT,\r
+                    JACIN,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,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 DONE:");\r
+                    System.out.println("EFFECTIVE STOPPING TOLERANCE: " + ESTOL[0]);\r
+                    if (ACCPT[0] && ESTOL <= GSUPE) {\r
+                        GACPT=true;\r
+                    }\r
+                    else {\r
+                        GACPT=false;\r
+                    }\r
+               \r
                      if (GACPT) {\r
-               C\r
-               C****   OUTPUT RESULTS\r
-               C\r
-                       IF (OULVL .LT. 4) THEN\r
+               \r
+                         // 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
-                       ELSE \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
-                       ENDIF\r
-                       WRITE(OCH,*) 'EFFECTIVE STOPPING TOLERANCE :',ESTOL\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
@@ -3166,8 +3247,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        CLOSE(OCH)\r
                      ENDIF\r
                  \r
-                     CALL WRTAIL(1,0,IER)*/\r
-               \r
+                     CALL WRTAIL(1,0,IER)\r
+               */\r
     } // private void JAPHYC\r
 \r
     private void ANGLE7(double BE[], int NA, boolean IN) {\r
@@ -5516,6 +5597,88 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        C\r
     }*/\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
+       // 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
+       // LOGICAL INTER\r
+        \r
+       // LOCAL VARIABLES\r
+       \r
+       int I;\r
+       double CAP;\r
+       // CHARACTER QTEXT(0:6)*22,LINE*72\r
+       String QTEXT[] = new String[7];\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(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
+       /*      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
+    } // private void RSLT72\r
+\r
+    \r
     private void LNSY11(double MATRX[][][], double RGTSD[], int MNEQN, int NCOLL, int ORDSG, boolean REFLN,\r
         int NQPTS, int TNSUA, int JATYP[], int PARNT[], int DGPOL[], int LOSUB[], int HISUB[], int NQUAD[],\r
         int LOQSB[], double TOLNR, double MIDPT[], double HALEN[], double H0VAL[], double COLSC[], double ACOEF[],\r
@@ -5948,9 +6111,489 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         return result;\r
     } // private double FNVAL\r
 \r
+    private void BCFVTF(double BCFSN[], double VTARG[], int DGPOL[], int JATYP[], int LOSUB[], int PARNT[],\r
+        int RFARC, int TNSUA, double H0VAL[], double JACIN[], double RFARG, double SOLUN[]) {\r
+       // INTEGER RFARC,TNSUA\r
+       // INTEGER DGPOL(*),JATYP(*),LOSUB(*),PARNT(*)\r
+       // REAL RFARG\r
+       // REAL BCFSN(*),H0VAL(*),JACIN(*),SOLUN(*),VTARG(*)\r
+       \r
+       //     PERFORMS VARIOUS PRELIMINARY TASKS TO PREPARE FOR THE POST-\r
+       //     PROCESSING QUADRATURE CALCULATIONS.\r
+       \r
+       //     SETS UP THE ARRAY OF COEFFICIENTS BCFSN NEEDED FOR THE \r
+       //     CALCULATION OF THE BOUNDARY CORRESPONDENCE FUNCTIONS FOR THE MAP\r
+       //     PHYSICAL --> CANONICAL; THESE ARE SIMPLY RELATED TO THE SOLUTION \r
+       //     COEFFICIENT ARRAY SOLUN.\r
+       \r
+       //     CALCULATES THE ARGUMENTS VTARG OF THE IMAGES ON THE UNIT CIRCLE\r
+       //     OF THE END POINTS OF ALL SUBARCS, ENFORCING THE USER'S ROTATION\r
+       //     CONDITION THAT PARFUN(RFARC,(-1.0,0.0)) BE MAPPED TO THE POINT\r
+       //     WITH ARGUMENT RFARG.\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       int AJT,DEG,I,J,J1,JT,LO;\r
+       double B1,RTH0,THET0,TUPI;\r
+       \r
+        TUPI=2.0*Math.PI;\r
+       VTARG[0]=0.0;\r
+       \r
+       for (I=1; I <= TNSUA; I++) {\r
+           JT=JATYP[I-1];\r
+           AJT=Math.abs(JT);\r
+           RTH0=Math.sqrt(H0VAL[AJT-1]);\r
+           B1=JACIN[AJT-1]+1.0;\r
+           DEG=DGPOL[I-1];\r
+           LO=LOSUB[I-1];\r
+           VTARG[I]=VTARG[I-1]+SOLUN[LO-1]*TUPI*RTH0;\r
+       \r
+           BCFSN[LO-1]=TUPI*SOLUN[LO-1]/(B1*RTH0);\r
+           for (J=1; J <= DEG; J++) {\r
+               J1=LO+J;\r
+               BCFSN[J1-1]=TUPI*SOLUN[J1-1]/Math.sqrt(J*(J+B1));\r
+           } // for (J=1; J <= DEG; J++)\r
+       \r
+           if (JT < 0) {\r
+               for (J=1; J <= DEG; J += 2) {\r
+                   J1=LO+J;\r
+                   BCFSN[J1-1]=-BCFSN[J1-1];\r
+               } // for (J=1; J <= DEG; J += 2)\r
+           } // if (JT < 0)\r
+       \r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+       I=1;\r
+       while (PARNT[I-1] != RFARC) {\r
+               I = I+1;\r
+       }\r
+       THET0=RFARG-VTARG[I-1];\r
+       \r
+       for (I = 1; I <= TNSUA; I++) {\r
+           VTARG[I-1]=VTARG[I-1]+THET0;\r
+       }\r
+       VTARG[TNSUA]=VTARG[0]+TUPI;\r
+       \r
+    } // private void BCFVTF\r
+\r
+    private void RECON(int ORDSG, boolean REFLN, int NCOLL,int TNSUA, int LOSUB[],\r
+               int HISUB[], double SOLUN[]) {\r
+        // INTEGER ORDSG,NCOLL,TNSUA\r
+        // INTEGER LOSUB(*),HISUB(*)\r
+        // REAL SOLUN(*)\r
+        // LOGICAL REFLN\r
 \r
+        // TO RECONSTITUTE THE FULL SOLUTION VECTOR FOR THE WHOLE BOUNDARY\r
+        // FROM THE SOLUTION ON THE FUNDAMENTAL BOUNDARY SECTION.\r
+\r
+        // LOCAL VARIABLES\r
+\r
+        int IA,IS,J,J1,J2,J3,LOM,NCFBS,ORDRG,TSFBS;\r
+        double SG;\r
+\r
+        if (ORDSG == 1) return;\r
+\r
+        TSFBS=TNSUA/ORDSG;\r
+        NCFBS=NCOLL/ORDSG;\r
+\r
+        SOLUN[NCOLL]=SOLUN[NCFBS];\r
+\r
+        if (REFLN) {\r
+            for (IA=1; IA <= TSFBS; IA++) {\r
+                J1=LOSUB[IA-1];\r
+                J2=HISUB[IA-1];\r
+                LOM=LOSUB[2*TSFBS-IA]-J1;\r
+                SG=-1.0;\r
+                for (J=J1; J <= J2; J++) {\r
+                    SG=-SG;\r
+                    J3=LOM+J;\r
+                    SOLUN[J3-1]=SG*SOLUN[J-1];\r
+                } // for (J=J1; J <= J2; J++)\r
+            } // for (IA=1; IA <= TSFBS; IA++)\r
+            ORDRG=ORDSG/2;\r
+            for (IS=2; IS <= ORDRG; IS++) {\r
+                LOM=(IS-1)*NCFBS*2;\r
+                for (J=1; J <= NCFBS*2; J++) {\r
+                    J1=LOM+J;\r
+                    SOLUN[J1-1]=SOLUN[J-1];\r
+                } // for (J=1; J <= NCFBS*2; J++)\r
+            } // for (IS=2; IS <= ORDRG; IS++)\r
+        } // if (REFLN)\r
+        else {\r
+            for (IS=2; IS <= ORDSG; IS++) {\r
+                LOM=(IS-1)*NCFBS;\r
+                for (J=1; J <= NCFBS; J++) {\r
+                    J1=LOM+J;\r
+                    SOLUN[J1-1]=SOLUN[J-1];\r
+                } // for (J=1; J <= NCFBS; J++)\r
+            } // for (IS=2; IS <= ORDSG; IS++)\r
+        } // else\r
+\r
+    } // private void RECON\r
+\r
+    private void SETIGL(double AIGLL[], int HISUB[], int JATYP[], int LOSUB[],\r
+               int NQPTS, double RIGLL[], int TNSUA) {\r
+        // INTEGER NQPTS,TNSUA\r
+        // INTEGER HISUB(*),JATYP(*),LOSUB(*)\r
+        // REAL AIGLL(*),RIGLL(*)\r
+\r
+        // **** COPY THE REFERENCE IGNORE LEVELS *RIGLL* INTO THE ACTUAL IGNORE\r
+        // **** IGNORE LEVEL ARRAY *AIGLL*\r
+\r
+        // LOCAL VARIABLES\r
+\r
+        int AJT,HIM,I,IA,LOD,LOM;\r
+\r
+        for (IA=1; IA <= TNSUA; IA++) {\r
+            AJT=Math.abs(JATYP[IA-1]);\r
+            LOD=(AJT-1)*NQPTS+1;\r
+            LOM=LOSUB[IA-1];\r
+            HIM=HISUB[IA-1];\r
+            for (I=LOM; I <= HIM; I++) {\r
+                AIGLL[I-1]=RIGLL[LOD+I-LOM-1];\r
+            } // for (I=LOM; I <= HIM; I++)\r
+        } // for (IA=1; IA <= TNSUA; IA++)\r
+\r
+    } // private void SETIGL\r
+\r
+    private void AXION1(int AXION[], int NEWDG[], double SOLUN[], int MDGPO,int TNSUA,\r
+        int DGPOL[], int LOSUB[], int HISUB[], double RIGLL[], double LGTOL, boolean ACCPT[],\r
+        double JACIN[], int JATYP[], int NJIND, double NEWHL[], double ESTOL[], int IER[]) {\r
+       // INTEGER AXION(*),NEWDG(*),MDGPO,TNSUA,DGPOL(*),LOSUB(*),HISUB(*),\r
+       //+JATYP(*),NJIND,IER\r
+       //REAL SOLUN(*),RIGLL(*),LGTOL,JACIN(*),NEWHL(*),ESTOL\r
+       //LOGICAL ACCPT\r
+       \r
+       //     TO DETERMINE THE ARRAY "AXION" WHICH SPECIFIES THE ACTIONS THAT \r
+       //     ARE TO TAKEN ON EACH SUBARC, AS FOLLOWS:\r
+       //        AXION(I)=-1 - REDUCE THE DEGREE ON ARC I\r
+       //        AXION(I)=0  - LEAVE THE DEGREE ON ARC I UNCHANGED\r
+       //        AXION(I)=1  - INCREASE THE DEGREE ON ARC I\r
+       //        AXION(I)=2  - SUBDIVIDE ARC I.\r
+       \r
+       //     IN CASE ABS(AXION(I))=1, NEWDG(I) RECORDS THE NEW DEGREE TO BE\r
+       //     USED ON ARC I.\r
+       \r
+       //     IN CASE AXION(I)<=0 FOR ALL I=1,TNSUA, THE SOLUTION IS DEEMED TO\r
+       //     BE ACCEPTABLE TO THE REQUIRED ACCURACY AND WE SET ACCPT=.TRUE.;\r
+       //     ACCPT=.FALSE. OTHERWISE.\r
+       \r
+       //     ALSO TO DETERMINE THE EFFECTIVE STOPPING TOLERANCE ESTOL; IF THE\r
+       //     USER HAD INPUT ESTOL AS THE VALUE FOR THE ARGUMENT MAXER IN\r
+       //     JAPHYC, THEN THE CURRENT SOLUTION WOULD BE ACCEPTED.  \r
+           \r
+       //     IER=0  - NORMAL EXIT\r
+       //     IER=54 - LOCAL PARAMETER MXCO MUST BE INCREASED\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MINDG = 0;\r
+       final int MXCO = 32;\r
+       int IFL[] = new int[1];\r
+       int AJT,CI,DG,I,IA,J,LOD,LOM,PNDG;\r
+       final double SAFEF = 5.0;\r
+       double AA[] = new double[1];\r
+       double POW[] = new double[1];\r
+       double EXA[] = new double[1];\r
+       double CONF[] = new double[1];\r
+       double BETA,LAM,TERM,THSLN,XX,VAR;\r
+       double COVAR[][] = new double[2][2];\r
+       final boolean CONSV = false;\r
+       // INTEGER CRITCO\r
+       // EXTERNAL CRITCO,STATS1\r
+       double POSCO[] = new double[MXCO];\r
+       \r
+       ESTOL[0]=0.0;\r
+       for (IA=1; IA <= TNSUA; IA++) {\r
+           DG=DGPOL[IA-1];\r
+           if (DG+1 > MXCO) {\r
+               IER[0]=54; \r
+               return;\r
+           }\r
+           LOM=LOSUB[IA-1];\r
+           AJT=Math.abs(JATYP[IA-1]);\r
+           LOD=(AJT-1)*(MDGPO+1)+1;\r
+           for (I=0; I <= DG; I++) {\r
+               J=LOM+I;\r
+               POSCO[I]=Math.abs(SOLUN[J-1])*RIGLL[LOD+I-1]/LGTOL;\r
+           } // for (I= 0; I <= DG; I++)   \r
+           CI=CRITCO(DG+1,POSCO)-1;  \r
+           PNDG=CI+1+MINDG;\r
+           if (CI == -1) {\r
+       \r
+               // ALL COEFFICIENTS ARE ACCEPTABLE.\r
+       \r
+               if (DG == MINDG) {\r
+                   AXION[IA-1]=0;\r
+                   NEWDG[IA-1]=DG;\r
+               } // if (DG == MINDG)\r
+               else {\r
+       \r
+                   // PROBABLY DECREASE THE DEGREE, BUT ONLY IGNORE THOSE\r
+                   // COEFFICIENTS WHICH ARE 'SAFELY' BELOW THE TOLERANCE\r
+                   // LIMIT.\r
+       \r
+                   for (I=0; I <= DG; I++) {\r
+                       POSCO[I]=POSCO[I]*SAFEF;\r
+                   } // for (I=0; I <= DG; I++)\r
+                   PNDG=CRITCO(DG+1,POSCO)+MINDG;\r
+                   if (PNDG >= DG) {\r
+                       AXION[IA-1]=0;\r
+                       NEWDG[IA-1]=DG;\r
+                   }\r
+                   else if (PNDG <= MINDG) {\r
+                       AXION[IA-1]=-1;\r
+                       NEWDG[IA-1]=MINDG;\r
+                   }\r
+                   else {\r
+                       AXION[IA-1]=-1;\r
+                       NEWDG[IA-1]=PNDG;\r
+                   }\r
+               } // else\r
+           } // if (CI == -1)\r
+           else if (PNDG == DG) {\r
+               AXION[IA-1]=0;\r
+               NEWDG[IA-1]=DG;\r
+           } // else if (PNDG == DG)\r
+           else if (PNDG < DG) {\r
+       \r
+               // PROBABLY DECREASE THE DEGREE, BUT ONLY IGNORE THOSE\r
+               // COEFFICIENTS WHICH ARE 'SAFELY' BELOW THE TOLERANCE\r
+               // LIMIT.\r
+               for (I=0; I <= DG; I++) {\r
+                   POSCO[I]=POSCO[I]*SAFEF;\r
+               }\r
+               PNDG=CRITCO(DG+1,POSCO)+MINDG;\r
+               if (PNDG >= DG) {\r
+                   AXION[IA-1]=0;\r
+                   NEWDG[IA-1]=DG;\r
+               }\r
+               else if (PNDG <= MINDG) {\r
+                   AXION[IA-1]=-1;\r
+                   NEWDG[IA-1]=MINDG;\r
+               }\r
+               else {\r
+                   AXION[IA-1]=-1;\r
+                   NEWDG[IA-1]=PNDG;\r
+               }\r
+           } // else if (PNDG < DG)\r
+           else if (DG == MDGPO) {\r
+       \r
+               // ARC SUBDIVISION IS REQUIRED AND ASSIGN NEW HALF-LENGTH.\r
+       \r
+               AXION[IA-1]=2;\r
+               BETA=JACIN[AJT-1];\r
+               LAM=Math.min(1.0,1.0+BETA);\r
+               if (AJT != NJIND) {\r
+                   NEWHL[IA-1]=Math.pow((1.0/POSCO[CI]),(1.0/(1.0+LAM+BETA)));\r
+               }\r
+               else {\r
+                   NEWHL[IA-1]=0.5;\r
+               }\r
+           } // else if (DG == MDGPO)\r
+           else {\r
+       \r
+               // MUST DECIDE BETWEEN ARC SUBDIVISION ARC DEGREE INCREASE.\r
+               // FIRST MAKE CONSERVATIVE ESTIMATES FOR THE COEFFICIENTS\r
+               // FOR DEGREES DG+1 TO MDGPO. \r
+       \r
+               double SN[] = new double[Math.min(32,DG+1)];\r
+               for (I = 0; I < Math.min(32, DG+1); I++) {\r
+                       SN[I] = SOLUN[LOM+I-1];\r
+               }\r
+               STATS1(SN,DG+1,AA,POW,EXA,COVAR,CONF,IFL);\r
+               if (IFL[0] == 1) {\r
+       \r
+                    // NOT ENOUGH DATA TO MAKE ESTIMATES, WHICH PRESUMES DG\r
+                   // IS RATHER SMALL; THEREFORE INCREASE THE DEGREE.\r
+       \r
+                   AXION[IA-1]=1;\r
+                   NEWDG[IA-1]=Math.min(DG+2,MDGPO);\r
+               } // if (IFL == 1)\r
+               else {\r
+                   for (I=DG+1; I <= MDGPO; I++) {\r
+                       if (CONSV) {\r
+                           XX=Math.log(1.0+I);\r
+                           VAR=COVAR[0][0]+XX*XX*COVAR[1][1]+2.0*XX*COVAR[0][1];\r
+                           THSLN=EXA[0]*Math.pow((I+1),POW[0])*Math.exp(CONF[0]*Math.sqrt(VAR));\r
+                       } // if (CONSV)\r
+                       else {\r
+                           THSLN=EXA[0]*Math.pow((I+1),POW[0]);\r
+                       }\r
+                       POSCO[I]=Math.abs(THSLN)*RIGLL[LOD+I-1]/LGTOL;\r
+                   } // for (I=DG+1; I <= MDGPO; I++)\r
+                   PNDG=CRITCO(MDGPO+1,POSCO)+MINDG;\r
+                   if (PNDG <= MDGPO) {\r
+       \r
+                       // INCREASE DEGREE\r
+       \r
+                       AXION[IA-1]=1;\r
+                       NEWDG[IA-1]=PNDG;\r
+                   }\r
+                   else {\r
+       \r
+                       // SUBDIVIDE ARC AND ASSIGN NEW HALF-LENGTH.\r
+       \r
+                       AXION[IA-1]=2;\r
+                       BETA=JACIN[AJT-1];\r
+                       LAM=Math.min(1.0,1.0+BETA);\r
+                       if (AJT != NJIND) {\r
+                           NEWHL[IA-1]=Math.pow((1.0/POSCO[CI]),(1.0/(1.0+LAM+BETA)));\r
+                       }\r
+                       else {\r
+                           NEWHL[IA-1]= 0.5;\r
+                       }\r
+                   }\r
+               }\r
+           } // else \r
+       \r
+           // NOW UPDATE THE EFFECTIVE STOPPING TOLERANCE\r
+       \r
+           J=HISUB[IA-1]-MINDG;\r
+           for (I=J; I <= HISUB[IA-1]; I++) {\r
+               TERM=Math.abs(SOLUN[I-1])*RIGLL[LOD+I-LOM-1];\r
+               ESTOL[0]=Math.max(ESTOL[0],Math.exp(TERM)-1.0);\r
+           } // for (I=J; I <= HISUB[IA-1]; I++)\r
+       } // for (IA=1; IA <= TNSUA; IA++)\r
+       \r
+       ACCPT[0]=true;\r
+       for (I=1; I <= TNSUA; I++) {\r
+           if (AXION[I-1] > 0) {\r
+               ACCPT[0]=false;\r
+               IER[0] = 0;\r
+               return;\r
+           }\r
+       } // for (I=1; I <= TNSUA; I++)\r
+       \r
+        IER[0]=0;\r
+       \r
+    } // private void AXION1\r
+    \r
+    private int CRITCO(int N, double POSCO[]) {\r
+        // INTEGER N\r
+        // REAL POSCO(*)\r
+\r
+        // GIVEN THE NON-NEGATIVE NUMBERS POSCO(I), I=1,...,N,  TO FIND THE\r
+        // INDEX CRITCO SUCH THAT POSCO(CRITCO) > 1 AND POSCO(I) <=1 FOR\r
+        // I=CRITCO+1,...,N.  IN CASE POSCO(I) <=1 FOR ALL I=1,...,N THEN\r
+        // POSCO=0.\r
+\r
+        // **** LOCAL VARIABLE\r
+      \r
+        int I;\r
+\r
+        I=N;\r
+        while (true) {\r
+            if (I == 0) {\r
+                return 0;\r
+            }\r
+            else if (POSCO[I-1] > 1.0) {\r
+               return I;\r
+            }\r
+            else {\r
+                I=I-1;\r
+                continue;\r
+            }\r
+        } // while true\r
+\r
+    } // private int CRITCO\r
+\r
+    private void STATS1(double SN[], int M, double A[], double B[], double EA[], double COV[][],\r
+               double CONF[], int IER[]) {\r
+        // INTEGER M,IER\r
+        // REAL A,B,CONF,EA\r
+        // REAL SN(*),COV(2,2)\r
+\r
+        // TO FIND THE LEAST SQUARES ESTIMATES A, B IN THE RELATIONSHIP\r
+\r
+        // LOG(SN(I)) = A + B*LOG(I) + ERROR, I=1,..,M, (M<=NDAT)\r
+\r
+        // AND ALSO TO GIVE EA=EXP(A) AND THE COVARIANCE MATRIX COV FOR\r
+        // ESTIMATES A,B.\r
+\r
+        // IER IS SET TO 1 IF LESS THAN 3 DATA PAIRS ARE COMPATIBLE WITH THE\r
+        // ABOVE MODEL AND NOTHING IS DONE.\r
+\r
+        // LOCAL VARIABLES..\r
+\r
+        final int NDAT = 32;\r
+       int I,N,N1;\r
+        double SX1,SX2,SY1,SXY,SE2,DET,DIFF,VAR;\r
+        double X[] = new double[NDAT];\r
+        double Y[] = new double[NDAT];\r
+        double par;\r
+\r
+        if (M < 3) {\r
+            IER[0]=1;\r
+            return;\r
+        }\r
+\r
+        N1=Math.min(M,NDAT);\r
+        N=N1;\r
+\r
+        for (I=1; I <= N1; I++) {\r
+            SY1=Math.abs(SN[I-1]);\r
+            if (SY1 > 0.0) {\r
+                Y[I-1]=Math.log(SY1);\r
+                X[I-1]=Math.log((double)(I));\r
+            }\r
+            else {\r
+                N=N-1;\r
+            }\r
+        } // for (I=1; I <= N1; I++)\r
+\r
+        if (N < 3) {\r
+            IER[0]=1;\r
+            return;\r
+        }\r
+\r
+        SX1=0.0;\r
+        SX2=0.0;\r
+        SY1=0.0;\r
+        SXY=0.0;\r
+        SE2=0.0;\r
+\r
+        for (I=1; I <= N; I++) {\r
+            SX1=SX1+X[I-1];\r
+            SX2=SX2+X[I-1]*X[I-1];\r
+            SY1=SY1+Y[I-1];\r
+            SXY=SXY+X[I-1]*Y[I-1];\r
+        } // for (I=1; I <= N; I++)\r
+\r
+        DET=N*SX2-SX1*SX1;\r
+        A[0]=(SX2*SY1-SX1*SXY)/DET;\r
+        EA[0]=Math.exp(A[0]);\r
+        B[0]=(N*SXY-SX1*SY1)/DET;\r
+\r
+        for (I=1; I <= N; I++) {\r
+               par = (Y[I-1]-A[0]-B[0]*X[I-1]);\r
+            SE2=SE2+par*par;\r
+        } // for (I=1; I <= N; I++)\r
+\r
+        SE2=SE2/(N-2)/DET;\r
+        COV[0][0]=SX2*SE2;\r
+        COV[0][1]=-SX1*SE2;\r
+        COV[1][0]=COV[0][1];\r
+        COV[1][1]=N*SE2;\r
+\r
+        CONF[0]=0.0;\r
+        I=N;\r
+           while (true) {\r
+               DIFF=Y[I-1]-A[0]-B[0]*X[I-1];\r
+               if (DIFF > 0.0) {\r
+                   VAR=COV[0][0]+X[I-1]*X[I-1]*COV[1][1]+2.0*X[I-1]*COV[0][1];\r
+                   CONF[0]=DIFF/Math.sqrt(VAR);\r
+                   break;\r
+               }\r
+               else {\r
+                   I=I-1;\r
+                   continue;\r
+               }\r
+           } // while (true)\r
+        IER[0]=0;\r
+\r
+    } // private void STATS1\r
 \r
\r
       /**\r
        * zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.\r
        * \r