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

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

index 2481726..ba88cb1 100644 (file)
@@ -158,7 +158,14 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        private boolean LCOPY[] = new boolean[IBNDS[0]];\r
        private boolean PNEWQ[] = new boolean[IBNDS[0]];\r
        private boolean LNSEG[] = new boolean[IBNDS[0]];\r
-    \r
+       \r
+       //COMMON /FNDEF/\r
+       private double BETA;\r
+       private double A1;\r
+       private double B1;\r
+       private double P0VAL;\r
+       private double SCALE;\r
+       private int TYPE;\r
        public SymmsIntegralMapping() {\r
                \r
        }\r
@@ -2637,9 +2644,10 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                     \r
                //     LOCAL VARAIBLES\r
                \r
-               int I,IMXER,INDEG,J,L,MDGPO,MNJXS,MNQUA,MNSUA,MQIN1,NCOLL,\r
+          int TNSUA[] = new int[1];\r
+          int I,IMXER,INDEG,J,L,MDGPO,MNJXS,MNQUA,MNSUA,MQIN1,NCOLL,\r
                     NEFF,NEQNS,NROWS,NTEST,\r
-                    TNSUA,ORDSG;\r
+                    ORDSG;\r
                int SOLCO = 0;\r
                int QIERC[] = new int[7];\r
                int QIERR[] = new int[7];\r
@@ -2648,8 +2656,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                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,MCQER,MQERR,MXERM,PI,R1MACH,RCOND,RQTOL,SSUPE,\r
+                    MCHEP,MQERR,MXERM,PI,R1MACH,RCOND,RQTOL,SSUPE,\r
                     TGTER,TOLNR;\r
                \r
                double ZMXER[] = new double[2];\r
@@ -2848,7 +2857,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                //**** HISUB NEEDED TO ACCESS COLPR AND ZCOLL CORRECTLY.  INITIALISE\r
                //**** DGPOL AND UPDATE LNSEG FOR ARC HALVING. \r
                \r
-               /*CPJAC3(NARCS,NQPTS,INDEG,DGPOL,JACIN,\r
+               CPJAC3(NARCS,NQPTS,INDEG,DGPOL,JACIN,\r
                     ACOEF,BCOEF,DIAG,SDIAG,TNSUA,\r
                     LOSUB,HISUB,JATYP,PARNT,MIDPT,\r
                     HALEN,COLPR,ZCOLL,LNSEG,LOOLD,\r
@@ -2857,7 +2866,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        WRTAIL(1,0,IER[0],null);\r
                    return; \r
                }\r
-               NCOLL=HISUB[TNSUA-1];\r
+               NCOLL=HISUB[TNSUA[0]-1];\r
                NEQNS=NCOLL+1;\r
                NROWS=NCOLL/ORDSG+1;\r
                if (NEQNS > MNEQN) {\r
@@ -2877,37 +2886,37 @@ 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
-                     IF (IER .GT. 0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-                     WRITE(*,1) 'COMPOSITE GAUSSIAN RULES DONE:'\r
-               C\r
-               C**** SET UP LINEAR ALGEBRAIC SYSTEM.\r
-               C\r
-               23    CONTINUE\r
-                     SOLCO=SOLCO+1\r
-                     WRITE(*,24) '********SOLUTION',SOLCO,'********',NROWS,'EQUATIONS'\r
-               24    FORMAT(/,T18,A,1X,I2,1X,A,/,T25,I3,1X,A)\r
-               C\r
-                     CALL LNSY11(MATRX,RSNPH(SOLUN),MNEQN,NCOLL,ORDSG,REFLN,NQPTS,\r
-                    +TNSUA,ISNPH(JATYP),IGEOM(PARNT),ISNPH(DGPOL),ISNPH(LOSUB),\r
-                    +IWORK(HISUB),IWORK(NQUAD),IWORK(LOQSB),TOLNR,RGEOM(MIDPT),\r
-                    +RGEOM(HALEN),RSNPH(H0VAL),RWORK(COLSC),RSNPH(ACOEF),RSNPH(BCOEF),\r
-                    +RWORK(COLPR),RWORK(QCOMX),RWORK(QCOMW),CENTR,ZWORK(ZCOLL),INTER,\r
-                    +LWORK(LNSEG),RWORK(WORK),QIERR,MQERR,RSNPH(JACIN),RWORK(A1COF),\r
-                    +RWORK(B1COF),AQTOL,RQTOL,RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),\r
-                    +IWORK(LOOLD),IWORK(HIOLD))\r
-               C\r
-                     DO 25 I=0,6\r
-                       QIERC(I)=QIERC(I)+QIERR(I)\r
-               25    CONTINUE\r
-                     WRITE(*,1) 'LINEAR SYSTEM SET UP DONE:'\r
-               C\r
-               C**** SOLVE LINEAR SYSTEM BY GAUSSIAN ELIMINATION USING LINPACK\r
-               C\r
-                     CALL SGECO(MATRX(1,1,2),MNEQN,NROWS,IWORK(IPIVT),RCOND,\r
-                    +RWORK(WORK2))\r
+               NUQTL=false;\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(1,0,IER[0],null);\r
+                   return;    \r
+               }\r
+               System.out.println("COMPOSITE GAUSSIAN RULES DONE:");\r
+               \r
+               //**** SET UP LINEAR ALGEBRAIC SYSTEM.\r
+               \r
+           /*while (true) {\r
+                   SOLCO=SOLCO+1;\r
+                   System.out.prinln("********SOLUTION "+ SOLCO+ " ******** " + NROWS + " EQUATIONS");\r
+               \r
+                   LNSY11(MATRX,SOLUN,MNEQN,NCOLL,ORDSG,REFLN,NQPTS,\r
+                    TNSUA[0],JATYP,PARNT,DGPOL,LOSUB,\r
+                    HISUB,NQUAD,LOQSB,TOLNR,MIDPT,\r
+                    HALEN,H0VAL,COLSC,ACOEF,BCOEF,\r
+                    COLPR,QCOMX,QCOMW,CENTR,ZCOLL,INTER,\r
+                    LNSEG,WORK,QIERR,MQERR,JACIN,A1COF,\r
+                    B1COF,AQTOL,RQTOL,AQCOF,BQCOF,CQCOF,\r
+                    LOOLD,HIOLD);\r
+               \r
+                   for (I=0; I <= 6; I++) {\r
+                       QIERC[I]=QIERC[I]+QIERR[I];\r
+                   }\r
+                   System.out.println("LINEAR SYSTEM SET UP DONE:");\r
+               \r
+                   //**** 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
@@ -2947,7 +2956,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        GACPT=.FALSE.\r
                      ENDIF\r
                C\r
-                     IF (GACPT) THEN\r
+                     if (GACPT) {\r
                C\r
                C****   OUTPUT RESULTS\r
                C\r
@@ -2980,7 +2989,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        WRITE(OCH,*) '****THE SOLUTION IS ACCEPTED****'\r
                        WRITE(OCH,*) 'EFFECTIVE SIZE OF ALL SYSTEMS : ',NEFF\r
                        WRITE(OCH,*) \r
-                     ELSE\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
@@ -3072,8 +3083,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                          GOTO 999\r
                        ENDIF\r
                        WRITE(*,1) 'QUADRATURE UPDATES DONE:'\r
-                       GOTO 23\r
-                     ENDIF \r
+                       continue;\r
+                     } // else !GACPT \r
+           } // while (true)\r
                C\r
                      IF (OULVL.EQ.2 .OR. OULVL.EQ.3) THEN\r
                        CALL RSLT71(QIERC,RCOND,RSNPH(SOLUN),NEQNS,ISNPH(LOSUB),\r
@@ -4177,7 +4189,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 \r
     private void CPJAC3(int NARCS,int NQPTS,int INDEG, int DGPOL[], double JACIN[],\r
                double ACOEF[], double BCOEF[], double DIAG[],\r
-            double SDIAG[],int TNSUA, int LOSUB[], int HISUB[], int JATYP[], int PARNT[],\r
+            double SDIAG[],int TNSUA[], int LOSUB[], int HISUB[], int JATYP[], int PARNT[],\r
             double MIDPT[], double HALEN[], double COLPR[], double ZCOLL[][], boolean LNSEG[],\r
             int LOOLD[], int HIOLD[], double EPS,int IER[], boolean INIBT) {\r
         // INTEGER NARCS,NQPTS,INDEG,TNSUA,IER\r
@@ -4212,7 +4224,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         // COMPLEX PARFUN\r
         // EXTERNAL PARFUN,IMTQLH,MDNBT\r
  \r
-        TNSUA=2*NARCS;\r
+        TNSUA[0]=2*NARCS;\r
         for (I=1; I <= NARCS; I++) {\r
             BETA=JACIN[I-1];\r
             if (I == NARCS) {\r
@@ -4250,23 +4262,23 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
             LNSEG[J-2]=LNSEG[I-1];\r
         } // for (I=NARCS; I >= 1; I--)\r
 \r
-        for (I=1; I <= TNSUA; I++) {\r
+        for (I=1; I <= TNSUA[0]; I++) {\r
             DGPOL[I-1]=INDEG;\r
         }\r
 \r
         LOSUB[0]=1;\r
         HISUB[0]=1+DGPOL[0];\r
-        for (I=2; I <= TNSUA; I++) {\r
+        for (I=2; I <= TNSUA[0]; 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
+        } // for (I=2; I <= TNSUA[0]; I++)\r
 \r
-        for (I=1; I <= TNSUA; I++) {\r
+        for (I=1; I <= TNSUA[0]; I++) {\r
             LOOLD[I-1]=0;\r
             HIOLD[I-1]=-1;\r
-        } // for (I=1; I <= TNSUA; I++)\r
+        } // for (I=1; I <= TNSUA[0]; I++)\r
 \r
-        for (I=1; I <= TNSUA; I++) {\r
+        for (I=1; I <= TNSUA[0]; I++) {\r
             J=JATYP[I-1];\r
             P=PARNT[I-1];\r
             D=DGPOL[I-1];\r
@@ -4304,7 +4316,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                 PIN[1] = 0.0;\r
                 ZCOLL[K2-1]=PARFUN(P,PIN);\r
             } // for (K=1; K <= D1; K++)\r
-        } // for (I=1; I <= TNSUA; I++)\r
+        } // for (I=1; I <= TNSUA[0]; I++)\r
 \r
         // NORMAL EXIT\r
 \r
@@ -4410,7 +4422,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     \r
     private void ICOQR1(int NARCS,int NJIND,int NQPTS,int MDGPO,int MQIN1, double AQTOL, double QUPTS[], double QUWTS[],\r
         double JACIN[], double MIDPT[], double HALEN[], double ACOEF[], double BCOEF[], double H0VAL[], double COLSC[],\r
-        int NQUAD[], int LOQSB[], double QCOMX[], double QCOMW[], int MNQUA, double TOLOU[], double MCQER, double XENPT[],\r
+        int NQUAD[], int LOQSB[], double QCOMX[], double QCOMW[], int MNQUA, double TOLOU[], double MCQER[], double XENPT[],\r
         double XIVAL[][], double XIDST[], int IER[]) {\r
         //INTEGER NARCS,NJIND,NQPTS,MDGPO,MQIN1,IER,NQUAD(*),LOQSB(*),\r
        //+MNQUA\r
@@ -4476,7 +4488,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        double ci[] = new double[1];\r
        \r
        HI=0;\r
-       /*for (JI=1; JI <= NARCS; JI++) {\r
+       for (JI=1; JI <= NARCS; JI++) {\r
        \r
            // AT THE JI'TH CORNER, THE ANALYTIC ARC IN THE BACKWARDS DIRECTION\r
            // IS THE JI0'TH, IN THE FORWARDS DIRECTION IT IS THE JI'TH AND THE\r
@@ -4640,73 +4652,83 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        AJAC[I] = ACOEF[LO-1+I];\r
                        BJAC[I] = BCOEF[LO-1+I];\r
                }\r
+               double CSCAL[] = new double[MDGPO+1];\r
+               for (I = 0; I < MDGPO+1; I++) {\r
+                       CSCAL[I] = COLSC[LO-1+I];\r
+               }\r
+               double TOLU[] =new double[1];\r
+               TOLU[0] = TOLOU[JI-1];\r
                SUBIN7(ZZ,2,BETA,MDGPO,NQPTS,AJAC,BJAC,\r
-                   H0VAL[JI-1],COLSC(LO),AQTOL,TOLOU(JI),XENPT,QINTS,MQIN1,IER);\r
-               IF (IER .GT. 0) THEN\r
-                 RETURN\r
-               ENDIF\r
-       C\r
-       C       SET UP THE COMPOSITE RULE ABSCISSAE AND WEIGHTS FOR THIS\r
-       C       JACOBI INDEX.\r
-       C\r
-               NQUAD(JI)=QINTS*NQPTS\r
-               LOQSB(JI)=HI+1\r
-               IF (HI+NQUAD(JI) .GT. MNQUA) THEN\r
-                   IER=9\r
-                   RETURN\r
-               ENDIF\r
-               SUM1=BETA+1E+0\r
-               K=HI\r
-               DO 40 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 20 J=1,NQPTS\r
-                     K=K+1\r
-                     QCOMX(K)=MEAN+RR*QUPTS(LO+J)\r
-                     QCOMW(K)=RRB*QUWTS(LO+J)\r
-       20          CONTINUE\r
-                 ELSE\r
-                   LO=NARCS*NQPTS\r
-                   DO 30 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
-       30          CONTINUE\r
-                 ENDIF\r
-       40      CONTINUE\r
-               HI=HI+NQUAD(JI)\r
+                   H0VAL[JI-1],CSCAL,AQTOL,TOLU,XENPT,QINTS,MQIN1,IER);\r
+               TOLOU[JI-1] = TOLU[0];\r
+               if (IER[0] > 0) {\r
+                   return;\r
+               }\r
+       \r
+                // SET UP THE COMPOSITE RULE ABSCISSAE AND WEIGHTS FOR THIS\r
+               // JACOBI INDEX.\r
+       \r
+               NQUAD[JI-1]=QINTS[0]*NQPTS;\r
+               LOQSB[JI-1]=HI+1;\r
+               if (HI+NQUAD[JI-1] > MNQUA) {\r
+                   IER[0]=9;\r
+                   return;\r
+               }\r
+               SUM1=BETA+1.0;\r
+               K=HI;\r
+               for (I=1; I <= QINTS[0]; I++) {\r
+                   RR=(XENPT[I]-XENPT[I-1])*0.5;\r
+                   MEAN=(XENPT[I]+XENPT[I-1])*0.5; \r
+                   if (I == 1) {\r
+                       RRB=Math.pow(RR,SUM1);\r
+                       LO=LO-1;\r
+                       for (J=1; J <= NQPTS; J++) {\r
+                           K=K+1;\r
+                           QCOMX[K-1]=MEAN+RR*QUPTS[LO+J-1];\r
+                           QCOMW[K-1]=RRB*QUWTS[LO+J-1];\r
+                       } // for (J=1; J <= NQPTS; J++)\r
+                   } // if (I == 1)\r
+                   else {\r
+                       LO=NARCS*NQPTS;\r
+                       for (J=1; J <= NQPTS; J++) {\r
+                           K=K+1;\r
+                           QCOMX[K-1]=MEAN+RR*QUPTS[LO+J-1];\r
+                           QCOMW[K-1]=RR*QUWTS[LO+J-1]*Math.pow((1.0+QCOMX[K-1]),BETA);\r
+                       } // for (J=1; J <= NQPTS; J++)\r
+                   }\r
+               } // for (I=1; I <= QINTS[0]; I++)\r
+               HI=HI+NQUAD[JI-1];\r
        } // for (JI=1; JI <= NARCS; JI++)\r
-       C\r
-       C     ASSIGN INITIAL DATA FOR LEGENDRE ARCS \r
-       C\r
-             I=2*NJIND\r
-             RR=R1MACH(2)\r
-             XIDST(I)=RR\r
-             XIDST(I-1)=RR\r
-             XIVAL(I)=CMPLX(RR)\r
-             XIVAL(I-1)=CMPLX(RR)\r
-             LOQSB(NJIND)=HI+1\r
-             NQUAD(NJIND)=0\r
-       C\r
-       C     FIND THE MAXIMUM OF THE TOLOU ELEMENTS\r
-       C\r
-             MCQER=0E+0\r
-             DO 60 I=1,NARCS\r
-               MCQER=MAX(MCQER,TOLOU(I))\r
-       60    CONTINUE \r
-       C\r
-       C     NORMAL TERMINATION\r
-       C\r
-             IER=0\r
-       C*/\r
+       \r
+       // ASSIGN INITIAL DATA FOR LEGENDRE ARCS \r
+       \r
+       I=2*NJIND;\r
+       RR=Double.MAX_VALUE;\r
+       XIDST[I-1]=RR;\r
+       XIDST[I-2]=RR;\r
+       XIVAL[I-1][0]=RR;\r
+       XIVAL[I-1][1] = 0.0;\r
+       XIVAL[I-2][0]=RR;\r
+       XIVAL[I-2][1] = 0.0;\r
+       LOQSB[NJIND-1]=HI+1;\r
+       NQUAD[NJIND-1]=0;\r
+       \r
+       // FIND THE MAXIMUM OF THE TOLOU ELEMENTS\r
+       \r
+       MCQER[0]=0.0;\r
+        for (I=1; I <= NARCS; I++) {\r
+           MCQER[0]=Math.max(MCQER[0],TOLOU[I-1]);\r
+        } \r
+       \r
+       // NORMAL TERMINATION\r
+       \r
+       IER[0]=0;\r
+       \r
     } // private void ICOQR1\r
 \r
 \r
     private void SUBIN7(double ZZ[][], int NZZ, double BETA, int MAXDG, int NQUAD, double AJAC[], double BJAC[],\r
-        double H0JAC, double CSCAL[], double TOLIN, double TOLOU, double XENPT[], int QINTS[], int MQIN1,int IER[]) {\r
+        double H0JAC, double CSCAL[], double TOLIN, double TOLOU[], double XENPT[], int QINTS[], int MQIN1,int IER[]) {\r
        //INTEGER MAXDG,NQUAD,QINTS,NZZ,IER,MQIN1\r
        //REAL BETA,AJAC(*),BJAC(*),H0JAC,CSCAL(*),TOLIN,TOLOU,XENPT(*)\r
        //COMPLEX ZZ(*)\r
@@ -4745,7 +4767,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        final int NMAX = 100;\r
        double TAU[] = new double[1];\r
        double MAXRM[] = new double[1];\r
-       double TOL,RIGHT;\r
+       double TOL;\r
+       double RIGHT[] = new double[1];\r
        boolean T1FXD;\r
        \r
        double REMND[][] = new double[NMAX][2];\r
@@ -4765,37 +4788,40 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        if (IER[0] > 0) {\r
            return;\r
        }\r
-       TOLOU=MAXRM[0];\r
+       TOLOU[0]=MAXRM[0];\r
        XENPT[1]=TAU[0];\r
        \r
        if (XENPT[1] < 1.0) {\r
            QINTS[0]=2;\r
            T1FXD=false;\r
            TAU[0]=1.0;\r
-           RIGHT=-1.0;\r
-           DELEG7(ZZ,NZZ,BETA,RIGHT,TAU[0],T1FXD,MAXDG,NQUAD,AJAC,\r
+           RIGHT[0]=-1.0;\r
+           DELEG7(ZZ,NZZ,BETA,RIGHT,TAU,T1FXD,MAXDG,NQUAD,AJAC,\r
                    BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER);\r
            if (IER[0] > 0) {\r
                return;\r
            }\r
-           TOLOU=TOLOU+MAXRM[0];\r
+           TOLOU[0]=TOLOU[0]+MAXRM[0];\r
            T1FXD=true;\r
        \r
            while (true) {\r
        \r
-               if (XENPT[QINTS[0]-1] > RIGHT) {\r
-                   XENPT[QINTS[0]-1]=0.5*(XENPT[QINTS[0]-1]+RIGHT);\r
+               if (XENPT[QINTS[0]-1] > RIGHT[0]) {\r
+                   XENPT[QINTS[0]-1]=0.5*(XENPT[QINTS[0]-1]+RIGHT[0]);\r
                    XENPT[QINTS[0]]=1.0;\r
                    break;\r
                }\r
                else {\r
                  TAU[0]=1.0;\r
-                 DELEG7(ZZ,NZZ,BETA,XENPT[QINTS[0]-1],TAU[0],T1FXD,MAXDG,\r
+                 double TAU1[] = new double[1];\r
+                 TAU1[0] = XENPT[QINTS[0]-1];\r
+                 DELEG7(ZZ,NZZ,BETA,TAU1,TAU,T1FXD,MAXDG,\r
                           NQUAD,AJAC,BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER);\r
+                 XENPT[QINTS[0]-1] = TAU1[0];\r
                  if (IER[0] > 0) {\r
                      return;\r
                  }\r
-                 TOLOU=TOLOU+MAXRM[0];\r
+                 TOLOU[0]=TOLOU[0]+MAXRM[0];\r
                  QINTS[0]=QINTS[0]+1;\r
                  if (QINTS[0] >= MQIN1) {\r
                      IER[0]=11;\r
@@ -5084,7 +5110,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         return result;\r
     } // private double LGGAM\r
 \r
-    private void DELEG7(double ZZ[][], int NZZ, double BETA, double TAU1, double TAU2, boolean T1FXD,\r
+    private void DELEG7(double ZZ[][], int NZZ, double BETA, double TAU1[], double TAU2[], boolean T1FXD,\r
         int MAXDG, int NQUAD, double ACOEF[], double BCOEF[], double H0VAL, double REMND[][], \r
         double CSCAL[], double TOL, double MAXRM[], int IER[]) {\r
         //INTEGER MAXDG,NQUAD,IER,NZZ\r
@@ -5267,15 +5293,15 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
         // NOW COME THE FACTORS DEPENDENT ON TAU1 AND TAU2.\r
        \r
-       LOWER=TAU1;\r
-       UPPER=TAU2;\r
+       LOWER=TAU1[0];\r
+       UPPER=TAU2[0];\r
        FIRST= true;\r
        \r
        while (true) {\r
        \r
            HTOL=0.5*TOL;\r
-           RR=(TAU2-TAU1)*0.5;\r
-           MEAN=(TAU1+TAU2)*0.5;\r
+           RR=(TAU2[0]-TAU1[0])*0.5;\r
+           MEAN=(TAU1[0]+TAU2[0])*0.5;\r
            BB=(1.0+MEAN)/RR;\r
            BB1=BB+Math.sqrt(BB*BB-1.0);\r
        \r
@@ -5351,14 +5377,14 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                      // TAU2<1) OR TAU1 NEED INCREASING OTHERWISE (BUT THIS IS ONLY\r
                       // POSSIBLE IF TAU1>-1)\r
        \r
-                     if (T1FXD && TAU2 < 1.0) {\r
-                         LOWER=TAU2;\r
-                         TAU2=0.5*(LOWER+UPPER);\r
+                     if (T1FXD && TAU2[0] < 1.0) {\r
+                         LOWER=TAU2[0];\r
+                         TAU2[0]=0.5*(LOWER+UPPER);\r
                          continue;\r
                      }\r
-                     else if (!T1FXD && TAU1 > -1.0) {\r
-                         UPPER=TAU1;\r
-                         TAU1=0.5*(LOWER+UPPER);\r
+                     else if (!T1FXD && TAU1[0] > -1.0) {\r
+                         UPPER=TAU1[0];\r
+                         TAU1[0]=0.5*(LOWER+UPPER);\r
                          continue;\r
                      }\r
                      else {\r
@@ -5379,12 +5405,12 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                      FIRST=false;\r
                  }\r
                  if (T1FXD) {\r
-                     UPPER=TAU2;\r
-                     TAU2=0.5*(LOWER+UPPER);\r
+                     UPPER=TAU2[0];\r
+                     TAU2[0]=0.5*(LOWER+UPPER);\r
                  }\r
                  else {\r
-                     LOWER=TAU1;\r
-                     TAU1=0.5*(LOWER+UPPER);\r
+                     LOWER=TAU1[0];\r
+                     TAU1[0]=0.5*(LOWER+UPPER);\r
                  }\r
                  continue;\r
              } // else MAXRM[0] >= TOL\r
@@ -5489,6 +5515,440 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        1000  FORMAT(A,1X,I5)\r
        C\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
+        double BCOEF[], double COLPR[], double QCOMX[], double QCOMW[], double CENTR[], double ZCOLL[][],\r
+        boolean INTER, boolean LNSEG[], double WORK[], int QIERR[], double MQERR[], double JACIN[], double A1COF[],\r
+        double B1COF[], double AQTOL, double RQTOL, double AQCOF[], double BQCOF[], double CQCOF[],\r
+        int LOOLD[], int HIOLD[]) {\r
+       \r
+        // INTEGER MNEQN,NCOLL,ORDSG,NQPTS,TNSUA,JATYP(*),PARNT(*),\r
+       //+DGPOL(*),LOSUB(*),HISUB(*),NQUAD(*),LOQSB(*),QIERR(0:6),LOOLD(*),\r
+       //+HIOLD(*)\r
+       \r
+       //REAL MATRX(MNEQN,MNEQN,2),RGTSD(*),TOLNR,MIDPT(*),HALEN(*),\r
+       //+H0VAL(*),COLSC(*),ACOEF(*),BCOEF(*),COLPR(*),QCOMX(*),QCOMW(*),\r
+       //+MQERR,WORK(*),A1COF(*),B1COF(*),AQTOL,\r
+       //+JACIN(*),RQTOL,AQCOF(*),BQCOF(*),CQCOF(*)\r
+       \r
+       //COMPLEX CENTR,ZCOLL(*)\r
+       \r
+       //LOGICAL INTER,LNSEG(*),REFLN\r
+       \r
+       // TO SET UP THE INITIAL COLLOCATION MATRIX (MATRX) AND THE RIGHT\r
+       // HAND SIDE VECTOR (RGTSD).  WE ASSUME COLLOCATION DOES NOT TAKE \r
+       // PLACE AT END POINTS OF ARCS.\r
+       \r
+       //.......................................................................\r
+       //     THIS ROUTINE IS ADAPTED FROM LNSY10 AND IS DESIGNED TO FORM ONLY\r
+       //     THOSE LINEAR EQUATIONS THAT ARISE FROM COLLOCATING AT POINTS ON\r
+       //     THE FUNDAMENTAL PART OF THE BOUNDARY\r
+       //.......................................................................\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+        int IER[] = new int[1];\r
+       int I,J,I1,J1,J2,J3,JO,IA,IQ,IS,JI,AJI,PT,DG,LOM,HIM,LOQ,LOD,\r
+           NQ,NCFBS,NROWS,NCOLS,ORDRG,ISAMAX,NVAL,IC,LOC,HIC,TSFBS;\r
+       double MD,HL,RH,TQ,WQ,DD,LDD,SG,SJI,TC,C0,\r
+              FNVAL,U0,U1,CURR,PREV,NEXT,ABER0,ABER1;\r
+       double GTQ[] = new double[2];\r
+       double ZQ[] = new double[2];\r
+       double ZC[] = new double[2];\r
+       // COMPLEX GTQ,ZQ,ZC,PARFUN,DPARFN\r
+       boolean OLDIA;\r
+       //COMMON /FNDEF/BETA,A1,B1,P0VAL,SCALE,TYPE\r
+       // EXTERNAL DPARFN,FNVAL,ISAMAX,JAPAR7,PARFUN,QAWS\r
+       qaws qmod;\r
+       int routine = Integration2.DQAWSE;\r
+       double DOUT[];\r
+       double DIN[] = new double[2];\r
+       /**\r
+         * Gives an upper bound on the number of subintervals in the partition of lower, upper.\r
+         */\r
+        int limit = 100;\r
+       \r
+       NCFBS=NCOLL/ORDSG;\r
+       TSFBS=TNSUA/ORDSG;\r
+       NROWS=NCFBS+1;\r
+       NCOLS=NCOLL+1;\r
+       \r
+       for (J=1; J <= NCOLS; J++) {\r
+           for (I=1; I <= NROWS; I++) {\r
+               MATRX[I-1][J-1][1]=0.0;\r
+           }\r
+       } // for (J=1; J <= NCOLS; J++)\r
+       for (I=0; I <= 6; I++) {\r
+           QIERR[I]=0;\r
+       }\r
+       MQERR[0]=0.0;\r
+       \r
+       // NOW SELECT THE INTEGRATION ARC\r
+       \r
+       for (IA=1; IA <= TNSUA; IA++) {\r
+\r
+       \r
+           // INITIALISE DATA FOR THIS ARC\r
+       \r
+           PT=PARNT[IA-1];\r
+           MD=MIDPT[IA-1];\r
+           HL=HALEN[IA-1];\r
+           DG=DGPOL[IA-1];\r
+           LOM=LOSUB[IA-1];\r
+           HIM=HISUB[IA-1];\r
+           JI=JATYP[IA-1];\r
+           if (JI < 0) {\r
+               SJI=-1.0;\r
+           }\r
+           else {\r
+               SJI=1.0;\r
+           }\r
+           AJI=Math.abs(JI);\r
+           LOD=(AJI-1)*NQPTS+1;\r
+           RH=Math.sqrt(H0VAL[AJI-1]);\r
+           BETA=JACIN[AJI-1];\r
+           A1=A1COF[AJI-1];\r
+           B1=B1COF[AJI-1];\r
+           P0VAL=1.0/RH;\r
+       \r
+           // ASSIGN THE NON-ZERO ELEMENT FOR THE SIDE EQUATION FOR THIS ARC\r
+       \r
+           MATRX[NROWS-1][LOM-1][1]=RH*COLSC[LOD-1];\r
+       \r
+           // SET UP THE DIAGONAL BLOCK OF SCALED PRINCIPAL SINGULAR INTEGRALS\r
+           // IF THIS DOESN'T ALREADY EXIST.\r
+       \r
+           JO=LOOLD[IA-1];\r
+           OLDIA=(JO > 0);\r
+           if (!OLDIA && HIM <= NCFBS) {\r
+               for (I1=LOM; I1 <= HIM; I1++) {\r
+                   TC=SJI*COLPR[I1-1];\r
+                   SCALE=COLSC[LOD-1];\r
+                   TYPE=1;\r
+                   // Use dqawse which is the same as (quadpack/dqaws) but provides more\r
+                   // information and control\r
+                   qmod = new qaws(-1.0,TC,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                   qmod.driver();\r
+                   U0 = qmod.getIntegral();\r
+                   ABER0 = qmod.getAbserr();\r
+                   NVAL = qmod.getNeval();\r
+                   IER[0] = qmod.getErrorStatus();\r
+                   QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                   TYPE=2;\r
+                   qmod = new qaws(TC,1.0,routine,0.0,0.0,2,AQTOL,RQTOL,limit);\r
+                   qmod.driver();\r
+                   U1 = qmod.getIntegral();\r
+                   ABER1 = qmod.getAbserr();\r
+                   NVAL = qmod.getNeval();\r
+                   IER[0] = qmod.getErrorStatus();\r
+                   QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                   WORK[0]=U0+U1;\r
+                   MQERR[0]=Math.max(MQERR[0],ABER0+ABER1);\r
+                   if (DG > 0) {\r
+                       SCALE=COLSC[LOD];\r
+                       TYPE=3;\r
+                       qmod = new qaws(-1.0,TC,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                       qmod.driver();\r
+                           U0 = qmod.getIntegral();\r
+                           ABER0 = qmod.getAbserr();\r
+                           NVAL = qmod.getNeval();\r
+                           IER[0] = qmod.getErrorStatus();\r
+                           QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                       TYPE=4;\r
+                       qmod = new qaws(TC,1.0,routine,0.0,0.0,2,AQTOL,RQTOL,limit);\r
+                       qmod.driver();\r
+                           U1 = qmod.getIntegral();\r
+                           ABER1 = qmod.getAbserr();\r
+                           NVAL = qmod.getNeval();\r
+                           IER[0] = qmod.getErrorStatus();\r
+                           QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                       WORK[1]=U0+U1;\r
+                       MQERR[0]=Math.max(MQERR[0],ABER0+ABER1);\r
+       \r
+                       // NOW USE THE (WEAKLY) STABLE FORWARD RECURRENCE SCHEME  \r
+                       // FOR WORK(I),I=3,NQPTS (WITH SCALE FACTOR FOR WORK(2))\r
+       \r
+                       CURR=WORK[1];\r
+                       PREV=SCALE;\r
+                       for (I=1; I <= DG-1; I++) {\r
+                           J=LOD+I-1;\r
+                           NEXT=(AQCOF[J-1]*TC-BQCOF[J-1])*CURR-CQCOF[J-1]*PREV;\r
+                           WORK[I+1]=NEXT;\r
+                           PREV=CURR;\r
+                           CURR=NEXT;\r
+                       } // for (I=1; I <= DG-1; I++)\r
+       \r
+                       // ASSIGN CORRECT SCALE FACTORS.\r
+       \r
+                       for (I=3; I <= DG+1; I++) {\r
+                           J=LOD+I-1;\r
+                           WORK[I-1]=WORK[I-1]*COLSC[J-1]/SCALE;\r
+                       } // for (I=3; I <= DG+1; I++)\r
+                   } // if (DG > 0)\r
+       \r
+                   SG=1.0;\r
+                   for (J=LOM; J <= HIM; J++) {\r
+                       MATRX[I1-1][J-1][1]=MATRX[I1-1][J-1][1]+SG*WORK[J-LOM];\r
+                       SG=SG*SJI;\r
+                   } // for (J=LOM; J <= HIM; J++)\r
+       \r
+               } // for (I1=LOM; I1 <= HIM; I1++)\r
+           } // if (!OLDIA && HIM <= NCFBS)\r
+\r
+       \r
+           // INITIALISE SOME DATA FOR THE NON-SINGULAR INTEGRALS\r
+       \r
+           WORK[0]=1.0/RH;\r
+           NQ=NQUAD[AJI-1];\r
+           LOQ=LOQSB[AJI-1];\r
+       \r
+           for (IQ=1; IQ <= NQ; IQ++) {\r
+               I=LOQ+IQ-1;\r
+               TQ=QCOMX[I-1];\r
+               WQ=QCOMW[I-1];\r
+               double AA[] = new double[DG];\r
+               double BB[] = new double[DG];\r
+               for (I = 0; I < DG; I++) {\r
+                       AA[I] = ACOEF[LOD-1+I];\r
+                       BB[I] = BCOEF[LOD-1+I];\r
+               }\r
+               JAPAR7(WORK,TQ,AA,BB,DG);\r
+               if (JI < 0) {\r
+                   TQ=-TQ;\r
+               }\r
+               GTQ[0] = MD+HL*TQ;\r
+               GTQ[1] = 0.0;\r
+               ZQ=PARFUN(PT,GTQ);\r
+       \r
+               // ACCUMULATE ALL NEW ELEMENTS NOT ON THE DIAGONAL BLOCK\r
+       \r
+               for (IC=1; IC <= TSFBS; IC++) {\r
+                   J2=LOOLD[IC-1];\r
+                   if (IC != IA && (J2 == 0 || !OLDIA)) {\r
+                       LOC=LOSUB[IC-1];\r
+                       HIC=HISUB[IC-1];\r
+                       for (I1=LOC; I1 <= HIC; I1++) {\r
+                           ZC[0]=ZCOLL[I1-1][0];\r
+                           ZC[1]=ZCOLL[I1-1][1];\r
+                           DD=zabs(ZC[0]-ZQ[0],ZC[1]-ZQ[1]);\r
+                           LDD=Math.log(DD)*WQ;\r
+                           SG=1.0;\r
+                           for (J1=LOM; J1 <= HIM; J1++) {\r
+                               J=J1-LOM+1;\r
+                               I=J1-LOM+LOD;\r
+                               MATRX[I1-1][J1-1][1]=MATRX[I1-1][J1-1][1]+SG*WORK[J-1]*LDD*COLSC[I-1];\r
+                               SG=SG*SJI;\r
+                           } // for (J1=LOM; J1 <= HIM; J1++)\r
+                       } // for (I1=LOC; I1 <= HIC; I1++)\r
+                   } // if (IC != IA && (J2 == 0 || !OLDIA))\r
+               } // for (IC=1; IC <= TSFBS; IC++)\r
+       \r
+               // ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
+               // ANY NEW DIAGONAL BLOCK FOR THE NON-LINE-SEGMENT CASE.\r
+       \r
+               if (!LNSEG[IA-1] && !OLDIA && HIM <= NCFBS) {\r
+                   for (I1=LOM; I1 <= HIM; I1++) {\r
+                       TC=COLPR[I1-1];\r
+                       DD=Math.abs(TC-TQ);\r
+                       if (DD <= TOLNR) {\r
+                               DOUT = DPARFN(PT,GTQ);\r
+                               DD = zabs(DOUT[0],DOUT[1])*HL;\r
+                       }\r
+                       else {\r
+                           ZC[0]=ZCOLL[I1-1][0];\r
+                           ZC[1]=ZCOLL[I1-1][1];\r
+                           DD=zabs(ZC[0]-ZQ[0],ZC[1]-ZQ[1])/DD;\r
+                           if (DD < TOLNR) {\r
+                               DOUT = DPARFN(PT,GTQ);\r
+                                       DD = zabs(DOUT[0],DOUT[1])*HL;\r
+                           }\r
+                       }\r
+                       LDD=Math.log(DD)*WQ;\r
+                       SG=1.0;\r
+                       for (J1=LOM; J1 <= HIM; J1++) {\r
+                           J=J1-LOM+1;\r
+                           I=J1-LOM+LOD;\r
+                           MATRX[I1-1][J1-1][1]=MATRX[I1-1][J1-1][1]+SG*WORK[J-1]*LDD*COLSC[I-1];\r
+                           SG=SG*SJI;\r
+                       } // for (J1=LOM; J1 <= HIM; J1++)\r
+                   } // for (I1=LOM; I1 <= HIM; I1++)\r
+               } // if (!LNSEG[IA-1] && !OLDIA && HIM <= NCFBS)\r
+           } // for (IQ=1; IQ <= NQ; IQ++)\r
+       \r
+           // ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
+           // ANY NEW DIAGONAL BLOCK FOR THE LINE-SEGMENT CASE.\r
+       \r
+           if (LNSEG[IA-1] && !OLDIA && HIM <= NCFBS) {\r
+                 DIN[0] = MD;\r
+                 DIN[1] = 0.0;\r
+                 DOUT = DPARFN(PT,DIN);\r
+                 ZC[0] = DOUT[0]*HL;\r
+                 ZC[1] = DOUT[1]*HL;\r
+                 C0=zabs(ZC[0],ZC[1]);\r
+                 C0=RH*Math.log(C0)*COLSC[LOD-1];\r
+                 for (I1=LOM; I1 <= HIM; I1++) {\r
+                     MATRX[I1-1][LOM-1][1]=MATRX[I1-1][LOM-1][1]+C0;\r
+                 }\r
+           } // if (LNSEG[IA-1] && !OLDIA && HIM <= NCFBS)\r
+       \r
+           // NOW COPY OVER ANY RELEVANT ELEMENTS ALREADY AVAILABLE IN\r
+           // MATRX(*,*,1)\r
+       \r
+           if (OLDIA) {\r
+               for (IC=1; IC <= TSFBS; IC++) {\r
+                   J2=LOOLD[IC-1];\r
+                   if (J2 > 0) {\r
+                       LOC=LOSUB[IC-1];\r
+                       HIC=HISUB[IC-1];\r
+                       J2=J2-1;\r
+                       for (I1=LOC; I1 <= HIC; I1++) {\r
+                           J2=J2+1;\r
+                           for (J1=LOM; J1 <= HIM; J1++) {\r
+                               J=J1+JO-LOM;\r
+                               MATRX[I1-1][J1-1][1]=MATRX[J2-1][J-1][0];\r
+                           } // for (J1=LOM; J1 <= HIM; J1++)\r
+                       } // for (I1=LOC; I1 <= HIC; I1++)\r
+                   } // if (J2 > 0)\r
+               } // for (IC=1; IC <= TSFBS; IC++)\r
+           } // if (OLDIA)\r
+       \r
+       } // for (IA=1; IA <= TNSUA; IA++)\r
+       \r
+       // SET UP THE LAST COLUMN\r
+       \r
+       for (I=1; I <= NCFBS; I++) {\r
+           MATRX[I-1][NCOLS-1][1]=1.0;\r
+       }\r
+       \r
+        // FINALLY SET UP THE RIGHT HAND SIDE VECTOR\r
+       \r
+       if (INTER) {\r
+           for (I=1; I <= NCFBS; I++) {\r
+               ZC[0]=ZCOLL[I-1][0]-CENTR[0];\r
+               ZC[1]=ZCOLL[I-1][1]-CENTR[1];\r
+               RGTSD[I-1]=Math.log(zabs(ZC[0],ZC[1]));\r
+           }\r
+       }\r
+       else {\r
+           for (I=1; I <= NCFBS; I++) {\r
+               RGTSD[I-1]=0.0;\r
+           }\r
+       }\r
+       RGTSD[NROWS-1]=1.0;\r
+       \r
+       // COPY MATRX(*,*,2) ONTO MATRX(*,*,1)\r
+       \r
+       for (J=1; J <= NCOLS; J++) {\r
+           for (I=1; I <= NROWS; I++) {\r
+               MATRX[I-1][J-1][0]=MATRX[I-1][J-1][1];\r
+           }\r
+       } // for (J=1; J <= NCOLS; J++)\r
+       \r
+       // COMBINE COLUMNS OF MATRX(*,*,2) TOGETHER ACCORDING TO SYMMETRY\r
+       // SPECIFICATIONS\r
+       \r
+       if (ORDSG > 1) {\r
+           if (REFLN) {\r
+               ORDRG=ORDSG/2;\r
+               for (IS=2; IS <= ORDRG; IS++) {\r
+                   LOM=(IS-1)*NCFBS*2;\r
+                   for (I=1; I <= NROWS; I++) {\r
+                       for (J=1; J <= NCFBS*2; J++) {\r
+                           J1=LOM+J;\r
+                           MATRX[I-1][J-1][1]=MATRX[I-1][J-1][1]+MATRX[I-1][J1-1][1];\r
+                       } // for (J=1; J <= NCFBS*2; J++)\r
+                   } // for (I=1; I <= NROWS; I++)\r
+               } // for (IS=2; IS <= ORDRG; IS++)\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
+                       for (I=1; I <= NROWS; I++) {\r
+                           MATRX[I-1][J-1][1]=MATRX[I-1][J-1][1]+SG*MATRX[I-1][J3-1][1];\r
+                       }\r
+                   } // for (J=J1; J <= J2; J++)\r
+               } // for (IA=1; IA <= TSFBS; IA++)\r
+           } // if (REFLN)\r
+           else {\r
+               for (IS=2; IS <= ORDSG; IS++) {\r
+                   LOM=(IS-1)*NCFBS;\r
+                   for (I=1; I <= NROWS; I++) {\r
+                       for (J=1; J <= NCFBS; J++) {\r
+                           J1=LOM+J;\r
+                           MATRX[I-1][J-1][1]=MATRX[I-1][J-1][1]+MATRX[I-1][J1-1][1];\r
+                       } // for (J=1; J <= NCFBS; J++)\r
+                   } // for (I=1; I <= NROWS; I++)\r
+               } // for (IS=2; IS <= ORDSG; IS++)\r
+           } // else\r
+       \r
+           for (I=1; I <= NROWS; I++) {\r
+               MATRX[I-1][NROWS-1][1]=MATRX[I-1][NCOLS-1][1];\r
+           }\r
+       } // if (ORDSG > 1)\r
+       \r
+    } // private void LNSY11\r
+    \r
+    class qaws extends Integration2 {\r
+       \r
+       public qaws(double lower, double upper, int routine, double alfa,\r
+                double beta, int integr, double epsabs, double epsrel, int limit) {\r
+               super(lower, upper, routine, alfa, beta, integr, epsabs, epsrel, limit);\r
+       }\r
+       \r
+       public void driver() {\r
+                       super.driver();\r
+               }\r
+       \r
+       public double intFunc(double X) {\r
+               double result;\r
+\r
+            if (TYPE == 1) {\r
+                result=P0VAL*SCALE;\r
+            }\r
+            else if (TYPE == 2) {\r
+                result=Math.pow((1.0+X),BETA)*P0VAL*SCALE;\r
+            }\r
+            else if (TYPE == 3) {\r
+                result=(X-B1)*P0VAL*SCALE/A1;\r
+            }\r
+            else {\r
+                result=Math.pow((1.0+X),BETA)*(X-B1)*P0VAL*SCALE/A1;\r
+            }\r
+            return result;     \r
+       }\r
+    }\r
+\r
+    private double FNVAL(double X) {\r
+\r
+        // LOCAL VARIABLES.\r
+\r
+        // INTEGER TYPE\r
+        // REAL BETA,A1,B1,P0VAL,SCALE\r
+        // COMMON /FNDEF/BETA,A1,B1,P0VAL,SCALE,TYPE\r
+       double result;\r
+\r
+        if (TYPE == 1) {\r
+            result=P0VAL*SCALE;\r
+        }\r
+        else if (TYPE == 2) {\r
+            result=Math.pow((1.0+X),BETA)*P0VAL*SCALE;\r
+        }\r
+        else if (TYPE == 3) {\r
+            result=(X-B1)*P0VAL*SCALE/A1;\r
+        }\r
+        else {\r
+            result=Math.pow((1.0+X),BETA)*(X-B1)*P0VAL*SCALE/A1;\r
+        }\r
+        return result;\r
+    } // private double FNVAL\r
+\r
+\r
 \r
  \r
       /**\r