Porting continuues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 20 Nov 2017 23:12:47 +0000 (23:12 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 20 Nov 2017 23:12:47 +0000 (23:12 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15267 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 5c371d1..116aebb 100644 (file)
@@ -36,11 +36,44 @@ public class DoublyConnectedSC extends AlgorithmBase {
        private final int EXAMPLE_GIVEN_FOR_CHECKING_INVERSE_MAP = 6;\r
        private final int UPPER_HALF_PLANE_WITH_HORIZONTAL_SLIT = 7;\r
        \r
+       private final int dscfun = 1;\r
+       \r
+       // Below are the 5 variables in the common block PARAM1\r
+       private double W02[][] = new double[30][2];\r
+       private double W12[][] = new double[30][2];\r
+       private double Z02[][] = new double[30][2];\r
+       private double Z12[][] = new double[30][2];\r
+       private double C2[] = new double[2];\r
+       \r
+       // Below are the 6 variables in the common block PARAM2\r
+       private double U2[] = new double[1];\r
+       private double PHI02[] = new double[30];\r
+       private double PHI12[] = new double[30];\r
+       private double ALFA02[] = new double[30];\r
+       private double ALFA12[] = new double[30];\r
+       private double QWORK2[] = new double[1660];\r
+       \r
+       // Below are the 7 variables in the common block PARAM3\r
+       private int M2, N2, NPTQ2, ISHAPE2, LINEARC2,  NSHAPE;\r
+       private int IND[] = new int[20];\r
+       \r
        // Below are the 4 variables in the common block PARAM4\r
-       private double DLAM;\r
-       private int IU;\r
        private double UARY[] = new double[8];\r
        private double VARY[] = new double[3];\r
+       private double DLAM;\r
+       private int IU;\r
+       \r
+       // Below are the 2 variables in the common block PARAM5:\r
+       // Screen display of the residual of the system as the iteration goes on, 1 for "YES", 2 for "NO"\r
+       private int ISPRT;\r
+       private int ICOUNT;\r
+       \r
+       // Common blocks ..\r
+    // COMMON /PARAM1/W02,W12,Z02,Z12,C2\r
+    // COMMON /PARAM2/U2,PHI02,PHI12,ALFA02,ALFA12,QWORK2\r
+    // COMMON /PARAM3/M2,N2,NPTQ2,ISHAPE2,LINEARC2,NSHAPE,IND\r
+    // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+    // COMMON /PARAM5/ISPRT,ICOUNT\r
        \r
        private SchwarzChristoffelMapping scm = new SchwarzChristoffelMapping();\r
        \r
@@ -58,10 +91,11 @@ public class DoublyConnectedSC extends AlgorithmBase {
                \r
        }\r
        \r
-       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV) {\r
+       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV, int ISPRT) {\r
                this.IPOLY = IPOLY;\r
                this.NPTQ = NPTQ;\r
                this.ISOLV = ISOLV;\r
+               this.ISPRT = ISPRT;\r
                testRoutine = true;\r
        }\r
        \r
@@ -450,6 +484,168 @@ double neweps;
 \r
        }\r
        \r
+       private double[] WTHETA(double U, double W[]) {\r
+           //  -------------------------\r
+           //    EVALUATES THETA-FUNCTION AT W,WHERE U IS THE\r
+           //    INNER RADIUS OF THE ANNULUS.THE DEFINITION OF\r
+           //    THETA-FUNCTION CAN BE FOUND IN REFERENCE 3.\r
+       \r
+           //     .. Scalar Arguments ..\r
+           //  DOUBLE COMPLEX W\r
+\r
+           //     .. Scalars in Common ..\r
+           // DOUBLE PRECISION DLAM\r
+           // INTEGER IU\r
+\r
+           //    .. Arrays in Common ..\r
+           // DOUBLE PRECISION UARY(8),VARY(3)\r
+\r
+           //     .. Local Scalars ..\r
+               double result[] = new double[2];\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+               double WT[] = new double[2];\r
+               double WWN[] = new double[2];\r
+               double WWN0[] = new double[2];\r
+               double WWP[] = new double[2];\r
+               double WWP0[] = new double[2];\r
+               double ZI[] = new double[2];\r
+           // DOUBLE COMPLEX WT,WWN,WWN0,WWP,WWP0,ZI\r
+           double PI;\r
+           int K;\r
+       \r
+           //     .. Common blocks ..\r
+           // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+       \r
+           WWP[0] = 1.0;\r
+           WWP[1] = 0.0;\r
+           WWN[0] = 1.0;\r
+           WWN[1] = 0.0;\r
+           result[0] = 1.0;\r
+           result[1] = 0.0;\r
+           if (U < 0.63) {\r
+               WWP0[0] = -W[0];\r
+               WWP0[1] = -W[1];\r
+               scm.zdiv(-1.0, 0.0, W[0], W[1], cr, ci);\r
+               WWN0[0] = cr[0];\r
+               WWN0[1] = ci[0];\r
+               for (K = 1; K <= IU; K++) {\r
+                       scm.zmlt(WWP[0], WWP[1], WWP0[0], WWP0[1], cr, ci);\r
+                       WWP[0] = cr[0];\r
+                       WWP[1] = ci[0];\r
+                   scm.zmlt(WWN[0], WWN[1], WWN0[0], WWN0[1], cr, ci);\r
+                   WWN[0] = cr[0];\r
+                   WWN[1] = ci[0];\r
+                   result[0] = result[0] + UARY[K-1]* (WWP[0]+WWN[0]);\r
+                   result[1] = result[1] + UARY[K-1]* (WWP[1]+WWN[1]);\r
+               } // for (K = 1; K <= IU; K++)\r
+               return result;\r
+          } // if (U < 0.63)\r
+\r
+          PI = Math.PI;\r
+          ZI[0] = 0.0;\r
+          ZI[1] = 1.0;\r
+          double logr = Math.log(scm.zabs(-W[0], -W[1]));\r
+          double logi = Math.atan2(-W[1], -W[0]);\r
+          scm.zmlt(-ZI[0], -ZI[1], logr, logi, cr, ci);\r
+          WT[0] = cr[0];\r
+          WT[1] = ci[0];\r
+          if (U < 0.94) {\r
+                 double base = Math.exp(WT[0]/DLAM);\r
+                 double arg = WT[1]/DLAM;\r
+                 WWP[0] = base * Math.cos(arg);\r
+                 WWP[1] = base * Math.sin(arg);\r
+                 scm.zdiv(1.0, 0.0, WWP[0], WWP[1], cr, ci);\r
+             WWN[0] = cr[0];\r
+             WWN[1] = ci[0];\r
+             result[0] = result[0] + VARY[0]* (WWP[0]+WWN[0]); \r
+             result[1] = result[1] + VARY[0]* (WWP[1]+WWN[1]);   \r
+          } // if (U < 0.94)\r
+          scm.zmlt(WT[0], WT[1], WT[0], WT[1], cr, ci);\r
+          double argr = -cr[0]/(4.0 * PI * DLAM);\r
+          double argi = -ci[0]/(4.0 * PI * DLAM);\r
+          double base = Math.exp(argr);\r
+          double expr = base * Math.cos(argi);\r
+          double expi = base * Math.sin(argi);\r
+          scm.zmlt(expr, expi, result[0], result[1], cr, ci);\r
+          result[0] = cr[0]/Math.sqrt(DLAM);\r
+          result[1] = ci[0]/Math.sqrt(DLAM);\r
+          return result;\r
+       }\r
+\r
+\r
+       \r
+       private double[] WPROD(double W[], int M, int N, double U,\r
+                       double W0[][], double W1[][], double ALFA0[], double ALFA1[]) {\r
+       //   --------------------------------------------\r
+       //    COMPUTES THE PRODUCT (D-SC INTEGRAND):\r
+       // M                                   N\r
+       //PROD THETA(W/U*W0(K))**(ALFA0(K)-1)*PROD THETA(U*W/W1(K))**(ALFA1(K)-1)\r
+       // K=1                                 K=1\r
+       //    THE CALLING SEQUENCE IS EXPLAINED IN THE ABOVE FORMULA.\r
+       \r
+       //     .. Scalar Arguments ..\r
+       //      DOUBLE COMPLEX W\r
+            \r
+       //     .. Array Arguments ..\r
+       //      DOUBLE COMPLEX W0(M),W1(N)\r
+       //      DOUBLE PRECISION ALFA0(M),ALFA1(N)\r
+       \r
+       //     .. Scalars in Common ..\r
+       //      DOUBLE PRECISION DLAM\r
+       //      INTEGER IU\r
+       \r
+       //     .. Arrays in Common ..\r
+       //      DOUBLE PRECISION UARY(8),VARY(3)\r
+       \r
+       //     .. Local Scalars ..\r
+               double result[] = new double[2];\r
+               double WSUM[] = new double[2];\r
+               double WTH[] = new double[2];\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+               double win[] = new double[2];\r
+               double wout[];\r
+               double base;\r
+       //      DOUBLE COMPLEX WSUM,WTH\r
+             int K;\r
+\r
+       //     .. External Functions ..\r
+       //      DOUBLE COMPLEX WTHETA\r
+       //      EXTERNAL WTHETA\r
+       \r
+       //     .. Common blocks ..\r
+       //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+             WSUM[0] = 0.0;\r
+             WSUM[1] = 0.0;\r
+             for (K = 1; K <= M; K++) {\r
+                 scm.zdiv(W[0], W[1], U*W0[K][0], U*W0[K][1], cr, ci);\r
+                 win[0] = cr[0];\r
+                 win[1] = ci[0];\r
+                 wout = WTHETA(U, win);\r
+                 WTH[0] = Math.log(scm.zabs(wout[0], wout[1]));\r
+                 WTH[1] = Math.atan2(wout[1], wout[0]);\r
+                 WSUM[0] = WSUM[0] + (ALFA0[K-1]-1.0)*WTH[0];\r
+                 WSUM[1] = WSUM[1] + (ALFA0[K-1]-1.0)*WTH[1];\r
+             } // for (K = 1; K <= M; K++)\r
+             for (K = 1; K <= N; K++) {\r
+                 scm.zdiv(U*W[0], U*W[1], W1[K][0], W1[K][1], cr, ci);\r
+                 win[0] = cr[0];\r
+                 win[1] = ci[0];\r
+                 wout = WTHETA(U, win);\r
+                 WTH[0] = Math.log(scm.zabs(wout[0], wout[1]));\r
+                 WTH[1] = Math.atan2(wout[1], wout[0]);\r
+                 WSUM[0] = WSUM[0] + (ALFA1[K-1]-1.0)*WTH[0];\r
+                 WSUM[1] = WSUM[1] + (ALFA1[K-1]-1.0)*WTH[1];\r
+             } // for (K = 1; K <= N; K++)\r
+             base = Math.exp(WSUM[0]);\r
+             result[0] = base * Math.cos(WSUM[1]);\r
+             result[1] = base * Math.sin(WSUM[1]);\r
+             return result;\r
+       }\r
+\r
+\r
+       \r
        private double[] WQSUM(double WA[], double PHIA,int KWA, int IC, double WB[],\r
                        double PHIB, double RADIUS,int M,int N, double U, double W0[][], double W1[][],\r
                        double ALFA0[], double ALFA1[], int NPTQ,double QWORK[], int LINEARC) {\r
@@ -475,53 +671,76 @@ double neweps;
                \r
                //     .. Local Scalars ..\r
                double result[] = new double[2];\r
-               /*      DOUBLE COMPLEX W,WC,WH,ZI\r
-                     DOUBLE PRECISION PWC,PWH\r
-                     INTEGER I,IOFFST,IWT1,IWT2\r
-               C     ..\r
-               C     .. External Functions ..\r
-                     DOUBLE COMPLEX WPROD\r
-                     EXTERNAL WPROD\r
-               C     ..\r
-               C     .. Intrinsic Functions ..\r
-                     INTRINSIC EXP\r
-               C     ..\r
+               double W[] = new double[2];\r
+               double WC[] = new double[2];\r
+               double WH[] = new double[2];\r
+               double ZI[] = new double[2];\r
+               //     DOUBLE COMPLEX W,WC,WH,ZI\r
+               double PWC,PWH;\r
+               int I,IOFFST,IWT1,IWT2;\r
+               double prod[] = null;\r
+               double base;\r
+               double arg;\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+\r
+               //     .. External Functions ..\r
+               //      DOUBLE COMPLEX WPROD\r
+               //      EXTERNAL WPROD\r
+               \r
                //     .. Common blocks ..\r
                //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
-               C     ..\r
-                     WQSUM = (0.D0,0.D0)\r
-               C\r
-               C   INDEX ARRANGEMENT:\r
-                     IWT1 = NPTQ* (IC*M+KWA-1) + 1\r
-                     IF (KWA.EQ.0) IWT1 = NPTQ* (M+N) + 1\r
-                     IWT2 = IWT1 + NPTQ - 1\r
-                     IOFFST = NPTQ* (M+N+1)\r
-               C\r
-               C   COMPUTE GAUSS-JACOBI SUM(W(J)*PROD(X(J))):\r
-                     IF (LINEARC.EQ.1) GO TO 20\r
-               C\r
-               C   INTEGRATE ALONG A LINE SEGMENT:\r
-                     WH = (WB-WA)/2.D0\r
-                     WC = (WA+WB)/2.D0\r
-                     DO 10 I = IWT1,IWT2\r
-                         W = WC + WH*QWORK(I)\r
-                         WQSUM = WQSUM + QWORK(IOFFST+I)*\r
-                    +            WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1)\r
-                  10 CONTINUE\r
-                     WQSUM = WQSUM*WH\r
-                     RETURN\r
-               C\r
-               C   INTEGRATE ALONG A CIRCULAR ARC:\r
-                  20 ZI = (0.D0,1.D0)\r
-                     PWH = (PHIB-PHIA)/2.D0\r
-                     PWC = (PHIB+PHIA)/2.D0\r
-                     DO 30 I = IWT1,IWT2\r
-                         W = RADIUS*EXP(ZI* (PWC+PWH*QWORK(I)))\r
-                         WQSUM = WQSUM + QWORK(IOFFST+I)*W*\r
-                    +            WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1)\r
-                  30 CONTINUE\r
-                     WQSUM = WQSUM*PWH*ZI*/\r
-                  return result;\r
+           result[0] = 0.0;\r
+               result[1] = 0.0;\r
+               \r
+               //   INDEX ARRANGEMENT:\r
+               IWT1 = NPTQ* (IC*M+KWA-1) + 1;\r
+               if (KWA == 0) {\r
+                   IWT1 = NPTQ* (M+N) + 1;\r
+               }\r
+           IWT2 = IWT1 + NPTQ - 1;\r
+               IOFFST = NPTQ* (M+N+1);\r
+               \r
+               //   COMPUTE GAUSS-JACOBI SUM(W(J)*PROD(X(J))):\r
+               if (LINEARC != 1) {\r
+               \r
+                   //   INTEGRATE ALONG A LINE SEGMENT:\r
+                   WH[0] = (WB[0]-WA[0])/2.0;\r
+                   WH[1] = (WB[1]-WA[1])/2.0;\r
+                   WC[0] = (WA[0]+WB[0])/2.0;\r
+                   WC[1] = (WA[1]+WB[1])/2.0;\r
+                   for (I = IWT1; I <= IWT2; I++) {\r
+                         W[0] = WC[0] + WH[0]*QWORK[I-1];\r
+                         W[1] = WC[1] + WH[1]*QWORK[I-1];\r
+                         prod = WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1);\r
+                         result[0] = result[0] + QWORK[IOFFST + I-1]*prod[0];\r
+                         result[1] = result[1] + QWORK[IOFFST + I-1]*prod[1];\r
+                   } // for (I = IWT1; I <= IWT2; I++)\r
+                   scm.zmlt(result[0], result[1], WH[0], WH[1], cr, ci);\r
+                   result[0] = cr[0];\r
+                   result[1] = ci[0];\r
+                   return result;\r
+           } if (LINEARC != 1)\r
+               \r
+               //   INTEGRATE ALONG A CIRCULAR ARC:\r
+           ZI[0] = 0.0;\r
+           ZI[1] = 1.0;\r
+               PWH = (PHIB-PHIA)/2.0;\r
+               PWC = (PHIB+PHIA)/2.0;\r
+               for (I = IWT1; I <= IWT2; I++) {\r
+                       base = RADIUS*Math.exp(ZI[0]* (PWC+PWH*QWORK[I-1]));\r
+                       arg = ZI[1]* (PWC+PWH*QWORK[I-1]);\r
+                       W[0] = base * Math.cos(arg);\r
+                       W[1] = base * Math.sin(arg);\r
+                   prod = WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1);\r
+                   scm.zmlt(W[0], W[1], prod[0], prod[1], cr, ci);\r
+                   result[0] = result[0] + QWORK[IOFFST+I-1] * cr[0];\r
+                   result[1] = result[1] + QWORK[IOFFST+I-1] * ci[0];\r
+               } // for (I = IWT1; I <= IWT2; I++)\r
+               scm.zmlt(result[0], result[1], PWH*ZI[0], PWH*ZI[1], cr, ci);\r
+               result[0] = cr[0];\r
+               result[1] = ci[0];\r
+               return result;\r
        }\r
 \r
        \r
@@ -801,6 +1020,204 @@ double neweps;
              } // for (I = 1; I <= M - 2; I++)\r
                  return;\r
        }\r
+       \r
+       private void DSCFUN(int NDIM, double X[], double FVAL[],int IFLAG[]) {\r
+           //    --------------------------------------\r
+           //  FORMS THE NONLINEAR SYSTEM SATISFIED BY D-SC PARAMETERS.THE\r
+           //  SUBROUTINE WILL BE CALLED BY NONLINEAR SYSTEM SOLVER HYBRD.\r
+           //  SEE ROUTINE DSCSOLV FOR THE VARIABLES IN THE COMMON BLOCKS.\r
+           //  NDIM: THE DIMENSION OF THE SYSTEM.\r
+           //  X:    UNKNOWNS\r
+           //  FVAL: THE VECTOR DEFINED BY THE SYSTEM.\r
+           //  IFLAG:(ON OUTPUT)=1,THE ITERATION WAS SUCCESSFULLY COMPLETED.\r
+           //         =2,3,OR 4, UNSUCCESSFUL TERMINATION OF THE ITERATION.\r
+           //         =5 MAY INDICATE THE TOLERANCE GIVEN IN THE CALLING\r
+           //         SEQUENCE OF HYBRD IS TOO SMALL.\r
+       \r
+       \r
+           //  TRANSFORM UNCONSTRAINED X(K) TO ACTUAL D-SC PARAMETERS\r
+           //  AND COMPUTE DATA TO BE USED IN WTHETA:\r
+       \r
+           //     .. Scalar Arguments ..\r
+           //  INTEGER IFLAG,NDIM\r
+           //     ..\r
+           //     .. Array Arguments ..\r
+           //  DOUBLE PRECISION FVAL(NDIM),X(NDIM)\r
+       \r
+           \r
+            //.. Local Scalars ..\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+               double WA[] = new double[2];\r
+               double WAI[] = new double[2];\r
+               double WARC[] = new double[2];\r
+               double WB[] = new double[2];\r
+               double WBI[] = new double[2];\r
+               double WCIRCLE[] = new double[2];\r
+               double WIN1[] = new double[2];\r
+               double WIN2[] = new double[2];\r
+               double WIN3[] = new double[2];\r
+               double WIN4[] = new double[2];\r
+               double WINT1[] = new double[2];\r
+               double WINT2[] = new double[2];\r
+               double WINT3[] = new double[2];\r
+               double WLINE[] = new double[2];\r
+               double WLINE1[] = new double[2];\r
+               double WLINE2[] = new double[2];\r
+               double WX[] = new double[2];\r
+               double ZI[] = new double[2];\r
+           //  DOUBLE COMPLEX WA,WAI,WARC,WB,WBI,WCIRCLE,WIN1,WIN2,WIN3,WIN4,\r
+           //               WINT1,WINT2,WINT3,WLINE,WLINE1,WLINE2,WX,ZI\r
+           double FMAXN,PHIA,PHIB,RADIUS,TEST1;\r
+           int I,J,K;\r
+       \r
+           //     .. External Functions ..\r
+           // DOUBLE COMPLEX WQUAD\r
+           // DOUBLE PRECISION FMAX\r
+           // EXTERNAL WQUAD,FMAX\r
+       \r
+           //     .. External Subroutines ..\r
+           // EXTERNAL THDATA,XWTRAN\r
+       \r
+           //     .. Common blocks ..\r
+           // COMMON /PARAM1/W02,W12,Z02,Z12,C2\r
+           // COMMON /PARAM2/U2,PHI02,PHI12,ALFA02,ALFA12,QWORK2\r
+           // COMMON /PARAM3/M2,N2,NPTQ2,ISHAPE2,LINEARC2,NSHAPE,IND\r
+           // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+           // COMMON /PARAM5/ISPRT,ICOUNT\r
+\r
+             XWTRAN(M2,N2,X,U2,C2,W02,W12,PHI02,PHI12);\r
+             THDATA(U2);\r
+             ZI[0] = 0.0;\r
+             ZI[1] = 1.0;\r
+             ICOUNT = ICOUNT + 1;\r
+       \r
+             //  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE INNER POLYGON:\r
+             double wout[] = WQUAD(W12[N2-1],PHI12[N2-1],N2,1,W12[0],PHI12[0],1,\r
+                            1,U2[0],M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
+             scm.zmlt(C2[0], C2[1], wout[0], wout[1], cr, ci);\r
+             WIN1[0] = Z12[0][0] - Z12[N2-1][0] - cr[0];\r
+             WIN1[1] = Z12[0][1] - Z12[N2-1][1] - ci[0];\r
+             scm.zmlt(ZI[0], ZI[1], WIN1[0], WIN1[1], cr, ci);\r
+             FVAL[0] = ci[0];\r
+             FVAL[1] = WIN1[1];\r
+       \r
+             //  N-1 SIDE LENGTH CONDITIONS FOR THE INNER POLYGON:\r
+             for (I = 1; I <= N2 - 1; I++) {\r
+                 WINT1 = WQUAD(W12[I-1],PHI12[I-1],I,1,W12[I],PHI12[I],I+1,1,U2[0],M2,N2,\r
+                        U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
+                 scm.zmlt(C2[0], C2[1], WINT1[0], WINT1[1], cr, ci);\r
+                 FVAL[I+1] = scm.zabs(Z12[I][0]-Z12[I-1][0],Z12[I][1]-Z12[I-1][1]) - scm.zabs(cr[0], ci[0]);\r
+             } // for (I = 1; I <= N2 - 1; I++)\r
+       \r
+             //  TWO EQUATIONS TO FIX THE RELATIVE POSITION OF THE INNER POLYGON:\r
+             TEST1 = Math.cos(PHI12[N2-1]);\r
+            /* if (TEST1 < U2[0]) {\r
+       \r
+                 // IF THE LINE PATH FROM W02[M2-1] TO W12[N2-1] IS OUT OF DOMAIN,THE\r
+              //  COMBINATION OF TWO DIFFERENT PATHS WILL BE USED INSTEAD:\r
+             WX = DCMPLX(U2,0.D0)\r
+             WLINE = WQUAD(W02(M2),0.D0,M2,0,WX,0.D0,0,2,0.D0,M2,N2,U2,W02,W12,ALFA02,\r
+            +        ALFA12,NPTQ2,QWORK2,0,2)\r
+             IF (PHI12(N2).LE.0.D0) THEN\r
+                 WARC = WQUAD(W12(N2),PHI12(N2),N2,1,WX,0.D0,0,2,U2,M2,N2,U2,W02,W12,\r
+            +           ALFA02,ALFA12,NPTQ2,QWORK2,1,2)\r
+                 WIN2 = WLINE - WARC\r
+\r
+             ELSE\r
+                 WARC = WQUAD(WX,0.D0,0,2,W12(N2),PHI12(N2),N2,1,U2,M2,N2,U2,W02,W12,\r
+            +           ALFA02,ALFA12,NPTQ2,QWORK2,1,2)\r
+                 WIN2 = WLINE + WARC\r
+             END IF\r
+\r
+             GO TO 30\r
+             } // if (TEST1 < U2[0])\r
+             else {\r
+\r
+                 WIN2 = WQUAD(W02(M2),0.D0,M2,0,W12(N2),0.D0,N2,1,0.D0,M2,N2,U2,W02,W12,ALFA02,\r
+            +       ALFA12,NPTQ2,QWORK2,0,2)\r
+             }\r
+             FVAL(N2+2) = DIMAG(ZI* (Z12(N2)-Z02(M2)-C2*WIN2))\r
+             FVAL(N2+3) = DIMAG(Z12(N2)-Z02(M2)-C2*WIN2)\r
+       C\r
+       C  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE OUTER POLYGON:\r
+             WIN3 = Z02(1) - Z02(M2) - C2*WQUAD(W02(M2),PHI02(M2),M2,0,W02(1),PHI02(1),1,\r
+            +       0,1.D0,M2,N2,U2,W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2)\r
+             FVAL(N2+4) = DIMAG(ZI*WIN3)\r
+             FVAL(N2+5) = DIMAG(WIN3)\r
+       C\r
+             IF (M2.EQ.3) THEN\r
+       C\r
+       C  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
+                 IF (ISPRT.NE.1) GO TO 40\r
+                 FMAXN = FMAX(NDIM,FVAL)\r
+                 WRITE (6,FMT=9000) ICOUNT,FMAXN\r
+          40     CONTINUE\r
+                 RETURN\r
+\r
+             END IF\r
+\r
+             IF (ISHAPE2.EQ.1) GO TO 70\r
+       C\r
+       C  M-3 SIDE LENGTH CONDITIONS OF THE OUTER POLYGON:\r
+             DO 50 J = 1,M2 - 3\r
+                 WINT2 = WQUAD(W02(J),PHI02(J),J,0,W02(J+1),PHI02(J+1),J+1,0,1.D0,\r
+            +            M2,N2,U2,W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2)\r
+                 FVAL(N2+5+J) = ABS(Z02(J+1)-Z02(J)) - ABS(C2*WINT2)\r
+          50 CONTINUE\r
+       C\r
+       C  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
+             IF (ISPRT.NE.1) GO TO 60\r
+             FMAXN = FMAX(NDIM,FVAL)\r
+             WRITE (6,FMT=9010) ICOUNT,FMAXN\r
+          60 CONTINUE\r
+             RETURN\r
+       C\r
+       C  OUTER POLYGON CONTAINS SOME INFINITE VERTICES & FOR EACH OF THEM\r
+       C  TWO LENGTH CONDITIONS WILL BE REPLACED BY A COMPLEX INTEGRAL:\r
+          70 DO 100 K = 1,NSHAPE - 1\r
+                 IF (IND(K+1).EQ.2 .OR. IND(K).GE.IND(K+1)-2) GO TO 90\r
+                 DO 80 J = IND(K) + 1,IND(K+1) - 2\r
+                     WINT3 = WQUAD(W02(J),PHI02(J),J,0,W02(J+1),PHI02(J+1),J+1,0,\r
+            +                1.D0,M2,N2,U2,W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2)\r
+                     FVAL(N+5+J) = ABS(Z02(J+1)-Z02(J)) - ABS(C2*WINT3)\r
+          80     CONTINUE\r
+          90     IF (K.EQ.NSHAPE-1 .OR. IND(K+1).EQ.M2-1) GO TO 100\r
+       C\r
+       C  THE COMBINATION  OF THREE DIFFERENT PATHS  IS USED TO INTEGRATE\r
+       C  FROM WA TO WB TO AVOID DOMAIN PROBLEM.THE BY-PRODUCT OF THIS IS\r
+       C  THAT IT IS NUMERICALLY  MORE EFFICIENT THAN A SINGLE LINE PATH:\r
+                 WA = W02(IND(K+1)-1)\r
+                 WB = W02(IND(K+1)+1)\r
+                 PHIA = PHI02(IND(K+1)-1)\r
+                 PHIB = PHI02(IND(K+1)+1)\r
+                 RADIUS = (1.D0+U2)/2.D0\r
+                 WAI = RADIUS*EXP(ZI*PHIA)\r
+                 WBI = RADIUS*EXP(ZI*PHIB)\r
+                 WLINE1 = WQUAD(WA,0.D0,IND(K+1)-1,0,WAI,0.D0,0,2,0.D0,M2,N2,U2,\r
+            +             W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2)\r
+                 WLINE2 = WQUAD(WB,0.D0,IND(K+1)+1,0,WBI,0.D0,0,2,0.D0,M2,N2,U2,\r
+            +             W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2)\r
+                 WCIRCLE = WQUAD(WAI,PHIA,0,2,WBI,PHIB,0,2,RADIUS,M2,N2,U2,W02,W12,\r
+            +              ALFA02,ALFA12,NPTQ2,QWORK2,1,2)\r
+                 WIN4 = C2* (WLINE1+WCIRCLE-WLINE2)\r
+                 FVAL(N2+5+IND(K+1)-1) = DIMAG(ZI*\r
+            +                           (Z02(IND(K+1)+1)-Z02(IND(K+1)-1)-WIN4))\r
+                 FVAL(N2+5+IND(K+1)) = DIMAG(Z02(IND(K+1)+1)-Z02(IND(K+1)-1)-WIN4)\r
+         100 CONTINUE\r
+       C\r
+       C  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
+             IF (ISPRT.NE.1) GO TO 110\r
+             FMAXN = FMAX(NDIM,FVAL)\r
+             WRITE (6,FMT=9020) ICOUNT,FMAXN\r
+         110 CONTINUE\r
+             RETURN\r
+\r
+        9000 FORMAT (2X,I5,6X,E10.4)\r
+        9010 FORMAT (2X,I5,6X,E10.4)\r
+        9020 FORMAT (2X,I5,6X,E10.4)*/\r
+       }\r
+\r
 \r
        \r
        private void DSCSOLV(double TOL,int IGUESS,int M,int N, double U[],\r
@@ -847,29 +1264,16 @@ double neweps;
        //      DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N),\r
        //                      QWORK(NPTQ* (2* (N+M)+3))\r
                \r
-               // Scalars in Common\r
-               double C2[] = new double[2];\r
-               double U2;\r
-               int ICOUNT, ISHAPE2, ISPRT, LINEARC2, M2, N2, NPTQ2, NSHAPE;\r
-               \r
-               // Arrays in Common\r
-               double W02[][] = new double[30][2];\r
-               double W12[][] = new double[30][2];\r
-               double Z02[][] = new double[30][2];\r
-               double Z12[][] = new double[30][2];\r
-               double ALFA02[] = new double[30];\r
-               double ALFA12[] = new double[30];\r
-               double PHI02[] = new double[30];\r
-               double PHI12[] = new double[30];\r
-               double QWORK2[] = new double[1660];\r
-               int IND[] = new int[20];\r
-               \r
                // Local scalars\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
                double C1[] = new double[2];\r
-               double WINT[] = new double[2];\r
+               double WINT[];\r
                double ZI[] = new double[2];\r
                double AVE, BOTM, DSTEP, FACTOR, PI, TOP;\r
-               int I, INFO, K, KM, KN, MAXFUN, NFEV, NM, NWDIM;\r
+               int INFO[] = new int[1];\r
+               int NFEV[] = new int[1];\r
+               int I, K, KM, KN, MAXFUN, NM, NWDIM;\r
                \r
                // Local arrays\r
                double DIAG[] = new double[42];\r
@@ -978,7 +1382,65 @@ double neweps;
         THDATA(U);\r
         WINT = WQUAD(W0[M-1],0.0,M,0,W1[N-1],0.0,N,1,0.0,M,N,U[0],W0,W1,ALFA0,\r
                            ALFA1,NPTQ,QWORK,0,2);\r
-\r
+        scm.zdiv(Z1[N-1][0] - Z0[M-1][0], Z1[N-1][1] - Z0[M-1][1], WINT[0], WINT[1], cr, ci);\r
+        C1[0] = cr[0];\r
+        C1[1] = ci[0];\r
+        scm.zmlt(ZI[0], ZI[1], C1[0], C1[1], cr, ci);\r
+        X[1] = ci[0];\r
+        X[2] = C1[1];\r
+        \r
+        //  HYBRD CONTROL PARAMETERS:\r
+        DSTEP = 1.0E-7;\r
+        MAXFUN = 200* (M+N);\r
+        FACTOR = 2.0;\r
+        NM = M + N + 2;\r
+      \r
+        //  COPY RELEVANT DATA TO THE COMMON BLOCK IN DSCFUN:\r
+        M2 = M;\r
+        N2 = N;\r
+        ISHAPE2 = ISHAPE;\r
+        LINEARC2 = LINEARC;\r
+        NPTQ2 = NPTQ;\r
+        for (K = 0; K < M; K++) {\r
+            W02[K][0] = W0[K][0];\r
+            W02[K][1] = W0[K][1];\r
+            PHI02[K] = PHI0[K];\r
+            Z02[K][0] = Z0[K][0];\r
+            Z02[K][1] = Z0[K][1];\r
+            ALFA02[K] = ALFA0[K];\r
+        } // for (K = 0; K < M; K++) \r
+        for (K = 0; K < N; K++) {\r
+            W12[K][0] = W1[K][0];\r
+            W12[K][1] = W1[K][1];\r
+            PHI12[K] = PHI1[K];\r
+            Z12[K][0] = Z1[K][0];\r
+            Z12[K][1] = Z1[K][1];\r
+            ALFA12[K] = ALFA1[K];\r
+        } // for (K = 0; K < N; K++)\r
+        NWDIM = NPTQ* (2* (M+N)+3);\r
+        for (I = 0; I < NWDIM; I++) {\r
+            QWORK2[I] = QWORK[I];\r
+        } // for (I = 0; I < NWDIM; I++)\r
+        \r
+        if (ISPRT != 1) {\r
+            System.out.println("Time for obtaining a converged result may range");\r
+            System.out.println("from several seconds to several minutes or longer");\r
+            System.out.println("depending on the complexity of the geometry of the");\r
+            System.out.println("region and the local computing system");\r
+            System.out.println("so be patient");\r
+        }\r
+        \r
+        // Solve nonlinear system with hybrid\r
+        double QTF[] = new double[NM];\r
+        double WA1[] = new double[NM];\r
+        double WA2[] = new double[NM];\r
+        double WA3[] = new double[NM];\r
+        double WA4[] = new double[NM];\r
+        HYBRD(dscfun, NM, X, FVAL, TOL, MAXFUN, NM, NM, DSTEP, DIAG, 1, FACTOR,\r
+                 -1, INFO, NFEV, FJAC, 42, QW, 903, QTF, WA1, WA2, WA3, WA4);\r
+        for (I = 0; I < NM; I++) {\r
+               QW[I+903] = QTF[I];\r
+        }\r
        }\r
 \r
        \r
@@ -1133,6 +1595,520 @@ double neweps;
                return;\r
        }\r
        \r
+       private void HYBRD(int FCN,int N, double X[], double FVEC[], double XTOL, int MAXFEV,\r
+                       int ML, int MU, double EPSFCN, double DIAG[], int MODE,\r
+                   double FACTOR, int NPRINT,int INFO[], int NFEV[], double FJAC[][],\r
+                   int LDFJAC, double R[],int LR, double QTF[], double WA1[],\r
+                   double WA2[], double WA3[], double WA4[]) {\r
+               //     ***********\r
+               //\r
+               //     SUBROUTINE HYBRD\r
+               //\r
+               //     THE PURPOSE OF HYBRD IS TO FIND A ZERO OF A SYSTEM OF\r
+               //     N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION\r
+               //     OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A\r
+               //     SUBROUTINE WHICH CALCULATES THE FUNCTIONS. THE JACOBIAN IS\r
+               //     THEN CALCULATED BY A FORWARD-DIFFERENCE APPROXIMATION.\r
+               //\r
+               //     THE SUBROUTINE STATEMENT IS\r
+               //\r
+               //       SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,\r
+               //                        DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,\r
+               //                        LDFJAC,R,LR,QTF,WA1,WA2,WA3,WA4)\r
+               //\r
+               //     WHERE\r
+           //      dscfun = 1 for DSCFUN\r
+               //       FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH\r
+               //         CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED\r
+               //         IN AN EXTERNAL STATEMENT IN THE USER CALLING\r
+               //         PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.\r
+               //\r
+               //         SUBROUTINE FCN(N,X,FVEC,IFLAG)\r
+               //         INTEGER N,IFLAG\r
+               //         DOUBLE PRECISION X(N),FVEC(N)\r
+               //         ----------\r
+               //         CALCULATE THE FUNCTIONS AT X AND\r
+               //         RETURN THIS VECTOR IN FVEC.\r
+               //         ---------\r
+               //         RETURN\r
+               //         END\r
+               //\r
+               //         THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS\r
+               //         THE USER WANTS TO TERMINATE EXECUTION OF HYBRD.\r
+               //         IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.\r
+               //\r
+               //       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
+               //         OF FUNCTIONS AND VARIABLES.\r
+               //\r
+               //       X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN\r
+               //         AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X\r
+               //         CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.\r
+               //\r
+               //       FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS\r
+               //         THE FUNCTIONS EVALUATED AT THE OUTPUT X.\r
+               //\r
+               //       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION\r
+               //         OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE\r
+               //         ITERATES IS AT MOST XTOL.\r
+               //\r
+               //       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION\r
+               //         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV\r
+               //         BY THE END OF AN ITERATION.\r
+               //\r
+               //       ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES\r
+               //         THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE\r
+               //         JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET\r
+               //         ML TO AT LEAST N - 1.\r
+               //\r
+               //       MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES\r
+               //         THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE\r
+               //         JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET\r
+               //         MU TO AT LEAST N - 1.\r
+               //\r
+               //       EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE\r
+               //         STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS\r
+               //         APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE\r
+               //         FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS\r
+               //         THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE\r
+               //         ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE\r
+               //         PRECISION.\r
+               //\r
+               //       DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE\r
+               //         BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG\r
+               //         MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS\r
+               //         MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES.\r
+               //\r
+               //       MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE\r
+               //         VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2,\r
+               //         THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER\r
+               //         VALUES OF MODE ARE EQUIVALENT TO MODE = 1.\r
+               //\r
+               //       FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE\r
+               //         INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF\r
+               //         FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE\r
+               //         TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE\r
+               //         INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE.\r
+               //\r
+               //       NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED\r
+               //         PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE,\r
+               //         FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST\r
+               //         ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND\r
+               //         IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE\r
+               //         FOR PRINTING. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS\r
+               //         OF FCN WITH IFLAG = 0 ARE MADE.\r
+               //\r
+               //       INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS\r
+               //         TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)\r
+               //         VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,\r
+               //         INFO IS SET AS FOLLOWS.\r
+               //\r
+               //         INFO = 0   IMPROPER INPUT PARAMETERS.\r
+               //\r
+               //         INFO = 1   RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES\r
+               //                    IS AT MOST XTOL.\r
+               //\r
+               //         INFO = 2   NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED\r
+               //                    MAXFEV.\r
+               //\r
+               //         INFO = 3   XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN\r
+               //                    THE APPROXIMATE SOLUTION X IS POSSIBLE.\r
+               //\r
+               //         INFO = 4   ITERATION IS NOT MAKING GOOD PROGRESS, AS\r
+               //                    MEASURED BY THE IMPROVEMENT FROM THE LAST\r
+               //                    FIVE JACOBIAN EVALUATIONS.\r
+               //\r
+               //         INFO = 5   ITERATION IS NOT MAKING GOOD PROGRESS, AS\r
+               //                    MEASURED BY THE IMPROVEMENT FROM THE LAST\r
+               //                    TEN ITERATIONS.\r
+               //\r
+               //       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF\r
+               //         CALLS TO FCN.\r
+               //\r
+               //       FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE\r
+               //         ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION\r
+               //         OF THE FINAL APPROXIMATE JACOBIAN.\r
+               //\r
+               //       LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N\r
+               //         WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.\r
+               //\r
+               //       R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE\r
+               //         UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION\r
+               //         OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE.\r
+               //\r
+               //       LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN\r
+               //         (N*(N+1))/2.\r
+               //\r
+               //       QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS\r
+               //         THE VECTOR (Q TRANSPOSE)*FVEC.\r
+               //\r
+               //       WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N.\r
+               //\r
+               //     SUBPROGRAMS CALLED\r
+               //\r
+               //       USER-SUPPLIED ...... FCN\r
+               //\r
+               //       MINPACK-SUPPLIED ... DOGLEG,DPMPAR,ENORM,FDJAC1,\r
+               //                            QFORM,QRFAC,R1MPYQ,R1UPDT\r
+               //\r
+               //       FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,MIN0,MOD\r
+               //\r
+               //     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. MARCH 1980.\r
+               //     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE\r
+               //\r
+               //     **********\r
+               //     .. Scalar Arguments ..\r
+               //      DOUBLE PRECISION EPSFCN,FACTOR,XTOL\r
+               //      INTEGER INFO,LDFJAC,LR,MAXFEV,ML,MODE,MU,N,NFEV,NPRINT\r
+               //     ..\r
+               //     .. Array Arguments ..\r
+               //      DOUBLE PRECISION DIAG(N),FJAC(LDFJAC,N),FVEC(N),QTF(N),R(LR),\r
+               //                      WA1(N),WA2(N),WA3(N),WA4(N),X(N)\r
+               //     ..\r
+               //     .. Subroutine Arguments ..\r
+               //      EXTERNAL FCN\r
+               //     ..\r
+               //     .. Local Scalars ..\r
+                     double ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,P0001,P001,\r
+                                     P1,P5,PNORM,PRERED,RATIO,SUM,TEMP,XNORM,ZERO;\r
+                     int IFLAG[] = new int[1];\r
+                     int I,ITER,J,JM1,L,MSUM,NCFAIL,NCSUC,NSLOW1,NSLOW2;\r
+                     boolean JEVAL,SING;\r
+               //     ..\r
+           //     .. Local Arrays ..\r
+                     int IWA[] = new int[1];\r
+               //     ..\r
+               //     .. External Functions ..\r
+               //      DOUBLE PRECISION DPMPAR,ENORM\r
+               //      EXTERNAL DPMPAR,ENORM\r
+               //     ..\r
+               //     .. External Subroutines ..\r
+               //      EXTERNAL DOGLEG,FDJAC1,QFORM,QRFAC,R1MPYQ,R1UPDT\r
+               //     ..\r
+               //     .. Intrinsic Functions ..\r
+               //      INTRINSIC DABS,DMAX1,DMIN1,MIN0,MOD\r
+               //     ..\r
+           //     .. Data statements ..\r
+               //\r
+               ONE = 1.0;\r
+               P1 = 1.0E-1;\r
+               P5 = 5.0E-1;\r
+               P001 = 1.0E-3;\r
+               P0001 = 1.0E-4;\r
+               ZERO = 0.0;\r
+               //\r
+               //     ..\r
+               //\r
+               //     EPSMCH IS THE MACHINE PRECISION.\r
+               //\r
+               EPSMCH = MACHEP;\r
+               //\r
+               INFO[0] = 0;\r
+               IFLAG[0] = 0;\r
+               NFEV[0] = 0;\r
+               //\r
+               //     CHECK THE INPUT PARAMETERS FOR ERRORS.\r
+               //\r
+               if (N <= 0 || XTOL < ZERO || MAXFEV <= 0 || ML < 0 ||\r
+                        MU < 0 || FACTOR <= ZERO || LDFJAC < N ||\r
+                        LR < (N* (N+1))/2) {\r
+                       if (NPRINT > 0) {\r
+                               if (FCN == dscfun) {\r
+                                       DSCFUN(N, X, FVEC, IFLAG);\r
+                               }\r
+                               return;\r
+                       } // if (NPRINT > 0)\r
+                       return;\r
+               }\r
+               if (MODE == 2) {\r
+                   for (J = 0; J < N; J++) {\r
+                       if (DIAG[J] <= ZERO) {\r
+                               if (NPRINT > 0) {\r
+                                               if (FCN == dscfun) {\r
+                                                       DSCFUN(N, X, FVEC, IFLAG);\r
+                                               }\r
+                                       } // if (NPRINT > 0)\r
+                               return;\r
+                       } // if (DIAG[J] <= ZERO)\r
+                   } // for (J = 0; J < N; J++)\r
+               } // if (mode == 2)  \r
+               //\r
+               //     EVALUATE THE FUNCTION AT THE STARTING POINT\r
+               //     AND CALCULATE ITS NORM.\r
+               //\r
+               IFLAG[0] = 1;\r
+               if (FCN == dscfun) {\r
+                   DSCFUN(N,X,FVEC,IFLAG);\r
+               }\r
+               NFEV[0] = 1;\r
+               if (IFLAG[0] < 0) {\r
+                   INFO[0] = IFLAG[0];\r
+                   IFLAG[0] = 0;\r
+                   if (NPRINT > 0) {\r
+                               if (FCN == dscfun) {\r
+                                       DSCFUN(N, X, FVEC, IFLAG);\r
+                               }\r
+                       } // if (NPRINT > 0)\r
+               return;\r
+               } // if (IFLAG[0] < 0)\r
+               /*      FNORM = ENORM(N,FVEC)\r
+               C\r
+               C     DETERMINE THE NUMBER OF CALLS TO FCN NEEDED TO COMPUTE\r
+               C     THE JACOBIAN MATRIX.\r
+               C\r
+                     MSUM = MIN0(ML+MU+1,N)\r
+               C\r
+               C     INITIALIZE ITERATION COUNTER AND MONITORS.\r
+               C\r
+                     ITER = 1\r
+                     NCSUC = 0\r
+                     NCFAIL = 0\r
+                     NSLOW1 = 0\r
+                     NSLOW2 = 0\r
+               C\r
+               C     BEGINNING OF THE OUTER LOOP.\r
+               C\r
+                  30 CONTINUE\r
+                     JEVAL = .TRUE.\r
+               C\r
+               C        CALCULATE THE JACOBIAN MATRIX.\r
+               C\r
+                     IFLAG = 2\r
+                     CALL FDJAC1(FCN,N,X,FVEC,FJAC,LDFJAC,IFLAG,ML,MU,EPSFCN,WA1,WA2)\r
+                     NFEV = NFEV + MSUM\r
+                     IF (IFLAG.LT.0) GO TO 300\r
+               C\r
+               C        COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.\r
+               C\r
+                     CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)\r
+               C\r
+               C        ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING\r
+               C        TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.\r
+               C\r
+                     IF (ITER.NE.1) GO TO 70\r
+                     IF (MODE.EQ.2) GO TO 50\r
+                     DO 40 J = 1,N\r
+                         DIAG(J) = WA2(J)\r
+                         IF (WA2(J).EQ.ZERO) DIAG(J) = ONE\r
+                  40 CONTINUE\r
+                  50 CONTINUE\r
+               C\r
+               C        ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X\r
+               C        AND INITIALIZE THE STEP BOUND DELTA.\r
+               C\r
+                     DO 60 J = 1,N\r
+                         WA3(J) = DIAG(J)*X(J)\r
+                  60 CONTINUE\r
+                     XNORM = ENORM(N,WA3)\r
+                     DELTA = FACTOR*XNORM\r
+                     IF (DELTA.EQ.ZERO) DELTA = FACTOR\r
+                  70 CONTINUE\r
+               C\r
+               C        FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.\r
+               C\r
+                     DO 80 I = 1,N\r
+                         QTF(I) = FVEC(I)\r
+                  80 CONTINUE\r
+                     DO 120 J = 1,N\r
+                         IF (FJAC(J,J).EQ.ZERO) GO TO 110\r
+                         SUM = ZERO\r
+                         DO 90 I = J,N\r
+                             SUM = SUM + FJAC(I,J)*QTF(I)\r
+                  90     CONTINUE\r
+                         TEMP = -SUM/FJAC(J,J)\r
+                         DO 100 I = J,N\r
+                             QTF(I) = QTF(I) + FJAC(I,J)*TEMP\r
+                 100     CONTINUE\r
+                 110     CONTINUE\r
+                 120 CONTINUE\r
+               C\r
+               C        COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.\r
+               C\r
+                     SING = .FALSE.\r
+                     DO 150 J = 1,N\r
+                         L = J\r
+                         JM1 = J - 1\r
+                         IF (JM1.LT.1) GO TO 140\r
+                         DO 130 I = 1,JM1\r
+                             R(L) = FJAC(I,J)\r
+                             L = L + N - I\r
+                 130     CONTINUE\r
+                 140     CONTINUE\r
+                         R(L) = WA1(J)\r
+                         IF (WA1(J).EQ.ZERO) SING = .TRUE.\r
+                 150 CONTINUE\r
+               C\r
+               C        ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.\r
+               C\r
+                     CALL QFORM(N,N,FJAC,LDFJAC,WA1)\r
+               C\r
+               C        RESCALE IF NECESSARY.\r
+               C\r
+                     IF (MODE.EQ.2) GO TO 170\r
+                     DO 160 J = 1,N\r
+                         DIAG(J) = DMAX1(DIAG(J),WA2(J))\r
+                 160 CONTINUE\r
+                 170 CONTINUE\r
+               C\r
+               C        BEGINNING OF THE INNER LOOP.\r
+               C\r
+                 180 CONTINUE\r
+               C\r
+               C           IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.\r
+               C\r
+                     IF (NPRINT.LE.0) GO TO 190\r
+                     IFLAG = 0\r
+                     IF (MOD(ITER-1,NPRINT).EQ.0) CALL FCN(N,X,FVEC,IFLAG)\r
+                     IF (IFLAG.LT.0) GO TO 300\r
+                 190 CONTINUE\r
+               C\r
+               C           DETERMINE THE DIRECTION P.\r
+               C\r
+                     CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)\r
+               C\r
+               C           STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.\r
+               C\r
+                     DO 200 J = 1,N\r
+                         WA1(J) = -WA1(J)\r
+                         WA2(J) = X(J) + WA1(J)\r
+                         WA3(J) = DIAG(J)*WA1(J)\r
+                 200 CONTINUE\r
+                     PNORM = ENORM(N,WA3)\r
+               C\r
+               C           ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.\r
+               C\r
+                     IF (ITER.EQ.1) DELTA = DMIN1(DELTA,PNORM)\r
+               C\r
+               C           EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.\r
+               C\r
+                     IFLAG = 1\r
+                     CALL FCN(N,WA2,WA4,IFLAG)\r
+                     NFEV = NFEV + 1\r
+                     IF (IFLAG.LT.0) GO TO 300\r
+                     FNORM1 = ENORM(N,WA4)\r
+               C\r
+               C           COMPUTE THE SCALED ACTUAL REDUCTION.\r
+               C\r
+                     ACTRED = -ONE\r
+                     IF (FNORM1.LT.FNORM) ACTRED = ONE - (FNORM1/FNORM)**2\r
+               C\r
+               C           COMPUTE THE SCALED PREDICTED REDUCTION.\r
+               C\r
+                     L = 1\r
+                     DO 220 I = 1,N\r
+                         SUM = ZERO\r
+                         DO 210 J = I,N\r
+                             SUM = SUM + R(L)*WA1(J)\r
+                             L = L + 1\r
+                 210     CONTINUE\r
+                         WA3(I) = QTF(I) + SUM\r
+                 220 CONTINUE\r
+                     TEMP = ENORM(N,WA3)\r
+                     PRERED = ZERO\r
+                     IF (TEMP.LT.FNORM) PRERED = ONE - (TEMP/FNORM)**2\r
+               C\r
+               C           COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED\r
+               C           REDUCTION.\r
+               C\r
+                     RATIO = ZERO\r
+                     IF (PRERED.GT.ZERO) RATIO = ACTRED/PRERED\r
+               C\r
+               C           UPDATE THE STEP BOUND.\r
+               C\r
+                     IF (RATIO.GE.P1) GO TO 230\r
+                     NCSUC = 0\r
+                     NCFAIL = NCFAIL + 1\r
+                     DELTA = P5*DELTA\r
+                     GO TO 240\r
+\r
+                 230 CONTINUE\r
+                     NCFAIL = 0\r
+                     NCSUC = NCSUC + 1\r
+                     IF (RATIO.GE.P5 .OR. NCSUC.GT.1) DELTA = DMAX1(DELTA,PNORM/P5)\r
+                     IF (DABS(RATIO-ONE).LE.P1) DELTA = PNORM/P5\r
+                 240 CONTINUE\r
+               C\r
+               C           TEST FOR SUCCESSFUL ITERATION.\r
+               C\r
+                     IF (RATIO.LT.P0001) GO TO 260\r
+               C\r
+               C           SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.\r
+               C\r
+                     DO 250 J = 1,N\r
+                         X(J) = WA2(J)\r
+                         WA2(J) = DIAG(J)*X(J)\r
+                         FVEC(J) = WA4(J)\r
+                 250 CONTINUE\r
+                     XNORM = ENORM(N,WA2)\r
+                     FNORM = FNORM1\r
+                     ITER = ITER + 1\r
+                 260 CONTINUE\r
+               C\r
+               C           DETERMINE THE PROGRESS OF THE ITERATION.\r
+               C\r
+                     NSLOW1 = NSLOW1 + 1\r
+                     IF (ACTRED.GE.P001) NSLOW1 = 0\r
+                     IF (JEVAL) NSLOW2 = NSLOW2 + 1\r
+                     IF (ACTRED.GE.P1) NSLOW2 = 0\r
+               C\r
+               C           TEST FOR CONVERGENCE.\r
+               C\r
+                     IF (DELTA.LE.XTOL*XNORM .OR. FNORM.EQ.ZERO) INFO = 1\r
+                     IF (INFO.NE.0) GO TO 300\r
+               C\r
+               C           TESTS FOR TERMINATION AND STRINGENT TOLERANCES.\r
+               C\r
+                     IF (NFEV.GE.MAXFEV) INFO = 2\r
+                     IF (P1*DMAX1(P1*DELTA,PNORM).LE.EPSMCH*XNORM) INFO = 3\r
+                     IF (NSLOW2.EQ.5) INFO = 4\r
+                     IF (NSLOW1.EQ.10) INFO = 5\r
+                     IF (INFO.NE.0) GO TO 300\r
+               C\r
+               C           CRITERION FOR RECALCULATING JACOBIAN APPROXIMATION\r
+               C           BY FORWARD DIFFERENCES.\r
+               C\r
+                     IF (NCFAIL.EQ.2) GO TO 290\r
+               C\r
+               C           CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN\r
+               C           AND UPDATE QTF IF NECESSARY.\r
+               C\r
+                     DO 280 J = 1,N\r
+                         SUM = ZERO\r
+                         DO 270 I = 1,N\r
+                             SUM = SUM + FJAC(I,J)*WA4(I)\r
+                 270     CONTINUE\r
+                         WA2(J) = (SUM-WA3(J))/PNORM\r
+                         WA1(J) = DIAG(J)* ((DIAG(J)*WA1(J))/PNORM)\r
+                         IF (RATIO.GE.P0001) QTF(J) = SUM\r
+                 280 CONTINUE\r
+               C\r
+               C           COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.\r
+               C\r
+                     CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)\r
+                     CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)\r
+                     CALL R1MPYQ(1,N,QTF,1,WA2,WA3)\r
+               C\r
+               C           END OF THE INNER LOOP.\r
+               C\r
+                     JEVAL = .FALSE.\r
+                     GO TO 180\r
+\r
+                 290 CONTINUE\r
+               C\r
+               C        END OF THE OUTER LOOP.\r
+               C\r
+                     GO TO 30\r
+\r
+                 300 CONTINUE\r
+               C\r
+               C     TERMINATION, EITHER NORMAL OR USER IMPOSED.\r
+               C\r
+                     IF (IFLAG.LT.0) INFO = IFLAG\r
+                     IFLAG = 0\r
+                     IF (NPRINT.GT.0) CALL FCN(N,X,FVEC,IFLAG)\r
+               return;*/\r
+      }\r
+\r
+       \r
        private void QINIT(int M, int N, double ALFA0[], double ALFA1[], int NPTQ, double QWORK[]) {\r
            // Computes the Gauss-Jacobi nodes & weights for Gauss-Jacobi quadrature.  Work array\r
                // QWORK must be dimensioned at least NPTQ*(2(M+N) + 3).  It is divided up into\r