In WPROD changed W0[K] to W0[K-1] and W1[K] to W1[K-1].
[mipav.git] / mipav / src / gov / nih / mipav / model / algorithms / DoublyConnectedSC.java
1 package gov.nih.mipav.model.algorithms;\r
2 \r
3 import gov.nih.mipav.view.MipavUtil;\r
4 \r
5 public class DoublyConnectedSC extends AlgorithmBase {\r
6         // This is a port of the FORTRAN ACM TOMS Algorithm 785, collected\r
7         // algorithms form ACM.  The original ACM work was published in Transactions\r
8         // on Mathematical Software, Vol. 24, No. 3, September, 1998, pp. 317-333.\r
9         // The original FORTRAN code was written by Dr. Chenglie Hu as described in\r
10         // his article "A Software Package for Computing Schwarz-Christoffel Conformal\r
11         // Transformation for Doubly Connected Polyhgonal Regions."\r
12         \r
13         //************* DSCPACK *************************************************\r
14         // THIS IS A COLLECTION OF SUBROUTINES TO SOLVE PARAMETER PROBLEM OF   *\r
15         // THE SCHWARZ-CHRISTOFFEL TRANSFORMATION FOR DOUBLY CONNECTED REGIONS *\r
16     // A DRIVER IS NEEDED FOR A/AN PARTICULAR PROBLEM OR APPLICATION.      *\r
17         // AUTHOR:                                                             *\r
18         // CHENGLIE HU   APRIL,1994 (REVISED JULY,1995) AT WICHITA STATE UNIV. *\r
19         //               A FURTHER REVISION WAS DONE AT FHSU (JULY, 1997)      *\r
20         // REFERENCES:1. H.DAEPPEN, DIE SCHWARZ-CRISTOFFEL ABBILDUNG FUER      *\r
21         //               ZWEIFACH ZUSAMMENHAENGENDE GEBIETE MIT ANWENDUNGEN.   *\r
22         //               PH.D. DISSERTATION, ETH ZUERICH.                      *\r
23         //            2. L.N.TREFETHEN, SCPACK USER'S GUIDE(MIT REPORT 1989)   *\r
24         //            3. HENRICI,APPLIED & COMPUTATIONAL COMPLEX ANALYSIS,VOL.3*\r
25         //            4. C.HU, APPLICATION OF COMPUTATIONAL COMPLEX ANALYSIS TO*\r
26         //               SOME FREE BOUNDARY AND VORTEX FLOWS.PH.D. DISS. 1995  *\r
27         //**********************************************************************\r
28 \r
29         \r
30         // geometries of the polygon region in the 7 test routines\r
31         private final int SQUARE_SYMMETRIC_REGION = 1;\r
32         private final int MILDLY_CROWDED_INFINITE_REGION = 2;\r
33         private final int HEAVILY_CROWDED_REGION = 3;\r
34         private final int CHINESE_CHARACTER_STRUCTURED_REGION = 4;\r
35         private final int FOUR_DIRECTION_INFINITE_REGION = 5;\r
36         private final int EXAMPLE_GIVEN_FOR_CHECKING_INVERSE_MAP = 6;\r
37         private final int UPPER_HALF_PLANE_WITH_HORIZONTAL_SLIT = 7;\r
38         \r
39         private final int dscfun = 1;\r
40         \r
41         // Below are the 5 variables in the common block PARAM1\r
42         private double W02[][] = new double[30][2];\r
43         private double W12[][] = new double[30][2];\r
44         private double Z02[][] = new double[30][2];\r
45         private double Z12[][] = new double[30][2];\r
46         private double C2[] = new double[2];\r
47         \r
48         // Below are the 6 variables in the common block PARAM2\r
49         private double U2[] = new double[1];\r
50         private double PHI02[] = new double[30];\r
51         private double PHI12[] = new double[30];\r
52         private double ALFA02[] = new double[30];\r
53         private double ALFA12[] = new double[30];\r
54         private double QWORK2[] = new double[1660];\r
55         \r
56         // Below are the 7 variables in the common block PARAM3\r
57         private int M2, N2, NPTQ2, ISHAPE2, LINEARC2,  NSHAPE;\r
58         private int IND[] = new int[20];\r
59         \r
60         // Below are the 4 variables in the common block PARAM4\r
61         private double UARY[] = new double[8];\r
62         private double VARY[] = new double[3];\r
63         private double DLAM;\r
64         private int IU;\r
65         \r
66         // Below are the 2 variables in the common block PARAM5:\r
67         // Screen display of the residual of the system as the iteration goes on, 1 for "YES", 2 for "NO"\r
68         private int ISPRT;\r
69         private int ICOUNT;\r
70         \r
71         // Common blocks ..\r
72     // COMMON /PARAM1/W02,W12,Z02,Z12,C2\r
73     // COMMON /PARAM2/U2,PHI02,PHI12,ALFA02,ALFA12,QWORK2\r
74     // COMMON /PARAM3/M2,N2,NPTQ2,ISHAPE2,LINEARC2,NSHAPE,IND\r
75     // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
76     // COMMON /PARAM5/ISPRT,ICOUNT\r
77         \r
78         private SchwarzChristoffelMapping scm = new SchwarzChristoffelMapping();\r
79         \r
80         // Specifies the geometry of the polygon region for 7 test examples\r
81         private int IPOLY;\r
82         // The number of Gauss-Jacobi points\r
83         // Recommended values for NPTQ are 2-8.\r
84         private int NPTQ;\r
85         // 1 for solving the nonlinear system, 2 for not solving the nonlinear system\r
86         private int ISOLV;\r
87         // Inverse points\r
88         private double INVERSE_POINTS[][] = null; ;\r
89         private boolean testRoutine = false;\r
90         private double MACHEP = 2.2204460E-16;\r
91         \r
92         public DoublyConnectedSC() {\r
93                 \r
94         }\r
95         \r
96         public DoublyConnectedSC(int IPOLY, int NPTQ, int ISOLV, int ISPRT, double INVERSE_POINTS[][]) {\r
97                 this.IPOLY = IPOLY;\r
98                 this.NPTQ = NPTQ;\r
99                 this.ISOLV = ISOLV;\r
100                 this.ISPRT = ISPRT;\r
101                 this.INVERSE_POINTS = INVERSE_POINTS;\r
102                 testRoutine = true;\r
103         }\r
104         \r
105         public void runAlgorithm() {\r
106                 int I;\r
107                 int M[] = new int[1];\r
108                 int N[] = new int[1];\r
109                 double Z0[][] = new double[30][2];\r
110                 double Z1[][] = new double[30][2];\r
111                 double ALFA0[] = new double[30];\r
112                 double ALFA1[] = new double[30];\r
113                 double QWORK[] = new double[1660];\r
114                 int ISHAPE;\r
115                 int IGUESS;\r
116                 int LINEARC;\r
117                 double TOL;\r
118                 double U[] = new double[1];\r
119                 double C[] = new double[2];\r
120                 double W0[][] = new double[30][2];\r
121                 double W1[][] = new double[30][2];\r
122                 double PHI0[] = new double[30];\r
123                 double PHI1[] = new double[30];\r
124                 double ZZ[] = new double[2];\r
125                 double EPS;\r
126                 double WW[];\r
127                 double ZZ0[];\r
128                 \r
129 double neweps;\r
130                 \r
131                 // MACHEP is a machine dependent parameter specifying the relative precision of\r
132                 // floating point arithmetic.\r
133                  // epsilon = D1MACH(4)\r
134         // Machine epsilon is the smallest positive epsilon such that\r
135         // (1.0 + epsilon) != 1.0.\r
136         // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
137         // epsilon = 2.2204460e-16\r
138         // epsilon is called the largest relative spacing\r
139         MACHEP = 1.0;\r
140         neweps = 1.0;\r
141 \r
142         while (true) {\r
143 \r
144             if (1.0 == (1.0 + neweps)) {\r
145                 break;\r
146             } else {\r
147                 MACHEP = neweps;\r
148                 neweps = neweps / 2.0;\r
149             }\r
150         } // while(true)\r
151                 \r
152                 if (testRoutine) {\r
153                         DSCDATA(IPOLY, M, N, Z0, Z1, ALFA0, ALFA1);\r
154                         ANGLES(N[0], Z1, ALFA1, 1);\r
155                         if ((IPOLY == 1) || (IPOLY == 3) || (IPOLY == 4) || (IPOLY == 6)) {\r
156                                 ANGLES(M[0], Z0, ALFA0, 0);\r
157                         }\r
158                         \r
159                         // Generate the Gauss-Jacobi weights & nodes and check the input\r
160                         QINIT(M[0],N[0],ALFA0, ALFA1, NPTQ, QWORK);\r
161                         ISHAPE = 0;\r
162                         if ((IPOLY == 2) || (IPOLY == 5) || (IPOLY == 7)) {\r
163                                 ISHAPE = 1;\r
164                         }\r
165                         CHECK(ALFA0, ALFA1, M[0], N[0], ISHAPE);\r
166                         \r
167                         // Specify some parameters of the calling sequence of DSCSOLV:\r
168                         IGUESS = 1;\r
169                         LINEARC = 1;\r
170                         TOL = 1.0E-10;\r
171                         \r
172                         // Solve the accessory parameter problem\r
173                         if (ISOLV == 1) {\r
174                                 DSCSOLV(TOL, IGUESS, M[0], N[0], U, C, W0, W1, PHI0, PHI1,\r
175                                                 Z0, Z1, ALFA0, ALFA1, NPTQ, QWORK, ISHAPE, LINEARC);\r
176                         }\r
177                         \r
178                         // Output will be arranged in DSCSOLV\r
179                         \r
180                         // Compute data for theta-function and test the accuracy\r
181                         THDATA(U);\r
182                     DSCTEST(M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,NPTQ,QWORK);\r
183                     \r
184             // Test inverse evaluations\r
185                     if ((INVERSE_POINTS != null) && (INVERSE_POINTS.length != 0)) {\r
186                         for (I = 0; I < INVERSE_POINTS.length; I++) {\r
187                             ZZ[0] = INVERSE_POINTS[I][0];       \r
188                             ZZ[1] = INVERSE_POINTS[I][1];\r
189                             EPS = 1.0E-6;\r
190                             WW = WDSC(ZZ,M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,\r
191                                                   NPTQ,QWORK,EPS,1);\r
192                     System.out.println("THE PREIMAGE OF ZZ = (" + WW[0] + ", " + WW[1] + ")");\r
193                     if (scm.zabs(WW[0], WW[1]) <= 1.0E-12) {\r
194                         continue;\r
195                     }\r
196                     System.out.println("CHECK BY MAPPING THE PREIMAGE BACK");\r
197                     ZZ0 = ZDSC(WW,0,2,M[0],N[0],U[0],C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
198                          PHI1,NPTQ,QWORK,1);\r
199                     System.out.println("THE POINT ENTERED = (" + ZZ0[0] + ", " + ZZ0[1] + ")");\r
200 \r
201                         } // for (I = 0; I < INVERSE_POINTS.length; I++)\r
202                     } // if ((INVERSE_POINTS != null) && (INVERSE_POINTS.length != 0))\r
203                 } // if (testRoutine)\r
204                 \r
205                 \r
206         }\r
207         \r
208         private void DSCDATA(int IPOLY, int M[], int N[], double Z0[][], double Z1[][], double ALFA0[], double ALFA1[]) {\r
209                 // DSCDATA generates data.\r
210                 double Q;\r
211                 double S;\r
212                 if (IPOLY == 1) {\r
213                         M[0] = 4;\r
214                         N[0] = 4;\r
215                         Q = Math.sqrt(2.0);\r
216                         Z0[0][0] = 1.0 + Q;\r
217                         Z0[0][1] = 1.0 + Q;\r
218                         Z0[1][0] = -1.0 - Q;\r
219                         Z0[1][1] = 1.0 + Q;\r
220                         Z0[2][0] = -1.0 - Q;\r
221                         Z0[2][1] = -1.0 - Q;\r
222                         Z0[3][0] = 1.0 + Q;\r
223                         Z0[3][1] = -1.0 - Q;\r
224                         Z1[0][0] = Q;\r
225                         Z1[0][1] = 0.0;\r
226                         Z1[1][0] = 0.0;\r
227                         Z1[1][1] = Q;\r
228                         Z1[2][0] = -Q;\r
229                         Z1[2][1] = 0.0;\r
230                         Z1[3][0] = 0.0;\r
231                         Z1[3][1] = -Q;\r
232                 } // if (IPOLY == 1)\r
233                 else if (IPOLY == 2) {\r
234                         M[0] = 12;\r
235                         N[0] = 6;\r
236                         Z0[0][0] = 0.6875;\r
237                         Z0[0][1] = 0.875;\r
238                         Z0[1][0] = 0.0;\r
239                         Z0[1][1] = 0.0;\r
240                         Z0[2][0] = 1.0;\r
241                         Z0[2][1] = 0.0;\r
242                         Z0[3][0] = 0.875;\r
243                         Z0[3][1] = 0.875;\r
244                         Z0[4][0] = 1.125;\r
245                         Z0[4][1] = 1.375;\r
246                         Z0[5][0] = 2.0;\r
247                         Z0[5][1] = 1.375;\r
248                         Z0[6][0] = 1.25;\r
249                         Z0[6][1] = 2.0;\r
250                         Z0[7][0] = 2.25;\r
251                         Z0[7][1] = 2.0;\r
252                         Z0[8][0] = 2.375;\r
253                         Z0[8][1] = 2.75;\r
254                         Z0[9][0] = 1.625;\r
255                         Z0[9][1] = 2.25;\r
256                         Z0[10][0] = 1.125;\r
257                     Z0[10][1] = 2.625;\r
258                     Z0[11][0] = -0.5;\r
259                     Z0[11][1] = 2.75;\r
260                     Z1[0][0] = 0.375;\r
261                     Z1[0][1] = 1.875;\r
262                     Z1[1][0] = 0.5;\r
263                     Z1[1][1] = 2.0;\r
264                     Z1[2][0] = 1.0;\r
265                     Z1[2][1] = 1.5;\r
266                     Z1[3][0] = 0.5;\r
267                     Z1[3][1] = 2.1875;\r
268                     Z1[4][0] = 0.5;\r
269                     Z1[4][1] = 2.5;\r
270                     Z1[5][0] = Z1[3][0];\r
271                     Z1[5][1] = Z1[3][1];\r
272                     ALFA0[0] = 1.39169261159339475;\r
273             ALFA0[1] = 0.28801540784794967;\r
274             ALFA0[2] = 0.454832764699133488;\r
275             ALFA0[3] = 1.19275085295129979;\r
276             ALFA0[4] = 1.35241638234956651;\r
277             ALFA0[5] = 0.0;\r
278             ALFA0[6] = 2.0;\r
279             ALFA0[7] = 0.552568456711253445;\r
280             ALFA0[8] = 0.260264501477747753;\r
281             ALFA0[9] = 1.39199980651013222;\r
282             ALFA0[10] = 0.819604487273064009;\r
283             ALFA0[11] = 0.295854728586457991;\r
284             INVERSE_POINTS = new double[1][2];\r
285             INVERSE_POINTS[0][0] = 0.5;\r
286             INVERSE_POINTS[0][1] = 0.5;\r
287             // Gives:\r
288             // THE PREIMAGE OF ZZ=          (-0.623005768209842,0.782172159414472)\r
289                 } // else if (IPOLY == 2)\r
290                 else if (IPOLY == 3) {\r
291                         M[0] = 11;\r
292                         N[0] = 6;\r
293                         Z0[0][0] = 0.5;\r
294                         Z0[0][1] = 2.5;\r
295                         Z0[1][0] = 0.5;\r
296                         Z0[1][1] = 0.5;\r
297                         Z0[2][0] = 1.0;\r
298                         Z0[2][1] = 0.5;\r
299                         Z0[3][0] = 1.0;\r
300                         Z0[3][1] = 1.0;\r
301                         Z0[4][0] = 1.0;\r
302                         Z0[4][1] = 0.5;\r
303                         Z0[5][0] = 0.5;\r
304                         Z0[5][1] = 0.5;\r
305                         Z0[6][0] = 0.5;\r
306                         Z0[6][1] = 2.5;\r
307                         Z0[7][0] = 0.0;\r
308                         Z0[7][1] = 2.5;\r
309                         Z0[8][0] = 0.0;\r
310                         Z0[8][1] = 0.0;\r
311                         Z0[9][0] = 2.0;\r
312                         Z0[9][1] = 0.0;\r
313                         Z0[10][0] = 2.0;\r
314                         Z0[10][1] = 2.5;\r
315                         Z1[0][0] = 1.0;\r
316                         Z1[0][1] = 2.0;\r
317                         Z1[1][0] = 1.0;\r
318                         Z1[1][1] = 1.5;\r
319                         Z1[2][0] = 1.0;\r
320                         Z1[2][1] = 2.0;\r
321                         Z1[3][0] = 1.5;\r
322                         Z1[3][1] = 2.0;\r
323                         Z1[4][0] = 1.5;\r
324                         Z1[4][1] = 0.5;\r
325                         Z1[5][0] = 1.5;\r
326                         Z1[5][0] = 2.0;\r
327                         ALFA0[0] = 1.0/2.0;\r
328             ALFA0[1] = 1.0/2.0;\r
329             ALFA0[2] = 1.0/2.0;\r
330             ALFA0[3] = 2.0;\r
331             ALFA0[4] = 3.0/2.0;\r
332             ALFA0[5] = 3.0/2.0;\r
333             ALFA0[6] = 1.0/2.0;\r
334             ALFA0[7] = 1.0/2.0;\r
335             ALFA0[8] = 1.0/2.0;\r
336             ALFA0[9] = 1.0/2.0;\r
337             ALFA0[10] = 1.0/2.0;\r
338             ALFA1[0] = 3.0/2.0;\r
339             ALFA1[1] = 2.0;\r
340             ALFA1[2] = 1.0/2.0;\r
341             ALFA1[3] = 1.0/2.0;\r
342             ALFA1[4] = 2.0;\r
343             ALFA1[5] = 3.0/2.0;\r
344                 } // else if (IPOLY == 3)\r
345                 else if (IPOLY == 4) {\r
346                     M[0] = 4;\r
347                     N[0] = 17;\r
348                     Z0[0][0] = -1.0;\r
349                     Z0[0][1] = -1.0;\r
350                     Z0[1][0] = 1.0;\r
351                     Z0[1][1] = -1.0;\r
352                     Z0[2][0] = 1.0;\r
353                     Z0[2][1] = 1.0;\r
354                     Z0[3][0] = -1.0;\r
355                     Z0[3][1] = 1.0;\r
356                     Z1[0][0] = 0.0;\r
357                     Z1[0][1] = 0.5;\r
358                     Z1[1][0] = 0.0;\r
359                     Z1[1][1] = 0.0;\r
360                     Z1[2][0] = -0.5;\r
361                     Z1[2][1] = 0.0;\r
362                     Z1[3][0] = Z1[1][0];\r
363                     Z1[3][1] = Z1[1][1];\r
364                     Z1[4][0] = 0.0;\r
365                     Z1[4][1] = -0.5;\r
366                     Z1[5][0] = -0.5;\r
367                     Z1[5][1] = -0.5;\r
368                     Z1[6][0] = 0.5;\r
369                     Z1[6][1] = -0.5;\r
370                     Z1[7][0] = 0.25;\r
371                     Z1[7][1] = -0.5;\r
372                     Z1[8][0] = 0.25;\r
373                     Z1[8][1] = -0.25;\r
374                     Z1[9][0] = Z1[7][0];\r
375                     Z1[9][1] = Z1[7][1];\r
376                     Z1[10][0] = Z1[4][0];\r
377                     Z1[10][1] = Z1[4][1];\r
378                     Z1[11][0] = Z1[1][0];\r
379                     Z1[11][1] = Z1[1][1];\r
380                     Z1[12][0] = 0.5;\r
381                     Z1[12][1] = 0.0;\r
382                     Z1[13][0] = Z1[1][0];\r
383                     Z1[13][1] = Z1[1][1];\r
384                     Z1[14][0] = Z1[0][0];\r
385                     Z1[14][1] = Z1[0][1];\r
386                     Z1[15][0] = 0.5;\r
387                     Z1[15][1] = 0.5;\r
388                     Z1[16][0] = -0.5;\r
389                     Z1[16][1] = 0.5;\r
390                 } // else if (IPOLY == 4)\r
391                 else if (IPOLY == 5) {\r
392                         M[0] = 12;\r
393                         N[0] = 4;\r
394                         S = 1.0/6.0;\r
395                         Z0[0][0] = 0.0;\r
396                         Z0[0][1] = -8.0*S;\r
397                         Z0[1][0] = 0.0;\r
398                         Z0[1][1] = 0.0;\r
399                         Z0[2][0] = 1.0;\r
400                         Z0[2][1] = -1.5;\r
401                         Z0[3][0] = 1.0;\r
402                         Z0[3][1] = -0.5;\r
403                         Z0[4][0] = 0.0;\r
404                         Z0[4][1] = 0.0;\r
405                         Z0[5][0] = 0.0;\r
406                         Z0[5][1] = 2.0*S;\r
407                         Z0[6][0] = 0.0;\r
408                         Z0[6][1] = 0.0;\r
409                         Z0[7][0] = -1.0;\r
410                         Z0[7][1] = 0.0;\r
411                         Z0[8][0] = 0.0;\r
412                         Z0[8][1] = 0.0;\r
413                         Z0[9][0] = -7.0*S;\r
414                         Z0[9][1] = -1.0;\r
415                         Z0[10][0] = 0.0;\r
416                         Z0[10][1] = 0.0;\r
417                         Z0[11][0] = -2.0*S;\r
418                         Z0[11][1] = -10.0*S;\r
419                         ALFA0[0] = 1.75;\r
420             ALFA0[1] = 0.0;\r
421             ALFA0[2] = 1.0;\r
422             ALFA0[3] = 1.0;\r
423             ALFA0[4] = (Math.PI-3.5)/Math.PI;\r
424             ALFA0[5] = 3.5/Math.PI;\r
425             ALFA0[6] = 1.5;\r
426             ALFA0[7] = 1.0;\r
427             ALFA0[8] = 0.0;\r
428             ALFA0[9] = 1.75;\r
429             ALFA0[10] = 0.0;\r
430             ALFA0[11] = 1.0;\r
431             Z1[0][0] = 2.0*S;\r
432             Z1[0][1] = -0.5;\r
433             Z1[1][0] = -S;\r
434             Z1[1][1] = -2.0*S;\r
435             Z1[2][0] = -4.0*S;\r
436             Z1[2][1] = -5.0*S;\r
437             Z1[3][0] = 0.0;\r
438             Z1[3][1] = -1.0;\r
439                 } // else if (IPOLY == 5)\r
440                 else if (IPOLY == 6) {\r
441                     M[0] = 7;\r
442                     N[0] = 2;\r
443                     Z0[0][0] = -2.0;\r
444                     Z0[0][1] = -1.0;\r
445                     Z0[1][0] = 2.0;\r
446                     Z0[1][1] = -1.0;\r
447                     Z0[2][0] = 2.0;\r
448                     Z0[2][1] = 2.0;\r
449                     Z0[3][0] = -0.8;\r
450                     Z0[3][1] = 2.0;\r
451                     Z0[4][0] = 1.0;\r
452                     Z0[4][1] = 0.5;\r
453                     Z0[5][0] = -1.0;\r
454                     Z0[5][1] = 2.0;\r
455                     Z0[6][0] = -2.0;\r
456                     Z0[6][1] = 2.0;\r
457                     Z1[0][0] = 0.0;\r
458                     Z1[0][1] = 0.0;\r
459                     Z1[1][0] = -1.0;\r
460                     Z1[1][1] = 0.0;\r
461                 } // else if (IPOLY == 6)\r
462                 else {\r
463                         M[0] = 3;\r
464                         N[0] = 2;\r
465                         Z0[0][0] = 1.01;\r
466                         Z0[0][1] = 0.0;\r
467                         Z0[1][0] = 100.0;\r
468                         Z0[1][1] = 100.0;\r
469                         Z0[2][0] = -1.01;\r
470                         Z0[2][1] = 0.0;\r
471                         Z1[0][0] = 0.0;\r
472                         Z1[0][1] = 2.0;\r
473                         Z1[1][0] = 0.0;\r
474                         Z1[1][1] = 1.0;\r
475                         ALFA0[0] = 1.0;\r
476                         ALFA0[1] = -1.0;\r
477                         ALFA0[2] = 1.0;\r
478                 } // else               \r
479         }\r
480         \r
481         private void THDATA(double U[]) {\r
482             //    ----------------------\r
483             //    GENERATES DATA RELATED ONLY TO INNER RADIUS\r
484             //    U AND USED IN COMPUTING THE THETA-FUNCTION.\r
485 \r
486                 // .. Local Scalars ..\r
487             double PI;\r
488             int  K,N;\r
489         \r
490             //    .. Common blocks ..\r
491             // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
492         \r
493             PI = Math.PI;\r
494             if (U[0] >= 0.63) {\r
495                 VARY[0] = Math.exp(PI*PI/Math.log(U[0]));\r
496                 DLAM = -Math.log(U[0])/PI;\r
497                 return;\r
498             }\r
499             if (U[0] < 0.06) {\r
500                 IU = 3;\r
501             }\r
502             else if (U[0] < 0.19) {\r
503                 IU = 4;\r
504             }\r
505             else if (U[0] < 0.33) {\r
506                 IU = 5;\r
507             }\r
508             else if (U[0] < 0.45) {\r
509                 IU = 6;\r
510             }\r
511             else if (U[0] < 0.55) {\r
512                 IU = 7;\r
513             }\r
514             else {\r
515                 IU = 8;\r
516             }\r
517 \r
518             for (K = 1; K <= IU; K++) {\r
519                   N = K*K;\r
520                   UARY[K-1] = Math.pow(U[0], N);\r
521             } // for (K = 1; K <= IU; K++)\r
522             return;\r
523 \r
524         }\r
525         \r
526         private double[] WTHETA(double U, double W[]) {\r
527             //  -------------------------\r
528             //    EVALUATES THETA-FUNCTION AT W,WHERE U IS THE\r
529             //    INNER RADIUS OF THE ANNULUS.THE DEFINITION OF\r
530             //    THETA-FUNCTION CAN BE FOUND IN REFERENCE 3.\r
531         \r
532             //     .. Scalar Arguments ..\r
533             //  DOUBLE COMPLEX W\r
534 \r
535             //     .. Scalars in Common ..\r
536             // DOUBLE PRECISION DLAM\r
537             // INTEGER IU\r
538 \r
539             //    .. Arrays in Common ..\r
540             // DOUBLE PRECISION UARY(8),VARY(3)\r
541 \r
542             //     .. Local Scalars ..\r
543                 double result[] = new double[2];\r
544                 double cr[] = new double[1];\r
545                 double ci[] = new double[1];\r
546                 double WT[] = new double[2];\r
547                 double WWN[] = new double[2];\r
548                 double WWN0[] = new double[2];\r
549                 double WWP[] = new double[2];\r
550                 double WWP0[] = new double[2];\r
551                 double ZI[] = new double[2];\r
552             // DOUBLE COMPLEX WT,WWN,WWN0,WWP,WWP0,ZI\r
553             double PI;\r
554             int K;\r
555         \r
556             //     .. Common blocks ..\r
557             // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
558         \r
559             WWP[0] = 1.0;\r
560             WWP[1] = 0.0;\r
561             WWN[0] = 1.0;\r
562             WWN[1] = 0.0;\r
563             result[0] = 1.0;\r
564             result[1] = 0.0;\r
565             if (U < 0.63) {\r
566                 WWP0[0] = -W[0];\r
567                 WWP0[1] = -W[1];\r
568                 scm.zdiv(-1.0, 0.0, W[0], W[1], cr, ci);\r
569                 WWN0[0] = cr[0];\r
570                 WWN0[1] = ci[0];\r
571                 for (K = 1; K <= IU; K++) {\r
572                         scm.zmlt(WWP[0], WWP[1], WWP0[0], WWP0[1], cr, ci);\r
573                         WWP[0] = cr[0];\r
574                         WWP[1] = ci[0];\r
575                     scm.zmlt(WWN[0], WWN[1], WWN0[0], WWN0[1], cr, ci);\r
576                     WWN[0] = cr[0];\r
577                     WWN[1] = ci[0];\r
578                     result[0] = result[0] + UARY[K-1]* (WWP[0]+WWN[0]);\r
579                     result[1] = result[1] + UARY[K-1]* (WWP[1]+WWN[1]);\r
580                 } // for (K = 1; K <= IU; K++)\r
581                 return result;\r
582            } // if (U < 0.63)\r
583 \r
584            PI = Math.PI;\r
585            ZI[0] = 0.0;\r
586            ZI[1] = 1.0;\r
587            double logr = Math.log(scm.zabs(-W[0], -W[1]));\r
588            double logi = Math.atan2(-W[1], -W[0]);\r
589            scm.zmlt(-ZI[0], -ZI[1], logr, logi, cr, ci);\r
590            WT[0] = cr[0];\r
591            WT[1] = ci[0];\r
592            if (U < 0.94) {\r
593                   double base = Math.exp(WT[0]/DLAM);\r
594                   double arg = WT[1]/DLAM;\r
595                   WWP[0] = base * Math.cos(arg);\r
596                   WWP[1] = base * Math.sin(arg);\r
597                   scm.zdiv(1.0, 0.0, WWP[0], WWP[1], cr, ci);\r
598               WWN[0] = cr[0];\r
599               WWN[1] = ci[0];\r
600               result[0] = result[0] + VARY[0]* (WWP[0]+WWN[0]); \r
601               result[1] = result[1] + VARY[0]* (WWP[1]+WWN[1]);   \r
602            } // if (U < 0.94)\r
603            scm.zmlt(WT[0], WT[1], WT[0], WT[1], cr, ci);\r
604            double argr = -cr[0]/(4.0 * PI * DLAM);\r
605            double argi = -ci[0]/(4.0 * PI * DLAM);\r
606            double base = Math.exp(argr);\r
607            double expr = base * Math.cos(argi);\r
608            double expi = base * Math.sin(argi);\r
609            scm.zmlt(expr, expi, result[0], result[1], cr, ci);\r
610            result[0] = cr[0]/Math.sqrt(DLAM);\r
611            result[1] = ci[0]/Math.sqrt(DLAM);\r
612            return result;\r
613         }\r
614 \r
615 \r
616         \r
617         private double[] WPROD(double W[], int M, int N, double U,\r
618                         double W0[][], double W1[][], double ALFA0[], double ALFA1[]) {\r
619         //   --------------------------------------------\r
620         //    COMPUTES THE PRODUCT (D-SC INTEGRAND):\r
621         // M                                   N\r
622         //PROD THETA(W/U*W0(K))**(ALFA0(K)-1)*PROD THETA(U*W/W1(K))**(ALFA1(K)-1)\r
623         // K=1                                 K=1\r
624         //    THE CALLING SEQUENCE IS EXPLAINED IN THE ABOVE FORMULA.\r
625         \r
626         //     .. Scalar Arguments ..\r
627         //      DOUBLE COMPLEX W\r
628              \r
629         //     .. Array Arguments ..\r
630         //      DOUBLE COMPLEX W0(M),W1(N)\r
631         //      DOUBLE PRECISION ALFA0(M),ALFA1(N)\r
632         \r
633         //     .. Scalars in Common ..\r
634         //      DOUBLE PRECISION DLAM\r
635         //      INTEGER IU\r
636         \r
637         //     .. Arrays in Common ..\r
638         //      DOUBLE PRECISION UARY(8),VARY(3)\r
639         \r
640         //     .. Local Scalars ..\r
641                 double result[] = new double[2];\r
642                 double WSUM[] = new double[2];\r
643                 double WTH[] = new double[2];\r
644                 double cr[] = new double[1];\r
645                 double ci[] = new double[1];\r
646                 double win[] = new double[2];\r
647                 double wout[];\r
648                 double base;\r
649         //      DOUBLE COMPLEX WSUM,WTH\r
650               int K;\r
651 \r
652         //     .. External Functions ..\r
653         //      DOUBLE COMPLEX WTHETA\r
654         //      EXTERNAL WTHETA\r
655         \r
656         //     .. Common blocks ..\r
657         //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
658               WSUM[0] = 0.0;\r
659               WSUM[1] = 0.0;\r
660               for (K = 1; K <= M; K++) {\r
661                   scm.zdiv(W[0], W[1], U*W0[K-1][0], U*W0[K-1][1], cr, ci);\r
662                   win[0] = cr[0];\r
663                   win[1] = ci[0];\r
664                   wout = WTHETA(U, win);\r
665                   WTH[0] = Math.log(scm.zabs(wout[0], wout[1]));\r
666                   WTH[1] = Math.atan2(wout[1], wout[0]);\r
667                   WSUM[0] = WSUM[0] + (ALFA0[K-1]-1.0)*WTH[0];\r
668                   WSUM[1] = WSUM[1] + (ALFA0[K-1]-1.0)*WTH[1];\r
669               } // for (K = 1; K <= M; K++)\r
670               for (K = 1; K <= N; K++) {\r
671                   scm.zdiv(U*W[0], U*W[1], W1[K-1][0], W1[K-1][1], cr, ci);\r
672                   win[0] = cr[0];\r
673                   win[1] = ci[0];\r
674                   wout = WTHETA(U, win);\r
675                   WTH[0] = Math.log(scm.zabs(wout[0], wout[1]));\r
676                   WTH[1] = Math.atan2(wout[1], wout[0]);\r
677                   WSUM[0] = WSUM[0] + (ALFA1[K-1]-1.0)*WTH[0];\r
678                   WSUM[1] = WSUM[1] + (ALFA1[K-1]-1.0)*WTH[1];\r
679               } // for (K = 1; K <= N; K++)\r
680               base = Math.exp(WSUM[0]);\r
681               result[0] = base * Math.cos(WSUM[1]);\r
682               result[1] = base * Math.sin(WSUM[1]);\r
683               return result;\r
684         }\r
685 \r
686 \r
687         \r
688         private double[] WQSUM(double WA[], double PHIA,int KWA, int IC, double WB[],\r
689                         double PHIB, double RADIUS,int M,int N, double U, double W0[][], double W1[][],\r
690                         double ALFA0[], double ALFA1[], int NPTQ,double QWORK[], int LINEARC) {\r
691                 //   ------------------------------------------------------------\r
692                 //   CALCULATES THE  COMPLEX  INTEGRAL  FROM WA TO WB ALONG  A\r
693                 //   LINE SEGMENT (LINEARC=0)OR A CIRCULAR ARC (LINEARC=1)WITH\r
694                 //   POSSIBLE SINGULARITY AT WA, WHERE  KWA IS THE INDEX OF WA\r
695                 //   IN W0 (IC=0) OR IN W1 (IC=1). KWA=0 AND IC=2 ( NO OTHER\r
696                 //   VALUES PERMITTED ) IF WA IS NOT A PREVERTEX, WHICH THEN\r
697                 //   INDICATES THAT ONLY PURE GAUSSIAN QUARATURE IS NEEDED.\r
698                 //   PHIA & PHIB  ARE  ARGUMENTS OF  WA & WB  RESPECTIVELY. IF\r
699                 //   INTEGRATING ALONG A CIRCULAR ARC,RADIUS SHOULD BE ASSIGNED\r
700                 //   TO BE EITHER 1 OR U. ANY VALUE,HOWEVER,CAN BE ASSIGNED TO\r
701                 //   RADIUS IF INTEGRATING ALONG A LINE SEGMENT.SEE DOCUMENTATIONS\r
702                 //   OF WPROD AND QINIT FOR THE REST OF CALLING SEQUENCE.\r
703                 \r
704                 //     .. Scalar Arguments ..\r
705                 //    DOUBLE COMPLEX WA,WB\r
706                 \r
707                 //     .. Array Arguments ..\r
708                 //      DOUBLE COMPLEX W0(M),W1(N)\r
709                 //      DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3))\r
710                 \r
711                 //     .. Local Scalars ..\r
712                 double result[] = new double[2];\r
713                 double W[] = new double[2];\r
714                 double WC[] = new double[2];\r
715                 double WH[] = new double[2];\r
716                 double ZI[] = new double[2];\r
717                 //     DOUBLE COMPLEX W,WC,WH,ZI\r
718                 double PWC,PWH;\r
719                 int I,IOFFST,IWT1,IWT2;\r
720                 double prod[] = null;\r
721                 double base;\r
722                 double arg;\r
723                 double cr[] = new double[1];\r
724                 double ci[] = new double[1];\r
725 \r
726                 //     .. External Functions ..\r
727                 //      DOUBLE COMPLEX WPROD\r
728                 //      EXTERNAL WPROD\r
729                 \r
730                 //     .. Common blocks ..\r
731                 //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
732             result[0] = 0.0;\r
733                 result[1] = 0.0;\r
734                 \r
735                 //   INDEX ARRANGEMENT:\r
736                 IWT1 = NPTQ* (IC*M+KWA-1) + 1;\r
737                 if (KWA == 0) {\r
738                     IWT1 = NPTQ* (M+N) + 1;\r
739                 }\r
740             IWT2 = IWT1 + NPTQ - 1;\r
741                 IOFFST = NPTQ* (M+N+1);\r
742                 \r
743                 //   COMPUTE GAUSS-JACOBI SUM(W(J)*PROD(X(J))):\r
744                 if (LINEARC != 1) {\r
745                 \r
746                     //   INTEGRATE ALONG A LINE SEGMENT:\r
747                     WH[0] = (WB[0]-WA[0])/2.0;\r
748                     WH[1] = (WB[1]-WA[1])/2.0;\r
749                     WC[0] = (WA[0]+WB[0])/2.0;\r
750                     WC[1] = (WA[1]+WB[1])/2.0;\r
751                     for (I = IWT1; I <= IWT2; I++) {\r
752                           W[0] = WC[0] + WH[0]*QWORK[I-1];\r
753                           W[1] = WC[1] + WH[1]*QWORK[I-1];\r
754                           prod = WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1);\r
755                           result[0] = result[0] + QWORK[IOFFST + I-1]*prod[0];\r
756                           result[1] = result[1] + QWORK[IOFFST + I-1]*prod[1];\r
757                     } // for (I = IWT1; I <= IWT2; I++)\r
758                     scm.zmlt(result[0], result[1], WH[0], WH[1], cr, ci);\r
759                     result[0] = cr[0];\r
760                     result[1] = ci[0];\r
761                     return result;\r
762             } if (LINEARC != 1)\r
763                 \r
764                 //   INTEGRATE ALONG A CIRCULAR ARC:\r
765             ZI[0] = 0.0;\r
766             ZI[1] = 1.0;\r
767                 PWH = (PHIB-PHIA)/2.0;\r
768                 PWC = (PHIB+PHIA)/2.0;\r
769                 for (I = IWT1; I <= IWT2; I++) {\r
770                         base = RADIUS*Math.exp(ZI[0]* (PWC+PWH*QWORK[I-1]));\r
771                         arg = ZI[1]* (PWC+PWH*QWORK[I-1]);\r
772                         W[0] = base * Math.cos(arg);\r
773                         W[1] = base * Math.sin(arg);\r
774                     prod = WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1);\r
775                     scm.zmlt(W[0], W[1], prod[0], prod[1], cr, ci);\r
776                     result[0] = result[0] + QWORK[IOFFST+I-1] * cr[0];\r
777                     result[1] = result[1] + QWORK[IOFFST+I-1] * ci[0];\r
778                 } // for (I = IWT1; I <= IWT2; I++)\r
779                 scm.zmlt(result[0], result[1], PWH*ZI[0], PWH*ZI[1], cr, ci);\r
780                 result[0] = cr[0];\r
781                 result[1] = ci[0];\r
782                 return result;\r
783         }\r
784 \r
785         \r
786         private double[] WQUAD1(double WA[], double PHIA, int KWA,int IC, double WB[],\r
787                         double PHIB, double RADIUS,int M,int N, double U, double W0[][], double W1[][],\r
788                         double ALFA0[], double ALFA1[], int NPTQ, double QWORK[], int LINEARC) {\r
789                 //   -------------------------------------------------------------\r
790                 //   CALCULATES THE COMPLEX INTEGRAL OF WPROD FROM WA TO WB ALONG\r
791                 //   EITHER A CIRCULAR ARC OR A LINE-EGMENT.  COMPOUND ONE-SIDED\r
792                 //   GAUSS-JACOBI QUDRATURE IS USED.SEE SUBROUTINE WQSUM FOR THE\r
793                 //   CALLING SEQUENCE.\r
794                 \r
795                 //   CHECK FOR ZERO-LENGTH INTEGRAL:\r
796                 //     .. Scalar Arguments ..\r
797                 // DOUBLE COMPLEX WA,WB\r
798                      \r
799                 //     .. Array Arguments ..\r
800                 //      DOUBLE COMPLEX W0(M),W1(N)\r
801                 //      DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (N+M)+3))\r
802                 \r
803                 //     .. Local Scalars ..\r
804                 double WAA[] = new double[2];\r
805                 double WBB[] = new double[2];\r
806                 double ZI[] = new double[2];\r
807                 //    DOUBLE COMPLEX WAA,WBB,ZI\r
808                 double PHAA,PHBB,R;\r
809                 double result[] = new double[2];\r
810                 double result2[] = new double[2];\r
811                 double expb;\r
812                 double arg;\r
813                      \r
814                 //     .. External Functions ..\r
815                 //      DOUBLE COMPLEX WQSUM\r
816                 //      DOUBLE PRECISION DIST\r
817                 //      EXTERNAL WQSUM,DIST\r
818                 \r
819                 //     .. Common blocks ..\r
820                 //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
821         \r
822             if (scm.zabs(WA[0]-WB[0], WA[1]-WB[1]) <= 0.0) {\r
823                 result[0] = 0.0;\r
824                 result[1] = 0.0;\r
825                 return result;\r
826             }\r
827             \r
828             ZI[0] = 0.0;\r
829             ZI[1] = 1.0;\r
830             if (LINEARC != 1) {\r
831                 \r
832                     //   LINE SEGMENT INTEGRATION PATH IS CONCERNED BELOW:\r
833                     // STEP1:ONE-SIDED G-J QUADRATURE FOR LEFT ENDPT WA:\r
834                     R = Math.min(1.0,DIST(M,N,W0,W1,WA,KWA,IC)/scm.zabs(WB[0]-WA[0],WB[1]-WA[1]));\r
835                     WAA[0] = WA[0] + R* (WB[0]-WA[0]);\r
836                     WAA[1] = WA[1] + R* (WB[1]-WA[1]); \r
837                     result = WQSUM(WA,0.0,KWA,IC,WAA,0.0,0.0,M,N,U,W0,W1,ALFA0,\r
838                               ALFA1,NPTQ,QWORK,LINEARC);\r
839                 \r
840                     //   STEP2:ADJOIN INTERVALS OF PURE GAUSS QUADRATURE IF NECESSARY:\r
841                    while (R != 1.0) {\r
842                       R = Math.min(1.0,DIST(M,N,W0,W1,WAA,0,IC)/scm.zabs(WAA[0]-WB[0], WAA[1] - WB[1]));\r
843                       WBB[0] = WAA[0] + R* (WB[0]-WAA[0]);\r
844                       WBB[1] = WAA[1] + R* (WB[1]-WAA[1]);\r
845                       result2 = WQSUM(WAA,0.0,0,2,WBB,0.0,0.0,M,N,U,W0,W1,\r
846                               ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
847                       result[0] = result[0] + result2[0];\r
848                       result[1] = result[1] + result2[1];\r
849                       WAA[0] = WBB[0];\r
850                       WAA[1] = WBB[1];\r
851                    }\r
852                    return result;\r
853             } // if (LINEARC != 1)\r
854                 \r
855                 //   CIRCULAR ARC INTEGRATION PATH IS CONCERNED BELOW:\r
856                 //   STEP1:ONE-SIDED G-J QUADRATURE FOR LEFT ENDPT WA:\r
857                 R = Math.min(1.0,DIST(M,N,W0,W1,WA,KWA,IC)/scm.zabs(WB[0]-WA[0],WB[1]-WA[1]));\r
858                 PHAA = PHIA + R* (PHIB-PHIA);\r
859                 expb = RADIUS * Math.exp(ZI[0]*PHAA);\r
860                 arg = ZI[1] * PHAA;\r
861                 WAA[0] = expb * Math.cos(arg);\r
862                 WAA[1] = expb * Math.sin(arg);\r
863                 result = WQSUM(WA,PHIA,KWA,IC,WAA,PHAA,RADIUS,M,N,U,W0,W1,ALFA0,\r
864                               ALFA1,NPTQ,QWORK,LINEARC);\r
865                 \r
866                 //   STEP2:ADJOIN INTERVALS OF PURE GAUSS QUADRATURE IF NECESSARY:\r
867                 while (R != 1.0) {\r
868                       R = Math.min(1.0,DIST(M,N,W0,W1,WAA,0,IC)/scm.zabs(WAA[0]-WB[0],WAA[1]-WB[1]));\r
869                       PHBB = PHAA + R* (PHIB-PHAA);\r
870                       expb = RADIUS * Math.exp(ZI[0]*PHBB);\r
871                       arg = ZI[1] * PHBB;\r
872                       WBB[0] = expb * Math.cos(arg);\r
873                       WBB[1] = expb * Math.sin(arg);\r
874                       result2 = WQSUM(WAA,PHAA,0,2,WBB,PHBB,RADIUS,M,N,U,W0,W1,\r
875                              ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
876                       result[0] = result[0] + result2[0];\r
877                       result[1] = result[1] + result2[1];\r
878                       PHAA = PHBB;\r
879                       WAA[0] = WBB[0];\r
880                       WAA[1] = WBB[1];\r
881                 } //while (R != 1.0)\r
882                 return result;\r
883         }\r
884 \r
885         \r
886         private double[] WQUAD(double WA[], double PHIA,int KWA,int ICA, double WB[],\r
887                         double PHIB, int KWB,int ICB,double RADIUS,int M,int N, double U, double W0[][],\r
888                     double W1[][], double ALFA0[], double ALFA1[], int NPTQ, double QWORK[],\r
889                     int LINEARC,int IEVL) {\r
890                 //   ---------------------------------------------------------------\r
891                 //   CALCULATES THE COMPLEX INTEGRAL OF WPROD FROM WA TO WB\r
892                 //   ALONG A CIRCULAR ARC OR A LINE-SEGMENT.FUNCTION WQUAD1\r
893                 //   IS CALLED FOUR TIMES,ONE FOR EACH 1/4 OF THE INTERVAL.\r
894                 //   NOTE:  WQUAD1 ALLOWS  ONLY THE LEFT ENDPOINT  TO  BE A\r
895                 //   POSSIBLE SINGULARITY. SEE ROUTINE WQSUM FOR THE CALLING\r
896             //   SEQUENCE.\r
897                 \r
898                 //     .. Scalar Arguments ..\r
899                 //    DOUBLE COMPLEX WA,WB\r
900                      \r
901                 //     .. Array Arguments ..\r
902                 //     DOUBLE COMPLEX W0(M),W1(N)\r
903                 //     DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3))\r
904         \r
905                 //     .. Local Scalars \r
906                 double WMID[] = new double[2];\r
907                 double WMIDA[] = new double[2];\r
908                 double WMIDB[] = new double[2];\r
909                 double WQA[] = new double[2];\r
910                 double WQA1[] = new double[2];\r
911                 double WQA2[] = new double[2];\r
912                 double WQB[] = new double[2];\r
913                 double WQB1[] = new double[2];\r
914                 double WQB2[] = new double[2];\r
915                 double ZI[] = new double[2];\r
916                 // DOUBLE COMPLEX WMID,WMIDA,WMIDB,WQA,WQB,ZI\r
917                 double PHMID,PHMIDA,PHMIDB,PI;\r
918                 double expb;\r
919                 double arg;\r
920                 double result[] = new double[2];\r
921                 //     .. Common blocks ..\r
922                 //     COMMON /PARAM4/UARY,VARY,DLAM,IU\r
923                 \r
924                 PI = Math.PI;\r
925                 ZI[0] = 0.0;\r
926                 ZI[1] = 1.0;\r
927                 \r
928                 //   DETERMINE MIDPTS ON A LINE SEGMENT OR ON A CIRCULAR ARC:\r
929                 if (LINEARC == 0) {\r
930                         WMID[0] = (WA[0] + WB[0])/2.0;\r
931                         WMID[1] = (WA[1] + WB[1])/2.0;\r
932                     WMIDA[0] = (WA[0] + WMID[0])/2.0;\r
933                     WMIDA[1] = (WA[1] + WMID[1])/2.0;\r
934                     WMIDB[0] = (WB[0] + WMID[0])/2.0;\r
935                     WMIDB[1] = (WB[1] + WMID[1])/2.0;\r
936                     PHMID = 0.0;\r
937                     PHMIDA = 0.0;\r
938                     PHMIDB = 0.0;\r
939                 }\r
940                 else {\r
941                 if (IEVL != 1) {\r
942                         if (PHIB < PHIA) {\r
943                             PHIA = PHIA - 2.0*PI;\r
944                         }\r
945                 } // if (IEVL != 1)\r
946                    PHMID = (PHIA+PHIB)/2.0;\r
947                    expb = RADIUS * Math.exp(ZI[0] * PHMID);\r
948                    arg = ZI[1] * PHMID;\r
949                    WMID[0] = expb * Math.cos(arg);\r
950                    WMID[1] = expb * Math.sin(arg);\r
951                    PHMIDA = (PHIA+PHMID)/2.0;\r
952                    expb = RADIUS * Math.exp(ZI[0] * PHMIDA);\r
953                    arg = ZI[1] * PHMIDA;\r
954                    WMIDA[0] = expb * Math.cos(arg);\r
955                    WMIDA[1] = expb * Math.sin(arg);\r
956                    PHMIDB = (PHIB+PHMID)/2.0;\r
957                    expb = RADIUS * Math.exp(ZI[0]*PHMIDB);\r
958                    arg = ZI[1] * PHMIDB;\r
959                    WMIDB[0] = expb * Math.cos(arg);\r
960                    WMIDB[1] = expb * Math.sin(arg);\r
961                 }\r
962                 \r
963                 //   COMPOUND GAUSS-JACOBI PROCESS ACCORDING TO ONE-QUATER RULE:\r
964         WQA1 = WQUAD1(WA,PHIA,KWA,ICA,WMIDA,PHMIDA,RADIUS,M,N,U,W0,W1,\r
965            ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
966         WQA2 = WQUAD1(WMID,PHMID,0,2,WMIDA,PHMIDA,RADIUS,M,N,U,W0,W1,ALFA0,\r
967            ALFA1,NPTQ,QWORK,LINEARC);\r
968         WQA[0] = WQA1[0] - WQA2[0];\r
969         WQA[1] = WQA1[1] - WQA2[1];\r
970         WQB1 = WQUAD1(WB,PHIB,KWB,ICB,WMIDB,PHMIDB,RADIUS,M,N,U,W0,W1,\r
971            ALFA0,ALFA1,NPTQ,QWORK,LINEARC);\r
972         WQB2 = WQUAD1(WMID,PHMID,0,2,WMIDB,PHMIDB,RADIUS,M,N,U,W0,W1,ALFA0,\r
973            ALFA1,NPTQ,QWORK,LINEARC);\r
974         WQB[0] = WQB1[0] - WQB2[0];\r
975         WQB[1] = WQB1[1] - WQB2[1];\r
976         result[0] = WQA[0] - WQB[0];\r
977         result[1] = WQA[1] - WQB[1];\r
978                 return result;\r
979         }\r
980 \r
981 \r
982         \r
983         private void XWTRAN(int M,int N, double X[], double U[],\r
984                         double C[], double W0[][], double W1[][], double PHI0[],\r
985                         double PHI1[]) {\r
986                 //   -----------------------------------------------\r
987                 //  TRANSFORMS X[K-1](UNCONSTRAINED PARAMETERS) TO ACTUAL\r
988                 //  D-SC PARAMETERS:U,C,W0,W1.PHI0 & PHI1 ARE ARGUMENTS\r
989                 //  OF THE PREVERTICES CONTAINED IN W0 & W1.\r
990                 \r
991                 //     .. Scalar Arguments ..\r
992                 //      DOUBLE COMPLEX C\r
993                       \r
994                 //     .. Array Arguments ..\r
995                 //      DOUBLE COMPLEX W0(M),W1(N)\r
996                 //      DOUBLE PRECISION PHI0(M),PHI1(N),X(M+N+2)\r
997         \r
998                 //     .. Local Scalars ..\r
999               double DPH, PH, PHSUM, PI;\r
1000               int  I;\r
1001         \r
1002               PI = Math.PI;\r
1003               if (Math.abs(X[0]) <= 1.0E-14) {\r
1004                   U[0] = 0.5;\r
1005               }\r
1006               else {\r
1007                   U[0] = (X[0]-2.0-Math.sqrt(0.9216*X[0]*X[0]+4.0))/ (2.0*X[0]);\r
1008                   U[0] = (0.0196*X[0]-1.0)/ (U[0]*X[0]);\r
1009               }\r
1010 \r
1011               C[0] = X[1];\r
1012               C[1] = X[2];\r
1013               if (Math.abs(X[N+2]) <= 1.0E-14) {\r
1014                   PHI1[N-1] = 0.0;\r
1015               }\r
1016               else {\r
1017                   PH = (1.0+Math.sqrt(1.0+PI*PI*X[N+2]*X[N+2]))/X[N+2];\r
1018                   PHI1[N-1] = PI*PI/PH;\r
1019               }\r
1020 \r
1021               DPH = 1.0;\r
1022               PHSUM = DPH;\r
1023               \r
1024               for (I = 1; I <= N - 1; I++) {\r
1025                   DPH = DPH/Math.exp(X[2+I]);\r
1026                   PHSUM = PHSUM + DPH;\r
1027               } // for (I = 1; I <= N - 1; I++)\r
1028               DPH = 2.0*PI/PHSUM;\r
1029               PHI1[0] = PHI1[N-1] + DPH;\r
1030               W1[0][0] = U[0] * Math.cos(PHI1[0]);\r
1031               W1[0][1] = U[0] * Math.sin(PHI1[0]);\r
1032               W1[N-1][0] = U[0] * Math.cos(PHI1[N-1]);\r
1033               W1[N-1][1] = U[0] * Math.sin(PHI1[N-1]);\r
1034               PHSUM = PHI1[0];\r
1035               for (I = 1; I <= N - 2; I++) {\r
1036                   DPH = DPH/Math.exp(X[2+I]);\r
1037                   PHSUM = PHSUM + DPH;\r
1038                   PHI1[I] = PHSUM;\r
1039                   W1[I][0] = U[0] * Math.cos(PHSUM);\r
1040                   W1[I][1] = U[0] * Math.sin(PHSUM);\r
1041               } // for (I = 1; I <= N - 2; I++)\r
1042               DPH = 1.0;\r
1043               PHSUM = DPH;\r
1044               for (I = 1; I <= M - 1; I++) {\r
1045                   DPH = DPH/Math.exp(X[N+2+I]);\r
1046                   PHSUM = PHSUM + DPH;\r
1047               } // for (I = 1; I <= M - 1; I++)\r
1048               DPH = 2.0*PI/PHSUM;\r
1049               PHSUM = DPH;\r
1050               PHI0[0] = DPH;\r
1051               W0[0][0] = Math.cos(DPH);\r
1052               W0[0][1] = Math.sin(DPH);\r
1053               for (I = 1; I <= M - 2; I++) {\r
1054                   DPH = DPH/Math.exp(X[N+2+I]);\r
1055                   PHSUM = PHSUM + DPH;\r
1056                   PHI0[I] = PHSUM;\r
1057                   W0[I][0] = Math.cos(PHSUM);\r
1058                   W0[I][1] = Math.sin(PHSUM);\r
1059               } // for (I = 1; I <= M - 2; I++)\r
1060                   return;\r
1061         }\r
1062         \r
1063         private double FMAX(int MN, double FVAL[]) {\r
1064             //   -------------------------\r
1065             //    DETERMINES THE MAXIMUM-NORM OF THE VECTOR FVAL.\r
1066         \r
1067             //     .. Array Arguments ..\r
1068             // DOUBLE PRECISION FVAL(MN)\r
1069                 \r
1070             //      .. Local Scalars ..\r
1071                 double maxVal;\r
1072             double FV;\r
1073             int K;\r
1074         \r
1075             maxVal = 0.0;\r
1076             for (K = 0; K < MN; K++) {\r
1077                 FV = Math.abs(FVAL[K]);\r
1078                 if (FV > maxVal) {\r
1079                         maxVal = FV;\r
1080                 }\r
1081             } // for (K = 0; K < MN; K++)\r
1082             return maxVal;\r
1083         }\r
1084 \r
1085         \r
1086         private void DSCFUN(int NDIM, double X[], double FVAL[],int IFLAG[]) {\r
1087             //    --------------------------------------\r
1088             //  FORMS THE NONLINEAR SYSTEM SATISFIED BY D-SC PARAMETERS.THE\r
1089             //  SUBROUTINE WILL BE CALLED BY NONLINEAR SYSTEM SOLVER HYBRD.\r
1090             //  SEE ROUTINE DSCSOLV FOR THE VARIABLES IN THE COMMON BLOCKS.\r
1091             //  NDIM: THE DIMENSION OF THE SYSTEM.\r
1092             //  X:    UNKNOWNS\r
1093             //  FVAL: THE VECTOR DEFINED BY THE SYSTEM.\r
1094             //  IFLAG:(ON OUTPUT)=1,THE ITERATION WAS SUCCESSFULLY COMPLETED.\r
1095             //         =2,3,OR 4, UNSUCCESSFUL TERMINATION OF THE ITERATION.\r
1096             //         =5 MAY INDICATE THE TOLERANCE GIVEN IN THE CALLING\r
1097             //         SEQUENCE OF HYBRD IS TOO SMALL.\r
1098         \r
1099         \r
1100             //  TRANSFORM UNCONSTRAINED X(K) TO ACTUAL D-SC PARAMETERS\r
1101             //  AND COMPUTE DATA TO BE USED IN WTHETA:\r
1102         \r
1103             //     .. Scalar Arguments ..\r
1104             //  INTEGER IFLAG,NDIM\r
1105             //     ..\r
1106             //     .. Array Arguments ..\r
1107             //  DOUBLE PRECISION FVAL(NDIM),X(NDIM)\r
1108         \r
1109             \r
1110              //.. Local Scalars ..\r
1111                 double cr[] = new double[1];\r
1112                 double ci[] = new double[1];\r
1113             double wout[];\r
1114             double base;\r
1115             double diffR;\r
1116             double diffI;\r
1117                 double WA[] = new double[2];\r
1118                 double WAI[] = new double[2];\r
1119                 double WARC[] = new double[2];\r
1120                 double WB[] = new double[2];\r
1121                 double WBI[] = new double[2];\r
1122                 double WCIRCLE[] = new double[2];\r
1123                 double WIN1[] = new double[2];\r
1124                 double WIN2[] = new double[2];\r
1125                 double WIN3[] = new double[2];\r
1126                 double WIN4[] = new double[2];\r
1127                 double WINT1[] = new double[2];\r
1128                 double WINT2[] = new double[2];\r
1129                 double WINT3[] = new double[2];\r
1130                 double WLINE[] = new double[2];\r
1131                 double WLINE1[] = new double[2];\r
1132                 double WLINE2[] = new double[2];\r
1133                 double WX[] = new double[2];\r
1134                 double ZI[] = new double[2];\r
1135             //  DOUBLE COMPLEX WA,WAI,WARC,WB,WBI,WCIRCLE,WIN1,WIN2,WIN3,WIN4,\r
1136             //               WINT1,WINT2,WINT3,WLINE,WLINE1,WLINE2,WX,ZI\r
1137             double FMAXN,PHIA,PHIB,RADIUS,TEST1;\r
1138             int I,J,K;\r
1139         \r
1140             //     .. External Functions ..\r
1141             // DOUBLE COMPLEX WQUAD\r
1142             // DOUBLE PRECISION FMAX\r
1143             // EXTERNAL WQUAD,FMAX\r
1144         \r
1145             //     .. External Subroutines ..\r
1146             // EXTERNAL THDATA,XWTRAN\r
1147         \r
1148             //     .. Common blocks ..\r
1149             // COMMON /PARAM1/W02,W12,Z02,Z12,C2\r
1150             // COMMON /PARAM2/U2,PHI02,PHI12,ALFA02,ALFA12,QWORK2\r
1151             // COMMON /PARAM3/M2,N2,NPTQ2,ISHAPE2,LINEARC2,NSHAPE,IND\r
1152             // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
1153             // COMMON /PARAM5/ISPRT,ICOUNT\r
1154 \r
1155             XWTRAN(M2,N2,X,U2,C2,W02,W12,PHI02,PHI12);\r
1156               THDATA(U2);\r
1157               ZI[0] = 0.0;\r
1158               ZI[1] = 1.0;\r
1159               ICOUNT = ICOUNT + 1;\r
1160         \r
1161               //  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE INNER POLYGON:\r
1162               wout = WQUAD(W12[N2-1],PHI12[N2-1],N2,1,W12[0],PHI12[0],1,\r
1163                              1,U2[0],M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
1164               scm.zmlt(C2[0], C2[1], wout[0], wout[1], cr, ci);\r
1165               WIN1[0] = Z12[0][0] - Z12[N2-1][0] - cr[0];\r
1166               WIN1[1] = Z12[0][1] - Z12[N2-1][1] - ci[0];\r
1167               scm.zmlt(ZI[0], ZI[1], WIN1[0], WIN1[1], cr, ci);\r
1168               FVAL[0] = ci[0];\r
1169               FVAL[1] = WIN1[1];\r
1170         \r
1171               //  N-1 SIDE LENGTH CONDITIONS FOR THE INNER POLYGON:\r
1172               for (I = 1; I <= N2 - 1; I++) {\r
1173                   WINT1 = WQUAD(W12[I-1],PHI12[I-1],I,1,W12[I],PHI12[I],I+1,1,U2[0],M2,N2,\r
1174                          U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
1175                   scm.zmlt(C2[0], C2[1], WINT1[0], WINT1[1], cr, ci);\r
1176                   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
1177               } // for (I = 1; I <= N2 - 1; I++)\r
1178         \r
1179               //  TWO EQUATIONS TO FIX THE RELATIVE POSITION OF THE INNER POLYGON:\r
1180               TEST1 = Math.cos(PHI12[N2-1]);\r
1181               if (TEST1 < U2[0]) {\r
1182         \r
1183                   // IF THE LINE PATH FROM W02[M2-1] TO W12[N2-1] IS OUT OF DOMAIN,THE\r
1184               //  COMBINATION OF TWO DIFFERENT PATHS WILL BE USED INSTEAD:\r
1185                   WX[0] = U2[0];\r
1186                   WX[1] = 0.0;\r
1187                   WLINE = WQUAD(W02[M2-1],0.0,M2,0,WX,0.0,0,2,0.0,M2,N2,U2[0],W02,W12,ALFA02,\r
1188                      ALFA12,NPTQ2,QWORK2,0,2);\r
1189                   if (PHI12[N2-1] <= 0.0) {\r
1190                       WARC = WQUAD(W12[N2-1],PHI12[N2-1],N2,1,WX,0.0,0,2,U2[0],M2,N2,U2[0],W02,W12,\r
1191                         ALFA02,ALFA12,NPTQ2,QWORK2,1,2);\r
1192                       WIN2[0] = WLINE[0] - WARC[0];\r
1193                       WIN2[1] = WLINE[1] - WARC[1];\r
1194                   }\r
1195                   else {\r
1196                       WARC = WQUAD(WX,0.0,0,2,W12[N2-1],PHI12[N2-1],N2,1,U2[0],M2,N2,U2[0],W02,W12,\r
1197                         ALFA02,ALFA12,NPTQ2,QWORK2,1,2);\r
1198                       WIN2[0] = WLINE[0] + WARC[0];\r
1199                       WIN2[1] = WLINE[1] + WARC[1];\r
1200                   }\r
1201               } // if (TEST1 < U2[0])\r
1202               else {\r
1203 \r
1204                   WIN2 = WQUAD(W02[M2-1],0.0,M2,0,W12[N2-1],0.0,N2,1,0.0,M2,N2,U2[0],W02,W12,ALFA02,\r
1205                     ALFA12,NPTQ2,QWORK2,0,2);\r
1206               }\r
1207               scm.zmlt(C2[0], C2[1], WIN2[0], WIN2[1], cr, ci);\r
1208               diffR = Z12[N2-1][0] - Z02[M2-1][0] - cr[0];\r
1209               diffI = Z12[N2-1][1] - Z02[M2-1][1] - ci[0];\r
1210               scm.zmlt(ZI[0], ZI[1], diffR, diffI, cr, ci);\r
1211               FVAL[N2+1] = ci[0];\r
1212               FVAL[N2+2] = diffI;\r
1213         \r
1214               //  TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE OUTER POLYGON:\r
1215               wout = WQUAD(W02[M2-1],PHI02[M2-1],M2,0,W02[0],PHI02[0],1,\r
1216                                     0,1.0,M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
1217               scm.zmlt(C2[0], C2[1], wout[0], wout[1], cr, ci);\r
1218               WIN3[0] = Z02[0][0] - Z02[M2-1][0] - cr[0];\r
1219               WIN3[1] = Z02[0][1] - Z02[M2-1][1] - ci[0];\r
1220               scm.zmlt(ZI[0],ZI[1], WIN3[0], WIN3[1], cr, ci);\r
1221               FVAL[N2+3] = ci[0];\r
1222               FVAL[N2+4] = WIN3[1];\r
1223         \r
1224               if (M2 == 3) {\r
1225         \r
1226                   // CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
1227                   if (ISPRT == 1) {\r
1228                       FMAXN = FMAX(NDIM,FVAL);\r
1229                       System.out.println("# of iterations = " + ICOUNT);\r
1230                       System.out.println("Maximum norm of the residuals = " + FMAXN);\r
1231                   } // if (ISPRT == 1)\r
1232                   return;\r
1233 \r
1234               } // if (M2 == 3)\r
1235 \r
1236               if (ISHAPE2 != 1) {\r
1237         \r
1238                   //  M2-3 SIDE LENGTH CONDITIONS OF THE OUTER POLYGON:\r
1239                   for (J = 1; J <= M2 - 3; J++) {\r
1240                       WINT2 = WQUAD(W02[J-1],PHI02[J-1],J,0,W02[J],PHI02[J],J+1,0,1.0,\r
1241                           M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
1242                       scm.zmlt(C2[0], C2[1], WINT2[0], WINT2[1], cr, ci);\r
1243                       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
1244                   } // for (J = 1; J <= M2 - 3; J++)\r
1245         \r
1246                   //  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
1247                   if (ISPRT == 1) {\r
1248                       FMAXN = FMAX(NDIM,FVAL);\r
1249                       System.out.println("# of iterations = " + ICOUNT);\r
1250                       System.out.println("Maximum norm of the residuals = " + FMAXN);\r
1251                   } // if (ISPRT == 1)\r
1252                   return;\r
1253               } // if (ISHAPE2 != 1)\r
1254         \r
1255               //  OUTER POLYGON CONTAINS SOME INFINITE VERTICES & FOR EACH OF THEM\r
1256               //  TWO LENGTH CONDITIONS WILL BE REPLACED BY A COMPLEX INTEGRAL:\r
1257               for (K = 1; K <= NSHAPE - 1; K++) {\r
1258                   if (IND[K] != 2 && IND[K-1] < IND[K]-2) {\r
1259                       for (J = IND[K-1] + 1; J <= IND[K] - 2; J++) {\r
1260                           WINT3 = WQUAD(W02[J-1],PHI02[J-1],J,0,W02[J],PHI02[J],J+1,0,\r
1261                              1.0,M2,N2,U2[0],W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,LINEARC2,2);\r
1262                           scm.zmlt(C2[0], C2[1], WINT3[0], WINT3[1], cr, ci);\r
1263                           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
1264                       } // for (J = IND[K-1] + 2; J <= IND[K] - 1; J++) \r
1265                   } // if (IND[K] != 1 && IND[K-1] < IND[K]-2)\r
1266                   if (K == NSHAPE-1 || IND[K] == M2-1) {\r
1267                           continue;\r
1268                   }\r
1269         \r
1270                   //  THE COMBINATION  OF THREE DIFFERENT PATHS  IS USED TO INTEGRATE\r
1271                   //  FROM WA TO WB TO AVOID DOMAIN PROBLEM.THE BY-PRODUCT OF THIS IS\r
1272                   //  THAT IT IS NUMERICALLY  MORE EFFICIENT THAN A SINGLE LINE PATH:\r
1273                   WA[0] = W02[IND[K]-2][0];\r
1274                   WA[1] = W02[IND[K]-2][1];\r
1275                   WB[0] = W02[IND[K]][0];\r
1276                   WB[1] = W02[IND[K]][1];\r
1277                   PHIA = PHI02[IND[K]-2];\r
1278                   PHIB = PHI02[IND[K]];\r
1279                   RADIUS = (1.0+U2[0])/2.0;\r
1280                   scm.zmlt(ZI[0], ZI[1], PHIA, 0.0, cr, ci);\r
1281                   base = RADIUS * Math.exp(cr[0]);\r
1282                   WAI[0] = base * Math.cos(ci[0]);\r
1283                   WAI[1] = base * Math.sin(ci[0]);\r
1284                   scm.zmlt(ZI[0], ZI[1], PHIB, 0.0, cr, ci);\r
1285                   base = RADIUS * Math.exp(cr[0]);\r
1286                   WBI[0] = base * Math.cos(ci[0]);\r
1287                   WBI[1] = base * Math.sin(ci[0]);\r
1288                   WLINE1 = WQUAD(WA,0.0,IND[K]-1,0,WAI,0.0,0,2,0.0,M2,N2,U2[0],\r
1289                              W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2);\r
1290                   WLINE2 = WQUAD(WB,0.0,IND[K]+1,0,WBI,0.0,0,2,0.0,M2,N2,U2[0],\r
1291                              W02,W12,ALFA02,ALFA12,NPTQ2,QWORK2,0,2);\r
1292                   WCIRCLE = WQUAD(WAI,PHIA,0,2,WBI,PHIB,0,2,RADIUS,M2,N2,U2[0],W02,W12,\r
1293                            ALFA02,ALFA12,NPTQ2,QWORK2,1,2);\r
1294                   scm.zmlt(C2[0], C2[1], (WLINE1[0]+WCIRCLE[0]-WLINE2[0]),\r
1295                                   (WLINE1[1]+WCIRCLE[1]-WLINE2[1]), cr, ci);\r
1296                   WIN4[0] = cr[0];\r
1297                   WIN4[1] = ci[0];\r
1298                   diffR = Z02[IND[K]][0]-Z02[IND[K]-2][0]-WIN4[0];\r
1299                   diffI = Z02[IND[K]][1]-Z02[IND[K]-2][1]-WIN4[1];\r
1300                   scm.zmlt(ZI[0], ZI[1], diffR, diffI, cr, ci);\r
1301                   FVAL[N2+5+IND[K]-2] = ci[0];\r
1302                   FVAL[N2+5+IND[K]-1] = diffI;\r
1303               } // for (K = 1; K <= NSHAPE - 1; K++)\r
1304         \r
1305               //  CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL:\r
1306               if (ISPRT == 1) {\r
1307               FMAXN = FMAX(NDIM,FVAL);\r
1308               System.out.println("# of iterations = " + ICOUNT);\r
1309               System.out.println("Maximum norm of the residuals = " + FMAXN);\r
1310           } // if (ISPRT == 1)\r
1311           return;\r
1312         }\r
1313 \r
1314 \r
1315         \r
1316         private void DSCSOLV(double TOL,int IGUESS,int M,int N, double U[],\r
1317                         double C[], double W0[][], double W1[][], double PHI0[],\r
1318                         double PHI1[], double Z0[][],  double Z1[][], double ALFA0[],\r
1319                    double ALFA1[], int NPTQ, double QWORK[], int ISHAPE, int LINEARC) {\r
1320         //  -----------------------------------------------------------------\r
1321         //  SOLVES THE NONLINEAR SYSTEM FOR D-SC PARAMETERS.\r
1322         //  CALLING SEQUENCE:\r
1323         //  TOL       A TOLERANCE TO CONTROL THE CONVERGENCE IN HYBRD\r
1324         //  IGUESS    (=0 )A NON-EQUALLY SPACED INITIAL GUESS OR(=1)THE\r
1325         //            OTHER  EQUALLY-SPACED  INITIAL GUESS PROVIDED, OR\r
1326         //            (=2)USER-SUPPLIED INITIAL GUESS WHICH IS REQUIRED\r
1327         //            TO BE THE ARGUMENTS OF THE INITIAL PREVERTICES.\r
1328         //            ROUTINE ARGUM MAY BE USED FOR COMPUTING ARGUMENTS\r
1329         //            IF NEEDED. NOTE: C WILL BE COMPUTED IN THE ROUTINE\r
1330         //            (NOT SUPPLIED!)\r
1331         //  M,N,U,C,W0,W1,PHI0,PHI1,ALFA0,ALFA1,Z0,Z1\r
1332         //            CONSTITUTE THE GEOMETRY OF THE POLYGONAL REGION\r
1333         //            AND THE MAPPING FUNCTION. ON RETURN U,C,W0,& W1\r
1334         //            WILL  CONTAIN COMPUTED PARAMETERS. (PHI0 & PHI1\r
1335         //            WILL CONTAIN THE ARGUMENTS OF THE PREVERTICES.)\r
1336         //  QWORK     SEE CALLING SEQUENCE DOCUMENTATION IN QINIT.THE ARRAY\r
1337         //            MUST HAVE BEEN FILLED BY QINIT BEFORE CALLING DSCSOLV.\r
1338         //  NPTQ      THE NUMBER OF GAUSS-JACOBI POINTS USED\r
1339         //  ISHAPE    INDICATES  THAT  OUTER POLYGON  CONTAINS NO INFINITE\r
1340         //            VERTICES (ISHAPE=0) OR IT HAS SOME INFINITE VERTICES\r
1341         //            (ISHAPE=1).\r
1342         //  LINEARC   INTEGER VARIABLE TO CONTROL INTEGRATION PATH.IN PATICULAR:\r
1343         //            LINEARC=0:INTEGRATING ALONG LINE SEGMENT WHENEVER POSSIBLE\r
1344         //            LINEARC=1:INTEGRATING ALONG CIRCULAR ARC WHENEVER POSSIBLE\r
1345         //  THE DISCRIPTION OF ARRAY IND & VARIABLE NSHAPE NEEDED IN DSCFUN:\r
1346         //  IND       CONTAINS INDICES  CORRESPONDING TO  THOSE INFINITE\r
1347         //            VERTICES,BUT THE FIRST & THE LAST ENTRIS MUST BE -1\r
1348         //           & M-2(M IS THE # OF VERTICES ON THE OUTER POLYGON).\r
1349         //  NSHAPE    THE DIMENSION OF THE INTERGE ARRAY IND(NSHAPE)\r
1350         \r
1351         \r
1352         //     .. Scalar Arguments ..\r
1353         //      DOUBLE COMPLEX C\r
1354               \r
1355         //     .. Array Arguments ..\r
1356         //      DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
1357         //      DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N),\r
1358         //                      QWORK(NPTQ* (2* (N+M)+3))\r
1359                 \r
1360                 // Local scalars\r
1361                 double cr[] = new double[1];\r
1362                 double ci[] = new double[1];\r
1363                 double C1[] = new double[2];\r
1364                 double WINT[];\r
1365                 double ZI[] = new double[2];\r
1366                 double AVE, BOTM, DSTEP, FACTOR, PI, TOP;\r
1367                 int INFO[] = new int[1];\r
1368                 int NFEV[] = new int[1];\r
1369                 int I, K, KM, KN, MAXFUN, NM, NWDIM;\r
1370                 \r
1371                 // Local arrays\r
1372                 double DIAG[] = new double[42];\r
1373                 double FJAC[][] = new double[42][42];\r
1374                 double FVAL[] = new double[42];\r
1375                 double QW[] = new double[1114];\r
1376                 double X[] = new double[42];\r
1377                 \r
1378             // Common blocks ..\r
1379             // COMMON /PARAM1/W02,W12,Z02,Z12,C2\r
1380             // COMMON /PARAM2/U2,PHI02,PHI12,ALFA02,ALFA12,QWORK2\r
1381             // COMMON /PARAM3/M2,N2,NPTQ2,ISHAPE2,LINEARC2,NSHAPE,IND\r
1382             // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
1383             // COMMON /PARAM5/ISPRT,ICOUNT\r
1384 \r
1385                 ZI[0] = 0.0;\r
1386                 ZI[1] = 1.0;\r
1387                 PI = Math.PI;\r
1388                 ICOUNT = 0;\r
1389                 if (ISHAPE != 0) {\r
1390                         NSHAPE = 1;\r
1391                         for (I = 2; I <= M-1; I++) {\r
1392                             if (ALFA0[I-1] >= 0.0) {\r
1393                                 continue;\r
1394                             }\r
1395                             NSHAPE = NSHAPE + 1;\r
1396                             IND[NSHAPE-1] = I;\r
1397                         } // for (I = 2; I <= M-1; I++)\r
1398                         IND[0] = 0;\r
1399                         NSHAPE = NSHAPE + 1;\r
1400                         IND[NSHAPE-1] = M-1;\r
1401                 } // if (ISHAPE != 0)\r
1402                 \r
1403                 // Fix one prevertex\r
1404                 W0[M-1][0] = 1.0;\r
1405                 W0[M-1][1] = 0.0;\r
1406                 PHI0[M-1] = 0.0;\r
1407                 \r
1408                 // Following two value assignments are to satisfy the compiler WATFOR77:\r
1409                 X[1] = 0.0;\r
1410                 X[2] = 0.0;\r
1411                 // Initial guess (IGUESS = 0):\r
1412                 if (IGUESS == 0) {\r
1413                     X[0] = 1.0/0.5 - 1.0/0.46;\r
1414                     AVE = 2.0 * PI/(double)N;\r
1415                     for (I = 1; I <= N-2; I++) {\r
1416                         X[2 + I] = Math.log((AVE + 0.0001*I)/(AVE + 0.0001*(I+1)));\r
1417                     } // for (I = 1; I <= N-2; I++)\r
1418                     X[N+1] = Math.log((AVE+0.0001*(N-1))/(2.0*PI-(N-1)* (AVE+N*0.00005)));\r
1419                     X[N+2] = 1.0/ (4.0-0.1) - 1.0/ (4.0+0.1);\r
1420                     AVE = 2.0*PI/M;\r
1421                     for (I = 1; I <= M - 2; I++) {\r
1422                         X[N+2+I] = Math.log((AVE+0.0001*(I-1))/(AVE+0.0001*I));\r
1423                     } // for (I = 1; I <= M - 2; I++)\r
1424                     X[M+N+1] = Math.log((AVE+0.0001*(M-2))/(2.0*PI-(M-1)* (AVE+(M-2)*0.00005)));\r
1425                 } // if (IGUESS == 0)\r
1426                 else if (IGUESS == 1) {\r
1427                         //  INITIAL GUESS (IGUESS=1):\r
1428                     X[0] = 1.0/0.53 - 1.0/0.43;\r
1429                     for (I = 1; I <= N - 1; I++) {\r
1430                         X[2+I] = 0.0;\r
1431                     } // for (I = 1; I <= N - 1; I++)\r
1432                     X[N+2] = 1.0/ (4.0-0.1) - 1.0/ (4.0+0.1);\r
1433                     for (I = 1; I <= M - 1; I++) {\r
1434                         X[N+2+I] = 0.0;\r
1435                     } // for (I = 1; I <= M - 1; I++)\r
1436                 } // else if (IGUESS == 1)\r
1437                 else {\r
1438                         X[0] = 1.0/ (0.98-U[0]) - 1.0/ (U[0]-0.02);\r
1439                         for (K = 1; K <= N - 1; K++) {\r
1440                             KN = K - 1;\r
1441                             if (KN == 0) {\r
1442                                 KN = N;\r
1443                             }\r
1444                             TOP = PHI1[K-1] - PHI1[KN-1];\r
1445                             if (TOP < 0.0) {\r
1446                                 TOP = TOP + 2.0*PI;\r
1447                             }\r
1448                             BOTM = PHI1[K] - PHI1[K-1];\r
1449                             if (BOTM < 0.0) {\r
1450                                 BOTM = BOTM + 2.0*PI;\r
1451                             }\r
1452                             X[2+K] = Math.log(TOP) - Math.log(BOTM);\r
1453                         } // for (K = 1; K <= N - 1; K++)\r
1454                         X[N+2] = 1.0/ (PI-PHI1[N-1]) - 1.0/ (PI+PHI1[N-1]);\r
1455                         for (K = 1; K <= M - 1; K++) {\r
1456                             KM = K - 1;\r
1457                             if (KM == 0) {\r
1458                                 KM = M;\r
1459                             }\r
1460                             TOP = PHI0[K-1] - PHI0[KM-1];\r
1461                             if (TOP < 0.0) {\r
1462                                 TOP = TOP + 2.0*PI;\r
1463                             }\r
1464                             BOTM = PHI0[K] - PHI0[K-1];\r
1465                             if (BOTM < 0.0) {\r
1466                                 BOTM = BOTM + 2.0*PI;\r
1467                             }\r
1468                             X[N+2+K] = Math.log(TOP) - Math.log(BOTM);\r
1469                         } // for (K = 1; K <= M - 1; K++)\r
1470                 } // else\r
1471                 \r
1472                 // CALCULATE THE INITIAL GUESS X[1] & X[2] TO MATCH\r
1473                 // THE CHOICE FOR X[0],X[3],...,X[M+N+1]:\r
1474                 XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1);\r
1475         THDATA(U);\r
1476         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
1477                             ALFA1,NPTQ,QWORK,0,2);\r
1478         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
1479         C1[0] = cr[0];\r
1480         C1[1] = ci[0];\r
1481         scm.zmlt(ZI[0], ZI[1], C1[0], C1[1], cr, ci);\r
1482         X[1] = ci[0];\r
1483         X[2] = C1[1];\r
1484         \r
1485         //  HYBRD CONTROL PARAMETERS:\r
1486         DSTEP = 1.0E-7;\r
1487         MAXFUN = 200* (M+N);\r
1488         FACTOR = 2.0;\r
1489         NM = M + N + 2;\r
1490       \r
1491         //  COPY RELEVANT DATA TO THE COMMON BLOCK IN DSCFUN:\r
1492         M2 = M;\r
1493         N2 = N;\r
1494         ISHAPE2 = ISHAPE;\r
1495         LINEARC2 = LINEARC;\r
1496         NPTQ2 = NPTQ;\r
1497         for (K = 0; K < M; K++) {\r
1498             W02[K][0] = W0[K][0];\r
1499             W02[K][1] = W0[K][1];\r
1500             PHI02[K] = PHI0[K];\r
1501             Z02[K][0] = Z0[K][0];\r
1502             Z02[K][1] = Z0[K][1];\r
1503             ALFA02[K] = ALFA0[K];\r
1504         } // for (K = 0; K < M; K++) \r
1505         for (K = 0; K < N; K++) {\r
1506             W12[K][0] = W1[K][0];\r
1507             W12[K][1] = W1[K][1];\r
1508             PHI12[K] = PHI1[K];\r
1509             Z12[K][0] = Z1[K][0];\r
1510             Z12[K][1] = Z1[K][1];\r
1511             ALFA12[K] = ALFA1[K];\r
1512         } // for (K = 0; K < N; K++)\r
1513         NWDIM = NPTQ* (2* (M+N)+3);\r
1514         for (I = 0; I < NWDIM; I++) {\r
1515             QWORK2[I] = QWORK[I];\r
1516         } // for (I = 0; I < NWDIM; I++)\r
1517         \r
1518         if (ISPRT != 1) {\r
1519             System.out.println("Time for obtaining a converged result may range");\r
1520             System.out.println("from several seconds to several minutes or longer");\r
1521             System.out.println("depending on the complexity of the geometry of the");\r
1522             System.out.println("region and the local computing system");\r
1523             System.out.println("so be patient");\r
1524         }\r
1525         \r
1526         // Solve nonlinear system with hybrid\r
1527         double QTF[] = new double[NM];\r
1528         double WA1[] = new double[NM];\r
1529         double WA2[] = new double[NM];\r
1530         double WA3[] = new double[NM];\r
1531         double WA4[] = new double[NM];\r
1532         HYBRD(dscfun, NM, X, FVAL, TOL, MAXFUN, NM, NM, DSTEP, DIAG, 1, FACTOR,\r
1533                   -1, INFO, NFEV, FJAC, 42, QW, 903, QTF, WA1, WA2, WA3, WA4);\r
1534         for (I = 0; I < NM; I++) {\r
1535                 QW[I+903] = QTF[I];\r
1536         }\r
1537         \r
1538         //  CHECK ERROR INFORMATION. THE DESCRIPTION OF INFO CAN BE\r
1539         //  FOUND IN THE DOCUMENTATION OF HYBRD.(INFO=1: SUCCESSFUL EXIT)\r
1540         if (INFO[0] == 0) {\r
1541                 System.out.println("INFO[0] == 0, IMPROPER INPUT PARAMETERS.");\r
1542         }\r
1543         else if (INFO[0] == 1) {\r
1544                 System.out.println("INFO[0] == 1, SUCESSFUL EXIT");\r
1545                 System.out.println("RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES IS AT MOST XTOL.");\r
1546         }\r
1547         else if (INFO[0] == 2) {\r
1548                 System.out.println("INFO[0] == 2, UNSUCESSFUL EXIT");\r
1549                 System.out.println("NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED MAXFEV");\r
1550         }\r
1551         else if (INFO[0] == 3) {\r
1552                 System.out.println("INFO[0] == 3, UNSUCESSFUL EXIT");\r
1553                 System.out.println("XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN THE APPROXIMATE SOLUTION X IS POSSIBLE.");\r
1554         }\r
1555         else if (INFO[0] == 4) {\r
1556                 System.out.println("INFO[0] == 4, UNSUCESSFUL EXIT");\r
1557                 System.out.println("ITERATION IS NOT MAKING GOOD PROGRESS, ");\r
1558                 System.out.println("AS MEASURED BY THE IMPROVEMENT FROM THE LAST FIVE JACOBIAN EVALUATIONS.");\r
1559         }\r
1560         else if (INFO[0] == 5) {\r
1561                 System.out.println("INFO[0] == 5, UNSUCESSFUL EXIT");\r
1562                 System.out.println("ITERATION IS NOT MAKING GOOD PROGRESS, ");\r
1563             System.out.println("AS MEASURED BY THE IMPROVEMENT FROM THE LAST TEN ITERATIONS.");\r
1564         }\r
1565                 \r
1566         // COPY OUTPUT DATA FROM COMMON BLOCK AND PRINT THE RESULTS:\r
1567         XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1);\r
1568         DSCPRINT(M,N,C,U[0],W0,W1,PHI0,PHI1,TOL,NPTQ);\r
1569         return;\r
1570 \r
1571         }\r
1572         \r
1573         private double[] ZDSC(double WW[],int KWW,int IC,int M,int N, double U,\r
1574                         double C[], double W0[][], double W1[][], double Z0[][], double Z1[][],\r
1575                         double ALFA0[], double ALFA1[], double PHI0[], double PHI1[], int NPTQ,\r
1576                         double QWORK[], int IOPT) {\r
1577         //   ----------------------------------------------------------\r
1578                 //   COMPUTES THE FORWORD MAP Z(WW) BY INTEGRATING FROM WA TO\r
1579                 //   WW WHERE WA IS THE NEAREST PREVERTEX TO WW.KWW=K,IC=0 IF\r
1580                 //   WW=W0(K), OR KWW=K, IC=1 IF WW=W1(K),OR KWW=0 AND IC=2\r
1581                 //   OTHERWISE.\r
1582                 //   IOPT: =1 IS NORMALLY ASSUMED FOR ANY GEOMETRY.\r
1583                 //         =# OTHER THAN 1 ASSUMES THAT ONLY LINE SEGMENT PATH\r
1584                 //          IS USED.\r
1585                 //   SEE ROUTINE DSCSOLV FOR THE REST OF THE CALLING SEQUENCE.\r
1586                 \r
1587                 \r
1588                 //     .. Scalar Arguments ..\r
1589                 //    DOUBLE COMPLEX C,WW\r
1590                 //    DOUBLE PRECISION U\r
1591                 //    INTEGER IC,IOPT,KWW,M,N,NPTQ\r
1592                 \r
1593                 //     .. Array Arguments ..\r
1594                 //    DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
1595                 //    DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N),\r
1596                 //                     QWORK(NPTQ* (2* (M+N)+3))\r
1597         \r
1598                 //     .. Scalars in Common ..\r
1599                 //    DOUBLE PRECISION DLAM\r
1600                 //    INTEGER IU\r
1601                 \r
1602                 //     .. Arrays in Common ..\r
1603                 //    DOUBLE PRECISION UARY(8),VARY(3)\r
1604                 \r
1605                 //     .. Local Scalars ..\r
1606                 double WA[] = new double[2];\r
1607                 double WB[] = new double[2];\r
1608                 double WINT1[] = new double[2];\r
1609                 double WINT2[] = new double[2];\r
1610                 double WW0[] = new double[2];\r
1611                 double ZA[] = new double[2];\r
1612                 double ZI[] = new double[2];\r
1613                 double wout[];\r
1614                 double cr[] = new double[1];\r
1615                 double ci[] = new double[1];\r
1616                 double result[] = new double[2];\r
1617                 double base;\r
1618                 // DOUBLE COMPLEX WA,WB,WINT1,WINT2,WW0,ZA,ZI\r
1619                 double DWW0,PHIWB,PHIWW0,PI;\r
1620                 double factor;\r
1621                 int INEAR[] = new int[1];\r
1622                 int KNEAR[] = new int[1];\r
1623                 int IBD;\r
1624 \r
1625                 //     .. External Functions ..\r
1626                 //    DOUBLE COMPLEX WQUAD\r
1627                 //    DOUBLE PRECISION ARGUM\r
1628                 //    EXTERNAL WQUAD,ARGUM\r
1629         \r
1630                 //     .. External Subroutines ..\r
1631                 //    EXTERNAL NEARW\r
1632 \r
1633                 //     .. Common blocks ..\r
1634                 //    COMMON /PARAM4/UARY,VARY,DLAM,IU\r
1635                 \r
1636                 ZI[0] = 0.0;\r
1637                 ZI[1] = 1.0;\r
1638                 PI = Math.PI;\r
1639                 IBD = 1;\r
1640                 WW0[0] = WW[0];\r
1641                 WW0[1] = WW[1];\r
1642                 \r
1643                 if (IOPT == 1) {\r
1644                     if (Math.abs(scm.zabs(WW0[0],WW0[1])-1.0) <= 1.0E-11) {\r
1645                         factor = (1.0 + U)/2.0;\r
1646                         WW0[0] = factor * WW0[0];\r
1647                         WW0[1] = factor * WW0[1];\r
1648                         IBD = 2;\r
1649                     } // if (Math.abs(scm.zabs(WW0[0],WW0[1])-1.0) <= 1.0E-11)\r
1650                 } // if (IOPT == 1)\r
1651                 \r
1652                 // FIND THE NEAREST PREVERTEX TO THE PT WW0:\r
1653                 NEARW(M,N,W0,W1,ALFA0,WW0,KNEAR,INEAR);\r
1654                 if (INEAR[0] == 0) {\r
1655                     ZA[0] = Z0[KNEAR[0]-1][0];\r
1656                     ZA[1] = Z0[KNEAR[0]-1][1];\r
1657                     WA[0] = W0[KNEAR[0]-1][0];\r
1658                     WA[1] = W0[KNEAR[0]-1][1];\r
1659                 }\r
1660                 else {\r
1661                     ZA[0] = Z1[KNEAR[0]-1][0];\r
1662                     ZA[1] = Z1[KNEAR[0]-1][1];\r
1663                     WA[0] = W1[KNEAR[0]-1][0];\r
1664                     WA[1] = W1[KNEAR[0]-1][1];\r
1665                 }\r
1666 \r
1667                 if (IOPT != 1) {\r
1668                     wout = WQUAD(WA,0.0,KNEAR[0],INEAR[0],WW,0.0,KWW,IC,0.0,M,N,U,\r
1669                                            W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
1670                     scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
1671                     result[0] = ZA[0] + cr[0];\r
1672                     result[1] = ZA[1] + ci[0];\r
1673                     return result;\r
1674                 } // if (IOPT != 1)\r
1675                 \r
1676                 // DETERMINE THE POINT CLOSEST TO WA ON THE\r
1677                 // SAME CONCENTRIC CIRCLE AS WW0 IS ON:\r
1678                 PHIWW0 = ARGUM(scm.zabs(WW0[0], WW0[1]),WW0);\r
1679                 DWW0 = scm.zabs(WW0[0],WW0[1]);\r
1680                 if (INEAR[0] == 0) {\r
1681                     PHIWB = PHI0[KNEAR[0]-1];\r
1682                 }\r
1683                 else {\r
1684                     PHIWB = PHI1[KNEAR[0]-1];\r
1685                 }\r
1686 \r
1687                 scm.zmlt(ZI[0], ZI[1], PHIWB, 0.0, cr, ci);\r
1688                 base = DWW0 * Math.exp(cr[0]);\r
1689                 WB[0] = base * Math.cos(ci[0]);\r
1690                 WB[1] = base * Math.sin(ci[0]);\r
1691                 //\r
1692             //   INTEGRATION FROM WA TO WB ON A LINE SEGMENT:\r
1693                 if (scm.zabs(WB[0]-WA[0],WB[1]-WA[1]) <= 1.0E-11) {\r
1694                         WINT1[0] = 0.0;\r
1695                         WINT1[1] = 0.0;\r
1696                 }\r
1697                 else {\r
1698                       WINT1 = WQUAD(WA,0.0,KNEAR[0],INEAR[0],WB,0.0,0,2,0.0,M,N,U,W0,W1,\r
1699                             ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
1700                 }\r
1701                 \r
1702                 //   INTEGRATION FROM WB TO WW ON A CIRCULAR ARC:\r
1703                 if (scm.zabs(WB[0]-WW0[0],WB[1]-WW0[1]) <= 1.0E-11) {\r
1704                         WINT2[0] = 0.0;\r
1705                         WINT2[1] = 0.0;\r
1706             }\r
1707                 else {\r
1708                     if (Math.abs(PHIWB-2.0*PI-PHIWW0) < \r
1709                          Math.abs(PHIWB-PHIWW0)) PHIWB = PHIWB - 2.0*PI;\r
1710                     if (Math.abs(PHIWW0-2.0*PI-PHIWB) <\r
1711                          Math.abs(PHIWB-PHIWW0)) PHIWW0 = PHIWW0 - 2.0*PI;\r
1712                     if (scm.zabs(WB[0]-WA[0],WB[1]-WA[1]) <= 1.0E-11) {\r
1713                           WINT2 = WQUAD(WB,PHIWB,KNEAR[0],INEAR[0],WW0,PHIWW0,KWW,IC,DWW0,M,N,\r
1714                                  U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,1,1);\r
1715                     }\r
1716                     else {\r
1717                       WINT2 = WQUAD(WB,PHIWB,0,2,WW0,PHIWW0,0,2,DWW0,M,N,U,W0,W1,ALFA0,\r
1718                                     ALFA1,NPTQ,QWORK,1,1);\r
1719                     }\r
1720                 } // else\r
1721                 \r
1722                 // EVALUATE THE MAPPING FUNCTION. THE INTEGRATION PATH IS\r
1723                 // A COMBINATION OF A CIRCULAR ARC AND LINE SEGMENT(S):\r
1724                 scm.zmlt(C[0], C[1], WINT1[0] + WINT2[0], WINT1[1] + WINT2[1], cr, ci);\r
1725                 result[0] = ZA[0] + cr[0];\r
1726                 result[1] = ZA[1] + ci[0];\r
1727                 if (IBD == 2) {\r
1728                         wout = WQUAD(WW0,0.0,0,2,WW,0.0,KWW,IC,0.0,M,N,U,\r
1729                                                 W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
1730                         scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
1731                         result[0] = result[0] + cr[0];\r
1732                         result[1] = result[1] + ci[0];\r
1733                 } // if (IBD == 2)\r
1734 \r
1735                 return result;    \r
1736     }\r
1737         \r
1738         private double[] WDSC(double ZZ[], int M,int N, double U, double C[], double W0[][], double W1[][],\r
1739                         double Z0[][], double Z1[][], double ALFA0[], double ALFA1[], double PHI0[], double PHI1[],\r
1740                         int NPTQ, double QWORK[], double EPS,int IOPT) {\r
1741         //    ----------------------------------------------------\r
1742                 // COMPUTES THE INVERSE MAP AFTER  ALL ACCESSARY PARAMETERS HAVE\r
1743                 // BEEN DETERMINED.EULER'S SCHEME FOR SOLVING ODE IS USED TO GIVE\r
1744             // THE INITIAL GUESS WHICH IS THEN REFINED BY NEWTON'S ITERATION.\r
1745                 // ZZ IS THE POINT AT WHICH INVERSE MAP IS EVALUATED. EPS IS THE\r
1746                 // REQUIRED ACCURACY SUPPLIED BY THE USER.\r
1747                 // NOTE: ZZ IS NOT ALLOWED TO BE A VERTEX!\r
1748                 \r
1749                 \r
1750                 // IF THE INVERSE EVALUATION IS NOT SUCCESSFUL, THE NAME OF THE\r
1751                 // FUNCTION, WDSC, WILL BE ASSIGNED WITH:\r
1752                 //      .. Scalar Arguments ..\r
1753                 // DOUBLE COMPLEX C,ZZ\r
1754                 // DOUBLE PRECISION EPS,U\r
1755                 // INTEGER IOPT,M,N,NPTQ\r
1756         \r
1757                 // .. Array Arguments ..\r
1758                 // DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
1759                 // DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N),\r
1760                 //                  QWORK(NPTQ* (2* (M+N)+3))\r
1761                 \r
1762                 // .. Scalars in Common ..\r
1763                 // DOUBLE PRECISION DLAM\r
1764                 // INTEGER IU\r
1765                 \r
1766                 // .. Arrays in Common ..\r
1767                 // DOUBLE PRECISION UARY(8),VARY(3)\r
1768         \r
1769                 // .. Local Scalars ..\r
1770                 double result[] = new double[2];\r
1771                 double WFN[] = new double[2];\r
1772                 double WI[] = new double[2];\r
1773                 double ZS[] = new double[2];\r
1774                 double ZZ1[] = new double[2];\r
1775                 double denom;\r
1776                 // DOUBLE COMPLEX WFN,WI,ZS,ZZ1\r
1777                 int INZ[] = new int[1];\r
1778                 int KNZ[] = new int[1];\r
1779                 int I,IT,K;\r
1780         \r
1781                 // .. Local Arrays ..\r
1782                 double ZS0[][] = new double[50][2];\r
1783                 double ZS1[][] = new double[50][2];\r
1784                 double wout[];\r
1785                 double cr[] = new double[1];\r
1786                 double ci[] = new double[1];\r
1787                 double cr2[] = new double[1];\r
1788                 double ci2[] = new double[1];\r
1789                 double abswi;\r
1790                 double factor;\r
1791                 double zout[];\r
1792                 double absz;\r
1793                 // DOUBLE COMPLEX ZS0(50),ZS1(50)\r
1794                 \r
1795                 //     .. External Functions ..\r
1796                 //      DOUBLE COMPLEX WPROD,ZDSC\r
1797                 //      EXTERNAL WPROD,ZDSC\r
1798         \r
1799                 //     .. External Subroutines ..\r
1800                 //      EXTERNAL NEARZ\r
1801                 \r
1802                 //     .. Common blocks ..\r
1803                 //      COMMON /PARAM4/UARY,VARY,DLAM,IU\r
1804                 \r
1805                 result[0] = 0.0;\r
1806                 result[1] = 0.0;\r
1807                 \r
1808                 //  1.FORM WORK ARRAYS:\r
1809                 for (I = 0; I < M; I++) {\r
1810                     ZS0[I][0] = Z0[I][0];\r
1811                     ZS0[I][1] = Z0[I][1];\r
1812                 } // for (I = 0; I < M; I++)\r
1813                 for (I = 0; I < N; I++) {\r
1814                     ZS1[I][0] = Z1[I][0];\r
1815                     ZS1[I][1] = Z1[I][1];\r
1816                 } // for (I = 0; I < N; I++)\r
1817                 \r
1818                 // 2.GENERATE THE INITIAL VALUE FOR SOLVING ODE:\r
1819                 while (true) {\r
1820                     NEARZ(M,N,ZS0,ZS1,ALFA0,ZZ,KNZ,INZ);\r
1821                     if (INZ[0] == 2) {\r
1822                         System.err.println("THE INVERSE EVALUATION FAILED");    \r
1823                         System.err.println("AT POINT Z = (" + ZZ[0] + ", " + ZZ[1] + ")");\r
1824                         System.err.println("THE DESIGNED INVERSION PROCEDURE EXPERIENCED");\r
1825                         System.err.println("SOME ESSENTIAL DIFFICULATIES");\r
1826                         return result;\r
1827                     } // if (INZ[0] == 2)\r
1828                     if (INZ[0] == 0) {\r
1829                         if (KNZ[0] >= 2) {\r
1830                             WI[0] = (W0[KNZ[0]-1][0]+W0[KNZ[0]-2][0])/2.0;\r
1831                             WI[1] = (W0[KNZ[0]-1][1]+W0[KNZ[0]-2][1])/2.0;\r
1832                         }\r
1833                         else {\r
1834                             WI[0] = (W0[KNZ[0]-1][0]+W0[M-1][0])/2.0;\r
1835                     WI[1] = (W0[KNZ[0]-1][1]+W0[M-1][1])/2.0;\r
1836                         }\r
1837 \r
1838                     } // if (INZ[0] == 0)\r
1839                     else {\r
1840                         if (KNZ[0] >= 2) {\r
1841                               WI[0] = (W1[KNZ[0]-1][0]+W1[KNZ[0]-2][0])/2.0;\r
1842                               WI[1] = (W1[KNZ[0]-1][1]+W1[KNZ[0]-2][1])/2.0;\r
1843                         }\r
1844                         else {\r
1845                               WI[0] = (W0[KNZ[0]-1][0]+W1[N-1][0])/2.0;\r
1846                               WI[1] = (W0[KNZ[0]-1][1]+W1[N-1][1])/2.0;\r
1847                         }\r
1848 \r
1849                         denom = scm.zabs(WI[0], WI[1]);\r
1850                         WI[0] = U*WI[0]/denom;\r
1851                         WI[1] = U*WI[1]/denom;\r
1852                     } //else\r
1853 \r
1854                     ZS = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ,QWORK,1);\r
1855                 \r
1856                     // 3.SOLVE ODE INITIAL VALUE PROBLEM (ALONG LINE SEGMENT FROM\r
1857                     //   ZS TO ZZ) BY EULER'S SCHEME TO GENERATE THE INITIAL GUESS:\r
1858                     for (K = 1; K <= 20; K++) {\r
1859                           wout = WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1);\r
1860                           scm.zmlt(20.0*C[0], 20.0*C[1], wout[0], wout[1], cr, ci);\r
1861                           scm.zdiv(ZZ[0]-ZS[0], ZZ[1]-ZS[1], cr[0], ci[0], cr2, ci2);\r
1862                           WI[0] = WI[0] + cr2[0];\r
1863                           WI[1] = WI[1] + ci2[0];\r
1864                     } // for (K = 1; K <= 20; K++)\r
1865                     abswi = scm.zabs(WI[0], WI[1]);\r
1866                     if (abswi > 1.0) {\r
1867                         denom = abswi + Math.abs(1.0 - abswi);\r
1868                         WI[0] = WI[0]/denom;\r
1869                         WI[1] = WI[1]/denom;\r
1870                     }\r
1871                     if (abswi < U) {\r
1872                         factor = U*(abswi + Math.abs(U - abswi))/abswi;\r
1873                         WI[0] = factor * WI[0];\r
1874                         WI[1] = factor * WI[1];\r
1875                     }\r
1876                 \r
1877                     //  4.REFINE THE SOLUTION BY NEWTON'S ITERATION:\r
1878                     for (IT = 1; IT <= 15; IT++) {\r
1879                         zout = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,\r
1880                                              PHI1,NPTQ,QWORK,IOPT);\r
1881                         WFN[0] = ZZ[0] - zout[0];\r
1882                         WFN[1] = ZZ[1] - zout[1];\r
1883                         wout = WPROD(WI,M,N,U,W0,W1,ALFA0,ALFA1);\r
1884                         scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
1885                         scm.zdiv(WFN[0], WFN[1], cr[0], ci[0], cr2, ci2);\r
1886                         WI[0] = WI[0] + cr2[0];\r
1887                         WI[1] = WI[1] + ci2[0];\r
1888                         abswi = scm.zabs(WI[0], WI[1]);\r
1889                         if (abswi > 1.0) {\r
1890                                 denom = abswi + Math.abs(1.0 - abswi);\r
1891                                 WI[0] = WI[0]/denom;\r
1892                                 WI[1] = WI[1]/denom;\r
1893                             }\r
1894                             if (abswi < U) {\r
1895                                 factor = U*(abswi + Math.abs(U - abswi))/abswi;\r
1896                                 WI[0] = factor * WI[0];\r
1897                                 WI[1] = factor * WI[1];\r
1898                             }\r
1899                         if (scm.zabs(WFN[0],WFN[1]) < EPS) {\r
1900                                 result[0] = WI[0];\r
1901                                 result[1] = WI[1];\r
1902                                 ZZ1 = ZDSC(WI,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1,NPTQ, QWORK,1);\r
1903                                 break;\r
1904                         }\r
1905                         if (IT == 15) {\r
1906                                 // THE ITERATION FAILED TO MEET THE TOLERANCE IN 15 ITERATIONS.\r
1907                                 // TRY A DIFFERENT VERTEX AS A REFERENCE POINT:\r
1908                                 ZZ1[0] = 1.0 + ZZ[0];\r
1909                                 ZZ1[1] = 1.0 + ZZ[1];   \r
1910                         }\r
1911                     } // for (IT = 1; IT <= 15; IT++)\r
1912                     abswi = scm.zabs(WI[0], WI[1]);\r
1913                     absz = scm.zabs(ZZ[0]-ZZ1[0], ZZ[1]-ZZ1[1]);\r
1914                     if (abswi >= U && abswi <= 1.0 && absz <= 1.0E-3) {\r
1915                         return result;  \r
1916                     }\r
1917                     if (INZ[0] == 0) {\r
1918                         ZS0[KNZ[0]-1][0] = ZZ[0];\r
1919                 ZS0[KNZ[0]-1][1] = ZZ[1];\r
1920                     }\r
1921                     else {\r
1922                         ZS1[KNZ[0]-1][0] = ZZ[0];\r
1923                         ZS1[KNZ[0]-1][1] = ZZ[1];\r
1924                     }\r
1925 \r
1926                 } // while (true)\r
1927 \r
1928         }\r
1929 \r
1930 \r
1931         \r
1932         private void ANGLES(int MN, double Z01[][], double ALFA01[], int I01) {\r
1933             // Computes the interior angles of a doubly\r
1934                 // connected and bounded polygonal region.\r
1935                 // I01 := 0(outer polygon) := 1(inner polygon)\r
1936                 // MN: The number of vertices\r
1937                 // Z01: Array containing vertices\r
1938                 // ALFA01: Array containing interior turning angles on return\r
1939                 \r
1940                 int K, KL, KR;\r
1941                 double cr[] = new double[1];\r
1942                 double ci[] = new double[1];\r
1943                 \r
1944                 for (K = 1; K <= MN; K++) {\r
1945                         KL = (K+MN-2)%MN + 1;\r
1946                         KR = K%MN + 1;\r
1947                         scm.zdiv(Z01[KL-1][0] - Z01[K-1][0], Z01[KL-1][1] - Z01[K-1][1], \r
1948                                          Z01[KR-1][0] - Z01[K-1][0], Z01[KR-1][1] - Z01[K-1][1], cr, ci);\r
1949                         ALFA01[K-1] = Math.atan2(ci[0], cr[0])/Math.PI;\r
1950                         if (ALFA01[K-1] <= 0.0) {\r
1951                                 ALFA01[K-1] = ALFA01[K-1] + 2.0;\r
1952                         }\r
1953                 } // for (K = 1; K <= MN; K++)\r
1954                 if (I01 == 1) {\r
1955                         for (K = 1; K <= MN; K++) {\r
1956                             if (ALFA01[K-1] != 2.0) {\r
1957                                 ALFA01[K-1] = 2.0 - ALFA01[K-1];\r
1958                             }\r
1959                         } // for (K = 1; K <= MN; K++)\r
1960                 } // if (I01 == 1)\r
1961                 return;\r
1962         }\r
1963         \r
1964         private double ARGUM(double U1, double W01K[]) {\r
1965             // ----------------------------\r
1966             // COMPUTES THE ARGUMENT OF VECTOR W01K (I.E.,A COMPLEX NUMBER)\r
1967             // WHOSE LENGTH IS U1.\r
1968             // NOTE: BE SURE THAT THE SPECIFIED VALUE FOR\r
1969             // U1 MUST BE DOUBLE PRECISION.\r
1970         \r
1971         \r
1972             // .. Scalar Arguments ..\r
1973             // DOUBLE COMPLEX W01K\r
1974             // DOUBLE PRECISION U1\r
1975 \r
1976             // .. Local Scalars ..\r
1977             double PI,REW;\r
1978             double result;\r
1979         \r
1980             PI = Math.PI;\r
1981             REW = W01K[0];\r
1982             if (W01K[1] >= 0.0) {\r
1983                 result = Math.acos(REW/U1);\r
1984             }\r
1985             else {\r
1986                 result = Math.acos(-REW/U1) + PI;\r
1987             }\r
1988             return result;\r
1989         }\r
1990 \r
1991         \r
1992         private void NEARW(int M,int N, double W0[][], double W1[][],\r
1993                         double ALFA0[], double W[],int KNEAR[],int INEAR[]) {\r
1994             //   --------------------------------------------------\r
1995             // GIVEN PT W,THIS ROUTINE DETERMINES THE NEAREST PREVERTEX TO\r
1996         // THE PT W.ON RETURN INTEGER'INEAR'INDICATES THAT THE NEAREST\r
1997             // PT IS FOUND EITHER IN W0(INEAR=0) OR IN W1(INEAR=1).  KNEAR\r
1998             // CONTAINS THE CORRESPONDING INDEX IN W0 OR 1N W1.\r
1999             // NOTE: THE PREVERTICES CORRESPONDING TO INFINITE VERTICES\r
2000             //     WILL BE SKIPPED (OUTER CIRCLE ONLY).\r
2001             // SEE ROUTINE DSCSOLV FOR THE REST OF THE CALLING SEQUENCE.\r
2002         \r
2003         \r
2004             // .. Scalar Arguments ..\r
2005             // DOUBLE COMPLEX W\r
2006             // INTEGER INEAR,KNEAR,M,N\r
2007 \r
2008             // .. Array Arguments ..\r
2009             // DOUBLE COMPLEX W0(M),W1(N)\r
2010             // DOUBLE PRECISION ALFA0(M)\r
2011         \r
2012             // .. Local Scalars ..\r
2013             double D,DIST;\r
2014             int I;\r
2015         \r
2016             DIST = 2.0;\r
2017             for (I = 1; I <= M; I++) {\r
2018                 D = scm.zabs(W[0]-W0[I-1][0], W[1]-W0[I-1][1]);\r
2019                 if (D <= 1.0E-11 || D >= DIST ||\r
2020                     ALFA0[I-1] <= 0.0) {\r
2021                         continue;\r
2022                 }\r
2023                 KNEAR[0] = I;\r
2024                 DIST = D;\r
2025             } // for (I = 1; I <= M; I++)\r
2026             INEAR[0] = 0;\r
2027             for (I = 1; I <= N; I++) {\r
2028                 D = scm.zabs(W[0]-W1[I-1][0], W[1]-W1[I-1][1]);\r
2029                 if (D <= 1.0E-6 || D >= DIST) {\r
2030                         continue;\r
2031                 }\r
2032                 DIST = D;\r
2033                 KNEAR[0] = I;\r
2034                 INEAR[0] = 1;\r
2035             } // for (I = 1; I <= N; I++) \r
2036             return;\r
2037         }\r
2038 \r
2039         \r
2040         private void NEARZ(int M,int N, double Z0[][], double Z1[][],\r
2041                         double ALFA0[], double Z[], int KNZ[],int INZ[]) {\r
2042             //   --------------------------------------------------\r
2043             // GIVEN PT Z, THIS ROUTINE DETERMINES THE NEAREST VERTEX TO\r
2044             // THE PT Z.ON RETURN INTEGER'INZ'INDICATES THAT THE NEAREST\r
2045             // PT IS FOUND  EITHER IN Z0 (INZ=0 ) OR IN Z1 ( INZ=1). KNZ\r
2046             // CONTAINS THE CORRESPONDING INDEX IN Z0 OR 1N Z1.\r
2047             // (IF ON RETURN INZ=2, THAT MEANS NO APPROPRIATE VERTEX IS\r
2048             // FOUND, WHICH IS A DEVICE USED FOR INVERSE MAPPING ROUTINE.\r
2049             // NOTE: THE VERTICES CORRESPONDING TO INFINITE VERTICES\r
2050         //       WILL BE SKIPPED (OUTER POLYGON ONLY).\r
2051             // SEE ROUTINE DSCSOLV FOR THE REST OF CALLING SEQUENCE.\r
2052         \r
2053         \r
2054             // .. Scalar Arguments ..\r
2055             // DOUBLE COMPLEX Z\r
2056             // INTEGER INZ,KNZ,M,N\r
2057         \r
2058             // .. Array Arguments ..\r
2059             // DOUBLE COMPLEX Z0(M),Z1(N)\r
2060             // DOUBLE PRECISION ALFA0(M)\r
2061         \r
2062             // .. Local Scalars ..\r
2063             double D,DIST;\r
2064             int I;\r
2065         \r
2066             INZ[0] = 2;\r
2067             DIST = 99.0;\r
2068             for (I = 1; I <= M; I++) {\r
2069                 D = scm.zabs(Z[0]-Z0[I-1][0],Z[1]-Z0[I-1][1]);\r
2070                 if (D <= 1.0E-11 || D >= DIST ||\r
2071                      ALFA0[I-1] <= 0.0) {\r
2072                         continue;\r
2073                 }\r
2074                 KNZ[0] = I;\r
2075                 DIST = D;\r
2076                 INZ[0] = 0;\r
2077             } // for (I = 1; I <= M; I++)\r
2078             for (I = 1; I <= N; I++) {\r
2079                   D = scm.zabs(Z[0]-Z1[I-1][0],Z[1]-Z1[I-1][1]);\r
2080                   if (D <= 1.0E-11 || D >= DIST) {\r
2081                           continue;\r
2082                   }\r
2083                   DIST = D;\r
2084                   KNZ[0] = I;\r
2085                   INZ[0] = 1;\r
2086             } // for (I = 1; I <= N; I++)\r
2087             return;\r
2088         }\r
2089 \r
2090         \r
2091         private double DIST(int M,int N, double W0[][], double W1[][],\r
2092                         double W[], int KWA, int IC) {\r
2093         //   -----------------------------------\r
2094         //    DETERMINES THE DISTANCE FROM W TO THE NEAREST SINGULARITY\r
2095         //    OTHER THAN W ITSELF.( W COULD BE ONE OF THE PREVERTICES.)\r
2096         //    KWA IS THE INDEX OF W IN W0 (IF IC=0) OR W1 (IF IC=1), OR\r
2097         //    WK COULD BE 0 (IF ONE CHOOSES) WHEN W IS NOT A PREVERTEX.\r
2098     //\r
2099         //     .. Scalar Arguments ..\r
2100         //     DOUBLE COMPLEX W\r
2101         \r
2102         //     .. Array Arguments ..\r
2103         //      DOUBLE COMPLEX W0(M),W1(N)\r
2104                 \r
2105         //     .. Local Scalars ..\r
2106               double D;\r
2107               int I;\r
2108               double result;\r
2109         \r
2110               result = 2.0;\r
2111               for (I = 1; I <= M; I++) {\r
2112                   D = scm.zabs(W[0]-W0[I-1][0],W[1]-W0[I-1][1]);\r
2113                   if (I == KWA && IC == 0) {\r
2114                           continue;\r
2115                   }\r
2116                   result = Math.min(result,D);\r
2117               } // for (I = 1; I <= M; I++)\r
2118               for (I = 1; I <= N; I++) {\r
2119                   D = scm.zabs(W[0]-W1[I-1][0], W[1]-W1[I-1][1]);\r
2120                   if (I == KWA && IC == 1) {\r
2121                           continue;\r
2122                   }\r
2123                   result = Math.min(result,D);\r
2124               } // for (I = 1; I <= N; I++)\r
2125               return result;\r
2126 \r
2127         }\r
2128 \r
2129         \r
2130         private void CHECK(double ALFA0[], double ALFA1[], int M, int N, int ISHAPE) {\r
2131                 //   CHECKS IF THE INPUT DATA AND PARAMETERS ARE CORRECT.\r
2132                 //   NOTE1:    ANGLE-CHECKING MAKES SENSE ONLY IF ANGLES\r
2133                 //             ARE USER-PROVIDED.\r
2134                 //   NOTE2:    USERS ARE RESPONSIBLE FOR CHECKING THE FOLLOWING:\r
2135                 //   1. EACH COMPONENT OF THE OUTER POLYGON CONTAINS AT LEAST ONE\r
2136                 //      FINITE VERTEX.\r
2137                 //   2. Z1(N) IS THE CLOSEST (OR ALMOST CLOSEST)\r
2138                 //      INNER VERTEX TO Z1(M).\r
2139                 //   3. COUNTERCLOCKWISE ORIENTATION OF THE VERTICES.\r
2140                 \r
2141                 // ALFA0(M), ALFA1(N)\r
2142                 double EPS, SUM;\r
2143                 int K, MK;\r
2144                 \r
2145                 if (!((M >= 3) && (M <= 30) && (N <= 30) && (M+N <= 40))) {\r
2146                         System.err.println("M must be no less than 3");\r
2147                         System.err.println("M, N must be no greater than 30");\r
2148                         System.err.println("M+N must be no greater than 40");\r
2149                         System.exit(-1);\r
2150                 }\r
2151         \r
2152                 EPS = 0.00001;\r
2153                 SUM = 0.0;\r
2154                 MK = 0;\r
2155                 for (K = 0; K <= M; K++) {\r
2156                         if (ALFA0[K] > 0.0) {\r
2157                                 MK = MK + 1;\r
2158                         }\r
2159                         SUM = SUM + ALFA0[K];\r
2160                 } // for (K = 0; K <= M; K++)\r
2161                 if (!((MK < M && ISHAPE == 1) || (MK == M && ISHAPE == 0))) {\r
2162                         System.err.println("For finite regions ISHAPE must be 0");\r
2163                         System.err.println("and for infinite regions ISHAPE must be 1");\r
2164                         System.exit(-1);\r
2165                 }\r
2166                 \r
2167                 if (Math.abs(SUM - (M-2)) >= EPS) {\r
2168                         System.err.println("Some angles for outer polygon are wrong");\r
2169                         System.exit(-1);\r
2170                 }\r
2171                 \r
2172                 SUM = 0.0;\r
2173                 for (K = 0; K < N; K++) {\r
2174                         SUM = SUM + ALFA1[K];\r
2175                 }\r
2176                 if (Math.abs(SUM - (N+2)) >= EPS) {\r
2177                         System.err.println("Some angles for inner polygon are wrong");\r
2178                         System.exit(-1);\r
2179                 }\r
2180                 \r
2181                 if (ALFA0[0] <= 0.0) {\r
2182                         System.err.println("Z0[0] must be finite");\r
2183                         System.exit(-1);\r
2184                 }\r
2185                 \r
2186                 if (ALFA0[M-1] <= 0.0) {\r
2187                         System.err.println("Z0[M-1] must be finite");\r
2188                         System.exit(-1);\r
2189                 }\r
2190                 \r
2191                 if (ALFA0[M-3] <= 0.0) {\r
2192                         System.err.println("Z0[M-3] must be finite");\r
2193                         System.exit(-1);\r
2194                 }\r
2195                 \r
2196                 if (ALFA0[M-2] >= 2.0-EPS) {\r
2197                         System.err.println("Z0[M-2] must no be a re-entrant corner");\r
2198                         System.exit(-1);\r
2199                 }\r
2200                 \r
2201                 if ((ALFA0[M-2] >= 1.0-EPS) && (ALFA0[M-2] <= 1.0+EPS)) {\r
2202                         System.err.println("Z0[M-2] must not be an artificial corner");\r
2203                         System.exit(-1);\r
2204                 }\r
2205                 \r
2206                 System.out.println("Inputs have been checked with no error being found");\r
2207                 return;\r
2208         }\r
2209         \r
2210         private void DSCPRINT(int M,int N, double C[], double U, double W0[][], double W1[][],\r
2211                         double PHI0[], double PHI1[], double TOL, int NPTQ) {\r
2212             // ---------------------------------------------------------\r
2213             // PRINTS THE COMPUTED SC-PARAMETERS AND SOME\r
2214         // OTHER CONTROLING PARAMETERS FOR REFERENCE:\r
2215         \r
2216         \r
2217           // OPEN FILE DSCPACK FOR OUTPUT:\r
2218           //     .. Scalar Arguments ..\r
2219           // DOUBLE COMPLEX C\r
2220           // DOUBLE PRECISION TOL,U\r
2221           // INTEGER M,N,NPTQ\r
2222         \r
2223           // .. Array Arguments ..\r
2224           // DOUBLE COMPLEX W0(M),W1(N)\r
2225           // DOUBLE PRECISION PHI0(M),PHI1(N)\r
2226         \r
2227           //.. Local Scalars ..\r
2228           int K;\r
2229         \r
2230          System.out.println(" PARAMETERS DEFINING MAP:   (M = "+M+")   (N = "+N+")   (NPTQ = "+NPTQ+")   (TOL = "+TOL+")");\r
2231          System.out.println(" U = " + U);\r
2232          System.out.println(" C = (" + C[0] + ", " + C[1] + ")");\r
2233      for (K = 0; K < M; K++) {\r
2234          System.out.println("W0["+K+"] = (" + W0[K][0] + ", " + W0[K][1] + "   PHI0["+K+"] = " + PHI0[K]);\r
2235      }\r
2236      for (K = 0; K < N; K++) {\r
2237          System.out.println("W1["+K+"] = (" + W1[K][0] + ", " + W1[K][1] + "   PHI1["+K+"] = " + PHI1[K]);\r
2238      }\r
2239              \r
2240          return;\r
2241         \r
2242         }\r
2243         \r
2244         private void DSCTEST(int M,int N, double U, double C[], double W0[][], double W1[][],\r
2245                         double Z0[][], double Z1[][], double ALFA0[], double ALFA1[], int NPTQ, double QWORK[]) {\r
2246             //   ----------------------------------------------------------------\r
2247             // TESTS THE COMPUTED PARAMETERS FOR ACCURACY BY COMPUTING VERTICES\r
2248             // Z0(K) NUMERICALLY AND COMPARING WITH THE EXACT ONES. ON OUTPUT,\r
2249             // THE MAXIMUM AND MINIMUM ERRORS ACHIEVED WILL BE GIVEN TOGETHER\r
2250             // WITH THE CORRESPONDING INDICES KMAX AND KMIN RESPECTIVELY. SEE\r
2251             // ROUTINE DSCSOLV FOR THE CALLING SEQUENCE.\r
2252         \r
2253         \r
2254             // .. Scalar Arguments ..\r
2255             // DOUBLE COMPLEX C\r
2256             // DOUBLE PRECISION U\r
2257             // INTEGER M,N,NPTQ\r
2258         \r
2259             // .. Array Arguments ..\r
2260             // DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N)\r
2261             // DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3))\r
2262         \r
2263             // .. Scalars in Common ..\r
2264             // DOUBLE PRECISION DLAM\r
2265             // INTEGER IU\r
2266         \r
2267             // .. Arrays in Common ..\r
2268             // DOUBLE PRECISION UARY(8),VARY(3)\r
2269         \r
2270             // .. Local Scalars ..\r
2271                 double WA[] = new double[2];\r
2272                 double ZC[] = new double[2];\r
2273                 double ZTEST[] = new double[2];\r
2274                 double wout[];\r
2275                 double cr[] = new double[1];\r
2276                 double ci[] = new double[1];\r
2277             // DOUBLE COMPLEX WA,ZC,ZTEST\r
2278             double D,D1,DIST,ERRMAX,ERRMIN;\r
2279             int I,K;\r
2280             int KMAX = -1;\r
2281             int KMIN = -1;\r
2282             int KWA = 0;\r
2283         \r
2284             // .. External Functions ..\r
2285             // DOUBLE COMPLEX WQUAD\r
2286             // EXTERNAL WQUAD\r
2287         \r
2288             // .. Common blocks ..\r
2289             // COMMON /PARAM4/UARY,VARY,DLAM,IU\r
2290         \r
2291             ERRMAX = 0.0;\r
2292             ERRMIN = 99.0;\r
2293             for (K = 0; K < M - 1; K++) {\r
2294                 if (ALFA0[K] <= 0.0) {\r
2295                         continue;\r
2296                 }\r
2297                 DIST = 2.0;\r
2298                 for (I = 0; I < N; I++) {\r
2299                     D = scm.zabs(W0[K][0]-W1[I][0], W0[K][1]-W1[I][1]);\r
2300                     if (D >= DIST) {\r
2301                         continue;\r
2302                     }\r
2303                     DIST = D;\r
2304                     WA[0] = W1[I][0];\r
2305                     WA[1] = W1[I][1];\r
2306                     KWA = I;\r
2307                     ZC[0] = Z1[I][0];\r
2308                     ZC[1] = Z1[I][1];\r
2309                 } // for (I = 0; I < N; I++)\r
2310                 wout = WQUAD(WA,0.0,KWA,1,W0[K],0.0,K,0,0.0,M,N,U,\r
2311                                  W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,1);\r
2312                 scm.zmlt(C[0], C[1], wout[0], wout[1], cr, ci);\r
2313                 ZTEST[0] = ZC[0] + cr[0];\r
2314                 ZTEST[1] = ZC[1] + ci[0];\r
2315                 D1 = scm.zabs(Z0[K][0]-ZTEST[0],Z0[K][1]-ZTEST[1]);\r
2316                 if (D1 > ERRMAX) {\r
2317                     ERRMAX = D1;\r
2318                     KMAX = K;\r
2319                 } // if (D1 > ERRMAX)\r
2320 \r
2321                 if (D1 < ERRMIN) {\r
2322                     ERRMIN = D1;\r
2323                     KMIN = K;\r
2324                 } // if (D1 < ERRMIN)\r
2325 \r
2326             } // for (K = 0; K < M - 1; K++)\r
2327             System.out.println("ACCURACY TEST: ");\r
2328             System.out.println("MAXIMUM ERROR = " + ERRMAX + " ACHIEVED AT KMAX = " + KMAX);\r
2329             System.out.println("MINIMUM ERROR = " + ERRMIN + " ACHIEVED AT KMIN = " + KMIN);\r
2330             return;\r
2331     }\r
2332 \r
2333 \r
2334         \r
2335         private void HYBRD(int FCN,int N, double X[], double FVEC[], double XTOL, int MAXFEV,\r
2336                         int ML, int MU, double EPSFCN, double DIAG[], int MODE,\r
2337                     double FACTOR, int NPRINT,int INFO[], int NFEV[], double FJAC[][],\r
2338                     int LDFJAC, double R[],int LR, double QTF[], double WA1[],\r
2339                     double WA2[], double WA3[], double WA4[]) {\r
2340                 //     ***********\r
2341                 //\r
2342                 //     SUBROUTINE HYBRD\r
2343                 //\r
2344                 //     THE PURPOSE OF HYBRD IS TO FIND A ZERO OF A SYSTEM OF\r
2345                 //     N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION\r
2346                 //     OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A\r
2347                 //     SUBROUTINE WHICH CALCULATES THE FUNCTIONS. THE JACOBIAN IS\r
2348                 //     THEN CALCULATED BY A FORWARD-DIFFERENCE APPROXIMATION.\r
2349                 //\r
2350                 //     THE SUBROUTINE STATEMENT IS\r
2351                 //\r
2352                 //       SUBROUTINE HYBRD(FCN,N,X,FVEC,XTOL,MAXFEV,ML,MU,EPSFCN,\r
2353                 //                        DIAG,MODE,FACTOR,NPRINT,INFO,NFEV,FJAC,\r
2354                 //                        LDFJAC,R,LR,QTF,WA1,WA2,WA3,WA4)\r
2355                 //\r
2356                 //     WHERE\r
2357             //      dscfun = 1 for DSCFUN\r
2358                 //       FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH\r
2359                 //         CALCULATES THE FUNCTIONS. FCN MUST BE DECLARED\r
2360                 //         IN AN EXTERNAL STATEMENT IN THE USER CALLING\r
2361                 //         PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.\r
2362                 //\r
2363                 //         SUBROUTINE FCN(N,X,FVEC,IFLAG)\r
2364                 //         INTEGER N,IFLAG\r
2365                 //         DOUBLE PRECISION X(N),FVEC(N)\r
2366                 //         ----------\r
2367                 //         CALCULATE THE FUNCTIONS AT X AND\r
2368                 //         RETURN THIS VECTOR IN FVEC.\r
2369                 //         ---------\r
2370                 //         RETURN\r
2371                 //         END\r
2372                 //\r
2373                 //         THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS\r
2374                 //         THE USER WANTS TO TERMINATE EXECUTION OF HYBRD.\r
2375                 //         IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.\r
2376                 //\r
2377                 //       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER\r
2378                 //         OF FUNCTIONS AND VARIABLES.\r
2379                 //\r
2380                 //       X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN\r
2381                 //         AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X\r
2382                 //         CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.\r
2383                 //\r
2384                 //       FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS\r
2385                 //         THE FUNCTIONS EVALUATED AT THE OUTPUT X.\r
2386                 //\r
2387                 //       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION\r
2388                 //         OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE\r
2389                 //         ITERATES IS AT MOST XTOL.\r
2390                 //\r
2391                 //       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION\r
2392                 //         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST MAXFEV\r
2393                 //         BY THE END OF AN ITERATION.\r
2394                 //\r
2395                 //       ML IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES\r
2396                 //         THE NUMBER OF SUBDIAGONALS WITHIN THE BAND OF THE\r
2397                 //         JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET\r
2398                 //         ML TO AT LEAST N - 1.\r
2399                 //\r
2400                 //       MU IS A NONNEGATIVE INTEGER INPUT VARIABLE WHICH SPECIFIES\r
2401                 //         THE NUMBER OF SUPERDIAGONALS WITHIN THE BAND OF THE\r
2402                 //         JACOBIAN MATRIX. IF THE JACOBIAN IS NOT BANDED, SET\r
2403                 //         MU TO AT LEAST N - 1.\r
2404                 //\r
2405                 //       EPSFCN IS AN INPUT VARIABLE USED IN DETERMINING A SUITABLE\r
2406                 //         STEP LENGTH FOR THE FORWARD-DIFFERENCE APPROXIMATION. THIS\r
2407                 //         APPROXIMATION ASSUMES THAT THE RELATIVE ERRORS IN THE\r
2408                 //         FUNCTIONS ARE OF THE ORDER OF EPSFCN. IF EPSFCN IS LESS\r
2409                 //         THAN THE MACHINE PRECISION, IT IS ASSUMED THAT THE RELATIVE\r
2410                 //         ERRORS IN THE FUNCTIONS ARE OF THE ORDER OF THE MACHINE\r
2411                 //         PRECISION.\r
2412                 //\r
2413                 //       DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE\r
2414                 //         BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG\r
2415                 //         MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS\r
2416                 //         MULTIPLICATIVE SCALE FACTORS FOR THE VARIABLES.\r
2417                 //\r
2418                 //       MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE\r
2419                 //         VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2,\r
2420                 //         THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER\r
2421                 //         VALUES OF MODE ARE EQUIVALENT TO MODE = 1.\r
2422                 //\r
2423                 //       FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE\r
2424                 //         INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF\r
2425                 //         FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE\r
2426                 //         TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE\r
2427                 //         INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE.\r
2428                 //\r
2429                 //       NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED\r
2430                 //         PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE,\r
2431                 //         FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST\r
2432                 //         ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND\r
2433                 //         IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE\r
2434                 //         FOR PRINTING. IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS\r
2435                 //         OF FCN WITH IFLAG = 0 ARE MADE.\r
2436                 //\r
2437                 //       INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS\r
2438                 //         TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)\r
2439                 //         VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,\r
2440                 //         INFO IS SET AS FOLLOWS.\r
2441                 //\r
2442                 //         INFO = 0   IMPROPER INPUT PARAMETERS.\r
2443                 //\r
2444                 //         INFO = 1   RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES\r
2445                 //                    IS AT MOST XTOL.\r
2446                 //\r
2447                 //         INFO = 2   NUMBER OF CALLS TO FCN HAS REACHED OR EXCEEDED\r
2448                 //                    MAXFEV.\r
2449                 //\r
2450                 //         INFO = 3   XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN\r
2451                 //                    THE APPROXIMATE SOLUTION X IS POSSIBLE.\r
2452                 //\r
2453                 //         INFO = 4   ITERATION IS NOT MAKING GOOD PROGRESS, AS\r
2454                 //                    MEASURED BY THE IMPROVEMENT FROM THE LAST\r
2455                 //                    FIVE JACOBIAN EVALUATIONS.