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

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

index f10d6f7..2481726 100644 (file)
@@ -4634,8 +4634,14 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        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
+               double AJAC[] = new double[MDGPO];\r
+               double BJAC[] = new double[MDGPO];\r
+               for (I = 0; I < MDGPO; I++) {\r
+                       AJAC[I] = ACOEF[LO-1+I];\r
+                       BJAC[I] = BCOEF[LO-1+I];\r
+               }\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
@@ -4756,52 +4762,55 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        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
+       if (IER[0] > 0) {\r
+           return;\r
+       }\r
+       TOLOU=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
+                   BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER);\r
+           if (IER[0] > 0) {\r
+               return;\r
+           }\r
+           TOLOU=TOLOU+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
+                   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
+                          NQUAD,AJAC,BJAC,H0JAC,REMND,CSCAL,TOL,MAXRM,IER);\r
+                 if (IER[0] > 0) {\r
+                     return;\r
+                 }\r
+                 TOLOU=TOLOU+MAXRM[0];\r
+                 QINTS[0]=QINTS[0]+1;\r
+                 if (QINTS[0] >= MQIN1) {\r
+                     IER[0]=11;\r
+                     return;\r
+                 }\r
+                 XENPT[QINTS[0]-1]=TAU[0];\r
+                 continue;\r
+               } // else\r
+       } // while (true)\r
+    } // if (XENPT[1] < 1.0)\r
+       \r
+    // NORMAL TERMINATION\r
+       \r
+    IER[0]=0;\r
+       \r
     } // private void SUBIN7\r
     \r
     private void DEJAC7(double ZZ[][], int NZZ, double BETA, double TAU[], int MAXDG, int NQUAD,\r
@@ -4868,6 +4877,14 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        //EXTERNAL GAMMA,LGGAM\r
        double cr[] = new double[1];\r
        double ci[] = new double[1];\r
+       double r;\r
+       double theta;\r
+       double mag;\r
+       double arg;\r
+       double absXI1;\r
+       double var[] = new double[2];\r
+       double cr2[] = new double[1];\r
+       double ci2[] = new double[1];\r
        \r
        if (NZZ > NC ) {\r
            IER[0]=12;\r
@@ -4915,8 +4932,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        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 = zabs(1.0 + ZZ[I-1][0], ZZ[I-1][1]);\r
+               theta = Math.atan2(ZZ[I-1][1],1.0+ZZ[I-1][0]);\r
+               mag = Math.pow(r, BETA);\r
+               arg = BETA * theta;\r
+               GG[I-1][0] = mag * Math.cos(arg)*KK;\r
+               GG[I-1][1] = mag * Math.sin(arg)*KK;\r
        } // for (I=1; I <= NZZ; I++)\r
        \r
        // NOW GIVE THE JACOBI POLYNOMIALS THE SCALING CORRESPONDING TO\r
@@ -4932,78 +4954,117 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
            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
+               zmlt(ZZ[J-1][0] - BCOEF[0], ZZ[J-1][1], PRE[0]/ACOEF[0], PRE[1]/ACOEF[0], cr, ci);\r
+               CUR[0] = cr[0];\r
+               CUR[1] = ci[0];\r
+               zmlt(CUR[0], CUR[1], GG[J-1][0]*CSCAL[1], GG[J-1][1]*CSCAL[1], cr, ci);\r
+               CONST[1][J-1][0] = cr[0];\r
+               CONST[1][J-1][1] = ci[0];\r
+               for (I=2; I <= MAXDG; I++) {\r
+                       zmlt(ZZ[J-1][0] - BCOEF[I-1], ZZ[J-1][1], CUR[0], CUR[1], cr, ci);\r
+                       NXT[0] = (cr[0] - ACOEF[I-2]*PRE[0])/ACOEF[I-1];\r
+                       NXT[1] = (ci[0] - ACOEF[I-2]*PRE[1])/ACOEF[I-1];\r
+                   PRE[0] = CUR[0];\r
+                   PRE[1] = CUR[1];\r
+                   CUR[0]=NXT[0];\r
+                   CUR[1]=NXT[1];\r
+                   zmlt(CUR[0],CUR[1],GG[J-1][0]*CSCAL[I],GG[J-1][1]*CSCAL[I], cr, ci);\r
+                   CONST[I][J-1][0] = cr[0];\r
+                   CONST[I][J-1][1] = ci[0];\r
+               } // for (I=2; I <= MAXDG; I++)\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
+       \r
+       // NOW COME THE FACTORS DEPENDENT ON TAU\r
+       \r
+       TAU[0]=1.0;\r
+       LOWER=-1.0;\r
+       UPPER=1.0;\r
+       LIM=NZZ*(MAXDG+1);\r
+       \r
+       while (true) {\r
+       \r
+           HTOL=0.5*TOL;\r
+           K=0;\r
+           for (J=1; J <= NZZ; J++) {\r
+               XI[0]=(2.0*ZZ[J-1][0]+1.0-TAU[0])/(1.0+TAU[0]);\r
+               XI[1]= (2.0*ZZ[J-1][1])/(1.0 + TAU[0]);\r
+               zmlt(XI[0],XI[1],XI[0],XI[1],cr,ci);\r
+               r = zabs(cr[0]-1.0,ci[0]);\r
+               theta = Math.atan2(ci[0],cr[0]-1.0);\r
+               mag = Math.sqrt(r);\r
+               arg = 0.5*theta;\r
+               Z1[0] = mag * Math.cos(arg);\r
+               Z1[1] = mag * Math.sin(arg);\r
+               XI1[0]=XI[0]+Z1[0];\r
+               XI1[1]=XI[1]+Z1[1];\r
+               absXI1 = zabs(XI1[0],XI1[1]);\r
+               if (absXI1 < 1.0) {\r
+                   XI1[0]=XI[0]-Z1[0];\r
+                   XI1[1]=XI[1]-Z1[1];\r
+               }\r
+               r = zabs(XI1[0],XI1[1]);\r
+               theta = Math.atan2(XI1[1], XI1[0]);\r
+               mag = Math.pow(r, -TUK-1.0);\r
+               arg = theta * (-TUK-1.0);\r
+               var[0] = mag * Math.cos(arg);\r
+               var[1] = mag * Math.sin(arg);\r
+               zmlt(XI1[0],XI1[1],XI1[0],XI1[1],cr,ci);\r
+               zmlt(var[0],var[1],cr[0]-1.0,ci[0],cr2,ci2);\r
+               FF[0] = cr2[0]*(1.0+TAU[0])*0.5;\r
+               FF[1] = ci2[0]*(1.0+TAU[0])*0.5;\r
+               for (I=0; I <= MAXDG; I++) {\r
+                 K=K+1;\r
+                 zmlt(CONST[I][J-1][0],CONST[I][J-1][1],FF[0],FF[1],cr,ci);\r
+                 REMND[K-1][0]= cr[0];\r
+                 REMND[K-1][1] = ci[0];\r
+               } // for (I=0; I <= MAXDG; I++)\r
+           } // for (J=1; J <= NZZ; J++)\r
+       \r
+           MAXRM[0]=0.0;\r
+           for (I=1; I <= LIM; I++) {\r
+               TERM=Math.abs(REMND[I-1][0]);\r
+               MAXRM[0]=Math.max(MAXRM[0],TERM);\r
+           } // for (I=1; I <= LIM; I++)\r
+       \r
+           if (MAXRM[0] < TOL) {\r
+       \r
+               // ACCURACY IS ACHIEVED, BUT MAYBE TAU COULD BE INCREASED.\r
+       \r
+               if (MAXRM[0] < HTOL) {\r
+       \r
+                   // TAU NEEDS INCREASING, BUT THIS IS ONLY POSSIBLE IF TAU<1.\r
+       \r
+                   if (TAU[0] < 1.0) {\r
+                       LOWER=TAU[0];\r
+                       TAU[0]=0.5*(LOWER+UPPER);\r
+                       continue;\r
+                   } // if (TAU[0 < 1.0)\r
+                   else {\r
+                       break;\r
+                   }\r
+               } // if (MAXRM[0] < HTOL)\r
+               else {\r
+                       break;\r
+               }\r
+           }\r
+           else { // MAXRM[0] >= TOL\r
+       \r
+                // ACCURACY NOT ACHIEVED AND TAU NEEDS DECREASING.\r
+       \r
+               if (TAU[0] == 1.0) {\r
+                   TOL=HTOL;\r
+               }\r
+               UPPER=TAU[0];\r
+               TAU[0]=0.5*(LOWER+UPPER);\r
+               continue;\r
+           } // else MAXRM[0] >= TOL\r
+       } // while (true)\r
+       \r
+       // NORMAL TERMINATION\r
+       \r
+       IER[0]=0;\r
+       \r
     } // private void DEJAC7\r
     \r
 \r
@@ -5023,7 +5084,317 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         return result;\r
     } // private double LGGAM\r
 \r
-\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
+       //REAL BETA,TAU1,TAU2,H0VAL,TOL,MAXRM\r
+       //REAL ACOEF(*),BCOEF(*),CSCAL(*)\r
+       //LOGICAL T1FXD\r
+       //COMPLEX ZZ(*),REMND(*)\r
+       \r
+       // WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN\r
+       // USING AN NQUAD-POINT GAUSS-LEGENDRE RULE TO ESTIMATE THE INTEGRALS\r
+    \r
+       //   INTEGRAL  [(1+X)**BETA*P(X,I)*LOG(ZZ(J)-X)*dX], I=0,1,...,MAXDG\r
+       // TAU1<=X<=TAU2                                     J=1,2\r
+       \r
+       //     WHERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I\r
+       //     ASSOCIATED WITH THE WEIGHT (1+X)**BETA AND -1<TAU1<TAU2<=1.\r
+       //     THE REMAINDER CORRESPONDING TO P(.,I) AND ZZ(J) IS STORED IN \r
+       //     REMND(I+J+MAXDG*(J-1)), I=0,1,...,MAXDG, J=1,2.  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 TO DETERMINE A VALUE FOR EITHER\r
+       //     TAU2 (TAU1 REMAINS FIXED IF T1FXD IS "TRUE") OR TAU1 (TAU2 REMAINS\r
+       //     FIXED IF T1FXD IS "FALSE") SUCH THAT\r
+       \r
+       //         ABS( REAL(REMND(I)) )*CSCAL(I) < TOL , I=1,2*MAXDG+2\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
+       //     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
+       //\r
+       //     LOCAL VARIABLES..\r
+       \r
+       final int NC = 8;\r
+       final int NR = 30;\r
+       int I,J,K;\r
+       double KK,RI,TURI,RN,P0SCL,TUK,LOWER,HTOL,UFLOW,EXPON,\r
+            UPPER,TERM,HCO,PI,FORK,SUM,FAC1,FAC2,PVAL,RR,MEAN,BB,BB1,BB2,FFH;\r
+       double XI[] = new double[2];\r
+       double Z1[] = new double[2];\r
+       double XI1[] = new double[2];\r
+       double FFG[] = 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,FFG,PRE,CUR,NXT\r
+       boolean FIRST;\r
+       double GG[][] = new double[NC][2];\r
+       double HH[][] = new double[NC][2];\r
+       double CONGG[][][] = new double[NR][NC][2];\r
+       double CONHH[][][] = new double[NR][NC][2]; \r
+       // COMPLEX GG(NC),HH(NC),CONGG(NR,NC),CONHH(NR,NC)\r
+       // EXTERNAL GAMMA,LGGAM\r
+       double r;\r
+       double mag;\r
+       double theta;\r
+       double ang;\r
+       double cr[] = new double[1];\r
+       double ci[] = new double[1];\r
+       double absXI1;\r
+       double var[] = new double[2];\r
+       double cr2[] = new double[1];\r
+       double ci2[] = new double[1];\r
+       \r
+       // FIRST WE COMPUTE THE FACTORS WHICH ARE INDEPENDENT OF TAU1,TAU2\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
+       // **** SET THE LOGARITHMIC UNDERFLOW LIMIT\r
+       // B**(EMIN-1)\r
+       // EMIN = -1022;\r
+       UFLOW = Math.log(Math.pow(2.0, -1023.0));    \r
+       \r
+       PI=Math.PI;\r
+       KK=32.0/6.0;\r
+       for (I=2; I <= NQUAD; I++) {\r
+           RI=(double)(I);\r
+           TURI=2.0*RI;\r
+           KK=KK*4.0*RI/(TURI+1.0);\r
+           KK=KK*RI/(TURI-1.0);\r
+       } // for (I=2; I <= NQUAD; I++)\r
+       RN=(double)(NQUAD);\r
+       TUK=2.0*RN+1.0;\r
+       FORK=2.0*TUK;\r
+       KK=KK/FORK;\r
+       \r
+       if (BETA >= 20.0) {\r
+           EXPON=LGGAM(BETA+1.0)-BETA*Math.log(FORK);\r
+           if (EXPON <= UFLOW) {\r
+               HCO=0.0;\r
+           }\r
+           else {\r
+               HCO=Math.sin(PI*BETA)*Math.exp(EXPON)/PI;\r
+           }\r
+       } // if (BETA >= 20.0)\r
+       else {\r
+           HCO=Math.sin(PI*BETA)*GAMMA(BETA+1.0)/PI/Math.pow(FORK,BETA);\r
+       }\r
+       for (I=1; I <= NZZ; I++) {\r
+               r = zabs(1.0+ZZ[I-1][0],ZZ[I-1][1]);\r
+               mag = Math.pow(r, BETA);\r
+               theta = Math.atan2(ZZ[I-1][1], 1.0 + ZZ[I-1][0]);\r
+               ang = BETA * theta;\r
+               GG[I-1][0] = mag * Math.cos(ang) * KK;\r
+               GG[I-1][1] = mag * Math.sin(ang) * KK;\r
+               HH[I-1][0] = -Math.log(r)*HCO*KK;\r
+               HH[I-1][1] = -theta * HCO * 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]*CSCAL[0],GG[J-1][1]*CSCAL[0],cr,ci);\r
+           CONGG[0][J-1][0] = cr[0];\r
+           CONGG[0][J-1][1] = ci[0];\r
+           zmlt(CUR[0],CUR[1],HH[J-1][0]*CSCAL[0],HH[J-1][1]*CSCAL[0],cr,ci);\r
+           CONHH[0][J-1][0] = cr[0];\r
+           CONHH[0][J-1][1] = ci[0];\r
+           if (MAXDG >= 1) {\r
+               zmlt(ZZ[J-1][0] - BCOEF[0],ZZ[J-1][1],PRE[0]/ACOEF[0],PRE[1]/ACOEF[0],cr,ci);\r
+               CUR[0] = cr[0];\r
+               CUR[1] = ci[0];\r
+               zmlt(CUR[0],CUR[1],GG[J-1][0]*CSCAL[1],GG[J-1][1]*CSCAL[1],cr,ci);\r
+               CONGG[1][J-1][0] = cr[0];\r
+               CONGG[1][J-1][1] = ci[0];\r
+               for (I=2; I <= MAXDG; I++) {\r
+                       zmlt(ZZ[J-1][0]-BCOEF[I-1],ZZ[J-1][1],CUR[0],CUR[1],cr,ci);\r
+                       NXT[0] = (cr[0] - ACOEF[I-2]*PRE[0])/ACOEF[I-1];\r
+                       NXT[1] = (ci[0] - ACOEF[I-2]*PRE[1])/ACOEF[I-1];\r
+                   PRE[0]=CUR[0];\r
+                   PRE[1]=CUR[1];\r
+                   CUR[0]=NXT[0];\r
+                   CUR[1]=NXT[1];\r
+                   zmlt(CUR[0],CUR[1],GG[J-1][0]*CSCAL[I],GG[J-1][1]*CSCAL[I],cr,ci);\r
+                   CONGG[I][J-1][0] = cr[0];\r
+                   CONGG[I][J-1][1] = ci[0];\r
+               } // for (I=2; I <= MAXDG; I++)\r
+           } // if (MAXDG >= 1)\r
+       } // for (J=1; J <= NZZ; J++)\r
+       \r
+       // NOW COMPUTE THE POLYNOMIAL VALUES AT -1, SCALE AND ACCUMULATE INTO\r
+       // CONHH\r
+       \r
+       SUM=BETA+1.0;\r
+       FAC1=1.0/Math.sqrt(Math.pow(2.0,SUM));\r
+       FAC2=1.0;\r
+       for (I=1; I <= MAXDG; I++) {\r
+           SUM=SUM+2.0;\r
+           FAC1=-FAC1;\r
+           FAC2=(I+BETA)*FAC2/I;\r
+           PVAL=Math.sqrt(SUM)*FAC1*FAC2;\r
+           for (J=1; J <= 2; J++) {\r
+               CONHH[I][J-1][0]=PVAL*HH[J-1][0]*CSCAL[I];\r
+               CONHH[I][J-1][1]=PVAL*HH[J-1][1]*CSCAL[I];\r
+           } // for (J=1; J <= 2; J++)\r
+       } // for (I=1; I <= MAXDG; I++)\r
+       \r
+        // NOW COME THE FACTORS DEPENDENT ON TAU1 AND TAU2.\r
+       \r
+       LOWER=TAU1;\r
+       UPPER=TAU2;\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
+           BB=(1.0+MEAN)/RR;\r
+           BB1=BB+Math.sqrt(BB*BB-1.0);\r
+       \r
+           // **** NOW COMPUTE THE QUANTITY\r
+           // ****\r
+           //****      FFH=(RR*(BB1-1E+0/BB1))**(BETA+1E+0)/BB1**TUK\r
+           //****\r
+           // **** BUT CHECK FOR POSSIBLE UNDERFLOW\r
+       \r
+           BB2=BB1-1.0/BB1;\r
+           if (BB2 <= 0.0) {\r
+               FFH=0.0;\r
+           }\r
+           else {\r
+               EXPON=(BETA+1.0)*Math.log(RR*BB2)-TUK*Math.log(BB1);\r
+               if (EXPON <= UFLOW) {\r
+                   FFH=0.0;\r
+               }\r
+               else {\r
+                   FFH=Math.exp(EXPON);\r
+               }\r
+           }\r
+           K=0;\r
+           for (J=1; J <= NZZ; J++) {\r
+               XI[0]=(ZZ[J-1][0]-MEAN)/RR;\r
+               XI[1] = ZZ[J-1][1]/RR;\r
+               zmlt(XI[0],XI[0],XI[1],XI[1],cr,ci);\r
+               r = zabs(cr[0]-1.0,ci[0]);\r
+               mag = Math.sqrt(r);\r
+               theta = Math.atan2(ci[0], cr[0]-1.0);\r
+               ang = 0.5 * theta;\r
+               Z1[0] = mag;\r
+               Z1[1] = ang;\r
+               XI1[0]=XI[0]+Z1[0];\r
+               XI1[1]=XI[1]+Z1[1];\r
+               absXI1 = zabs(XI1[0],XI1[1]);\r
+               if (absXI1 < 1.0) {\r
+                 XI1[0]=XI[0]-Z1[0];\r
+                 XI1[1]=XI[1]-Z1[1];\r
+               }\r
+               r = zabs(XI1[0],XI1[1]);\r
+               mag = Math.pow(r, -TUK-1.0);\r
+               theta = Math.atan2(XI1[1],XI1[0]);\r
+               ang = theta *(-TUK - 1.0);\r
+               var[0] = mag * Math.cos(ang);\r
+               var[1] = mag * Math.sin(ang);\r
+               zmlt(XI1[0],XI1[1],XI1[0],XI1[1],cr,ci);\r
+               zmlt(var[0],var[1],cr[0]-1.0,ci[0],cr2,ci2);\r
+               FFG[0] = cr2[0] * RR;\r
+               FFG[1] = ci2[0] * RR;\r
+               for (I=0; I <= MAXDG; I++) {\r
+                   K=K+1;\r
+                   zmlt(CONGG[I][J-1][0],CONGG[I][J-1][1],FFG[0],FFG[1],cr,ci);\r
+                   REMND[K-1][0] = cr[0] + CONHH[I][J-1][0]*FFH;\r
+                   REMND[K-1][1] = cr[0] + CONHH[I][J-1][1]*FFH;\r
+               } // for (I=0; I <= MAXDG; I++)\r
+           } // for (J=1; J <= NZZ; J++)\r
+       \r
+             MAXRM[0]=0.0;\r
+             for (I=1; I <= 2*MAXDG+2; I++) {\r
+               TERM=Math.abs(REMND[I-1][0]);\r
+               MAXRM[0]=Math.max(MAXRM[0],TERM);\r
+             } // for (I=1; I <= 2*MAXDG+2; I++)\r
+       \r
+             if (MAXRM[0] < TOL) {\r
+       \r
+                 // ACCURACY IS ACHIEVED, BUT MAYBE TAU2 COULD BE INCREASED OR\r
+                 // TAU1 DECREASED\r
+       \r
+                 if (MAXRM[0] < HTOL) {\r
+       \r
+                     // TAU2 NEEDS INCREASING IF T1FXD (BUT THIS IS ONLY POSSIBLE IF \r
+                     // 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
+                         continue;\r
+                     }\r
+                     else if (!T1FXD && TAU1 > -1.0) {\r
+                         UPPER=TAU1;\r
+                         TAU1=0.5*(LOWER+UPPER);\r
+                         continue;\r
+                     }\r
+                     else {\r
+                         break;\r
+                     }\r
+                 } // if (MAXRM[0] < HTOL\r
+                 else {\r
+                         break;\r
+                 }\r
+             } // if (MAXRM[0] < TOL)\r
+             else { // MAXRM[0] >= TOL\r
+       \r
+                 // ACCURACY NOT ACHIEVED AND TAU2 NEEDS DECREASING OR TAU1 NEEDS\r
+                 // INCREASING.\r
+       \r
+                 if (FIRST) {\r
+                     TOL=HTOL;\r
+                     FIRST=false;\r
+                 }\r
+                 if (T1FXD) {\r
+                     UPPER=TAU2;\r
+                     TAU2=0.5*(LOWER+UPPER);\r
+                 }\r
+                 else {\r
+                     LOWER=TAU1;\r
+                     TAU1=0.5*(LOWER+UPPER);\r
+                 }\r
+                 continue;\r
+             } // else MAXRM[0] >= TOL\r
+       } // while (true)\r
+       \r
+       // NORMAL TERMINATION\r
+       \r
+       IER[0]=0;\r
+       \r
+    } // private void DELEG7\r
 \r
 \r
             /* SUBROUTINE RSLT71(QIERC,RCOND,SOLUN,NEQNS,LOSUB,HISUB,COLSC,\r