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

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

index bcb2924..b2e5827 100644 (file)
@@ -8654,68 +8654,658 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        //     INITIALISATION\r
        \r
        PI=Math.PI;\r
-             /*TAUI=1E+0\r
-             DP=DELTA*PI\r
-             DT=DELTA*2E+0\r
-             HDP=5E-1*DP\r
-             LEN=2E+0+DP\r
-             SINC=LEN/NINTS\r
-             STAR=SINC\r
+       TAUI=1.0;\r
+       DP=DELTA*PI;\r
+       DT=DELTA*2.0;\r
+       HDP=0.5*DP;\r
+       LEN=2.0+DP;\r
+       SINC=LEN/NINTS;\r
+       STAR=SINC;\r
+       \r
+       // START OF LOOP FOR INTERVAL HALVING\r
+       \r
+       while (true) {\r
+           NXI=0;\r
+           SS1=-STAR+SINC;\r
+           SS2=LEN;\r
+           for (SS=SS1; SS <= SS2; SS += SINC) {\r
+               if  (SS < HDP) {\r
+                   SS3=SS/DELTA;\r
+                   XI[0] = 1.0 + DELTA * Math.cos(SS3);\r
+                   XI[1] = DELTA * Math.sin(SS3);\r
+               }\r
+               else if (SS > 2.0+HDP) {\r
+                   SS3=(SS-2.0-HDP)/DELTA;\r
+                   XI[0] = -1.0 - DELTA * Math.sin(SS3);\r
+                   XI[1] = DELTA * Math.cos(SS3);\r
+               }\r
+               else {\r
+                   SS3=SS-HDP;\r
+                   XI[0] = 1.0-SS3;\r
+                   XI[1] = DELTA;\r
+               }\r
+               RXI=XI[0];\r
+               if (RXI < (TAUI+DT)) {\r
+                 NXI=NXI+1;\r
+                 if (NXI > MNXI) {\r
+                     IER[0]=25;\r
+                     return;\r
+                 }\r
+                 XIVAL[NXI-1][0]=XI[0];\r
+                 XIVAL[NXI-1][1]=XI[1];\r
+               } // \r
+           } // for (SS=SS1; SS <= SS2; SS += SINC)\r
+       \r
+           if (NXI == 0) {\r
+               SINC=STAR;\r
+               STAR=5E-1*STAR;\r
+               NINTS=NINTS*2;\r
+               continue;\r
+           } // if (NXI == 0)\r
+            \r
+           TAU[0]=TAUI;\r
+           DEPPJ9(XIVAL,NXI,BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL,SOLUN,TOL,MAXRM,IER);\r
+           if (IER[0] > 0) {\r
+               return;\r
+           }\r
+       \r
+           if (TAU[0] != TAUI) {\r
+               TAUI=TAU[0];\r
+               SINC=STAR;\r
+               STAR=0.5*STAR;\r
+               NINTS=NINTS*2;\r
+               continue;\r
+           }\r
+           break;\r
+       } // while (true)\r
+       \r
+        IER[0]=0;\r
+         \r
+    } // private void DEPPJ8\r
+    \r
+    private void DEPPJ9(double ZZ[][], int NZZ, double BETA, double TAU[], int NQUAD,\r
+        int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][],\r
+        double TOL, double MAXRM[], int IER[]) {\r
+        // INTEGER NQUAD,NZZ,IER,DGPOL\r
+        // REAL BETA,TAU,TOL,MAXRM,ACOEF(*),BCOEF(*),H0VAL\r
+        // COMPLEX SOLUN(*),ZZ(*)\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*FNPHI(X)*LOG(ZZ(I)-X)*dX], I=1,NZZ\r
+        //-1<=X<=TAU                                       \r
+\r
+        // WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY \r
+        // CORRESPONDENCE DERIVATIVE - JACOBI WEIGHT QUOTIENT AND ZZ IS\r
+        // A GIVEN ARRAY OF POINTS. \r
+\r
+        // THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO \r
+        // DEFINE FNPHI.\r
+\r
+        // THE MAXIMUM ABSOLUTE VALUE OF ALL THE REMAINDERS CORRESPONDING TO \r
+        // ZZ(I) , I=1,NZZ, IS STORED IN MAXRM.\r
\r
+        // THIS ROUTINE USES THE SIMPLEST POSSIBLE ESTIMATES; I.E. THE   \r
+        // LEADING TERM ONLY IN THE ASYMPTOTIC  EXPANSION AND THE WATSON-\r
+        // DOETSCH ESTIMATE FOR ANY INTEGRALS.\r
+\r
+        // THE PURPOSE OF THIS ROUTINE IS THEN TO DETERMINE A VALUE FOR TAU\r
+        // SUCH THAT\r
+\r
+        //                    MAXRM < TOL \r
+\r
+        // AND THAT, IF POSSIBLE, \r
+\r
+        //               0.5*TOL <= MAXRM < TOL.\r
+\r
+        // IER=0 - NORMAL EXIT\r
+        // IER=26- TOO MANY TEST POINTS ON DELTA CONTOUR; INCREASE \r
+        //         PARAMETER MAXNZ BELOW\r
+        // IER=14- BETA MAY CAUSE OVERFLOW IN GAMMA FUNCTION; AN ANGLE\r
+        //         ON THE BOUNDARY IS TOO SMALL\r
+\r
+        // LOCAL VARIABLES\r
+\r
+        final int MAXNZ = 200;\r
+       int I;\r
+        double S,KK,SUM1,RI,TURI,RN,TUK,LOWER,UPPER,TERM,HTOL,TAUI,OFLOW;\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 FNPHI[] = new double[2];\r
+        double REMND[] = new double[2];\r
+        // COMPLEX XI,Z1,XI1,FF,FNPHI,REMND\r
+        double GG[][] = new double[MAXNZ][2];\r
+        // COMPLEX GG(MAXNZ)\r
+        // REAL GAMMA, LGGAM\r
+        // COMPLEX CCJACS\r
+        // EXTERNAL GAMMA,CCJACS,LGGAM\r
+        double r;\r
+        double theta;\r
+        double mag;\r
+        double ang;\r
+        double realVar;\r
+        double imagVar;\r
+        double cr[] = new double[1];\r
+        double ci[] = new double[1];\r
+        double cr2[] = new double[1];\r
+        double ci2[] = new double[1];\r
+\r
+        if (NZZ > MAXNZ ) {\r
+\r
+            // SOME LOCAL ARRAY BOUNDS MUST BE INCREASED\r
+\r
+            IER[0]=26;\r
+            return;\r
+        } // if (NZZ > MAXNZ )\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
+            FNPHI=CCJACS(ZZ[I-1],DGPOL,ACOEF,BCOEF,H0VAL,SOLUN);\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
+            ang = BETA *theta;\r
+            realVar = mag * Math.cos(ang);\r
+            imagVar = mag * Math.sin(ang);\r
+            zmlt(realVar,imagVar,KK*FNPHI[0],KK*FNPHI[1],cr,ci);\r
+            GG[I-1][0] = cr[0];\r
+            GG[I-1][1] = ci[0];\r
+        } // for (I=1; I <= NZZ; I++)\r
+\r
+\r
+        // NOW COME THE FACTORS DEPENDENT ON TAU\r
+\r
+        LOWER=-1.0;\r
+        UPPER=TAU[0];\r
+        TAUI=TAU[0];\r
+\r
+        while (true) {\r
+\r
+            MAXRM[0]=0.0;\r
+            HTOL=0.5*TOL;\r
+            for (I=1; I <= NZZ; I++) {\r
+               XI[0] = (2.0*ZZ[I-1][0]+1.0-TAU[0])/(1.0+TAU[0]);\r
+               XI[1] = (2.0*ZZ[I-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
+                ang = 0.5 * theta;\r
+                Z1[0] = mag * Math.cos(ang);\r
+                Z1[1] = mag * Math.sin(ang);\r
+                XI1[0]=XI[0]+Z1[0];\r
+                XI1[1]=XI[1]+Z1[1];\r
+                if (zabs(XI1[0],XI1[1]) < 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
+                ang = theta * (-TUK-1.0);\r
+                realVar = mag * Math.cos(ang);\r
+                imagVar = mag * Math.sin(ang);\r
+                zmlt(XI1[0],XI1[1],XI1[0],XI1[1],cr,ci);\r
+                zmlt(realVar,imagVar,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
+                zmlt(GG[I-1][0],GG[I-1][1],FF[0],FF[1],cr,ci);\r
+                REMND[0] = cr[0];\r
+                REMND[1] = ci[0];\r
+                TERM = zabs(REMND[0],REMND[1]);\r
+                MAXRM[0]=Math.max(MAXRM[0],TERM);\r
+            } // for (I=1; I <= NZZ; 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<TAUI.\r
+\r
+                   if (TAU[0] < TAUI) {\r
+                       LOWER=TAU[0];\r
+                       TAU[0]=0.5*(LOWER+UPPER);\r
+                       continue;\r
+                   } // if (TAU[0] < TAUI)\r
+                   \r
+                } // if (MAXRM[0] < HTOL)\r
+                break;\r
+            } // if (MAXRM[0] < TOL)\r
+            else {\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\r
+        } // while (true)\r
+\r
+        // NORMAL EXIT\r
+\r
+        IER[0]=0;\r
+\r
+    } // private void DEPPJ9 \r
+         \r
+    //COMPLEX FUNCTION CCJACS(X,N,A,B,H,CO)\r
+    private double[] CCJACS(double X[], int N, double A[], double B[],\r
+        double H, double CO[][]) {\r
+        // INTEGER N\r
+        // REAL A(*),B(*),H\r
+        // COMPLEX CO(*),X\r
+\r
+        // ..TO CALCULATE SUMMATION{CO(K+1)*P(K,X)},K=0(1)N, WHERE P(K,X)\r
+        // ..DENOTES THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE K\r
+        // ..EVALUATED AT X, ARRAY CO STORES A GIVEN SET OF COEFFICIENTS,\r
+        // ..ARRAYS A,B STORE THE COEFFICIENTS IN THE THREE-TERM\r
+        // ..RECURRENCE FORMULA FOR THE JACOBI POLYNOMIALS (SEE ASONJ7)\r
+        // ..AND H IS THE SQUARED 2-NORM OF UNITY.\r
+\r
+        double PREV[] = new double[2];\r
+        double CURR[] = new double[2];\r
+        double NEXT[] = new double[2];\r
+        double result[] = new double[2];\r
+        double cr[] = new double[1];\r
+        double ci[] = new double[1];\r
+       // COMPLEX PREV,CURR,NEXT\r
+        int K;\r
+\r
+        if  (N == 0) {\r
+            result[0]=CO[0][0]/Math.sqrt(H);\r
+            result[1]=CO[0][1]/Math.sqrt(H);\r
+        }\r
+        else if (N > 0) {\r
+            PREV[0]=CO[N][0];\r
+            PREV[1]=CO[N][1];\r
+            zmlt(X[0]-B[N-1],X[1],PREV[0],PREV[1],cr,ci);\r
+            CURR[0] = CO[N-1][0] + cr[0]/A[N-1];\r
+            CURR[1] = CO[N-1][1] + ci[0]/A[N-1];\r
+            for (K=N-2; K >= 0; K--) {\r
+               zmlt(X[0]-B[K],X[1],CURR[0],CURR[1],cr,ci);\r
+               NEXT[0] = CO[K][0] + cr[0]/A[K] -A[K]*PREV[0]/A[K+1];\r
+               NEXT[1] = CO[K][1] + ci[0]/A[K] -A[K]*PREV[1]/A[K+1];\r
+                PREV[0]=CURR[0];\r
+                PREV[1]=CURR[1];\r
+                CURR[0]=NEXT[0];\r
+                CURR[1]=NEXT[1];\r
+            } // for (K=N-2; K >= 0; K--)\r
+            result[0]=CURR[0]/Math.sqrt(H);\r
+            result[1]=CURR[1]/Math.sqrt(H);\r
+        } // else if (N > 0)\r
+        else {\r
+               result[0] = 0.0;\r
+               result[1] = 0.0;\r
+        }\r
+        return result;\r
+    } // private double[] CCJACS\r
+\r
+    private void DEPPL8(double BETA, double TAU1[], double TAU2[], boolean T1FXD, int NQUAD,\r
+        int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][], double TOL,\r
+        double MAXRM[], int NINTS, double DELTA,int IER[]) {\r
+       // INTEGER NQUAD,IER,DGPOL,NINTS\r
+       // REAL BETA,TAU1,TAU2,TOL,MAXRM,H0VAL,DELTA\r
+       // REAL ACOEF(*),BCOEF(*)\r
+       // COMPLEX SOLUN(*)\r
+       // LOGICAL T1FXD\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*FNPHI(X)*LOG(ZZ-X)*dX]\r
+       //          TAU1<=X<=TAU2                                       \r
+       \r
+       //     WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY \r
+       //     CORRESPONDENCE DERIVATIVE - JACOBI WEIGHT QUOTIENT. \r
+        \r
+       //     ZZ IS ANY POINT ON A "DELTA-CONTOUR" IN THE UPPER HALF PLANE,\r
+       //     THIS CONTOUR BEING DEFINED BY THE PARAMETER DELTA.  TEST VALUES \r
+       //     FOR ZZ ARE ASSIGNED IN THE BODY OF THE ROUTINE AND THE LOCAL ARRAY\r
+       //     XIVAL IS USED FOR STORING THESE TEST VALUES.\r
+       \r
+       //     THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO \r
+       //     DEFINE FNPHI AND ARE PASSED TO DEPPL9 FOR THIS PURPOSE.\r
+       \r
+       //     MAXRM RECORDS THE MAXIMUM OF THE ABSOLUTE VALUES OF THE REMAINDER\r
+       //     ESTIMATES ASSOCIATED WITH THE DELTA-CONTOUR.\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
+       //                    MAXRM < TOL\r
+       \r
+       //     AND THAT, IF POSSIBLE, \r
+       \r
+       //                   0.5*TOL <= MAXRM < TOL.\r
+       \r
+       //     ON ENTRY NINTS IS THE NUMBER OF INITIAL TEST INTERVALS TO BE USED\r
+       //     ON THE UPPER HALF DELTA-CONTOUR; ON EXIT NINTS GIVES THE FINAL\r
+       //     NUMBEROF INTERVALS REQUIRED FOR CONVERGENCE TO TAU1 OR TAU2.\r
+       \r
+       //     IER=0 - NORMAL EXIT\r
+       //     IER=25- THE LOCAL ARRAY BOUND PARAMETER MNXI NEEDS INCREASING.\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MNXI = 100;\r
+       int NXI;\r
+       double LOMAX[] = new double[1];\r
+       double SS,PI,SS1,SS2,SS3,DP,LEN,TAU1I,TAU2I,SINC,STAR,RXI,DT,HDP;\r
+       double XI[] = new double[2];\r
+       // COMPLEX XI\r
+       boolean FIRST;\r
+       double XIVAL[][] = new double[MNXI][2];\r
+       // COMPLEX XIVAL(MNXI)\r
+       // EXTERNAL DEPPL9\r
+       \r
+       //     INITIALISATION\r
+       \r
+       PI=Math.PI;\r
+       TAU1I=TAU1[0];\r
+       TAU2I=TAU2[0];\r
+       DP=DELTA*PI;\r
+       DT=DELTA*2.0;\r
+       HDP=0.5*DP;\r
+       LEN=2.0+DP;\r
+       SINC=LEN/NINTS;\r
+       STAR=SINC;\r
+       FIRST=true;\r
+       MAXRM[0]=0.0;\r
+       \r
+       // START OF LOOP FOR INTERVAL HALVING\r
+       \r
+       while (true) {\r
+           NXI=0;\r
+           SS1=-STAR+SINC;\r
+           SS2=LEN;\r
+           for (SS = SS1; SS <= SS2; SS += SINC) {\r
+               if (SS < HDP) {\r
+                   SS3=SS/DELTA;\r
+                   XI[0] = 1.0 + DELTA*Math.cos(SS3);\r
+                   XI[1] = DELTA*Math.sin(SS3);\r
+               }\r
+               else if (SS > 2.0+HDP) {\r
+                   SS3=(SS-2.0-HDP)/DELTA;\r
+                   XI[0]= -1.0-DELTA*Math.sin(SS3);\r
+                   XI[1] = DELTA*Math.cos(SS3);\r
+               }\r
+               else {\r
+                   SS3=SS-HDP;\r
+                   XI[0] = 1.0-SS3;\r
+                   XI[1] = DELTA;\r
+               }\r
+               RXI=XI[0];\r
+               if (RXI < (TAU2I+DT) && RXI > (TAU1I-DT)) {\r
+                   NXI=NXI+1;\r
+                   if (NXI > MNXI) {\r
+                       IER[0]=25;\r
+                       return;\r
+                   }\r
+                 XIVAL[NXI-1][0]=XI[0];\r
+                 XIVAL[NXI-1][1]=XI[1];\r
+               } // if (RXI < (TAU2I+DT) && RXI > (TAU1I-DT))\r
+           } // for (SS = SS1; SS <= SS2; SS += SINC)\r
+       \r
+           if (NXI == 0) {\r
+               SINC=STAR;\r
+               STAR=0.5*STAR;\r
+               NINTS=NINTS*2;\r
+               continue;\r
+           }\r
+             \r
+           TAU1[0]=TAU1I;\r
+           TAU2[0]=TAU2I;\r
+           DEPPL9(XIVAL,NXI,BETA,TAU1,TAU2,T1FXD,NQUAD,DGPOL,ACOEF,\r
+            BCOEF,H0VAL,SOLUN,TOL,LOMAX,FIRST,IER); \r
+           if (IER[0] > 0) {\r
+               return;\r
+           }\r
+       \r
+           MAXRM[0]=Math.max(MAXRM[0],LOMAX[0]);\r
+           if (T1FXD && TAU2[0] != TAU2I) {\r
+               TAU2I=TAU2[0];\r
+               SINC=STAR;\r
+               STAR=0.5*STAR;\r
+               NINTS=NINTS*2;\r
+               continue;\r
+           }\r
+           else if (!T1FXD && TAU1[0] != TAU1I) {\r
+               TAU1I=TAU1[0];\r
+               SINC=STAR;\r
+               STAR=0.5*STAR;\r
+               NINTS=NINTS*2;\r
+               continue;\r
+           }\r
+           break;\r
+       } // while (true)\r
+       IER[0]=0;\r
+       \r
+    } // private void DEPPL8\r
+    \r
+    private void DEPPL9(double ZZ[][], int NZZ, double BETA, double TAU1[], double TAU2[], boolean T1FXD,\r
+        int NQUAD, int DGPOL, double ACOEF[], double BCOEF[], double H0VAL, double SOLUN[][], double TOL,\r
+        double MAXRM[], boolean FIRST, int IER[]) {\r
+       // INTEGER NQUAD,DGPOL,IER,NZZ\r
+       // REAL BETA,TAU1,TAU2,H0VAL,TOL,MAXRM\r
+       // REAL ACOEF(*),BCOEF(*)\r
+       // LOGICAL T1FXD,FIRST\r
+       // COMPLEX SOLUN(*),ZZ(*)\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*FNPHI(X)*LOG(ZZ(I)-X)*dX], I=1,NZZ\r
+       //     TAU1<=X<=TAU2                                     \r
+       \r
+       //     WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY \r
+       //     CORRESPONDENCE DERIVATIVE - JACOBI WEIGHT QUOTIENT AND ZZ IS\r
+       //     A GIVEN ARRAY OF POINTS.  \r
+       \r
+       //     THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO \r
+       //     DEFINE FNPHI.\r
+       \r
+       //     THE MAXIMUM ABSOLUTE VALUE OF ALL THE REMAINDERS CORRESPONDING TO \r
+       //     ZZ(I) , I=1,NZZ, IS STORED IN MAXRM. \r
+       \r
+       //     THIS ROUTINE USES THE SIMPLEST POSSIBLE ESTIMATES; I.E. THE  \r
+       //     LEADING TERM ONLY IN THE ASYMPTOTIC EXPANSION AND THE WATSON- \r
+       //     DOETSCH ESTIMATE FOR ANY 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
+       //                    MAXRM < TOL\r
+       \r
+       //     AND THAT, IF POSSIBLE, \r
+       \r
+       //                   0.5*TOL <= MAXRM < TOL.\r
+       \r
+       //     IER=0 - NORMAL EXIT\r
+       //     IER=26- TOO MANY TEST POINTS ON DELTA CONTOUR; INCREASE \r
+       //             PARAMETER MAXNZ BELOW\r
+       //     IER=14- BETA MAY CAUSE OVERFLOW IN GAMMA FUNCTION; AN ANGLE\r
+       //             ON THE BOUNDARY IS TOO SMALL\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MAXNZ = 200;\r
+       int I;\r
+       double KK,RI,TURI,RN,TUK,LOWER,HTOL,UPPER,HCO,PI,FORK,RR,\r
+            MEAN,BB,BB1,FFH,TERM,TAU1I,TAU2I,OFLOW;\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 FNPHI[] = new double[2];\r
+       double FPHI1[] = new double[2];\r
+       double REMND[] = new double[2];\r
+       // COMPLEX XI,Z1,XI1,FFG,FNPHI,FPHI1,REMND\r
+       // COMPLEX CCJACS\r
+       double GG[][] = new double[MAXNZ][2];\r
+       double HH[][] = new double[MAXNZ][2];\r
+       // COMPLEX GG(MAXNZ),HH(MAXNZ)\r
+       // EXTERNAL GAMMA,CCJACS,LGGAM\r
+       \r
+       if (NZZ > MAXNZ) {\r
+           IER[0]=26;\r
+           return;\r
+       }\r
+       \r
+       if (BETA > 24.0) {\r
+       \r
+           // TEST FOR POSSIBLE OVERFLOW IN GAMMA FUNCTION\r
+       \r
+           OFLOW=Math.log(Double.MAX_VALUE);\r
+           KK=LGGAM(BETA+1.0);\r
+           if (KK > OFLOW) {\r
+               IER[0]=14;\r
+               return;\r
+           }\r
+       } // if (BETA > 24.0)\r
+       \r
+       // FIRST WE COMPUTE THE FACTORS WHICH ARE INDEPENDENT OF TAU1,TAU2\r
+       \r
+       TAU1I=TAU1[0];\r
+       TAU2I=TAU2[0];\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
+       HCO=Math.sin(PI*BETA)*GAMMA(BETA+1.0)/PI/Math.pow(FORK,BETA);\r
+       /*      FPHI1=CCJACS((-1E+0,0E+0),DGPOL,ACOEF,BCOEF,H0VAL,SOLUN)\r
+             DO 125 I=1,NZZ\r
+               FNPHI=CCJACS(ZZ(I),DGPOL,ACOEF,BCOEF,H0VAL,SOLUN)\r
+               GG(I)=-(1E+0+ZZ(I))**BETA*KK*FNPHI\r
+               HH(I)=-CLOG(1E+0+ZZ(I))*HCO*KK*FPHI1\r
+       125   CONTINUE\r
+              \r
        C\r
-       C     START OF LOOP FOR INTERVAL HALVING\r
+       C     NOW COME THE FACTORS DEPENDENT ON TAU1 AND TAU2.\r
        C\r
-       10    CONTINUE\r
-             NXI=0\r
-             SS1=-STAR+SINC\r
-             SS2=LEN\r
-             DO 30 SS=SS1,SS2,SINC\r
-               IF (SS .LT. HDP) THEN\r
-                   SS3=SS/DELTA\r
-                   XI=1E+0+DELTA*CMPLX(COS(SS3),SIN(SS3))\r
-               ELSE IF (SS .GT. 2E+0+HDP) THEN\r
-                   SS3=(SS-2E+0-HDP)/DELTA\r
-                   XI=-1E+0+DELTA*(0E+0,1E+0)*CMPLX(COS(SS3),SIN(SS3))\r
-               ELSE\r
-                   SS3=SS-HDP\r
-                   XI=CMPLX(1E+0-SS3,DELTA)\r
+             LOWER=TAU1\r
+             UPPER=TAU2\r
+       C\r
+       250   CONTINUE\r
+       C\r
+             HTOL=5E-1*TOL\r
+             RR=(TAU2-TAU1)*5E-1\r
+             MEAN=(TAU1+TAU2)*5E-1\r
+             BB=(1E+0+MEAN)/RR\r
+             BB1=BB+SQRT(BB*BB-1E+0)\r
+             FFH=(RR*(BB1-1E+0/BB1))**(BETA+1E+0)/BB1**TUK\r
+             MAXRM=0E+0\r
+             DO 325 I=1,NZZ\r
+               XI=(ZZ(I)-MEAN)/RR\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
-               RXI=REAL(XI)\r
-               IF (RXI .LT. (TAUI+DT)) THEN\r
-                 NXI=NXI+1\r
-                 IF (NXI .GT. MNXI) THEN\r
-                   IER=25\r
-                   RETURN\r
+               FFG=XI1**(-TUK-1E+0)*(XI1*XI1-1E+0)*RR\r
+               REMND=GG(I)*FFG+HH(I)*FFH\r
+               TERM=ABS(REMND)\r
+               MAXRM=MAX(MAXRM,TERM)\r
+       325   CONTINUE\r
+       C\r
+             IF (MAXRM .LT. TOL) THEN\r
+       C\r
+       C       ACCURACY IS ACHIEVED, BUT MAYBE TAU2 COULD BE INCREASED OR\r
+       C       TAU1 DECREASED\r
+       C\r
+               IF (MAXRM .LT. HTOL) THEN\r
+       C\r
+       C         TAU2 NEEDS INCREASING IF T1FXD (BUT THIS IS ONLY POSSIBLE IF \r
+       C         TAU2<TAU2I) OR TAU1 NEED DECREASING OTHERWISE (BUT THIS IS \r
+       C         ONLY POSSIBLE IF TAU1>TAU1I)\r
+       C\r
+                 IF (T1FXD .AND. TAU2 .LT. TAU2I) THEN\r
+                   LOWER=TAU2\r
+                   TAU2=5E-1*(LOWER+UPPER)\r
+                   GOTO 250\r
+                 ELSE IF (.NOT.T1FXD .AND. TAU1 .GT. TAU1I) THEN\r
+                   UPPER=TAU1\r
+                   TAU1=5E-1*(LOWER+UPPER)\r
+                   GOTO 250\r
                  ENDIF\r
-                 XIVAL(NXI)=XI\r
                ENDIF\r
-       30    CONTINUE\r
+             ELSE\r
        C\r
-             IF (NXI .EQ. 0) THEN\r
-               SINC=STAR\r
-               STAR=5E-1*STAR\r
-               NINTS=NINTS*2\r
-               GOTO 10\r
-             ENDIF\r
-       C      \r
-             TAU=TAUI\r
-             CALL DEPPJ9(XIVAL,NXI,BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL,\r
-            +            SOLUN(1),TOL,MAXRM,IER) \r
-             IF (IER .GT. 0) THEN\r
-               RETURN\r
-             ENDIF\r
+       C       ACCURACY NOT ACHIEVED AND TAU2 NEEDS DECREASING OR TAU1 NEEDS\r
+       C       INCREASING.\r
        C\r
-             IF (TAU .NE. TAUI) THEN\r
-               TAUI=TAU\r
-               SINC=STAR\r
-               STAR=5E-1*STAR\r
-               NINTS=NINTS*2\r
-               GOTO 10\r
+               IF (FIRST) THEN\r
+                 TOL=HTOL\r
+                 FIRST=.FALSE.\r
+               ENDIF\r
+               IF (T1FXD) THEN\r
+                 UPPER=TAU2\r
+                 TAU2=5E-1*(LOWER+UPPER)\r
+               ELSE\r
+                 LOWER=TAU1\r
+                 TAU1=5E-1*(LOWER+UPPER)\r
+               ENDIF\r
+               GOTO 250\r
              ENDIF\r
        C\r
-             IER=0\r
+       C     NORMAL EXIT\r
        C\r
-             END*/\r
-    } // private void DEPPJ8\r
+             IER=0\r
+       C */\r
+    } // private void DEPPL9\r
+\r
 \r
 \r
       /**\r