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

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

index 98bd92e..5c371d1 100644 (file)
@@ -10,6 +10,23 @@ public class DoublyConnectedSC extends AlgorithmBase {
        // his article "A Software Package for Computing Schwarz-Christoffel Conformal\r
        // Transformation for Doubly Connected Polyhgonal Regions."\r
        \r
+       //************* DSCPACK *************************************************\r
+       // THIS IS A COLLECTION OF SUBROUTINES TO SOLVE PARAMETER PROBLEM OF   *\r
+       // THE SCHWARZ-CHRISTOFFEL TRANSFORMATION FOR DOUBLY CONNECTED REGIONS *\r
+    // A DRIVER IS NEEDED FOR A/AN PARTICULAR PROBLEM OR APPLICATION.      *\r
+       // AUTHOR:                                                             *\r
+       // CHENGLIE HU   APRIL,1994 (REVISED JULY,1995) AT WICHITA STATE UNIV. *\r
+       //               A FURTHER REVISION WAS DONE AT FHSU (JULY, 1997)      *\r
+       // REFERENCES:1. H.DAEPPEN, DIE SCHWARZ-CRISTOFFEL ABBILDUNG FUER      *\r
+       //               ZWEIFACH ZUSAMMENHAENGENDE GEBIETE MIT ANWENDUNGEN.   *\r
+       //               PH.D. DISSERTATION, ETH ZUERICH.                      *\r
+       //            2. L.N.TREFETHEN, SCPACK USER'S GUIDE(MIT REPORT 1989)   *\r
+       //            3. HENRICI,APPLIED & COMPUTATIONAL COMPLEX ANALYSIS,VOL.3*\r
+       //            4. C.HU, APPLICATION OF COMPUTATIONAL COMPLEX ANALYSIS TO*\r
+       //               SOME FREE BOUNDARY AND VORTEX FLOWS.PH.D. DISS. 1995  *\r
+       //**********************************************************************\r
+\r
+       \r
        // geometries of the polygon region in the 7 test routines\r
        private final int SQUARE_SYMMETRIC_REGION = 1;\r
        private final int MILDLY_CROWDED_INFINITE_REGION = 2;\r
@@ -32,6 +49,8 @@ public class DoublyConnectedSC extends AlgorithmBase {
        // The number of Gauss-Jacobi points\r
        // Recommended values for NPTQ are 2-8.\r
        private int NPTQ;\r
+       // 1 for solving the nonlinear system, 2 for not solving the nonlinear system\r
+       private int ISOLV;\r
        private boolean testRoutine = false;\r
        private double MACHEP = 2.2204460E-16;\r
        \r
@@ -39,9 +58,10 @@ public class DoublyConnectedSC extends AlgorithmBase {
                \r
        }\r
        \r
-       public DoublyConnectedSC(int IPOLY, int NPTQ) {\r
+       public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV) {\r
                this.IPOLY = IPOLY;\r
                this.NPTQ = NPTQ;\r
+               this.ISOLV = ISOLV;\r
                testRoutine = true;\r
        }\r
        \r
@@ -54,6 +74,15 @@ public class DoublyConnectedSC extends AlgorithmBase {
                double ALFA1[] = new double[30];\r
                double QWORK[] = new double[1660];\r
                int ISHAPE;\r
+               int IGUESS;\r
+               int LINEARC;\r
+               double TOL;\r
+               double U[] = new double[1];\r
+               double C[] = new double[2];\r
+               double W0[][] = new double[30][2];\r
+               double W1[][] = new double[30][2];\r
+               double PHI0[] = new double[30];\r
+               double PHI1[] = new double[30];\r
                \r
 double neweps;\r
                \r
@@ -92,6 +121,17 @@ double neweps;
                                ISHAPE = 1;\r
                        }\r
                        CHECK(ALFA0, ALFA1, M[0], N[0], ISHAPE);\r
+                       \r
+                       // Specify some parameters of the calling sequence of DSCSOLV:\r
+                       IGUESS = 1;\r
+                       LINEARC = 1;\r
+                       TOL = 1.0E-10;\r
+                       \r
+                       // Solve the accessory parameter problem\r
+                       if (ISOLV == 1) {\r
+                               DSCSOLV(TOL, IGUESS, M[0], N[0], U, C, W0, W1, PHI0, PHI1,\r
+                                               Z0, Z1, ALFA0, ALFA1, NPTQ, QWORK, ISHAPE, LINEARC);\r
+                       }\r
                } // if (testRoutine)\r
                \r
                \r
@@ -365,6 +405,583 @@ double neweps;
                } // else               \r
        }\r
        \r
+       private void THDATA(double U[]) {\r
+           //    ----------------------\r
+           //    GENERATES DATA RELATED ONLY TO INNER RADIUS\r
+           //    U AND USED IN COMPUTING THE THETA-FUNCTION.\r
+\r
+               // .. Local Scalars ..\r
+           double PI;\r
+           int  K,N;\r
+       \r
+           //    .. Common blocks ..\r
+           // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+       \r
+           PI = Math.PI;\r
+           if (U[0] >= 0.63) {\r
+               VARY[0] = Math.exp(PI*PI/Math.log(U[0]));\r
+               DLAM = -Math.log(U[0])/PI;\r
+               return;\r
+           }\r
+           if (U[0] < 0.06) {\r
+               IU = 3;\r
+           }\r
+           else if (U[0] < 0.19) {\r
+               IU = 4;\r
+           }\r
+           else if (U[0] < 0.33) {\r
+               IU = 5;\r
+           }\r
+           else if (U[0] < 0.45) {\r
+               IU = 6;\r
+           }\r
+           else if (U[0] < 0.55) {\r
+               IU = 7;\r
+           }\r
+           else {\r
+               IU = 8;\r
+           }\r
+\r
+           for (K = 1; K <= IU; K++) {\r
+                 N = K*K;\r
+                 UARY[K-1] = Math.pow(U[0], N);\r
+           } // for (K = 1; K <= IU; K++)\r
+           return;\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
+               //   ------------------------------------------------------------\r
+               //   CALCULATES THE  COMPLEX  INTEGRAL  FROM WA TO WB ALONG  A\r
+               //   LINE SEGMENT (LINEARC=0)OR A CIRCULAR ARC (LINEARC=1)WITH\r
+               //   POSSIBLE SINGULARITY AT WA, WHERE  KWA IS THE INDEX OF WA\r
+               //   IN W0 (IC=0) OR IN W1 (IC=1). KWA=0 AND IC=2 ( NO OTHER\r
+               //   VALUES PERMITTED ) IF WA IS NOT A PREVERTEX, WHICH THEN\r
+               //   INDICATES THAT ONLY PURE GAUSSIAN QUARATURE IS NEEDED.\r
+               //   PHIA & PHIB  ARE  ARGUMENTS OF  WA & WB  RESPECTIVELY. IF\r
+               //   INTEGRATING ALONG A CIRCULAR ARC,RADIUS SHOULD BE ASSIGNED\r
+               //   TO BE EITHER 1 OR U. ANY VALUE,HOWEVER,CAN BE ASSIGNED TO\r
+               //   RADIUS IF INTEGRATING ALONG A LINE SEGMENT.SEE DOCUMENTATIONS\r
+               //   OF WPROD AND QINIT FOR THE REST OF CALLING SEQUENCE.\r
+               \r
+               //     .. Scalar Arguments ..\r
+               //    DOUBLE COMPLEX WA,WB\r
+               \r
+               //     .. Array Arguments ..\r
+               //      DOUBLE COMPLEX W0(M),W1(N)\r
+               //      DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3))\r
+               \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
+               //     .. 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
+       }\r
+\r
+       \r
+       private double[] WQUAD1(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
+               //   -------------------------------------------------------------\r
+               //   CALCULATES THE COMPLEX INTEGRAL OF WPROD FROM WA TO WB ALONG\r
+               //   EITHER A CIRCULAR ARC OR A LINE-EGMENT.  COMPOUND ONE-SIDED\r
+               //   GAUSS-JACOBI QUDRATURE IS USED.SEE SUBROUTINE WQSUM FOR THE\r
+               //   CALLING SEQUENCE.\r
+               \r
+               //   CHECK FOR ZERO-LENGTH INTEGRAL:\r
+               //     .. Scalar Arguments ..\r
+               // DOUBLE COMPLEX WA,WB\r
+                    \r
+               //     .. Array Arguments ..\r
+               //      DOUBLE COMPLEX W0(M),W1(N)\r
+               //      DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (N+M)+3))\r
+               \r
+               //     .. Local Scalars ..\r
+               double WAA[] = new double[2];\r
+               double WBB[] = new double[2];\r
+               double ZI[] = new double[2];\r
+               //    DOUBLE COMPLEX WAA,WBB,ZI\r
+               double PHAA,PHBB,R;\r
+               double result[] = new double[2];\r
+               double result2[] = new double[2];\r
+               double expb;\r
+               double arg;\r
+                    \r
+               //     .. External Functions ..\r
+               //      DOUBLE COMPLEX WQSUM\r
+               //      DOUBLE PRECISION DIST\r
+               //      EXTERNAL WQSUM,DIST\r
+               \r
+               //     .. Common blocks ..\r
+               //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+       \r
+           if (scm.zabs(WA[0]-WB[0], WA[1]-WB[1]) <= 0.0) {\r
+               result[0] = 0.0;\r
+               result[1] = 0.0;\r
+               return result;\r
+           }\r
+           \r
+           ZI[0] = 0.0;\r
+           ZI[1] = 1.0;\r
+           if (LINEARC != 1) {\r
+               \r
+                   //   LINE SEGMENT INTEGRATION PATH IS CONCERNED BELOW:\r
+                   // STEP1:ONE-SIDED G-J QUADRATURE FOR LEFT ENDPT WA:\r
+                   R = Math.min(1.0,DIST(M,N,W0,W1,WA,KWA,IC)/scm.zabs(WB[0]-WA[0],WB[1]-WA[1]));\r
+                   WAA[0] = WA[0] + R* (WB[0]-WA[0]);\r
+                   WAA[1] = WA[1] + R* (WB[1]-WA[1]); \r
+                   result = WQSUM(WA,0.0,KWA,IC,WAA,0.0,0.0,M,N,U,W0,W1,ALFA0,\r
+                             ALFA1,NPTQ,QWORK,LINEARC);\r
+               \r
+                   //   STEP2:ADJOIN INTERVALS OF PURE GAUSS QUADRATURE IF NECESSARY:\r
+                  while (R != 1.0) {\r
+                     R = Math.min(1.0,DIST(M,N,W0,W1,WAA,0,IC)/scm.zabs(WAA[0]-WB[0], WAA[1] - WB[1]));\r
+                     WBB[0] = WAA[0] + R* (WB[0]-WAA[0]);\r
+                     WBB[1] = WAA[1] + R* (WB[1]-WAA[1]);\r
+                     result2 = WQSUM(WAA,0.0,0,2,WBB,0.0,0.0,M,N,U,W0,W1,\r
+                             ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
+                     result[0] = result[0] + result2[0];\r
+                     result[1] = result[1] + result2[1];\r
+                     WAA[0] = WBB[0];\r
+                     WAA[1] = WBB[1];\r
+                  }\r
+                  return result;\r
+           } // if (LINEARC != 1)\r
+               \r
+               //   CIRCULAR ARC INTEGRATION PATH IS CONCERNED BELOW:\r
+               //   STEP1:ONE-SIDED G-J QUADRATURE FOR LEFT ENDPT WA:\r
+               R = Math.min(1.0,DIST(M,N,W0,W1,WA,KWA,IC)/scm.zabs(WB[0]-WA[0],WB[1]-WA[1]));\r
+               PHAA = PHIA + R* (PHIB-PHIA);\r
+               expb = RADIUS * Math.exp(ZI[0]*PHAA);\r
+               arg = ZI[1] * PHAA;\r
+               WAA[0] = expb * Math.cos(arg);\r
+               WAA[1] = expb * Math.sin(arg);\r
+               result = WQSUM(WA,PHIA,KWA,IC,WAA,PHAA,RADIUS,M,N,U,W0,W1,ALFA0,\r
+                             ALFA1,NPTQ,QWORK,LINEARC);\r
+               \r
+               //   STEP2:ADJOIN INTERVALS OF PURE GAUSS QUADRATURE IF NECESSARY:\r
+               while (R != 1.0) {\r
+                     R = Math.min(1.0,DIST(M,N,W0,W1,WAA,0,IC)/scm.zabs(WAA[0]-WB[0],WAA[1]-WB[1]));\r
+                     PHBB = PHAA + R* (PHIB-PHAA);\r
+                     expb = RADIUS * Math.exp(ZI[0]*PHBB);\r
+                     arg = ZI[1] * PHBB;\r
+                     WBB[0] = expb * Math.cos(arg);\r
+                     WBB[1] = expb * Math.sin(arg);\r
+                     result2 = WQSUM(WAA,PHAA,0,2,WBB,PHBB,RADIUS,M,N,U,W0,W1,\r
+                            ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
+                     result[0] = result[0] + result2[0];\r
+                     result[1] = result[1] + result2[1];\r
+                     PHAA = PHBB;\r
+                     WAA[0] = WBB[0];\r
+                     WAA[1] = WBB[1];\r
+               } //while (R != 1.0)\r
+               return result;\r
+       }\r
+\r
+       \r
+       private double[] WQUAD(double WA[], double PHIA,int KWA,int ICA, double WB[],\r
+                       double PHIB, int KWB,int ICB,double RADIUS,int M,int N, double U, double W0[][],\r
+                   double W1[][], double ALFA0[], double ALFA1[], int NPTQ, double QWORK[],\r
+                   int LINEARC,int IEVL) {\r
+               //   ---------------------------------------------------------------\r
+               //   CALCULATES THE COMPLEX INTEGRAL OF WPROD FROM WA TO WB\r
+               //   ALONG A CIRCULAR ARC OR A LINE-SEGMENT.FUNCTION WQUAD1\r
+               //   IS CALLED FOUR TIMES,ONE FOR EACH 1/4 OF THE INTERVAL.\r
+               //   NOTE:  WQUAD1 ALLOWS  ONLY THE LEFT ENDPOINT  TO  BE A\r
+               //   POSSIBLE SINGULARITY. SEE ROUTINE WQSUM FOR THE CALLING\r
+           //   SEQUENCE.\r
+               \r
+               //     .. Scalar Arguments ..\r
+               //    DOUBLE COMPLEX WA,WB\r
+                    \r
+               //     .. Array Arguments ..\r
+               //     DOUBLE COMPLEX W0(M),W1(N)\r
+               //     DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3))\r
+       \r
+               //     .. Local Scalars \r
+               double WMID[] = new double[2];\r
+               double WMIDA[] = new double[2];\r
+               double WMIDB[] = new double[2];\r
+               double WQA[] = new double[2];\r
+               double WQA1[] = new double[2];\r
+               double WQA2[] = new double[2];\r
+               double WQB[] = new double[2];\r
+               double WQB1[] = new double[2];\r
+               double WQB2[] = new double[2];\r
+               double ZI[] = new double[2];\r
+               // DOUBLE COMPLEX WMID,WMIDA,WMIDB,WQA,WQB,ZI\r
+               double PHMID,PHMIDA,PHMIDB,PI;\r
+               double expb;\r
+               double arg;\r
+               double result[] = new double[2];\r
+               //     .. Common blocks ..\r
+               //     COMMON /PARAM4/UARY,VARY,DLAM,IU\r
+               \r
+               PI = Math.PI;\r
+               ZI[0] = 0.0;\r
+               ZI[1] = 1.0;\r
+               \r
+               //   DETERMINE MIDPTS ON A LINE SEGMENT OR ON A CIRCULAR ARC:\r
+               if (LINEARC == 0) {\r
+                       WMID[0] = (WA[0] + WB[0])/2.0;\r
+                       WMID[1] = (WA[1] + WB[1])/2.0;\r
+                   WMIDA[0] = (WA[0] + WMID[0])/2.0;\r
+                   WMIDA[1] = (WA[1] + WMID[1])/2.0;\r
+                   WMIDB[0] = (WB[0] + WMID[0])/2.0;\r
+                   WMIDB[1] = (WB[1] + WMID[1])/2.0;\r
+                   PHMID = 0.0;\r
+                   PHMIDA = 0.0;\r
+                   PHMIDB = 0.0;\r
+               }\r
+               else {\r
+               if (IEVL != 1) {\r
+                       if (PHIB < PHIA) {\r
+                           PHIA = PHIA - 2.0*PI;\r
+                       }\r
+               } // if (IEVL != 1)\r
+                  PHMID = (PHIA+PHIB)/2.0;\r
+                  expb = RADIUS * Math.exp(ZI[0] * PHMID);\r
+                  arg = ZI[1] * PHMID;\r
+                  WMID[0] = expb * Math.cos(arg);\r
+                  WMID[1] = expb * Math.sin(arg);\r
+                  PHMIDA = (PHIA+PHMID)/2.0;\r
+                  expb = RADIUS * Math.exp(ZI[0] * PHMIDA);\r
+                  arg = ZI[1] * PHMIDA;\r
+                  WMIDA[0] = expb * Math.cos(arg);\r
+                  WMIDA[1] = expb * Math.sin(arg);\r
+                  PHMIDB = (PHIB+PHMID)/2.0;\r
+                  expb = RADIUS * Math.exp(ZI[0]*PHMIDB);\r
+                  arg = ZI[1] * PHMIDB;\r
+                  WMIDB[0] = expb * Math.cos(arg);\r
+                  WMIDB[1] = expb * Math.sin(arg);\r
+               }\r
+               \r
+               //   COMPOUND GAUSS-JACOBI PROCESS ACCORDING TO ONE-QUATER RULE:\r
+        WQA1 = WQUAD1(WA,PHIA,KWA,ICA,WMIDA,PHMIDA,RADIUS,M,N,U,W0,W1,\r
+           ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
+        WQA2 = WQUAD1(WMID,PHMID,0,2,WMIDA,PHMIDA,RADIUS,M,N,U,W0,W1,ALFA0,\r
+           ALFA1,NPTQ,QWORK,LINEARC);\r
+        WQA[0] = WQA1[0] - WQA2[0];\r
+        WQA[1] = WQA1[1] - WQA2[1];\r
+        WQB1 = WQUAD1(WB,PHIB,KWB,ICB,WMIDB,PHMIDB,RADIUS,M,N,U,W0,W1,\r
+           ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
+        WQB2 = WQUAD1(WMID,PHMID,0,2,WMIDB,PHMIDB,RADIUS,M,N,U,W0,W1,ALFA0,\r
+           ALFA1,NPTQ,QWORK,LINEARC);\r
+        WQB[0] = WQB1[0] - WQB2[0];\r
+        WQB[1] = WQB1[1] - WQB2[1];\r
+        result[0] = WQA[0] - WQB[0];\r
+        result[1] = WQA[1] - WQB[1];\r
+               return result;\r
+       }\r
+\r
+\r
+       \r
+       private void XWTRAN(int M,int N, double X[], double U[],\r
+                       double C[], double W0[][], double W1[][], double PHI0[],\r
+                       double PHI1[]) {\r
+               //   -----------------------------------------------\r
+               //  TRANSFORMS X[K-1](UNCONSTRAINED PARAMETERS) TO ACTUAL\r
+               //  D-SC PARAMETERS:U,C,W0,W1.PHI0 & PHI1 ARE ARGUMENTS\r
+               //  OF THE PREVERTICES CONTAINED IN W0 & W1.\r
+               \r
+               //     .. Scalar Arguments ..\r
+               //      DOUBLE COMPLEX C\r
+                     \r
+               //     .. Array Arguments ..\r
+               //      DOUBLE COMPLEX W0(M),W1(N)\r
+               //      DOUBLE PRECISION PHI0(M),PHI1(N),X(M+N+2)\r
+       \r
+               //     .. Local Scalars ..\r
+             double DPH, PH, PHSUM, PI;\r
+             int  I;\r
+       \r
+             PI = Math.PI;\r
+             if (Math.abs(X[0]) <= 1.0E-14) {\r
+                 U[0] = 0.5;\r
+             }\r
+             else {\r
+                 U[0] = (X[0]-2.0-Math.sqrt(0.9216*X[0]*X[0]+4.0))/ (2.0*X[0]);\r
+                 U[0] = (0.0196*X[0]-1.0)/ (U[0]*X[0]);\r
+             }\r
+\r
+             C[0] = X[1];\r
+             C[1] = X[2];\r
+             if (Math.abs(X[N+2]) <= 1.0E-14) {\r
+                 PHI1[N-1] = 0.0;\r
+             }\r
+             else {\r
+                 PH = (1.0+Math.sqrt(1.0+PI*PI*X[N+2]*X[N+2]))/X[N+2];\r
+                 PHI1[N-1] = PI*PI/PH;\r
+             }\r
+\r
+             DPH = 1.0;\r
+             PHSUM = DPH;\r
+             \r
+             for (I = 1; I <= N - 1; I++) {\r
+                 DPH = DPH/Math.exp(X[2+I]);\r
+                 PHSUM = PHSUM + DPH;\r
+             } // for (I = 1; I <= N - 1; I++)\r
+             DPH = 2.0*PI/PHSUM;\r
+             PHI1[0] = PHI1[N-1] + DPH;\r
+             W1[0][0] = U[0] * Math.cos(PHI1[0]);\r
+             W1[0][1] = U[0] * Math.sin(PHI1[0]);\r
+             W1[N-1][0] = U[0] * Math.cos(PHI1[N-1]);\r
+             W1[N-1][1] = U[0] * Math.sin(PHI1[N-1]);\r
+             PHSUM = PHI1[0];\r
+             for (I = 1; I <= N - 2; I++) {\r
+                 DPH = DPH/Math.exp(X[2+I]);\r
+                 PHSUM = PHSUM + DPH;\r
+                 PHI1[I] = PHSUM;\r
+                 W1[I][0] = U[0] * Math.cos(PHSUM);\r
+                 W1[I][1] = U[0] * Math.sin(PHSUM);\r
+             } // for (I = 1; I <= N - 2; I++)\r
+             DPH = 1.0;\r
+             PHSUM = DPH;\r
+             for (I = 1; I <= M - 1; I++) {\r
+                 DPH = DPH/Math.exp(X[N+2+I]);\r
+                 PHSUM = PHSUM + DPH;\r
+             } // for (I = 1; I <= M - 1; I++)\r
+             DPH = 2.0*PI/PHSUM;\r
+             PHSUM = DPH;\r
+             PHI0[0] = DPH;\r
+             W0[0][0] = Math.cos(DPH);\r
+             W0[0][1] = Math.sin(DPH);\r
+             for (I = 1; I <= M - 2; I++) {\r
+                 DPH = DPH/Math.exp(X[N+2+I]);\r
+                 PHSUM = PHSUM + DPH;\r
+                 PHI0[I] = PHSUM;\r
+                 W0[I][0] = Math.cos(PHSUM);\r
+                 W0[I][1] = Math.sin(PHSUM);\r
+             } // for (I = 1; I <= M - 2; I++)\r
+                 return;\r
+       }\r
+\r
+       \r
+       private void DSCSOLV(double TOL,int IGUESS,int M,int N, double U[],\r
+                       double C[], double W0[][], double W1[][], double PHI0[],\r
+                       double PHI1[], double Z0[][],  double Z1[][], double ALFA0[],\r
+                  double ALFA1[], int NPTQ, double QWORK[], int ISHAPE, int LINEARC) {\r
+       //  -----------------------------------------------------------------\r
+       //  SOLVES THE NONLINEAR SYSTEM FOR D-SC PARAMETERS.\r
+       //  CALLING SEQUENCE:\r
+       //  TOL       A TOLERANCE TO CONTROL THE CONVERGENCE IN HYBRD\r
+       //  IGUESS    (=0 )A NON-EQUALLY SPACED INITIAL GUESS OR(=1)THE\r
+       //            OTHER  EQUALLY-SPACED  INITIAL GUESS PROVIDED, OR\r
+       //            (=2)USER-SUPPLIED INITIAL GUESS WHICH IS REQUIRED\r
+       //            TO BE THE ARGUMENTS OF THE INITIAL PREVERTICES.\r
+       //            ROUTINE ARGUM MAY BE USED FOR COMPUTING ARGUMENTS\r
+       //            IF NEEDED. NOTE: C WILL BE COMPUTED IN THE ROUTINE\r
+       //            (NOT SUPPLIED!)\r
+       //  M,N,U,C,W0,W1,PHI0,PHI1,ALFA0,ALFA1,Z0,Z1\r
+       //            CONSTITUTE THE GEOMETRY OF THE POLYGONAL REGION\r
+       //            AND THE MAPPING FUNCTION. ON RETURN U,C,W0,& W1\r
+       //            WILL  CONTAIN COMPUTED PARAMETERS. (PHI0 & PHI1\r
+       //            WILL CONTAIN THE ARGUMENTS OF THE PREVERTICES.)\r
+       //  QWORK     SEE CALLING SEQUENCE DOCUMENTATION IN QINIT.THE ARRAY\r
+       //            MUST HAVE BEEN FILLED BY QINIT BEFORE CALLING DSCSOLV.\r
+       //  NPTQ      THE NUMBER OF GAUSS-JACOBI POINTS USED\r
+       //  ISHAPE    INDICATES  THAT  OUTER POLYGON  CONTAINS NO INFINITE\r
+       //            VERTICES (ISHAPE=0) OR IT HAS SOME INFINITE VERTICES\r
+       //            (ISHAPE=1).\r
+       //  LINEARC   INTEGER VARIABLE TO CONTROL INTEGRATION PATH.IN PATICULAR:\r
+       //            LINEARC=0:INTEGRATING ALONG LINE SEGMENT WHENEVER POSSIBLE\r
+       //            LINEARC=1:INTEGRATING ALONG CIRCULAR ARC WHENEVER POSSIBLE\r
+       //  THE DISCRIPTION OF ARRAY IND & VARIABLE NSHAPE NEEDED IN DSCFUN:\r
+       //  IND       CONTAINS INDICES  CORRESPONDING TO  THOSE INFINITE\r
+       //            VERTICES,BUT THE FIRST & THE LAST ENTRIS MUST BE -1\r
+       //           & M-2(M IS THE # OF VERTICES ON THE OUTER POLYGON).\r
+       //  NSHAPE    THE DIMENSION OF THE INTERGE ARRAY IND(NSHAPE)\r
+       \r
+       \r
+       //     .. Scalar Arguments ..\r
+       //      DOUBLE COMPLEX C\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* (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 C1[] = new double[2];\r
+               double WINT[] = new double[2];\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
+               \r
+               // Local arrays\r
+               double DIAG[] = new double[42];\r
+               double FJAC[][] = new double[42][42];\r
+               double FVAL[] = new double[42];\r
+               double QW[] = new double[1114];\r
+               double X[] = new double[42];\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
+               ZI[0] = 0.0;\r
+               ZI[1] = 1.0;\r
+               PI = Math.PI;\r
+               ICOUNT = 0;\r
+               if (ISHAPE != 0) {\r
+                       NSHAPE = 1;\r
+                       for (I = 2; I <= M-1; I++) {\r
+                           if (ALFA0[I-1] >= 0.0) {\r
+                               continue;\r
+                           }\r
+                           NSHAPE = NSHAPE + 1;\r
+                           IND[NSHAPE-1] = I-1;\r
+                       } // for (I = 2; I <= M-1; I++)\r
+                       IND[0] = -1;\r
+                       NSHAPE = NSHAPE + 1;\r
+                       IND[NSHAPE-1] = M-2;\r
+               } // if (ISHAPE != 0)\r
+               \r
+               // Fix one prevertex\r
+               W0[M-1][0] = 1.0;\r
+               W0[M-1][1] = 0.0;\r
+               PHI0[M-1] = 0.0;\r
+               \r
+               // Following two value assignments are to satisfy the compiler WATFOR77:\r
+               X[1] = 0.0;\r
+               X[2] = 0.0;\r
+               // Initial guess (IGUESS = 0):\r
+               if (IGUESS == 0) {\r
+                   X[0] = 1.0/0.5 - 1.0/0.46;\r
+                   AVE = 2.0 * PI/(double)N;\r
+                   for (I = 1; I <= N-2; I++) {\r
+                       X[2 + I] = Math.log((AVE + 0.0001*I)/(AVE + 0.0001*(I+1)));\r
+                   } // for (I = 1; I <= N-2; I++)\r
+                   X[N+1] = Math.log((AVE+0.0001*(N-1))/(2.0*PI-(N-1)* (AVE+N*0.00005)));\r
+                   X[N+2] = 1.0/ (4.0-0.1) - 1.0/ (4.0+0.1);\r
+                   AVE = 2.0*PI/M;\r
+                   for (I = 1; I <= M - 2; I++) {\r
+                       X[N+2+I] = Math.log((AVE+0.0001*(I-1))/(AVE+0.0001*I));\r
+                   } // for (I = 1; I <= M - 2; I++)\r
+                   X[M+N+1] = Math.log((AVE+0.0001*(M-2))/(2.0*PI-(M-1)* (AVE+(M-2)*0.00005)));\r
+               } // if (IGUESS == 0)\r
+               else if (IGUESS == 1) {\r
+                       //  INITIAL GUESS (IGUESS=1):\r
+                   X[0] = 1.0/0.53 - 1.0/0.43;\r
+                   for (I = 1; I <= N - 1; I++) {\r
+                       X[2+I] = 0.0;\r
+                   } // for (I = 1; I <= N - 1; I++)\r
+                   X[N+2] = 1.0/ (4.0-0.1) - 1.0/ (4.0+0.1);\r
+                   for (I = 1; I <= M - 1; I++) {\r
+                       X[N+2+I] = 0.0;\r
+                   } // for (I = 1; I <= M - 1; I++)\r
+               } // else if (IGUESS == 1)\r
+               else {\r
+                       X[0] = 1.0/ (0.98-U[0]) - 1.0/ (U[0]-0.02);\r
+                       for (K = 1; K <= N - 1; K++) {\r
+                           KN = K - 1;\r
+                           if (KN == 0) {\r
+                               KN = N;\r
+                           }\r
+                           TOP = PHI1[K-1] - PHI1[KN-1];\r
+                           if (TOP < 0.0) {\r
+                               TOP = TOP + 2.0*PI;\r
+                           }\r
+                           BOTM = PHI1[K] - PHI1[K-1];\r
+                           if (BOTM < 0.0) {\r
+                               BOTM = BOTM + 2.0*PI;\r
+                           }\r
+                           X[2+K] = Math.log(TOP) - Math.log(BOTM);\r
+                       } // for (K = 1; K <= N - 1; K++)\r
+                       X[N+2] = 1.0/ (PI-PHI1[N-1]) - 1.0/ (PI+PHI1[N-1]);\r
+                       for (K = 1; K <= M - 1; K++) {\r
+                           KM = K - 1;\r
+                           if (KM == 0) {\r
+                               KM = M;\r
+                           }\r
+                           TOP = PHI0[K-1] - PHI0[KM-1];\r
+                           if (TOP < 0.0) {\r
+                               TOP = TOP + 2.0*PI;\r
+                           }\r
+                           BOTM = PHI0[K] - PHI0[K-1];\r
+                           if (BOTM < 0.0) {\r
+                               BOTM = BOTM + 2.0*PI;\r
+                           }\r
+                           X[N+2+K] = Math.log(TOP) - Math.log(BOTM);\r
+                       } // for (K = 1; K <= M - 1; K++)\r
+               } // else\r
+               \r
+               // CALCULATE THE INITIAL GUESS X[1] & X[2] TO MATCH\r
+               // THE CHOICE FOR X[0],X[3],...,X[M+N+1]:\r
+               XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1);\r
+        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
+       }\r
+\r
+       \r
        private void ANGLES(int MN, double Z01[][], double ALFA01[], int I01) {\r
            // Computes the interior angles of a doubly\r
                // connected and bounded polygonal region.\r
@@ -397,6 +1014,45 @@ double neweps;
                return;\r
        }\r
        \r
+       private double DIST(int M,int N, double W0[][], double W1[][],\r
+                       double W[], int KWA, int IC) {\r
+       //   -----------------------------------\r
+       //    DETERMINES THE DISTANCE FROM W TO THE NEAREST SINGULARITY\r
+       //    OTHER THAN W ITSELF.( W COULD BE ONE OF THE PREVERTICES.)\r
+       //    KWA IS THE INDEX OF W IN W0 (IF IC=0) OR W1 (IF IC=1), OR\r
+       //    WK COULD BE 0 (IF ONE CHOOSES) WHEN W IS NOT A PREVERTEX.\r
+    //\r
+       //     .. Scalar Arguments ..\r
+       //     DOUBLE COMPLEX W\r
+       \r
+       //     .. Array Arguments ..\r
+       //      DOUBLE COMPLEX W0(M),W1(N)\r
+               \r
+       //     .. Local Scalars ..\r
+             double D;\r
+             int I;\r
+             double result;\r
+       \r
+             result = 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 (I == KWA && IC == 0) {\r
+                         continue;\r
+                 }\r
+                 result = Math.min(result,D);\r
+             } // for (I = 1; I <= M; I++)\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 (I == KWA && IC == 1) {\r
+                         continue;\r
+                 }\r
+                 result = Math.min(result,D);\r
+             } // for (I = 1; I <= N; I++)\r
+             return result;\r
+\r
+       }\r
+\r
+       \r
        private void CHECK(double ALFA0[], double ALFA1[], int M, int N, int ISHAPE) {\r
                //   CHECKS IF THE INPUT DATA AND PARAMETERS ARE CORRECT.\r
                //   NOTE1:    ANGLE-CHECKING MAKES SENSE ONLY IF ANGLES\r
@@ -409,6 +1065,8 @@ double neweps;
                //   3. COUNTERCLOCKWISE ORIENTATION OF THE VERTICES.\r
                \r
                // ALFA0(M), ALFA1(N)\r
+               double EPS, SUM;\r
+               int K, MK;\r
                \r
                if (!((M >= 3) && (M <= 30) && (N <= 30) && (M+N <= 40))) {\r
                        System.err.println("M must be no less than 3");\r
@@ -417,6 +1075,62 @@ double neweps;
                        System.exit(-1);\r
                }\r
        \r
+               EPS = 0.00001;\r
+               SUM = 0.0;\r
+               MK = 0;\r
+               for (K = 0; K <= M; K++) {\r
+                       if (ALFA0[K] > 0.0) {\r
+                               MK = MK + 1;\r
+                       }\r
+                       SUM = SUM + ALFA0[K];\r
+               } // for (K = 0; K <= M; K++)\r
+               if (!((MK < M && ISHAPE == 1) || (MK == M && ISHAPE == 0))) {\r
+                       System.err.println("For finite regions ISHAPE must be 0");\r
+                       System.err.println("and for infinite regions ISHAPE must be 1");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               if (Math.abs(SUM - (M-2)) >= EPS) {\r
+                       System.err.println("Some angles for outer polygon are wrong");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               SUM = 0.0;\r
+               for (K = 0; K < N; K++) {\r
+                       SUM = SUM + ALFA1[K];\r
+               }\r
+               if (Math.abs(SUM - (N+2)) >= EPS) {\r
+                       System.err.println("Some angles for inner polygon are wrong");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               if (ALFA0[0] <= 0.0) {\r
+                       System.err.println("Z0[0] must be finite");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               if (ALFA0[M-1] <= 0.0) {\r
+                       System.err.println("Z0[M-1] must be finite");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               if (ALFA0[M-3] <= 0.0) {\r
+                       System.err.println("Z0[M-3] must be finite");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               if (ALFA0[M-2] >= 2.0-EPS) {\r
+                       System.err.println("Z0[M-2] must no be a re-entrant corner");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               if ((ALFA0[M-2] >= 1.0-EPS) && (ALFA0[M-2] <= 1.0+EPS)) {\r
+                       System.err.println("Z0[M-2] must not be an artificial corner");\r
+                       System.exit(-1);\r
+               }\r
+               \r
+               System.out.println("Inputs have been checked with no error being found");\r
+               return;\r
        }\r
        \r
        private void QINIT(int M, int N, double ALFA0[], double ALFA1[], int NPTQ, double QWORK[]) {\r