Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 27 Nov 2017 21:58:53 +0000 (21:58 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 27 Nov 2017 21:58:53 +0000 (21:58 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15276 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 58d8ef4..ae04b9c 100644 (file)
@@ -124,6 +124,7 @@ public class DoublyConnectedSC extends AlgorithmBase {
                double ZZ[] = new double[2];\r
                double EPS;\r
                double WW[];\r
+               double ZZ0[];\r
                \r
 double neweps;\r
                \r
@@ -188,6 +189,14 @@ double neweps;
                            EPS = 1.0E-6;\r
                            WW = WDSC(ZZ,M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,\r
                                                  NPTQ,QWORK,EPS,1);\r
+                    System.out.println("THE PREIMAGE OF ZZ = (" + WW[0] + ", " + WW[1] + ")");\r
+                    if (scm.zabs(WW[0], WW[1]) <= 1.0E-12) {\r
+                       continue;\r
+                    }\r
+                    System.out.println("CHECK BY MAPPING THE PREIMAGE BACK");\r
+                    ZZ0 = ZDSC(WW,0,2,M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
+                         PHI1,NPTQ,QWORK,1);\r
+                    System.out.println("THE POINT ENTERED = (" + ZZ0[0] + ", " + ZZ0[1] + ")");\r
 \r
                        } // for (I = 0; I < INVERSE_POINTS.length; I++)\r
                    } // if ((INVERSE_POINTS != null) && (INVERSE_POINTS.length != 0))\r
@@ -1143,7 +1152,7 @@ double neweps;
            // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
            // COMMON /PARAM5/ISPRT,ICOUNT\r
 \r
-             XWTRAN(M2,N2,X,U2,C2,W02,W12,PHI02,PHI12);\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
@@ -1246,27 +1255,27 @@ double neweps;
              //  OUTER POLYGON CONTAINS SOME INFINITE VERTICES & FOR EACH OF THEM\r
              //  TWO LENGTH CONDITIONS WILL BE REPLACED BY A COMPLEX INTEGRAL:\r
              for (K = 1; K <= NSHAPE - 1; K++) {\r
-                 if (IND[K] != 1 && IND[K-1] < IND[K]-2) {\r
-                     for (J = IND[K-1] + 2; J <= IND[K] - 1; J++) {\r
+                 if (IND[K] != 2 && IND[K-1] < IND[K]-2) {\r
+                     for (J = IND[K-1] + 1; J <= IND[K] - 2; J++) {\r
                          WINT3 = WQUAD(W02[J-1],PHI02[J-1],J,0,W02[J],PHI02[J],J+1,0,\r
                             1.0,M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
                          scm.zmlt(C2[0], C2[1], WINT3[0], WINT3[1], cr, ci);\r
                          FVAL[N2+4+J] = scm.zabs(Z02[J][0]-Z02[J-1][0], Z02[J][1]-Z02[J-1][1]) - scm.zabs(cr[0], ci[0]);\r
                      } // for (J = IND[K-1] + 2; J <= IND[K] - 1; J++) \r
                  } // if (IND[K] != 1 && IND[K-1] < IND[K]-2)\r
-                 if (K == NSHAPE-1 || IND[K] == M2-2) {\r
+                 if (K == NSHAPE-1 || IND[K] == M2-1) {\r
                          continue;\r
                  }\r
        \r
                  //  THE COMBINATION  OF THREE DIFFERENT PATHS  IS USED TO INTEGRATE\r
                  //  FROM WA TO WB TO AVOID DOMAIN PROBLEM.THE BY-PRODUCT OF THIS IS\r
                  //  THAT IT IS NUMERICALLY  MORE EFFICIENT THAN A SINGLE LINE PATH:\r
-                 WA[0] = W02[IND[K]-1][0];\r
-                 WA[1] = W02[IND[K]-1][1];\r
-                 WB[0] = W02[IND[K]+1][0];\r
-                 WB[1] = W02[IND[K]+1][1];\r
-                 PHIA = PHI02[IND[K]-1];\r
-                 PHIB = PHI02[IND[K]+1];\r
+                 WA[0] = W02[IND[K]-2][0];\r
+                 WA[1] = W02[IND[K]-2][1];\r
+                 WB[0] = W02[IND[K]][0];\r
+                 WB[1] = W02[IND[K]][1];\r
+                 PHIA = PHI02[IND[K]-2];\r
+                 PHIB = PHI02[IND[K]];\r
                  RADIUS = (1.0+U2[0])/2.0;\r
                  scm.zmlt(ZI[0], ZI[1], PHIA, 0.0, cr, ci);\r
                  base = RADIUS * Math.exp(cr[0]);\r
@@ -1286,11 +1295,11 @@ double neweps;
                                  (WLINE1[1]+WCIRCLE[1]-WLINE2[1]), cr, ci);\r
                  WIN4[0] = cr[0];\r
                  WIN4[1] = ci[0];\r
-                 diffR = Z02[IND[K]+1][0]-Z02[IND[K]-1][0]-WIN4[0];\r
-                 diffI = Z02[IND[K]+1][1]-Z02[IND[K]-1][1]-WIN4[1];\r
+                 diffR = Z02[IND[K]][0]-Z02[IND[K]-2][0]-WIN4[0];\r
+                 diffI = Z02[IND[K]][1]-Z02[IND[K]-2][1]-WIN4[1];\r
                  scm.zmlt(ZI[0], ZI[1], diffR, diffI, cr, ci);\r
-                 FVAL[N2+5+IND[K]-1] = ci[0];\r
-                 FVAL[N2+5+IND[K]] = diffI;\r
+                 FVAL[N2+5+IND[K]-2] = ci[0];\r
+                 FVAL[N2+5+IND[K]-1] = diffI;\r
              } // for (K = 1; K <= NSHAPE - 1; K++)\r
        \r
              //  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
@@ -1384,11 +1393,11 @@ double neweps;
                                continue;\r
                            }\r
                            NSHAPE = NSHAPE + 1;\r
-                           IND[NSHAPE-1] = I-1;\r
+                           IND[NSHAPE-1] = I;\r
                        } // for (I = 2; I <= M-1; I++)\r
-                       IND[0] = -1;\r
+                       IND[0] = 0;\r
                        NSHAPE = NSHAPE + 1;\r
-                       IND[NSHAPE-1] = M-2;\r
+                       IND[NSHAPE-1] = M-1;\r
                } // if (ISHAPE != 0)\r
                \r
                // Fix one prevertex\r
@@ -1561,6 +1570,171 @@ double neweps;
 \r
        }\r
        \r
+       private double[] ZDSC(double WW[],int KWW,int IC,int M,int N, double U,\r
+                       double C[], double W0[][], double W1[][], double Z0[][], double Z1[][],\r
+                       double ALFA0[], double ALFA1[], double PHI0[], double PHI1[], int NPTQ,\r
+                       double QWORK[], int IOPT) {\r
+        //   ----------------------------------------------------------\r
+               //   COMPUTES THE FORWORD MAP Z(WW) BY INTEGRATING FROM WA TO\r
+               //   WW WHERE WA IS THE NEAREST PREVERTEX TO WW.KWW=K,IC=0 IF\r
+               //   WW=W0(K), OR KWW=K, IC=1 IF WW=W1(K),OR KWW=0 AND IC=2\r
+               //   OTHERWISE.\r
+               //   IOPT: =1 IS NORMALLY ASSUMED FOR ANY GEOMETRY.\r
+               //         =# OTHER THAN 1 ASSUMES THAT ONLY LINE SEGMENT PATH\r
+               //          IS USED.\r
+               //   SEE ROUTINE DSCSOLV FOR THE REST OF THE CALLING SEQUENCE.\r
+               \r
+               \r
+               //     .. Scalar Arguments ..\r
+               //    DOUBLE COMPLEX C,WW\r
+               //    DOUBLE PRECISION U\r
+               //    INTEGER IC,IOPT,KWW,M,N,NPTQ\r
+               \r
+               //     .. Array Arguments ..\r
+               //    DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
+               //    DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N),\r
+               //                     QWORK(NPTQ* (2* (M+N)+3))\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 WA[] = new double[2];\r
+               double WB[] = new double[2];\r
+               double WINT1[] = new double[2];\r
+               double WINT2[] = new double[2];\r
+               double WW0[] = new double[2];\r
+               double ZA[] = new double[2];\r
+               double ZI[] = new double[2];\r
+               double wout[];\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+               double result[] = new double[2];\r
+               double base;\r
+               // DOUBLE COMPLEX WA,WB,WINT1,WINT2,WW0,ZA,ZI\r
+               double DWW0,PHIWB,PHIWW0,PI;\r
+               double factor;\r
+               int INEAR[] = new int[1];\r
+               int KNEAR[] = new int[1];\r
+               int IBD;\r
+\r
+               //     .. External Functions ..\r
+               //    DOUBLE COMPLEX WQUAD\r
+               //    DOUBLE PRECISION ARGUM\r
+               //    EXTERNAL WQUAD,ARGUM\r
+       \r
+               //     .. External Subroutines ..\r
+               //    EXTERNAL NEARW\r
+\r
+               //     .. Common blocks ..\r
+               //    COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+               \r
+               ZI[0] = 0.0;\r
+               ZI[1] = 1.0;\r
+               PI = Math.PI;\r
+               IBD = 1;\r
+               WW0[0] = WW[0];\r
+               WW0[1] = WW[1];\r
+               \r
+               if (IOPT == 1) {\r
+                   if (Math.abs(scm.zabs(WW0[0],WW0[1])-1.0) <= 1.0E-11) {\r
+                       factor = (1.0 + U)/2.0;\r
+                       WW0[0] = factor * WW0[0];\r
+                       WW0[1] = factor * WW0[1];\r
+                       IBD = 2;\r
+                   } // if (Math.abs(scm.zabs(WW0[0],WW0[1])-1.0) <= 1.0E-11)\r
+               } // if (IOPT == 1)\r
+               \r
+               // FIND THE NEAREST PREVERTEX TO THE PT WW0:\r
+               NEARW(M,N,W0,W1,ALFA0,WW0,KNEAR,INEAR);\r
+               if (INEAR[0] == 0) {\r
+                   ZA[0] = Z0[KNEAR[0]-1][0];\r
+                   ZA[1] = Z0[KNEAR[0]-1][1];\r
+                   WA[0] = W0[KNEAR[0]-1][0];\r
+                   WA[1] = W0[KNEAR[0]-1][1];\r
+               }\r
+               else {\r
+                   ZA[0] = Z1[KNEAR[0]-1][0];\r
+                   ZA[1] = Z1[KNEAR[0]-1][1];\r
+                   WA[0] = W1[KNEAR[0]-1][0];\r
+                   WA[1] = W1[KNEAR[0]-1][1];\r
+               }\r
+\r
+               if (IOPT != 1) {\r
+                   wout = WQUAD(WA,0.0,KNEAR[0],INEAR[0],WW,0.0,KWW,IC,0.0,M,N,U,\r
+                                          W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
+                   scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
+                   result[0] = ZA[0] + cr[0];\r
+                   result[1] = ZA[1] + ci[0];\r
+                   return result;\r
+               } // if (IOPT != 1)\r
+               \r
+               // DETERMINE THE POINT CLOSEST TO WA ON THE\r
+               // SAME CONCENTRIC CIRCLE AS WW0 IS ON:\r
+               PHIWW0 = ARGUM(scm.zabs(WW0[0], WW0[1]),WW0);\r
+               DWW0 = scm.zabs(WW0[0],WW0[1]);\r
+               if (INEAR[0] == 0) {\r
+                   PHIWB = PHI0[KNEAR[0]-1];\r
+               }\r
+               else {\r
+                   PHIWB = PHI1[KNEAR[0]-1];\r
+               }\r
+\r
+               scm.zmlt(ZI[0], ZI[1], PHIWB, 0.0, cr, ci);\r
+               base = DWW0 * Math.exp(cr[0]);\r
+               WB[0] = base * Math.cos(ci[0]);\r
+               WB[1] = base * Math.sin(ci[0]);\r
+               //\r
+           //   INTEGRATION FROM WA TO WB ON A LINE SEGMENT:\r
+               if (scm.zabs(WB[0]-WA[0],WB[1]-WA[1]) <= 1.0E-11) {\r
+                       WINT1[0] = 0.0;\r
+                       WINT1[1] = 0.0;\r
+               }\r
+               else {\r
+                     WINT1 = WQUAD(WA,0.0,KNEAR[0],INEAR[0],WB,0.0,0,2,0.0,M,N,U,W0,W1,\r
+                           ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
+               }\r
+               \r
+               //   INTEGRATION FROM WB TO WW ON A CIRCULAR ARC:\r
+               if (scm.zabs(WB[0]-WW0[0],WB[1]-WW0[1]) <= 1.0E-11) {\r
+                       WINT2[0] = 0.0;\r
+                       WINT2[1] = 0.0;\r
+           }\r
+               else {\r
+                   if (Math.abs(PHIWB-2.0*PI-PHIWW0) < \r
+                        Math.abs(PHIWB-PHIWW0)) PHIWB = PHIWB - 2.0*PI;\r
+                   if (Math.abs(PHIWW0-2.0*PI-PHIWB) <\r
+                        Math.abs(PHIWB-PHIWW0)) PHIWW0 = PHIWW0 - 2.0*PI;\r
+                   if (scm.zabs(WB[0]-WA[0],WB[1]-WA[1]) <= 1.0E-11) {\r
+                         WINT2 = WQUAD(WB,PHIWB,KNEAR[0],INEAR[0],WW0,PHIWW0,KWW,IC,DWW0,M,N,\r
+                                U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,1,1);\r
+                   }\r
+                   else {\r
+                     WINT2 = WQUAD(WB,PHIWB,0,2,WW0,PHIWW0,0,2,DWW0,M,N,U,W0,W1,ALFA0,\r
+                                   ALFA1,NPTQ,QWORK,1,1);\r
+                   }\r
+               } // else\r
+               \r
+               // EVALUATE THE MAPPING FUNCTION. THE INTEGRATION PATH IS\r
+               // A COMBINATION OF A CIRCULAR ARC AND LINE SEGMENT(S):\r
+               scm.zmlt(C[0], C[1], WINT1[0] + WINT2[0], WINT1[1] + WINT2[1], cr, ci);\r
+               result[0] = ZA[0] + cr[0];\r
+               result[1] = ZA[1] + ci[0];\r
+               if (IBD == 2) {\r
+                       wout = WQUAD(WW0,0.0,0,2,WW,0.0,KWW,IC,0.0,M,N,U,\r
+                                               W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
+                       scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
+                       result[0] = result[0] + cr[0];\r
+                       result[1] = result[1] + ci[0];\r
+               } // if (IBD == 2)\r
+\r
+               return result;    \r
+    }\r
+       \r
        private double[] WDSC(double ZZ[], int M,int N, double U, double C[], double W0[][], double W1[][],\r
                        double Z0[][], double Z1[][], double ALFA0[], double ALFA1[], double PHI0[], double PHI1[],\r
                        int NPTQ, double QWORK[], double EPS,int IOPT) {\r
@@ -1598,6 +1772,7 @@ double neweps;
                double WI[] = new double[2];\r
                double ZS[] = new double[2];\r
                double ZZ1[] = new double[2];\r
+               double denom;\r
                // DOUBLE COMPLEX WFN,WI,ZS,ZZ1\r
                int INZ[] = new int[1];\r
                int KNZ[] = new int[1];\r
@@ -1606,6 +1781,15 @@ double neweps;
                // .. Local Arrays ..\r
                double ZS0[][] = new double[50][2];\r
                double ZS1[][] = new double[50][2];\r
+               double wout[];\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
+               double abswi;\r
+               double factor;\r
+               double zout[];\r
+               double absz;\r
                // DOUBLE COMPLEX ZS0(50),ZS1(50)\r
                \r
                //     .. External Functions ..\r
@@ -1622,7 +1806,7 @@ double neweps;
                result[1] = 0.0;\r
                \r
                //  1.FORM WORK ARRAYS:\r
-               /*for (I = 0; I < M; I++) {\r
+               for (I = 0; I < M; I++) {\r
                    ZS0[I][0] = Z0[I][0];\r
                    ZS0[I][1] = Z0[I][1];\r
                } // for (I = 0; I < M; I++)\r
@@ -1647,65 +1831,100 @@ double neweps;
                            WI[1] = (W0[KNZ[0]-1][1]+W0[KNZ[0]-2][1])/2.0;\r
                        }\r
                        else {\r
-                             WI[0] = (W0[KNZ[0]-1][0]+W0[M-1][0]))/2.0\r
-                      WI[1] = (W0[KNZ[0]-1][1]+W0[M-1][1]))/2.0\r
+                           WI[0] = (W0[KNZ[0]-1][0]+W0[M-1][0])/2.0;\r
+                    WI[1] = (W0[KNZ[0]-1][1]+W0[M-1][1])/2.0;\r
                        }\r
 \r
                    } // if (INZ[0] == 0)\r
                    else {\r
-                         IF (KNZ.GE.2) THEN\r
-                             WI = (W1(KNZ)+W1(KNZ-1))/2.D0\r
-\r
-                         ELSE\r
-                             WI = (W0(KNZ)+W1(N))/2.D0\r
-                         END IF\r
+                       if (KNZ[0] >= 2) {\r
+                             WI[0] = (W1[KNZ[0]-1][0]+W1[KNZ[0]-2][0])/2.0;\r
+                             WI[1] = (W1[KNZ[0]-1][1]+W1[KNZ[0]-2][1])/2.0;\r
+                       }\r
+                       else {\r
+                             WI[0] = (W0[KNZ[0]-1][0]+W1[N-1][0])/2.0;\r
+                             WI[1] = (W0[KNZ[0]-1][1]+W1[N-1][1])/2.0;\r
+                       }\r
 \r
-                         WI = U*WI/ABS(WI)\r
+                       denom = scm.zabs(WI[0], WI[1]);\r
+                       WI[0] = U*WI[0]/denom;\r
+                       WI[1] = U*WI[1]/denom;\r
                    } //else\r
 \r
-                     ZS = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ,\r
-                    +     QWORK,1)\r
-               C\r
-               C  3.SOLVE ODE INITIAL VALUE PROBLEM (ALONG LINE SEGMENT FROM\r
-               C    ZS TO ZZ) BY EULER'S SCHEME TO GENERATE THE INITIAL GUESS:\r
-                     DO 40 K = 1,20\r
-                         WI = WI + (ZZ-ZS)/ (20.D0*C*WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1))\r
-                  40 CONTINUE\r
-                     IF (ABS(WI).GT.1.D0) WI = WI/ (ABS(WI)+ABS(1.D0-ABS(WI)))\r
-                     IF (ABS(WI).LT.U) WI = U*WI* (ABS(WI)+ABS(U-ABS(WI)))/ABS(WI)\r
-               C\r
-               C  4.REFINE THE SOLUTION BY NEWTON'S ITERATION:\r
-                     DO 50 IT = 1,15\r
-                         WFN = ZZ - ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
-                    +          PHI1,NPTQ,QWORK,IOPT)\r
-                         WI = WI + WFN/ (C*WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1))\r
-                         IF (ABS(WI).GT.1.D0) WI = WI/ (ABS(WI)+ABS(1.D0-ABS(WI)))\r
-                         IF (ABS(WI).LT.U) WI = U*WI* (ABS(WI)+ABS(U-ABS(WI)))/ABS(WI)\r
-                         IF (ABS(WFN).LT.EPS) GO TO 60\r
-                  50 CONTINUE\r
-               C\r
-               C  THE ITERATION FAILED TO MEET THE TOLERANCE IN 15 ITERATIONS.\r
-               C  TRY A DIFFERENT VERTEX AS A REFERENCE POINT:\r
-                     ZZ1 = (1.D0,1.D0) + ZZ\r
-                     GO TO 70\r
-\r
-                  60 WDSC = WI\r
-               C\r
-               C  5.VERIFICATION:\r
-                     ZZ1 = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ,\r
-                    +      QWORK,1)\r
-                  70 IF (ABS(WI).GE.U .AND. ABS(WI).LE.1.D0 .AND.\r
-                    +    ABS(ZZ-ZZ1).LE.1.D-3) GO TO 90\r
-                     IF (INZ.EQ.0) THEN\r
-                         ZS0(KNZ) = ZZ\r
-\r
-                     ELSE\r
-                         ZS1(KNZ) = ZZ\r
-                     END IF\r
-\r
-               } // while (true)*/\r
+                   ZS = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ,QWORK,1);\r
+               \r
+                   // 3.SOLVE ODE INITIAL VALUE PROBLEM (ALONG LINE SEGMENT FROM\r
+                   //   ZS TO ZZ) BY EULER'S SCHEME TO GENERATE THE INITIAL GUESS:\r
+                   for (K = 1; K <= 20; K++) {\r
+                         wout = WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1);\r
+                         scm.zmlt(20.0*C[0], 20.0*C[1], wout[0], wout[1], cr, ci);\r
+                         scm.zdiv(ZZ[0]-ZS[0], ZZ[1]-ZS[1], cr[0], ci[0], cr2, ci2);\r
+                         WI[0] = WI[0] + cr2[0];\r
+                         WI[1] = WI[1] + ci2[0];\r
+                   } // for (K = 1; K <= 20; K++)\r
+                   abswi = scm.zabs(WI[0], WI[1]);\r
+                   if (abswi > 1.0) {\r
+                       denom = abswi + Math.abs(1.0 - abswi);\r
+                       WI[0] = WI[0]/denom;\r
+                       WI[1] = WI[1]/denom;\r
+                   }\r
+                   if (abswi < U) {\r
+                       factor = U*(abswi + Math.abs(U - abswi))/abswi;\r
+                       WI[0] = factor * WI[0];\r
+                       WI[1] = factor * WI[1];\r
+                   }\r
+               \r
+                   //  4.REFINE THE SOLUTION BY NEWTON'S ITERATION:\r
+                   for (IT = 1; IT <= 15; IT++) {\r
+                       zout = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
+                                            PHI1,NPTQ,QWORK,IOPT);\r
+                       WFN[0] = ZZ[0] - zout[0];\r
+                       WFN[1] = ZZ[1] - zout[1];\r
+                       wout = WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1);\r
+                       scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
+                       scm.zdiv(WFN[0], WFN[1], cr[0], ci[0], cr2, ci2);\r
+                       WI[0] = WI[0] + cr2[0];\r
+                       WI[1] = WI[1] + ci2[0];\r
+                       abswi = scm.zabs(WI[0], WI[1]);\r
+                       if (abswi > 1.0) {\r
+                               denom = abswi + Math.abs(1.0 - abswi);\r
+                               WI[0] = WI[0]/denom;\r
+                               WI[1] = WI[1]/denom;\r
+                           }\r
+                           if (abswi < U) {\r
+                               factor = U*(abswi + Math.abs(U - abswi))/abswi;\r
+                               WI[0] = factor * WI[0];\r
+                               WI[1] = factor * WI[1];\r
+                           }\r
+                       if (scm.zabs(WFN[0],WFN[1]) < EPS) {\r
+                               result[0] = WI[0];\r
+                               result[1] = WI[1];\r
+                               ZZ1 = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ, QWORK,1);\r
+                               break;\r
+                       }\r
+                       if (IT == 15) {\r
+                               // THE ITERATION FAILED TO MEET THE TOLERANCE IN 15 ITERATIONS.\r
+                               // TRY A DIFFERENT VERTEX AS A REFERENCE POINT:\r
+                               ZZ1[0] = 1.0 + ZZ[0];\r
+                               ZZ1[1] = 1.0 + ZZ[1];   \r
+                       }\r
+                   } // for (IT = 1; IT <= 15; IT++)\r
+                   abswi = scm.zabs(WI[0], WI[1]);\r
+                   absz = scm.zabs(ZZ[0]-ZZ1[0], ZZ[1]-ZZ1[1]);\r
+                   if (abswi >= U && abswi <= 1.0 && absz <= 1.0E-3) {\r
+                       return result;  \r
+                   }\r
+                   if (INZ[0] == 0) {\r
+                       ZS0[KNZ[0]-1][0] = ZZ[0];\r
+                ZS0[KNZ[0]-1][1] = ZZ[1];\r
+                   }\r
+                   else {\r
+                       ZS1[KNZ[0]-1][0] = ZZ[0];\r
+                       ZS1[KNZ[0]-1][1] = ZZ[1];\r
+                   }\r
+\r
+               } // while (true)\r
 \r
-               return result;\r
        }\r
 \r
 \r
@@ -1742,6 +1961,82 @@ double neweps;
                return;\r
        }\r
        \r
+       private double ARGUM(double U1, double W01K[]) {\r
+           // ----------------------------\r
+           // COMPUTES THE ARGUMENT OF VECTOR W01K (I.E.,A COMPLEX NUMBER)\r
+           // WHOSE LENGTH IS U1.\r
+           // NOTE: BE SURE THAT THE SPECIFIED VALUE FOR\r
+           // U1 MUST BE DOUBLE PRECISION.\r
+       \r
+       \r
+           // .. Scalar Arguments ..\r
+           // DOUBLE COMPLEX W01K\r
+           // DOUBLE PRECISION U1\r
+\r
+           // .. Local Scalars ..\r
+           double PI,REW;\r
+           double result;\r
+       \r
+           PI = Math.PI;\r
+           REW = W01K[0];\r
+           if (W01K[1] >= 0.0) {\r
+               result = Math.acos(REW/U1);\r
+           }\r
+           else {\r
+               result = Math.acos(-REW/U1) + PI;\r
+           }\r
+           return result;\r
+       }\r
+\r
+       \r
+       private void NEARW(int M,int N, double W0[][], double W1[][],\r
+                       double ALFA0[], double W[],int KNEAR[],int INEAR[]) {\r
+           //   --------------------------------------------------\r
+           // GIVEN PT W,THIS ROUTINE DETERMINES THE NEAREST PREVERTEX TO\r
+        // THE PT W.ON RETURN INTEGER'INEAR'INDICATES THAT THE NEAREST\r
+           // PT IS FOUND EITHER IN W0(INEAR=0) OR IN W1(INEAR=1).  KNEAR\r
+           // CONTAINS THE CORRESPONDING INDEX IN W0 OR 1N W1.\r
+           // NOTE: THE PREVERTICES CORRESPONDING TO INFINITE VERTICES\r
+           //     WILL BE SKIPPED (OUTER CIRCLE ONLY).\r
+           // SEE ROUTINE DSCSOLV FOR THE REST OF THE CALLING SEQUENCE.\r
+       \r
+       \r
+           // .. Scalar Arguments ..\r
+           // DOUBLE COMPLEX W\r
+           // INTEGER INEAR,KNEAR,M,N\r
+\r
+           // .. Array Arguments ..\r
+           // DOUBLE COMPLEX W0(M),W1(N)\r
+           // DOUBLE PRECISION ALFA0(M)\r
+       \r
+           // .. Local Scalars ..\r
+           double D,DIST;\r
+           int I;\r
+       \r
+           DIST = 2.0;\r
+           for (I = 1; I <= M; I++) {\r
+               D = scm.zabs(W[0]-W0[I-1][0], W[1]-W0[I-1][1]);\r
+               if (D <= 1.0E-11 || D >= DIST ||\r
+                   ALFA0[I-1] <= 0.0) {\r
+                       continue;\r
+               }\r
+               KNEAR[0] = I;\r
+               DIST = D;\r
+           } // for (I = 1; I <= M; I++)\r
+           INEAR[0] = 0;\r
+           for (I = 1; I <= N; I++) {\r
+               D = scm.zabs(W[0]-W1[I-1][0], W[1]-W1[I-1][1]);\r
+               if (D <= 1.0E-6 || D >= DIST) {\r
+                       continue;\r
+               }\r
+               DIST = D;\r
+               KNEAR[0] = I;\r
+               INEAR[0] = 1;\r
+           } // for (I = 1; I <= N; I++) \r
+           return;\r
+       }\r
+\r
+       \r
        private void NEARZ(int M,int N, double Z0[][], double Z1[][],\r
                        double ALFA0[], double Z[], int KNZ[],int INZ[]) {\r
            //   --------------------------------------------------\r