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

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

index de0966f..f10d6f7 100644 (file)
@@ -2838,42 +2838,45 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                // **** SET UP THE ARRAY RIGLL OF REFERENCE IGNORE LEVELS.\r
                \r
                IGNLVL(RIGLL,COLSC,ACOEF,BCOEF,H0VAL,JACIN,NJIND,NQPTS,IER);\r
-               /*      IF (IER .GT. 0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-               C\r
-               C**** SET UP THE ARRAY OF COLLOCATION POINTS PARAMETER VALUES, COLPR,\r
-               C**** THE ARRAY OF COLLOCATION POINTS ZCOLL AND THE ARRAYS LOSUB AND \r
-               C**** HISUB NEEDED TO ACCESS COLPR AND ZCOLL CORRECTLY.  INITIALISE\r
-               C**** DGPOL AND UPDATE LNSEG FOR ARC HALVING. \r
-               C\r
-                     CALL CPJAC3(NARCS,NQPTS,INDEG,ISNPH(DGPOL),RSNPH(JACIN),\r
-                    +RSNPH(ACOEF),RSNPH(BCOEF),RWORK(DIAG),RWORK(SDIAG),TNSUA,\r
-                    +ISNPH(LOSUB),IWORK(HISUB),ISNPH(JATYP),IGEOM(PARNT),RGEOM(MIDPT),\r
-                    +RGEOM(HALEN),RWORK(COLPR),ZWORK(ZCOLL),LWORK(LNSEG),IWORK(LOOLD),\r
-                    +IWORK(HIOLD),EPS,IER,INIBT)\r
-                     IF (IER .GT. 0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-                     NCOLL=IWORK(HISUB+TNSUA-1)\r
-                     NEQNS=NCOLL+1\r
-                     NROWS=NCOLL/ORDSG+1\r
-                     IF (NEQNS .GT. MNEQN) THEN\r
-                       IER=8\r
-                       GOTO 999\r
-                     ENDIF\r
-                     WRITE(*,1) 'COLLOCATION POINT CHOICE DONE:'\r
-               C\r
-               C**** SET UP THE COMPOSITE GAUSSIAN QUADRATURE RULES, STORING ABSCISSAE\r
-               C**** AND WEIGHTS IN QCOMX AND QCOMW.  SET UP ARRAYS NQUAD,LOQSB\r
-               C**** NEEDED TO ACCESS THESE DATA.  RECORD MAXIMUM QUADRATURE ERRORS\r
-               C**** FOR COLUMN SCALED INTEGRALS IN ARRAY TOLOU.\r
-               C\r
-                     CALL ICOQR1(NARCS,NJIND,NQPTS,MDGPO,MQIN1,AQTOL,RSNPH(QUPTS),\r
-                    +RSNPH(QUWTS),RSNPH(JACIN),RGEOM(MIDPT),RGEOM(HALEN),RSNPH(ACOEF),\r
-                    +RSNPH(BCOEF),RSNPH(H0VAL),RWORK(COLSC),IWORK(NQUAD),IWORK(LOQSB),\r
-                    +RWORK(QCOMX),RWORK(QCOMW),MNQUA,RWORK(TOLOU),MCQER,RWORK(XENPT),\r
-                    +ZWORK(XIVAL),RWORK(XIDST),IER)\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(1,0,IER[0],null);\r
+                   return; \r
+               }\r
+               \r
+               // **** SET UP THE ARRAY OF COLLOCATION POINTS PARAMETER VALUES, COLPR,\r
+               // **** THE ARRAY OF COLLOCATION POINTS ZCOLL AND THE ARRAYS LOSUB AND \r
+               //**** 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
+                    ACOEF,BCOEF,DIAG,SDIAG,TNSUA,\r
+                    LOSUB,HISUB,JATYP,PARNT,MIDPT,\r
+                    HALEN,COLPR,ZCOLL,LNSEG,LOOLD,\r
+                    HIOLD,EPS,IER,INIBT);\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(1,0,IER[0],null);\r
+                   return; \r
+               }\r
+               NCOLL=HISUB[TNSUA-1];\r
+               NEQNS=NCOLL+1;\r
+               NROWS=NCOLL/ORDSG+1;\r
+               if (NEQNS > MNEQN) {\r
+                   IER[0]=8;\r
+                   WRTAIL(1,0,IER[0],null);\r
+                   return;     \r
+               }\r
+               System.out.println("COLLOCATION POINT CHOICE DONE:");\r
+               \r
+               //**** SET UP THE COMPOSITE GAUSSIAN QUADRATURE RULES, STORING ABSCISSAE\r
+               //**** AND WEIGHTS IN QCOMX AND QCOMW.  SET UP ARRAYS NQUAD,LOQSB\r
+               //**** NEEDED TO ACCESS THESE DATA.  RECORD MAXIMUM QUADRATURE ERRORS\r
+               //**** FOR COLUMN SCALED INTEGRALS IN ARRAY TOLOU.\r
+               \r
+               ICOQR1(NARCS,NJIND,NQPTS,MDGPO,MQIN1,AQTOL,QUPTS,\r
+                    QUWTS,JACIN,MIDPT,HALEN,ACOEF,\r
+                    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
@@ -3165,7 +3168,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
 \r
         // **** LOCAL VARIABLES\r
         int K,B0,B1,B2;\r
-        double X,Y,ANG,PI,R1MACH,EPSA,XI,APP;\r
+        double X,Y,ANG,PI,EPSA,XI,APP;\r
         double U[] = new double[2];\r
         double V[] = new double[2];\r
         double DIN[] = new double[2];\r
@@ -3746,7 +3749,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        // LOCAL VARIABLES\r
        \r
         int I,J,J1,JI,K,K1,KMAX,LO,LO1,M,N;\r
-       double BETA,ROOTH,P0SCL,X,MAXVL,SUM1,SUM2;\r
+       double ROOTH,P0SCL,X,MAXVL,SUM1,SUM2;\r
        // EXTERNAL JAPAR7\r
        \r
        K1=0;\r
@@ -3756,7 +3759,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        }\r
        \r
        for (JI=1; JI <= NJIND; JI++) {\r
-           BETA=JACIN[JI-1];\r
+           //BETA=JACIN[JI-1];\r
            ROOTH=Math.sqrt(H0VAL[JI-1]);\r
            P0SCL=1.0/ROOTH;\r
            LO=(JI-1)*NQPTS;\r
@@ -3880,7 +3883,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
        final int MNQPT = 20;\r
        int JT,K,K1,LO,N;\r
-       double B1,B2,BETA,EXPON,H1VAL,LGTWO,OFLOW,PI,RH,RH1,SUP,T;\r
+       double H1VAL[] = new double[1];\r
+       double B1,B2,BETA,EXPON,LGTWO,OFLOW,PI,RH,RH1,SUP,T;\r
        boolean SYMM;\r
        double A1COF[] = new double[MNQPT];\r
        double B1COF[] = new double[MNQPT];\r
@@ -3905,7 +3909,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        PI= Math.PI;\r
        LGTWO=Math.log(2.0);\r
        OFLOW=Math.log(Double.MAX_VALUE);\r
-       /*for (JT=1; JT <= NJIND; JT++) {\r
+       for (JT=1; JT <= NJIND; JT++) {\r
            BETA=JACIN[JT-1];\r
            B1=BETA+1.0;\r
            B2=BETA+2.0;\r
@@ -3936,7 +3940,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                continue;\r
            }\r
        \r
-           RH1=Math.sqrt(H1VAL);\r
+           RH1=Math.sqrt(H1VAL[0]);\r
            K1=K1+1;\r
        \r
            // ****   COMPUTE THE QUANTITY\r
@@ -3973,32 +3977,33 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                } // for (K=1; K <= N; K++)\r
        \r
                IMTQLH(N,DIAG,SDIAG,IER);\r
-                 IF (IER.GT.0) THEN\r
-                   IER=6\r
-                   RETURN\r
-                 ENDIF\r
-       C\r
-                 SUP=0E+0\r
-                 DO 20 K=1,N\r
-                   T=DIAG(K)\r
-                   PP(1)=1E+0/RH1\r
-                   CALL JAPAR7(PP,T,A1COF,B1COF,N)\r
-                   T=(1E+0-T)*(1E+0+T)**B1*PP(N)\r
-                   SUP=MAX(SUP,ABS(T))\r
-       20        CONTINUE\r
-       C\r
-                 K1=K1+1\r
-                 IF (SYMM) THEN\r
-                   RIGLL(K1)=2E+0*PI*SUP*COLSC(K1)/SQRT(N*(N+B1))\r
-                 ELSE\r
-                   RIGLL(K1)=2E+0*PI*SUP/SQRT(N*(N+B1))\r
-                 ENDIF\r
+               if (IER[0] > 0) {\r
+                   IER[0]=6;\r
+                   return;\r
+               }\r
+       \r
+               SUP=0.0;\r
+               for (K=1; K <= N; K++) {\r
+                   T=DIAG[K-1];\r
+                   PP[0]=1.0/RH1;\r
+                   JAPAR7(PP,T,A1COF,B1COF,N);\r
+                   T=(1.0-T)*Math.pow((1.0+T),B1)*PP[N-1];\r
+                   SUP=Math.max(SUP,Math.abs(T));\r
+               } // for (K=1; K <= N; K++)\r
+       \r
+               K1=K1+1;\r
+               if (SYMM) {\r
+                   RIGLL[K1-1]=2.0*PI*SUP*COLSC[K1-1]/Math.sqrt(N*(N+B1));\r
+               }\r
+               else {\r
+                   RIGLL[K1-1]=2.0*PI*SUP/Math.sqrt(N*(N+B1));\r
+               }\r
            } // for (N=2; N <= NQPTS-1; N++)\r
        } // for (JT=1; JT <= NJIND; JT++)\r
        \r
        // NORMAL EXIT\r
        \r
-       IER[0]=0;*/\r
+       IER[0]=0;\r
        \r
     } // private void IGNLVL\r
     \r
@@ -4006,7 +4011,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     \r
         int I,J,L,M,II,MML;\r
         //REAL D(N),E(N)\r
-        double B,C,F,G,P,R,S,TST1,TST2,PYTHAG;\r
+        double B,C,F,G,P,R,S,TST1,TST2;\r
     \r
         // THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1,\r
         // NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,\r
@@ -4060,74 +4065,964 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         }\r
         E[N-1] = 0.0;\r
     \r
-        /*for (L = 1; L <= N; L++) {\r
+        for (L = 1; L <= N; L++) {\r
             J = 0;\r
             // .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........\r
-            while (true) {\r
+            loopW: while (true) {\r
                 for (M = L; M <= N; M++) {\r
-                IF (M .EQ. N) GO TO 120\r
-                TST1 = ABS(D(M)) + ABS(D(M+1))\r
-                TST2 = TST1 + ABS(E(M))\r
-                IF (TST2 .EQ. TST1) GO TO 120\r
+                    if (M == N) break;\r
+                    TST1 = Math.abs(D[M-1]) + Math.abs(D[M]);\r
+                    TST2 = TST1 + Math.abs(E[M-1]);\r
+                    if (TST2 == TST1) break;\r
                 } // for (M = L; M <= N; M++)\r
-    C\r
-      120    P = D(L)\r
-             IF (M .EQ. L) GO TO 215\r
-             IF (J .EQ. 30) GO TO 1000\r
-             J = J + 1\r
-    C     .......... FORM SHIFT ..........\r
-             G = (D(L+1) - P) / (2.0E0 * E(L))\r
-             R = PYTHAG(G,1.0E0)\r
-             G = D(M) - P + E(L) / (G + SIGN(R,G))\r
-             S = 1.0E0\r
-             C = 1.0E0\r
-             P = 0.0E0\r
-             MML = M - L\r
-    C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........\r
-             DO 200 II = 1, MML\r
-                I = M - II\r
-                F = S * E(I)\r
-                B = C * E(I)\r
-                R = PYTHAG(F,G)\r
-                E(I+1) = R\r
-                IF (R .EQ. 0.0E0) GO TO 210\r
-                S = F / R\r
-                C = G / R\r
-                G = D(I+1) - P\r
-                R = (D(I) - G) * S + 2.0E0 * C * B\r
-                P = S * R\r
-                D(I+1) = G + P\r
-                G = C * R - B\r
-      200    CONTINUE\r
-    C\r
-             D(L) = D(L) - P\r
-             E(L) = G\r
-             E(M) = 0.0E0\r
-             GO TO 105\r
-    C     .......... RECOVER FROM UNDERFLOW ..........\r
-      210    D(I+1) = D(I+1) - P\r
-             E(M) = 0.0E0\r
-             GO TO 105\r
-    C     .......... ORDER EIGENVALUES ..........\r
-      215    IF (L .EQ. 1) GO TO 250\r
-    C     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........\r
-             DO 230 II = 2, L\r
-                I = L + 2 - II\r
-                IF (P .GE. D(I-1)) GO TO 270\r
-                D(I) = D(I-1)\r
-      230    CONTINUE\r
-    C\r
-      250    I = 1\r
-      270    D(I) = P\r
-            } // while (true)\r
+    \r
+             P = D[L-1];\r
+             if (M != L) {\r
+                 if (J == 30) {\r
+                     // .......... SET ERROR -- NO CONVERGENCE TO AN\r
+                        //                EIGENVALUE AFTER 30 ITERATIONS ..........\r
+                        IERR[0] = L;\r
+                        return;\r
+                 } // if (J == 30)\r
+                 J = J + 1;\r
+                 // .......... FORM SHIFT ..........\r
+                 G = (D[L] - P) / (2.0 * E[L-1]);\r
+                 R = PYTHAG(G,1.0);\r
+                 double sigR;\r
+                 if (G >= 0) {\r
+                        sigR = Math.abs(R);\r
+                 }\r
+                 else {\r
+                        sigR = -Math.abs(R);\r
+                 }\r
+                 G = D[M-1] - P + E[L-1] / (G + sigR);\r
+                 S = 1.0;\r
+                 C = 1.0;\r
+                 P = 0.0;\r
+                 MML = M - L;\r
+                 sect: {\r
+                     // .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........\r
+                     for (II = 1; II <= MML; II++) {\r
+                         I = M - II;\r
+                         F = S * E[I-1];\r
+                         B = C * E[I-1];\r
+                         R = PYTHAG(F,G);\r
+                         E[I] = R;\r
+                         if (R == 0.0) break sect;\r
+                         S = F / R;\r
+                         C = G / R;\r
+                         G = D[I] - P;\r
+                         R = (D[I-1] - G) * S + 2.0 * C * B;\r
+                         P = S * R;\r
+                         D[I] = G + P;\r
+                         G = C * R - B;\r
+                     } // for (II = 1; II <= MML; II++)\r
+    \r
+                     D[L-1] = D[L-1] - P;\r
+                     E[L-1] = G;\r
+                     E[M-1] = 0.0;\r
+                     continue loopW;\r
+                 } // sect\r
+                 //.......... RECOVER FROM UNDERFLOW ..........\r
+                 D[I] = D[I] - P;\r
+                 E[M-1] = 0.0;\r
+                 continue loopW;\r
+             } // if (M != L)\r
+             else { // if (M == L)\r
+             // .......... ORDER EIGENVALUES ..........\r
+             secb: {\r
+             if (L != 1) {\r
+                 // .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........\r
+                 for (II = 2; II <= L; II++) {\r
+                     I = L + 2 - II;\r
+                     if (P >= D[I-2]) break secb;\r
+                     D[I-1] = D[I-2];\r
+                 } // for (II = 2; II <= L; II++)\r
+             } // if (L != 1)\r
+    \r
+             I = 1;\r
+             } // secb:\r
+             D[I-1] = P;\r
+             break loopW;\r
+             } // else if (M == L)\r
+            } // loopW: while (true)\r
         } // for (L = 1; L <= N; L++)\r
-    C\r
-          GO TO 1001\r
-    C     .......... SET ERROR -- NO CONVERGENCE TO AN\r
-    C                EIGENVALUE AFTER 30 ITERATIONS ..........\r
-     1000 IERR = L\r
-     1001 RETURN*/\r
+        return;\r
+   \r
+    } // private void IMTQLH\r
+    \r
+    private double PYTHAG(double A, double B) {\r
+\r
+        // FINDS SQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW\r
+\r
+        double P,R,S,T,U;\r
+        P = Math.max(Math.abs(A),Math.abs(B));\r
+        if (P == 0.0E0) {\r
+            return P;  \r
+        }\r
+        R = (Math.min(Math.abs(A),Math.abs(B))/P);\r
+        R = R*R;\r
+        while (true) {\r
+           T = 4.0 + R;\r
+           if (T == 4.0) {\r
+                  return P;\r
+           }\r
+           S = R/T;\r
+           U = 1.0 + 2.0*S;\r
+           P = U*P;\r
+           double div = (S/U);\r
+           R = div*div * R;\r
+        } // while (true)\r
+    } // private double PYTHAG\r
+\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 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
+        // INTEGER DGPOL(*),LOSUB(*),HISUB(*),JATYP(*),PARNT(*)\r
+        // REAL EPS,JACIN(*),ACOEF(*),BCOEF(*),DIAG(*),SDIAG(*),\r
+        // +MIDPT(*),HALEN(*),COLPR(*),LOOLD(*),HIOLD(*)\r
+       // JAPHYC call has INT LOOLD(*) and HIOLD(*).\r
+        // COMPLEX ZCOLL(*)\r
+        // LOGICAL LNSEG(*),INIBT\r
+\r
+        // **** TO MAKE THE INITIAL ASSIGNMENT OF THE COLLOCATION PARAMETERS \r
+        // **** (STORED IN COLPR), THE COLLOCATION POINTS ON THE PHYSICAL \r
+        // **** BOUNDARY (STORED IN ZCOLL) AND THE ARRAYS LOSUB AND HISUB\r
+        // **** NEEDED TO ACCESS THIS DATA CORRECTLY.  ALSO TO SET UP THE\r
+        // **** ARRAYS\r
+        // ****        JATYP - THE JACOBI INDEX TYPE OF EACH SUBARC\r
+        // ****        PARNT - THE ORIGINAL 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 INITIAL POLYNOMIAL DEGREE ON EACH SUBARC\r
+        // ****   LNSEG - THE INITIAL LINE SEGMENT TYPE OF EACH SUBARC.\r
+\r
+        // **** IER=0 - NORMAL EXIT\r
+        // **** IER=7 - FAILURE TO CONVERGE IN EIGENVALUE ROUTINE IMTQLH\r
+\r
+        // LOCAL VARIABLES\r
+\r
+        int IFAIL[] = new int[1];\r
+       int D,D1,FIRST,I,J,K,K1,K2,P,PREV;\r
+        double S,TC,ALFA,BETA,MED,MDNBT;\r
+        double PIN[] = new double[2];\r
+        // COMPLEX PARFUN\r
+        // EXTERNAL PARFUN,IMTQLH,MDNBT\r
\r
+        TNSUA=2*NARCS;\r
+        for (I=1; I <= NARCS; I++) {\r
+            BETA=JACIN[I-1];\r
+            if (I == NARCS) {\r
+                ALFA=JACIN[0];\r
+            }\r
+            else {\r
+                ALFA=JACIN[I];\r
+            }\r
+            if (INIBT) {\r
+                MED=MDNBT(ALFA,BETA);\r
+            }\r
+            else {\r
+                MED=0.0;\r
+            }\r
+            J=2*I-1;\r
+            MIDPT[J-1]=0.5*(MED-1.0);\r
+            HALEN[J-1]=0.5*(MED+1.0);\r
+            PARNT[J-1]=I;\r
+            JATYP[J-1]=I;\r
+            J=J+1;\r
+            MIDPT[J-1]=0.5*(MED+1.0);\r
+            HALEN[J-1]=0.5*(1.0-MED);\r
+            PARNT[J-1]=I;\r
+            if (I == NARCS) {\r
+                JATYP[J-1]=-1;\r
+            }\r
+            else {\r
+                JATYP[J-1]=-I-1;\r
+            }\r
+        } // for (I=1; I <= NARCS; I++)\r
+\r
+        for (I=NARCS; I >= 1; I--) {\r
+            J=2*I;\r
+            LNSEG[J-1]=LNSEG[I-1];\r
+            LNSEG[J-2]=LNSEG[I-1];\r
+        } // for (I=NARCS; I >= 1; I--)\r
+\r
+        for (I=1; I <= TNSUA; 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
+            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
+            LOOLD[I-1]=0;\r
+            HIOLD[I-1]=-1;\r
+        } // for (I=1; 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
+                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
+            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
+                ZCOLL[K2-1]=PARFUN(P,PIN);\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 CPJAC3\r
+\r
+    private double MDNBT(double ALFA, double BETA) {\r
+\r
+        // MDNBT IS THE MEDIAN OF THE BETA DISTRIBUTION DEFINED BY THE\r
+        // DENSITY (1-X)**ALFA * (1+X)**BETA ON (-1,1).\r
+\r
+        // LOCAL VARIABLES..\r
+\r
+        int N;\r
+        double TOL,CC,CONST,SOLD,SNEW,SMNAB,GAMMA,FF,AA,BB,SS;\r
+        // EXTERNAL SMNAB,GAMMA\r
+\r
+        TOL=100.0*EPS;\r
+\r
+        if (Math.abs(ALFA-BETA) <= TOL) {\r
+            return 0.0;\r
+        }\r
+\r
+        if (ALFA > BETA) {\r
+            CC=BETA+1.0;\r
+        }\r
+        else {\r
+            CC=ALFA+1.0;\r
+        }\r
+\r
+        FF=1.0;\r
+        AA=ALFA+1.0;\r
+        BB=BETA+1.0;\r
+        SS=ALFA+BETA+2.0;\r
+\r
+    while (AA > 30.0) {\r
+        FF=(AA-1.0)*FF/(SS-1.0);\r
+        AA=AA-1.0;\r
+        SS=SS-1.0;\r
+    }\r
+\r
+    while (BB > 30.0) {\r
+        FF=(BB-1.0)*FF/(SS-1.0);\r
+        BB=BB-1.0;\r
+        SS=SS-1.0;\r
+    }\r
+\r
+    CONST=0.5*FF*GAMMA(AA)/GAMMA(SS);\r
+    CONST=CONST*GAMMA(BB);\r
+\r
+    N=0;\r
+    SOLD=Math.pow((CC*CONST),(1.0/CC));\r
+\r
+    while (true) {\r
+        N=N+1;\r
+        if (ALFA > BETA) {\r
+            SNEW=Math.pow((CONST/SMNAB(N,ALFA,BETA,SOLD)),(1.0/CC));\r
+        }\r
+        else {\r
+            SNEW=Math.pow((CONST/SMNAB(N,BETA,ALFA,SOLD)),(1.0/CC));\r
+        }\r
+\r
+        if (Math.abs(1.0-SOLD/SNEW) >= TOL) {\r
+            SOLD=SNEW;\r
+            continue;\r
+        }\r
+        else {\r
+               break;\r
+        }\r
+    } // while (true)\r
+\r
+    double result;\r
+    if (ALFA > BETA) {\r
+        result=2.0*SNEW-1.0;\r
+    }\r
+    else {\r
+        result=1.0-2.0*SNEW;\r
     }\r
+    return result;\r
+    } // private double MDNBT\r
+    \r
+    private double SMNAB(int N, double AL, double BE, double X) {\r
+\r
+        // EVALUATES A SUM NEEDED TO DETERMINE THE MEDIAN OF THE BETA\r
+        // DISTRIBUTION.\r
+\r
+        // LOCAL VARIABLES..\r
+\r
+        int I;\r
+        double SUM,TERM;\r
+\r
+        TERM=1.0/(1.0+BE);\r
+        SUM=TERM;\r
+        for (I=1; I <= N; I++) {\r
+            TERM=-X*TERM*(AL-I+1)*(I+BE)/I/(I+BE+1);\r
+            SUM=SUM+TERM;\r
+        } // for (I=1; I <= N; I++)\r
+\r
+        return SUM;\r
+\r
+    } // private double SMNAB\r
+    \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
+        double XIVAL[][], double XIDST[], int IER[]) {\r
+        //INTEGER NARCS,NJIND,NQPTS,MDGPO,MQIN1,IER,NQUAD(*),LOQSB(*),\r
+       //+MNQUA\r
+       //REAL AQTOL,QUPTS(*),QUWTS(*),JACIN(*),MIDPT(*),HALEN(*),ACOEF(*),\r
+       //+BCOEF(*),H0VAL(*),COLSC(*),QCOMX(*),QCOMW(*),TOLOU(*),XENPT(*),\r
+       //+XIDST(*),MCQER\r
+       //COMPLEX XIVAL(*)\r
+       \r
+       // THE MAIN PURPOSE OF THIS ROUTINE IS TO SET UP 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
+       //     MCQRE - THE INFINITY NORM OF TOLOU.\r
+       //     IER   - IER=0 FOR NORMAL TERMINATION.\r
+       //             IER=9 THE REQUIRED TOTAL NUMBER OF COMPOSITE QUADRATURE\r
+       //                   POINTS EXCEEDS THE LIMIT MNQUA.\r
+       \r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       int QINTS[] = new int[1];\r
+       int I,J,K,J0,J1,J2,J3,JI,JI0,JI1,HI,LO;\r
+       double DST[] = new double[4];\r
+       double BETA,H1,H2,T0,T1,T2,T3,SUM1,RR,RRB,MEAN,RXI,IXI;\r
+       final double ONE[] = new double[]{1.0,0.0};\r
+       double Z0[] = new double[2];\r
+       double Z1[] = new double[2];\r
+       double Z2[] = new double[2];\r
+       double Z3[] = new double[2];\r
+       double XI[][] = new double[4][2];\r
+       double GT[] = new double[2];\r
+       //COMPLEX ONE,Z0,Z1,Z2,Z3,XI(4),PARFUN,DPARFN,GT\r
+       //PARAMETER (ONE=(1E+0,0E+0))\r
+       //EXTERNAL DPARFN,PARFUN,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
+       HI=0;\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
+           // ONE BEYOND THAT IS THE JI1'TH.  THE FOUR SUBARCS ON THE JI0'TH\r
+           // AND JI'TH ANALYTIC ARCS ARE THE J0'TH, J1'TH, J2'TH AND J3'TH,\r
+           // STARTING AT THE BEGINING OF ARC JI0 AND ENDING AT THE END OF ARC\r
+           // JI.\r
+       \r
+           BETA=JACIN[JI-1];\r
+           LO=(JI-1)*NQPTS+1;\r
+           J2=2*JI-1;\r
+           J3=J2+1;\r
+           if (JI == 1) {\r
+               JI0=NARCS;\r
+               J1=2*NARCS;\r
+           }\r
+           else {\r
+               JI0=JI-1;\r
+               J1=J2-1;\r
+           }\r
+           J0=J1-1;\r
+           if (JI == NARCS) {\r
+               JI1=1;\r
+           }\r
+           else {\r
+               JI1=JI+1;\r
+           }\r
+       \r
+           // NEXT WE FIX THE LOCAL PARAMETER VALUES OF THE COLLOCATION\r
+           // POINTS NEAREST TO SUBARCS J1 AND J2.  FOR SUBARC J1 THE NEAR \r
+           // POINTS HAVE LOCAL PARAMETER VALUES T0 ON J0 AND T2 ON J2.  FOR\r
+           // SUBARC J2 THE NEAR POINTS HAVE LOCAL PARAMETER VALUES T1 ON J1\r
+           // AND T3 ON J3.\r
+       \r
+           T0=QUPTS[JI0*NQPTS-1];\r
+           T2=QUPTS[LO-1];\r
+           T1=-T2;\r
+           T3=-QUPTS[JI1*NQPTS-1];\r
+       \r
+           // NOW CONVERT THESE LOCAL PARAMETER VALUES FOR THE SUBARCS TO\r
+           // GLOBAL PARMETER VALUES FOR THE MAIN ARCS JI0 AND JI.\r
+       \r
+           T0=MIDPT[J0-1]+HALEN[J0-1]*T0;\r
+           H1=HALEN[J1-1];\r
+           T1=MIDPT[J1-1]+H1*T1;\r
+           H2=HALEN[J2-1];\r
+           T2=MIDPT[J2-1]+H2*T2;\r
+           T3=MIDPT[J3-1]+HALEN[J3-1]*T3;\r
+       \r
+           // NOW COMPUTE THE POSITIONS OF THE FOUR NEAR POINTS ON THE\r
+           // PHYSICAL BOUNDARY.\r
+           PIN[0] = T0;\r
+           PIN[1] = 0;\r
+           Z0=PARFUN(JI0,PIN);\r
+           PIN[0] = T1;\r
+           PIN[1] = 0;\r
+           Z1=PARFUN(JI0,PIN);\r
+           PIN[0] = T2;\r
+           PIN[1] = 0;\r
+           Z2=PARFUN(JI,PIN);\r
+           PIN[0] = T3;\r
+           PIN[1] = 0;\r
+           Z3=PARFUN(JI,PIN);\r
+       \r
+            // FIND THE APPROXIMATE PARAMETRIC PREIMAGE OF Z0 WRT SUBARC J1.\r
+           GT[0] = MIDPT[J1-1]-H1;\r
+           GT[1] = 0.0;\r
+           POUT = PARFUN(JI0,GT);\r
+           DOUT = DPARFN(JI0,GT);\r
+           zdiv(Z0[0] - POUT[0],Z0[1] - POUT[1],DOUT[0],DOUT[1],cr,ci);\r
+           XI[0][0] = -1.0 + cr[0]/H1;\r
+           XI[0][1] = ci[0]/H1;\r
+       \r
+           // CONVERT TO  PARAMETRIC PREIMAGE WRT SUBARC J2.\r
+       \r
+           XI[0][0]=-XI[0][0]; \r
+           XI[0][1]=-XI[0][1];\r
+       \r
+           // FIND THE APPROXIMATE PARAMETRIC PREIMAGE OF Z1 WRT SUBARC J2.\r
+       \r
+           PIN[0] = -1.0;\r
+           PIN[1] = 0.0;\r
+           POUT = PARFUN(JI,PIN);\r
+           DOUT = DPARFN(JI,PIN);\r
+           zdiv(Z1[0] - POUT[0],Z1[1] - POUT[1], DOUT[0], DOUT[1], cr, ci);\r
+           XI[1][0] = -1.0 + cr[0]/H2;\r
+           XI[1][1] = ci[0]/H2;\r
+       \r
+           // FIND THE APPROXIMATE PARAMETRIC PREIMAGE OF Z2 WRT SUBARC J1.\r
+           POUT = PARFUN(JI0,ONE);\r
+           DOUT = DPARFN(JI0,ONE);\r
+           zdiv(Z2[0] - POUT[0],Z2[1] - POUT[1],DOUT[0],DOUT[1],cr,ci);\r
+           XI[2][0] = 1.0 + cr[0]/H1;\r
+           XI[2][1]= ci[0]/H1; \r
+       \r
+           // CONVERT TO  PARAMETRIC PREIMAGE WRT SUBARC J2.\r
+       \r
+           XI[2][0]=-XI[2][0]; \r
+           XI[2][1]=-XI[2][1]; \r
+       \r
+           // FIND THE APPROXIMATE PARAMETRIC PREIMAGE OF Z3 WRT SUBARC J2.\r
+       \r
+               GT[0]=MIDPT[J2-1]+H2;\r
+               GT[1] = 0.0;\r
+               POUT = PARFUN(JI,GT);\r
+               DOUT = DPARFN(JI,GT);\r
+               zdiv(Z3[0]-POUT[0],Z3[1]-POUT[1],DOUT[0],DOUT[1],cr,ci);\r
+               XI[3][0] = 1.0 + cr[0]/H2;\r
+               XI[3][1] = ci[0]/H2;\r
+       \r
+               // SELECT THE PREIMAGE NEAREST -1 AND THE ONE NEAREST +1.\r
+       \r
+               for (J=1; J <= 4; 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 <= 4; J++)\r
+       \r
+               if (DST[1] < DST[2]) {\r
+                   XIVAL[J2-1][0]=XI[1][0];\r
+                   XIVAL[J2-1][1]=XI[1][1];\r
+                   XIDST[J2-1]=DST[1];\r
+               }\r
+               else {\r
+                   XIVAL[J2-1][0]=XI[2][0];\r
+                   XIVAL[J2-1][1]=XI[2][1];\r
+                   XIDST[J2-1]=DST[2];\r
+               }\r
+       \r
+               if (DST[0] < DST[3]) {\r
+                   XIVAL[J3-1][0]=XI[0][0];\r
+                   XIVAL[J3-1][1]=XI[0][1];\r
+                   XIDST[J3-1]=DST[0];\r
+               }\r
+               else {\r
+                   XIVAL[J3-1][0]=XI[3][0];\r
+                   XIVAL[J3-1][1]=XI[3][1];\r
+                   XIDST[J3-1]=DST[3];\r
+               }\r
+       \r
+               // NOW DETERMINE THE NUMBER AND LOCATION OF THE  QUADRATURE \r
+               // SUBINTERVALS NEEDED TO MEET THE TOLERANCE AT XIVAL(J2) AND \r
+                // XIVAL(J3).\r
+       \r
+               double ZZ[][] = new double[2][2];\r
+               for (I = 0; I < 2; I++) {\r
+                       ZZ[I][0] = XIVAL[J2-1+I][0];\r
+                       ZZ[I][1] = XIVAL[J2-1+I][1];\r
+               }\r
+               SUBIN7(ZZ,2,BETA,MDGPO,NQPTS,ACOEF(LO),BCOEF(LO),\r
+                   H0VAL(JI),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
+       } // 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
+    } // 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
+       //INTEGER MAXDG,NQUAD,QINTS,NZZ,IER,MQIN1\r
+       //REAL BETA,AJAC(*),BJAC(*),H0JAC,CSCAL(*),TOLIN,TOLOU,XENPT(*)\r
+       //COMPLEX ZZ(*)\r
+       \r
+       //     CALCULATES THE NUMBER OF QUADRATURE INTERVALS (QINTS) REQUIRED\r
+       //     FOR THE COMPOSITE GAUSS-JACOBI/GAUSS-LEGENDRE ESTIMATION OF\r
+       \r
+       //       INTEGRAL  [(1+X)**BETA*P(X,I)*LOG|ZZ(J)-X|*dX], I=0,1,...,MAXDG\r
+       //      -1<=X<=1                                         J=1,NZZ\r
+       \r
+       //     WHERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I\r
+       //     ASSOCIATED WITH THE WEIGHT (1+X)**BETA AND ZZ(J),J=1,..,NZZ, ARE \r
+       //     GIVEN POINTS CLOSE TO [-1,1].  \r
+       \r
+       //     THE ENDPOINTS OF THESE INTERVALS ARE RETURNED IN VECTOR XENPT, \r
+       //     WITH XENPT(1)=-1<XENPT(2)<...<1=XENPT(QINTS+1).\r
+       \r
+       //     IF Q(I,J) DENOTES THE ABSOLUTE QUADRATURE ERROR FOR THE INTEGRAL\r
+       //     ASSOCIATED WITH P(.,I) AND ZZ(J) THEN WE REQUIRE THAT\r
+       \r
+       //                 MAX       Q(I,J)*CSCAL(I) < TOLIN,\r
+       //              I=0,MAXDG\r
+       //               J=1,NZZ\r
+       \r
+       //     WITH THE MAXIMUM ON THE LEFT BEING REASONABLY CLOSE TO TOLIN.\r
+       //     TOLOU RETURNS THE COMPUTED VALUE OF THE ABOVE MAXIMUM.\r
+        \r
+       //     IER=0 - NORMAL EXIT\r
+       //     IER=10- PARAMETER NMAX LOCAL TO THIS ROUTINE NEEDS INCREASING TO\r
+       //             BE AT LEAST NZZ*(MAXDG+1)\r
+       //     IER=11- REQUESTED NUMBER OF QUADRATURE PANELS EXCEEDS THAT DEFINED\r
+       //             BY MQIN1   \r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int NMAX = 100;\r
+       double TAU[] = new double[1];\r
+       double MAXRM[] = new double[1];\r
+       double TOL,RIGHT;\r
+       boolean T1FXD;\r
+       \r
+       double REMND[][] = new double[NMAX][2];\r
+       //COMPLEX REMND(NMAX)\r
+       \r
+       // EXTERNAL DEJAC7,DELEG7\r
+       \r
+       if (NZZ*(MAXDG+1) > NMAX) {\r
+           IER[0]=10;\r
+           return;\r
+       }\r
+       \r
+       QINTS[0]=1;\r
+       XENPT[0]=-1.0;\r
+       TOL=TOLIN;\r
+       DEJAC7(ZZ,NZZ,BETA,TAU,MAXDG,NQUAD,AJAC,BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER);\r
+       /*      IF (IER .GT. 0) THEN\r
+               RETURN\r
+             ENDIF\r
+             TOLOU=MAXRM\r
+             XENPT(2)=TAU\r
+       C\r
+             IF (XENPT(2) .LT. 1E+0) THEN\r
+               QINTS=2\r
+               T1FXD=.FALSE.\r
+               TAU=1E+0\r
+               RIGHT=-1E+0\r
+               CALL DELEG7(ZZ,NZZ,BETA,RIGHT,TAU,T1FXD,MAXDG,NQUAD,AJAC,\r
+            +              BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER)\r
+               IF (IER .GT. 0) THEN\r
+                 RETURN\r
+               ENDIF\r
+               TOLOU=TOLOU+MAXRM\r
+               T1FXD=.TRUE.\r
+       C\r
+       100     CONTINUE\r
+       C\r
+               IF (XENPT(QINTS) .GT. RIGHT) THEN\r
+                 XENPT(QINTS)=5E-1*(XENPT(QINTS)+RIGHT)\r
+                 XENPT(QINTS+1)=1E+0\r
+               ELSE\r
+                 TAU=1E+0\r
+                 CALL DELEG7(ZZ,NZZ,BETA,XENPT(QINTS),TAU,T1FXD,MAXDG,\r
+            +                NQUAD,AJAC,BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER)\r
+                 IF (IER .GT. 0) THEN\r
+                   RETURN\r
+                 ENDIF\r
+                 TOLOU=TOLOU+MAXRM\r
+                 QINTS=QINTS+1\r
+                 IF (QINTS .GE. MQIN1) THEN\r
+                   IER=11\r
+                   RETURN\r
+                 ENDIF\r
+                 XENPT(QINTS)=TAU\r
+                 GOTO 100\r
+               ENDIF\r
+             ENDIF\r
+       C\r
+       C     NORMAL TERMINATION\r
+       C\r
+             IER=0\r
+       C*/\r
+    } // private void SUBIN7\r
+    \r
+    private void DEJAC7(double ZZ[][], int NZZ, double BETA, double TAU[], int MAXDG, int NQUAD,\r
+        double ACOEF[], double BCOEF[], double H0VAL,\r
+       double REMND[][], double CSCAL[], double TOL, double MAXRM[],int IER[]) {\r
+       //INTEGER MAXDG,NQUAD,NZZ,IER\r
+       //REAL BETA,TAU,H0VAL,TOL,MAXRM\r
+       //REAL ACOEF(*),BCOEF(*),CSCAL(*)\r
+       //COMPLEX ZZ(*),REMND(*)\r
+       \r
+       //     WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN\r
+       //     USING AN NQUAD-POINT GAUSS-JACOBI RULE TO ESTIMATE THE INTEGRALS\r
+       \r
+       //       INTEGRAL  [(1+X)**BETA*P(X,I)*LOG(ZZ(J)-X)*dX], I=0,MAXDG\r
+       //      -1<=X<=TAU                                       J=1,NZZ\r
+       \r
+       //     WHERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I\r
+       //     ASSOCIATED WITH THE WEIGHT (1+X)**BETA.  THE REMAINDER \r
+       //     CORRESPONDING TO P(.,I) AND ZZ(J) IS STORED IN \r
+       //     REMND(I+J+MAXDG*(J-1)), I=0,MAXDG, J=1,NZZ.  THIS ROUTINE USES\r
+       //     THE SIMPLEST POSSIBLE ESTIMATES; I.E. THE LEADING TERM ONLY IN\r
+       //     THE ASYMPTOTIC EXPANSION AND THE WATSON-DOETSCH ESTIMATE FOR ANY\r
+       //     INTEGRALS.\r
+       \r
+       //     THE PURPOSE OF THIS ROUTINE IS THEN TO DETERMINE A VALUE FOR TAU\r
+       //     SUCH THAT\r
+       \r
+       //         ABS( REAL(REMND(I)) )*CSCAL(I) < TOL , I=1,NZZ*(MAXDG+1)\r
+       \r
+       //     AND THAT, IF POSSIBLE, \r
+       \r
+       //        0.5*TOL <= ABS( REAL(REMND(I)) )*CSCAL(I) < TOL\r
+       \r
+       //     FOR AT LEAST ONE VALUE OF I.\r
+       \r
+       //     IER=0 - NORMAL EXIT\r
+       //     IER=12- LOCAL PARAMETER NC NEEDS INCREASING TO AT LEAST NZZ\r
+       //             (THIS ERROR CAN'T ARISE IN THE PRESENT VERSION, SINCE\r
+       //              NZZ IS FIXED AT 2)\r
+       //     IER=13- LOCAL PARAMETER NR NEEDS INCREASING TO AT LEAST MAXDG\r
+       //             (AT PRESENT MAXDG=NQPTS-1) \r
+       //     IER=14- A JACOBI INDEX MAY BE LARGE ENOUGH TO CAUSE OVERFLOW IN\r
+       //             THE GAMMA FUNCTION (AN ANGLE ON THE PHYSICAL BOUNDARY\r
+       //             MUST BE LESS THAN ABOUT 6 DEGREES)\r
+       \r
+       //     LOCAL VARIABLES..\r
+       \r
+       final int NC = 8;\r
+       final int NR = 30;\r
+       int I,J,K,LIM;\r
+       double S,KK,SUM1,RI,TURI,RN,OFLOW,P0SCL,TUK,LOWER,UPPER,TERM,HTOL;\r
+       double XI[] = new double[2];\r
+       double Z1[] = new double[2];\r
+       double XI1[] = new double[2];\r
+       double FF[] = new double[2];\r
+       double PRE[] = new double[2];\r
+       double CUR[] = new double[2];\r
+       double NXT[] = new double[2];\r
+       //COMPLEX XI,Z1,XI1,FF,PRE,CUR,NXT\r
+       double GG[][] = new double[NC][2];\r
+       double CONST[][][] = new double[NR][NC][2];\r
+       //COMPLEX GG(NC),CONST(NR,NC)\r
+       //REAL GAMMA, LGGAM\r
+       //EXTERNAL GAMMA,LGGAM\r
+       double cr[] = new double[1];\r
+       double ci[] = new double[1];\r
+       \r
+       if (NZZ > NC ) {\r
+           IER[0]=12;\r
+           return;\r
+       }\r
+       \r
+       if (MAXDG >= NR) {\r
+           IER[0]=13;\r
+           return;\r
+       }\r
+       \r
+       S=BETA+4.0;\r
+       if (S > 20.0) {\r
+       \r
+           // TEST FOR POSSIBLE OVERFLOW IN GAMMA FUNCTION\r
+       \r
+           OFLOW=Math.log(Double.MAX_VALUE);\r
+           KK=LGGAM(S);\r
+           if (KK > OFLOW) {\r
+               IER[0]=14;\r
+               return;\r
+           }\r
+           else {\r
+               KK=Math.exp(-KK);\r
+           }\r
+       } // if (S > 20.0)\r
+       else {\r
+           KK=1.0/GAMMA(S);\r
+       }\r
+       \r
+       // FIRST WE COMPUTE THE FACTORS WHICH ARE INDEPENDENT OF TAU\r
+       \r
+       S=S-1.0;\r
+       KK=Math.pow(4.0,S)*KK*GAMMA(BETA+2.0)/(S-1.0);\r
+       SUM1=BETA+1.0;\r
+       for (I=2; I <= NQUAD; I++) {\r
+           RI=(double)(I);\r
+           TURI=2.0*RI;\r
+           KK=KK*16.0*(RI+BETA)/(TURI+SUM1);\r
+           KK=KK*RI/(TURI+BETA);\r
+           KK=KK*(RI+BETA)/(TURI+BETA);\r
+           KK=KK*RI/(TURI+BETA-1.0);\r
+       } // for (I=2; I <= NQUAD; I++)\r
+       RN=(double)(NQUAD);\r
+       TUK=2.0*RN+SUM1;\r
+       KK=-KK/TUK/2.0;\r
+       \r
+       /*for (I=1; I <= NZZ; I++) {\r
+           GG[I-1]=Math.pow((1.0+ZZ[I-1]),BETA)*KK;\r
+       } // for (I=1; I <= NZZ; I++)\r
+       \r
+       // NOW GIVE THE JACOBI POLYNOMIALS THE SCALING CORRESPONDING TO\r
+        // [-1,1] AS STANDARD INTERVAL\r
+       \r
+       P0SCL=1.0/Math.sqrt(H0VAL);\r
+       for (J=1; J <= NZZ; J++) {\r
+           PRE[0]=P0SCL;\r
+           PRE[1] = 0.0;\r
+           CUR[0] = PRE[0];\r
+           CUR[1] = PRE[1];\r
+           zmlt(CUR[0], CUR[1], GG[J-1][0],GG[J-1][1], cr, ci);\r
+           CONST[0][J-1][0] = cr[0] * CSCAL[0];\r
+           CONST[0][J-1][1] = ci[0] * CSCAL[0];\r
+           if (MAXDG >= 1) {\r
+                 CUR=(ZZ(J)-BCOEF(1))*PRE/ACOEF(1)\r
+                 CONST(2,J)=CUR*GG(J)*CSCAL(2)\r
+                 DO 200 I=2,MAXDG\r
+                   NXT=((ZZ(J)-BCOEF(I))*CUR-ACOEF(I-1)*PRE)/ACOEF(I)\r
+                   PRE=CUR\r
+                   CUR=NXT\r
+                   CONST(I+1,J)=CUR*GG(J)*CSCAL(I+1)\r
+       200       CONTINUE\r
+           } // if (MAXDG >= 1)\r
+       } // for (J=1; J <= NZZ; J++)\r
+       C\r
+       C     NOW COME THE FACTORS DEPENDENT ON TAU\r
+       C\r
+             TAU=1E+0\r
+             LOWER=-1E+0\r
+             UPPER=1E+0\r
+             LIM=NZZ*(MAXDG+1)\r
+       C\r
+       250   CONTINUE\r
+       C\r
+             HTOL=5E-1*TOL\r
+             K=0\r
+             DO 325 J=1,NZZ\r
+               XI=(2E+0*ZZ(J)+1E+0-TAU)/(1E+0+TAU)\r
+               Z1=SQRT(XI*XI-1E+0)\r
+               XI1=XI+Z1\r
+               IF (ABS(XI1) .LT. 1E+0) THEN\r
+                 XI1=XI-Z1\r
+               ENDIF\r
+               FF=XI1**(-TUK-1E+0)*(XI1*XI1-1E+0)*(1E+0+TAU)*5E-1\r
+               DO 300 I=0,MAXDG\r
+                 K=K+1\r
+                 REMND(K)=CONST(I+1,J)*FF\r
+       300     CONTINUE\r
+       325   CONTINUE\r
+       C\r
+             MAXRM=0E+0\r
+             DO 600 I=1,LIM\r
+               TERM=ABS(REAL(REMND(I)))\r
+               MAXRM=MAX(MAXRM,TERM)\r
+       600   CONTINUE\r
+       C\r
+             IF (MAXRM .LT. TOL) THEN\r
+       C\r
+       C       ACCURACY IS ACHIEVED, BUT MAYBE TAU COULD BE INCREASED.\r
+       C\r
+               IF (MAXRM .LT. HTOL) THEN\r
+       C\r
+       C         TAU NEEDS INCREASING, BUT THIS IS ONLY POSSIBLE IF TAU<1.\r
+       C\r
+                 IF (TAU .LT. 1E+0) THEN\r
+                   LOWER=TAU\r
+                   TAU=5E-1*(LOWER+UPPER)\r
+                   GOTO 250\r
+                 ENDIF\r
+               ENDIF\r
+             ELSE\r
+       C\r
+       C       ACCURACY NOT ACHIEVED AND TAU NEEDS DECREASING.\r
+       C\r
+               IF (TAU .EQ. 1E+0) THEN\r
+                 TOL=HTOL\r
+               ENDIF\r
+               UPPER=TAU\r
+               TAU=5E-1*(LOWER+UPPER)\r
+               GOTO 250\r
+             ENDIF\r
+       C\r
+       C     NORMAL TERMINATION\r
+       C\r
+             IER=0\r
+       C*/\r
+    } // private void DEJAC7\r
+    \r
+\r
+    private double LGGAM(double X) {\r
+\r
+        // **** TO ESTIMATE THE LOGARITHM OF THE GAMMA FUNCTION FOR  L A R G E\r
+        // **** POSITIVE VALUES OF X USING THE ASYMPTOTIC FORMULA FROM ABRAMOWITZ\r
+        // **** AND STEGUN, SECTION 6.1.41\r
+\r
+        double PI,W,result;\r
+\r
+        PI=Math.PI;\r
+        W=1.0/X/X;\r
+        result=\r
+                ((((W/9.9E+1-1E+0/1.4E+2)*W+1E+0/1.05E+2)*W-1E+0/3E+1)*W+1E+0)/\r
+                (1.2E+1*X) + 5E-1*Math.log(2E+0*PI) - X + (X-5E-1)*Math.log(X);\r
+        return result;\r
+    } // private double LGGAM\r
+\r
 \r
 \r
 \r