Started debugging PARGEN.
[mipav.git] / mipav / src / gov / nih / mipav / model / algorithms / SymmsIntegralMapping.java
1 package gov.nih.mipav.model.algorithms;\r
2 \r
3 import java.io.File;\r
4 import java.io.IOException;\r
5 import java.io.RandomAccessFile;\r
6 import java.util.Scanner;\r
7 import java.util.Vector;\r
8 \r
9 import gov.nih.mipav.model.structures.ModelImage;\r
10 import gov.nih.mipav.model.structures.jama.GeneralizedEigenvalue;\r
11 import gov.nih.mipav.model.structures.jama.LinearEquations2;\r
12 import gov.nih.mipav.view.MipavUtil;\r
13 import gov.nih.mipav.view.Preferences;\r
14 \r
15 public class SymmsIntegralMapping extends AlgorithmBase {\r
16         String fileDir;\r
17 \r
18         // EPS returns the distance from 1.0 to the next larger double-precision\r
19         // number, that is, eps = 2^-52.\r
20         private double EPS;\r
21         // Filename to receive FORTRAN output code\r
22         private String FORTFL = null;\r
23         // Does the domain have any symmetry?\r
24         private boolean SYMTY;\r
25         // If SYMTY is true, are there any reflectional symmetries?\r
26         private boolean REFLN;\r
27         // If SYMTY is true, what are the coordinates of the center of symmetry?\r
28         private final int MNARC = 100;\r
29         private double CENSY[] = new double[2];\r
30         // If SYMTY is true, number of arcs on the fundamental boundary section\r
31         // If SYMTY is false, number of arcs on the boundary\r
32         // NARCS <= MNARC-1 = 99\r
33         private int NARCS;\r
34         // NUMDER is of length MNARC\r
35         // NUMDER is initially false\r
36         private boolean NUMDER[] = new boolean[MNARC];\r
37         // ARCTY is of length MNARC\r
38         // Type of arc, 1 = LINE SEGMENT, 2 = CIRCULAR ARC SEGMENT, 3 = CARTESIAN\r
39         // PARAMETRIC FUNCTION\r
40         // 4 = POLAR FUNCTION\r
41         private int ARCTY[] = new int[MNARC];\r
42         // STAPT is of length MNARC,2\r
43         // Initial point on line, initial point on circle, or initial point on curve\r
44         // IF (SYMTY) THEN\r
45         // WRITE(*,*) 'COORDINATES OF FINAL POINT ON THIS ARC?'\r
46         // STAPT[NARCS]=CMPLX(FINAL X point, FINAL Y point)\r
47         // ELSE\r
48         // STAPT[NARCS]=STAPT[1]\r
49         // ENDIF\r
50         private double STAPT[][] = new double[MNARC][2];\r
51         // Start with GMCO = 0\r
52         // For types 2, 3, and 4 increment GMCO and put\r
53         // GMCO=GMCO+1\r
54         // PGM[IA-1]=GMCO where IA is the number of the curve from 1 to NARCS\r
55         // For type 2 circular arc:\r
56         // RGM[GMCO-1]= X center of circle\r
57         // GMCO=GMCO+1\r
58         // RGM[GMCO-1]= Y center of circle\r
59         // GMCO=GMCO+1\r
60         // RGM[GMCO-1]=ALPHA*PI where alpha is the signed angle subtended at center\r
61         // in units of PI\r
62         // ALPHA is positive for CCW and negative for CW.\r
63         // For type 3 CARTESIAN PARAMETRIC FUNCTION:\r
64         // RGM[GMCO-1]= initial parameter value\r
65         // GMCO=GMCO+1\r
66         // RGM[GMCO-1]= final parameter value\r
67         // For type 4 polar function:\r
68         // RGM[GMCO-1]= (initial polar angle in units of PI) *PI\r
69         // GMCO=GMCO+1\r
70         // RGM[GMCO-1]= (final polar angle in units of PI) *PI\r
71         // PGM is of length MNARC\r
72         private int PGM[] = new int[MNARC];\r
73         // RGM is of length 3*MNARC\r
74         private double RGM[] = new double[3 * MNARC];\r
75         // Start with TXCO = 0\r
76         // For types 3 and 4\r
77         // For each IA = 1 to NARCS do J = 1,2\r
78         // TXCO=TXCO+1\r
79         // PTX[IA-1+(J-1)*MNARC]=TXCO\r
80         // DEFN[TXCO-1][0]= TXT for real part\r
81         // DEFN[TXCO-1][1] = TXT for imaginary part\r
82         // where TXT = A COMPLEX EXPRESSION for TYPE 3 and a REAL EXPRESSION FOR\r
83         // TYPE 4\r
84         // for J = 1 and TYPE = 3 JAVA EXPRESSION FOR PARFUN\r
85         // for J = 2 and TYPE = 3 JAVA EXPRESSION FOR DPARFN\r
86         // for J = 1 and TYPE = 4 JAVA EXPRESSION FOR RADIUS\r
87         // for J = 2 and TYPE = 4 JAVA EXPRESSION FOR RADIUS DERIVATIVE\r
88         // PTX is of length 2*MNARC\r
89         private int PTX[] = new int[2 * MNARC];\r
90         private int NTX[] = new int[2 * MNARC];\r
91         // CHARACTER DEFN(MNARC*2)*72\r
92         // Holds text for real for types 3 and 4 and imaginary parts for type 3\r
93         // Start imaginary text with ui. All text following ui is imaginary.\r
94         String DEFN[] = new String[2 * MNARC];\r
95         private boolean traditionalInput = true;\r
96         Scanner input = new Scanner(System.in);\r
97         private double zzset[][] = new double[400][2];\r
98         private int IBNDS[] = new int[5];\r
99         private int IGEOM[] = new int[4]; // Written in original JAPHYC\r
100         private int PARNT[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
101         private double RGEOM[] = new double[2]; // Written in original JAPHYC\r
102         private double HALEN[] = new double[IBNDS[0]]; // Written in original JAPHYC\r
103         private double MIDPT[] = new double[IBNDS[0]]; // Written in original JAPHYC\r
104         private double VTARG[] = new double[IBNDS[0]]; // Written in original JPAHYC\r
105         private int ISNPH[] = new int[6]; // Written in original JAPHYC\r
106         private int DGPOL[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
107         private int JATYP[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
108         private int LOSUB[] = new int[IBNDS[0]]; // Written in original JAPHYC\r
109         // Parts of RSNPH\r
110         private int NQPTS;\r
111         private int NJIND = NARCS + 1;\r
112         private int TNGQP = NQPTS * NJIND;\r
113         private int MNEQN;\r
114         private double ACOEF[] = new double[TNGQP]; // Written in original JAPHYC\r
115         private double BCOEF[] = new double[TNGQP]; // Written in original JAPHYC\r
116         private double AICOF[] = new double[TNGQP]; // Written in original JAPHYC\r
117         private double BICOF[] = new double[TNGQP]; // Written in original JAPHYC\r
118         private double QUPTS[] = new double[TNGQP]; // Written in original JAPHYC\r
119         private double QUWTS[] = new double[TNGQP]; // Written in original JAPHYC\r
120         private double H0VAL[] = new double[NJIND]; // Written in original JAPHYC\r
121         private double HIVAL[] = new double[NJIND]; // Written in original JAPHYC\r
122         private double JACIN[] = new double[NJIND]; // Written in original JAPHYC\r
123         private double ERARC[] = new double[IBNDS[0]]; // Written in original JAPHYC\r
124         private double BCFSN[] = new double[MNEQN]; // Written in original JAPHYC\r
125         private double SOLUN[] = new double[MNEQN]; // Written in original JAPHYC\r
126         // Parts of IWORK\r
127         private int IPIVT[] = new int[MNEQN];\r
128         private int LOQSB[] = new int[IBNDS[1]];\r
129         private int NQUAD[] = new int[IBNDS[1]];\r
130         private int HISUB[] = new int[IBNDS[0]];\r
131         private int LOTES[] = new int[IBNDS[0]];\r
132         private int HITES[] = new int[IBNDS[0]];\r
133         private int AXION[] = new int[IBNDS[0]];\r
134         private int NEWDG[] = new int[IBNDS[0]];\r
135         private int ICOPY[] = new int[IBNDS[0]];\r
136         private int LOOLD[] = new int[IBNDS[0]];\r
137         private int HIOLD[] = new int[IBNDS[0]];\r
138         // Parts of RWORK\r
139         private double WORK2[] = new double[MNEQN];\r
140         private double COLPR[] = new double[MNEQN];\r
141         private double A1COF[] = new double[IBNDS[1]];\r
142         private double B1COF[] = new double[IBNDS[1]];\r
143         private double TOLOU[] = new double[IBNDS[1]];\r
144         private double XIDST[] = new double[2 * IBNDS[1]];\r
145         private double XENPT[] = new double[IBNDS[2]];\r
146         private double QCOMX[] = new double[IBNDS[3]];\r
147         private double QCOMW[] = new double[IBNDS[3]];\r
148         private double RCOPY[] = new double[IBNDS[0]];\r
149         private double NEWHL[] = new double[IBNDS[0]];\r
150         private double AQCOF[] = new double[TNGQP];\r
151         private double BQCOF[] = new double[TNGQP];\r
152         private double CQCOF[] = new double[TNGQP];\r
153         private double COLSC[] = new double[TNGQP];\r
154         private double RIGLL[] = new double[TNGQP];\r
155         private double WORK[] = new double[2 * NQPTS];\r
156         private double DIAG[] = new double[NQPTS];\r
157         private double SDIAG[] = new double[NQPTS];\r
158         private double WORKT[] = new double[2 * NQPTS * NQPTS];\r
159         private double WORKQ[] = new double[NQPTS * NQPTS];\r
160         // Parts of ZWORK\r
161         private double ZCOLL[][] = new double[MNEQN][2];\r
162         private double XIVAL[][] = new double[2 * IBNDS[1]][2];\r
163         // Parts of LWORK\r
164         private boolean NEWQU[] = new boolean[IBNDS[1]];\r
165         private boolean LCOPY[] = new boolean[IBNDS[0]];\r
166         private boolean PNEWQ[] = new boolean[IBNDS[0]];\r
167         private boolean LNSEG[] = new boolean[IBNDS[0]];\r
168 \r
169         // CENTR - COMPLEX\r
170         // THE POINT IN THE PHYSICAL PLANE THAT IS TO BE\r
171         // MAPPED TO THE CENTRE OF THE UNIT DISC. FOR\r
172         // EXTERIOR DOMAINS CENTR MUST BE SOME POINT IN THE\r
173         // COMPLEMENTARY INTERIOR PHYSICAL DOMAIN.\r
174         // IN CASE ABS(ISYGP).GT.1 THEN CENTR MUST ALSO BE\r
175         // A CENTRE OF SYMMETRY FOR THE PHYSICAL DOMAIN.\r
176         private double CENTR[] = new double[2]; // Written in original JAPHYC\r
177         // INTER - LOGICAL\r
178         // TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE\r
179         // OTHERWISE.\r
180         private boolean INTER; // Written in in original JAPHYC\r
181 \r
182         // COMMON /FNDEF/\r
183         private double BETA;\r
184         private double A1;\r
185         private double B1;\r
186         private double P0VAL;\r
187         private double SCALE;\r
188         private int TYPE;\r
189 \r
190         // From GQPHYC:\r
191         // IQUPH - INTEGER ARRAY\r
192         // AN INTEGER VECTOR OF SIZE AT LEAST 2*IBNDS(1) + 4,\r
193         // WHERE IBNDS(1) (=IGEOM(4)) IS THE VALUE PREVIOUSLY\r
194         // SUPPLIED TO JAPHYC; IQUPH STORES POINTERS TO\r
195         // ACCESS RQUPH AND ZQUPH.\r
196         private int IQUPH[] = new int[4];\r
197         private int LQSBF[] = new int[IBNDS[0]];\r
198         private int NPPQF[] = new int[IBNDS[0]];\r
199         // MQUPH - INTEGER\r
200         // THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
201         // IN THE FINAL GLOBAL RULE. (THE VALUE OF THIS\r
202         // ARGUMENT IS LINKED TO THOSE OF ARGUMENTS NQPTS\r
203         // AND IBNDS(1) PREVIOUSLY SUPPLIED TO JAPHYC VIA\r
204         // MQUPH <= (MQIN1-1)*NQPTS*IBNDS(1))\r
205         private int MQUPH;\r
206         // RQUPH - REAL ARRAY\r
207         // A REAL VECTOR OF SIZE AT LEAST 3*MQUPH + 1; STORES\r
208         // THE REAL QUADRATURE DATA.\r
209         private double RQUPH[] = new double[1];\r
210         private double TPPQF[] = new double[MQUPH];\r
211         private double TRRAD[] = new double[MQUPH];\r
212         private double WPPQF[] = new double[MQUPH];\r
213         // ZQUPH - COMPLEX ARRAY\r
214         // A COMPLEX VECTOR OF SIZE AT LEAST MQUPH + 1;\r
215         // STORES THE QUADRATURE POINTS ON THE PHYSICAL\r
216         // BOUNDARY.\r
217         private double FACTR[] = new double[2]; // At first location of ZQUPH\r
218         private double ZPPQF[][] = new double[MQUPH][2];\r
219         // From JACANP:\r
220         // ISNCA - INTEGER ARRAY\r
221         // AN INTEGER VECTOR OF SIZE AT LEAST\r
222         // 4*IBNDS(1) + 6 ;\r
223         // ISNCA MAINLY STORES POINTERS TO ACCESS RSNCA AND\r
224         // ZSNCA.\r
225         private int ISNCA[] = new int[6];\r
226         private int DGPOC[] = new int[IBNDS[0]];\r
227         private int JTYPC[] = new int[IBNDS[0]];\r
228         private int LSUBC[] = new int[IBNDS[0]];\r
229         private int PRNSA[] = new int[IBNDS[0]];\r
230         // RSNCA - REAL ARRAY\r
231         // A REAL VECTOR OF SIZE AT LEAST\r
232         // 2*IBNDS(1) + (4 + 6*NQPTS)*(NARCS + 1) + 2,\r
233         // WHERE NARCS, NQPTS ARE INPUT ARGUMENTS TO JAPHYC.\r
234         // (NOTE: NARCS=IGEOM(1), NQPTS=IGEOM(2))\r
235         // STORES DATA RELATING TO THREE-TERM RECURRENCE\r
236         // SCHEMES, ELEMENTARY GAUSS-JACOBI QUADRATURE RULES,\r
237         // AND THE ARGUMENTS OF SUB-ARC ENDPOINTS ON THE UNIT\r
238         // DISC.\r
239         private double RSNCA[] = new double[1];\r
240         private double ACOFC[] = new double[TNGQP];\r
241         private double BCOFC[] = new double[TNGQP];\r
242         private double AICOC[] = new double[TNGQP];\r
243         private double BICOC[] = new double[TNGQP];\r
244         private double QUPTC[] = new double[TNGQP];\r
245         private double QUWTC[] = new double[TNGQP];\r
246         private double H0VLC[] = new double[NJIND];\r
247         private double HIVLC[] = new double[NJIND];\r
248         private double JAINC[] = new double[NJIND];\r
249         private double COARG[] = new double[NJIND];\r
250         private double PHPAS[] = new double[IBNDS[0]];\r
251         private double VARGC[] = new double[IBNDS[0]];\r
252         // ZSNCA - COMPLEX ARRAY\r
253         // Also called CSNCA elsewhere\r
254         // A COMPLEX VECTOR OF SIZE AT LEAST 2*IBNDS(2) + 1;\r
255         // STORES THE JACOBI COEFFICIENTS FOR THE COMPLEX\r
256         // (INVERSE) BOUNDARY CORRESPONDENCE FUNCTION AND\r
257         // ITS DERIVATIVE.\r
258         private double ZSNCA[] = new double[2]; // At first location of ZSNCA\r
259         private double BFSNC[][] = new double[IBNDS[1]][2];\r
260         private double SOLNC[][] = new double[IBNDS[1]][2];\r
261         // From GQCANP:\r
262         // MQIN1 - INTEGER\r
263         // DEFINES THE NUMBER OF PANELS ALLOWED IN A\r
264         // COMPOSITE RULE. SPECIFICALLY, MQIN1 = 1 + (THE\r
265         // MAXIMUM NUMBER OF PANELS IN A COMPOSITE RULE FOR\r
266         // A SINGLE SUB-ARC ON THE BOUNDARY)\r
267         private int MQIN1;\r
268         // MQUCA - INTEGER\r
269         // THE MAXIMUM NUMBER OF QUADRATURE POINTS ALLOWED\r
270         // IN THE FINAL GLOBAL RULE. THE VALUE OF THIS\r
271         // ARGUMENT IS LINKED TO THOSE OF ARGUMENTS NQPTS\r
272         // AND IBNDS(1) PREVIOUSLY SUPPLIED TO JACANP VIA\r
273         // MQUCA <= (MQIN1-1)*NQPTS*IBNDS(1). (NOTE THAT\r
274         // NQPTS = ISNCA(2) 'JACANP'IBNDS(1) =ISNCA(5) )\r
275         private int MQUCA = (MQIN1 - 1) * NQPTS * IBNDS[0];\r
276         // IQUCA - INTEGER ARRAY\r
277         // AN INTEGER VECTOR OF SIZE AT LEAST 2*IBNDS(1) + 4,\r
278         // WHERE IBNDS(1) (=ISNCA(5)) IS THE VALUE PREVIOUSLY\r
279         // SUPPLIED TO JACANP; IQUCA MAINLY STORES POINTERS\r
280         // TO ACCESS ZQUCA.\r
281         private int IQUCA[] = new int[4];\r
282         private int LQSBG[] = new int[IBNDS[0]];\r
283         private int NPPQG[] = new int[IBNDS[0]];\r
284         // ZQUCA - COMPLEX ARRAY\r
285         // A COMPLEX VECTOR OF SIZE AT LEAST 2*MQUCA+1;\r
286         // STORES THE QUADRATURE POINTS AND WEIGHTS.\r
287         private double ZQUCA[] = new double[2]; // At first location of ZQUCA\r
288         private double WPPQG[][] = new double[MQUCA][2];\r
289         private double ZPPQG[][] = new double[MQUCA][2];\r
290         // COMMON /DSDTDA/PT,MD,HL\r
291         private int PT;\r
292         private double MD;\r
293         private double HL;\r
294         //From LEVCUR:\r
295         Vector<Double> Contour[]; // x y pairs\r
296         Vector<Double>Ray[]; // x y pairs\r
297         // From TSTPLT\r
298         Vector<Double>Boundary; // x y pairs\r
299 \r
300         public SymmsIntegralMapping() {\r
301 \r
302         }\r
303 \r
304         public SymmsIntegralMapping(ModelImage destImg, ModelImage srcImg, String FORTFL, boolean SYMTY, boolean REFLN,\r
305                         double CENSY[], int NARCS, boolean NUMDER[], int ARCTY[], double STAPT[][], int PGM[], double RGM[],\r
306                         int PTX[], int NTX[], String DEFN[]) {\r
307                 super(destImg, srcImg);\r
308                 this.FORTFL = FORTFL;\r
309                 this.SYMTY = SYMTY;\r
310                 this.REFLN = REFLN;\r
311                 this.CENSY = CENSY;\r
312                 this.NARCS = NARCS;\r
313                 this.NUMDER = NUMDER;\r
314                 this.ARCTY = ARCTY;\r
315                 this.STAPT = STAPT;\r
316                 this.PGM = PGM;\r
317                 this.RGM = RGM;\r
318                 this.PTX = PTX;\r
319                 this.NTX = NTX;\r
320                 this.DEFN = DEFN;\r
321         }\r
322 \r
323         public void runAlgorithm() {\r
324                 fileDir = srcImage.getFileInfo(0).getFileDirectory();\r
325 \r
326                 // eps returns the distance from 1.0 to the next larger double-precision\r
327                 // number, that is, eps = 2^-52.\r
328                 // epsilon = D1MACH(4)\r
329                 // Machine epsilon is the smallest positive epsilon such that\r
330                 // (1.0 + epsilon) != 1.0.\r
331                 // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
332                 // epsilon = 2.2204460e-16\r
333                 // epsilon is called the largest relative spacing\r
334                 EPS = 1.0;\r
335                 double neweps = 1.0;\r
336 \r
337                 while (true) {\r
338 \r
339                         if (1.0 == (1.0 + neweps)) {\r
340                                 break;\r
341                         } else {\r
342                                 EPS = neweps;\r
343                                 neweps = neweps / 2.0;\r
344                         }\r
345                 } // while(true)\r
346         }\r
347 \r
348         public void PARGEN() {\r
349                 // .......................................................................\r
350                 // AN INTERACTIVE PREPROCESSOR TO HELP THE C O N F P A C K USER\r
351                 // GENERATE THE FORTRAN CODE DEFINING THE BOUNDARY PARAMETRISATION\r
352                 // AND ITS DERIVATIVE.\r
353                 //\r
354                 // THE FOLLOWING CONVENTIONS ARE ASSUMED:\r
355                 //\r
356                 // ** IF A QUESTION REQUIRES A YES OR NO ANSWER, THE DETECTION OF 'Y' OR\r
357                 // 'y' IN THE FIRST 6 INPUT CHARACTERS IS TAKEN AS YES, ANYTHING ELSE\r
358                 // AS NO.\r
359                 //\r
360                 // ** WHEN ASKED FOR COORDINATES, THIS ALWAYS MEANS CARTESIAN COORDIN-\r
361                 // ATES AND THESE SHOULD BE SUPPLIED AS TWO REAL NUMBERS, AS EITHER\r
362                 //\r
363                 // X,Y OR X Y\r
364                 //\r
365                 // WITHOUT PARENTHESES.\r
366                 //\r
367                 // ** FOUR TYPES OF ARCS ARE CURRENTLY TREATED, WITH NUMERICAL CODES TO\r
368                 // DENOTE THE TYPE AS FOLLOWS:\r
369                 // 1:= LINE SEGMENT.\r
370                 // 2:= CIRCULAR ARC SEGMENT.\r
371                 // CONVENTIONS:\r
372                 // 1 - THE ANGLE SUBTENDED AT THE CENTRE IS POSITIVE FOR\r
373                 // ANTICLOCKWISE TRAVERSAL OF THE ARC AND NEGATIVE FOR\r
374                 // CLOCKWISE TRAVERSAL.\r
375                 // 3:= THE USER IS ASKED TO SUPPLY THE FORTRAN 77 ARITHMETIC\r
376                 // EXPRESSIONS WHICH DEFINE THE CARTESIAN PARAMETRIC FUNCTION\r
377                 // AND THE DERIVATIVE OF THIS FUNCTION.\r
378                 // CONVENTIONS:\r
379                 // 1 - THE PARAMETER MUST BE DENOTED BY T.\r
380                 // 2 - THE REAL CONSTANT PI=3.14159.. AND THE COMPLEX CONSTANT\r
381                 // UI=(0.0,1.0) MAY BE USED IN THE ARITHMETIC EXPRESSIONS.\r
382                 // 4:= THE USER IS ASKED TO SUPPLY THE FORTRAN 77 ARITHMETIC\r
383                 // EXPRESSIONS WHICH DEFINE THE POLAR COORDINATE AS A FUNCTION\r
384                 // OF POLAR ANGLE AND THE DERIVATIVE OF THIS FUNCTION.\r
385                 // CONVENTIONS:\r
386                 // 1 - THE POLAR ANGLE MUST BE DENOTED BY T.\r
387                 // 2 - THE REAL CONSTANT PI=3.14159.. MAY BE USED IN THE ARITH-\r
388                 // METIC EXPRESSIONS.\r
389                 // 3 - PARGEN ASSIGNS THE EXPRESSION FOR THE RADIUS TO THE\r
390                 // COMPLEX VARIABLE ZRAD; IF REQUIRED THE USER MAY THERE-\r
391                 // FORE USE THE VARIABLE ZRAD IN THE EXPRESSION FOR THE\r
392                 // DERIVATIVE OF THE RADIUS WRT POLAR ANGLE\r
393                 // IN ADDITION, FOR TYPES 3 AND 4, THE FOLLOWING CONVENTIONS HOLD:\r
394                 // 1 - ONLY USE UP TO 66 CHARACTERS PER LINE.\r
395                 // 2 - IF THE EXPRESSION OCCUPIES MORE THAN ONE LINE THEN THERE\r
396                 // IS NO NEED TO SUPPLY ANY CONTINUATION CHARACTER.\r
397                 // 3 - ONLY USE THOSE FORTRAN 77 INTRINSIC MATHS FUNCTIONS\r
398                 // WHICH ACCEPT COMPLEX ARGUMENTS AND ARE ANALYTIC; I.E.,\r
399                 // IN STANDARD FORTRAN,\r
400                 // SQRT, EXP, LOG, SIN, COS\r
401                 // 4 - THE WHOLE EXPRESSION SHOULD BE TERMINATED WITH A\r
402                 // REPEATED DIVISION SIGN, I.E. //.\r
403                 //\r
404                 // ** THE CODE ISN'T VERY ROBUST IN THAT THERE IS LITTLE PROVISION FOR\r
405                 // INTERACTIVE CORRECTION OF ERRORS. HOWEVER, ALL THE USER'S INPUT\r
406                 // IS AUTOMATICALLY OUTPUT TO A FILE NAMED pgenin. IF THE USER\r
407                 // REALISES THAT AN ITEM HAS BEEN INPUT INCORRECTLY, THE BEST POLICY\r
408                 // IS TO CARRY ON TO THE END OF THE INPUT PHASE AND THEN TERMINATE\r
409                 // THE EXECUTION; THE FILE pgenin CAN BE EDITED, RENAMED AND\r
410                 // SUBMITTED AS STANDARD INPUT FOR A SECOND RUN OF PARGEN.\r
411                 // (THE SUGGESTION TO RENAME IS TO AVOID ANY POSSIBLE DIFFICULTY\r
412                 // ARISING FROM READING THE FILE AS STANDARD INPUT WHILST ALSO\r
413                 // WRITING OUTPUT TO THE SAME FILE.)\r
414                 // IF THE USER REALISES ONLY AT A LATER STAGE (E.G. IN PLOTTING THE\r
415                 // BOUNDARY) THAT AN ITEM HAS BEEN INPUT INCORRECTLY THEN pgenin\r
416                 // MAY STILL BE AVAILABLE FOR RE-USE AS ABOVE.\r
417                 //\r
418                 // SUBROUTINES OR FUNCTIONS NEEDED\r
419                 // - THE CONFPACK LIBRARY.\r
420                 // - THE REAL FUNCTION R1MACH.\r
421                 //\r
422                 // .......................................................................\r
423                 // AUTHOR: DAVID HOUGH, COVENTRY POLYTECHNIC, UK\r
424                 // LAST UPDATE: 15 FEB 1991\r
425                 // .......................................................................\r
426                 //\r
427                 // LOCAL VARIABLES\r
428                 //\r
429                 int GMCO, IA, I, J, L, SW, TXCO;\r
430                 int TYPE = 1;\r
431                 int IER[] = new int[1];\r
432                 int ORDRG[] = new int[1];\r
433                 int ORDSG[] = new int[1];\r
434                 double X = 0.0;\r
435                 double Y = 0.0;\r
436                 double ALPHA, PI;\r
437                 // COMPLEX CENSY,RTUNI,U2\r
438                 double RTUNI[] = new double[2];\r
439                 double U2[] = new double[2];\r
440                 // CHARACTER TXT*72,TABC*6,FORTFL*72,CH*2,SIG(10)*2,WID(10)*2,REDD*6,\r
441                 // +FMT1*8,FMT2*9\r
442                 String TXT;\r
443                 String CH;\r
444                 String SIG[] = new String[] { "7", "8", "9", "10", "11", "12", "13", "14", "15", "16" };\r
445                 String WID[] = new String[] { "15", "16", "17", "18", "19", "20", "21", "22", "23", "24" };\r
446                 String REDD;\r
447                 // String FMT1;\r
448                 // String FMT2;\r
449 \r
450                 // PARAMETER (MNARC=100,TABC=' +',CHNL=20,CHIN=21)\r
451                 final String TABC = "     +";\r
452                 final int CHNL = 20;\r
453                 // final int CHIN = 21;\r
454                 String line;\r
455                 String tokens[];\r
456 \r
457                 File file;\r
458                 RandomAccessFile raFile = null;\r
459                 boolean validInput;\r
460 \r
461                 // EXTERNAL CHRIN,HEADER,R1MACH,SYINF1,WRFUN1,WRFUN2,WRHEAD,WRSYM1,\r
462                 // +WRSYM2,WRSYM3,WRTAIL\r
463 \r
464                 WRHEAD(6, 0, null);\r
465 \r
466                 PI = Math.PI;\r
467 \r
468                 // **** DETERMINE NUMBER OF SIGNIFICANT FIGURES REQUIRED TO MATCH\r
469                 // MACHINE\r
470                 // **** PRECISION AND SET UP POINTER SW TO SIG AND WID\r
471 \r
472                 SW = (int) (-Math.log10(EPS)) + 2;\r
473                 if (SW <= 7) {\r
474                         SW = 1;\r
475                 } else if (SW >= 16) {\r
476                         SW = 10;\r
477                 } else {\r
478                         SW = SW - 6;\r
479                 }\r
480 \r
481                 // **** SET UP THE EDIT DESCRIPTOR AND FORMAT SPECIFICATION FOR FLOATING\r
482                 // **** POINT OUTPUT\r
483 \r
484                 REDD = "E" + WID[SW - 1] + "." + SIG[SW - 1];\r
485                 // FMT1="("+REDD+")";\r
486                 // FMT2="(2"+REDD+")";\r
487 \r
488                 if (traditionalInput) {\r
489                         System.out.println("ENTER FILENAME TO RECEIVE OUTPUT JAVA CODE");\r
490                         FORTFL = input.next();\r
491                 } // if (traditionalInput)\r
492 \r
493                 // **** WRITE THE SOURCE CODE FOR PARFUN\r
494                 fileDir = "C:\\conformal mapping\\CONFPACK\\";\r
495                 file = new File(fileDir + FORTFL);\r
496                 \r
497                 try {\r
498                         raFile = new RandomAccessFile(file, "rw");\r
499                 } catch (IOException e) {\r
500                         MipavUtil.displayError("IOException " + e + " on raFile = new RandomAccessFile(file, rw)");\r
501                         System.exit(-1);\r
502                 }\r
503                 // Necessary so that if this is an overwritten file there isn't any\r
504                 // junk at the end\r
505                 try {\r
506                         raFile.setLength(0);\r
507                 } catch (IOException e) {\r
508                         MipavUtil.displayError("IOException " + e + " on raFile.setLength(0)");\r
509                         System.exit(-1);\r
510                 }\r
511 \r
512                 // OPEN(CHNL,FILE=FORTFL)\r
513 \r
514                 if (traditionalInput) {\r
515                         validInput = false;\r
516                         while (!validInput) {\r
517                                 System.out.println("DOES THE DOMAIN HAVE ANY SYMMETRY (Y/N)?");\r
518                                 String sym = input.next();\r
519                                 String firstSym = sym.substring(0, 1);\r
520                                 if (firstSym.equalsIgnoreCase("Y")) {\r
521                                         SYMTY = true;\r
522                                         validInput = true;\r
523                                 } else if (firstSym.equalsIgnoreCase("N")) {\r
524                                         SYMTY = false;\r
525                                         validInput = true;\r
526                                 } else {\r
527                                         System.out.println(sym + " is not a valid response");\r
528                                 }\r
529                         } // while (!validInput)\r
530 \r
531                         if (SYMTY) {\r
532                                 validInput = false;\r
533                                 while (!validInput) {\r
534                                         System.out.println("ARE THERE ANY REFLECTIONAL SYMMETRIES (Y/N)?");\r
535                                         String ref = input.next();\r
536                                         String firstRef = ref.substring(0, 1);\r
537                                         if (firstRef.equalsIgnoreCase("Y")) {\r
538                                                 REFLN = true;\r
539                                                 validInput = true;\r
540                                         } else if (firstRef.equalsIgnoreCase("N")) {\r
541                                                 REFLN = false;\r
542                                                 validInput = true;\r
543                                         } else {\r
544                                                 System.out.println(ref + " is not a valid response");\r
545                                         }\r
546                                 } // while (!validInput)\r
547 \r
548                                 validInput = false;\r
549                                 while (!validInput) {\r
550                                     System.out.println("What are the coordinates of the center of symmetry (X,Y)?");\r
551                                     try {\r
552                                         line = input.next();\r
553                                         tokens = line.split(",");\r
554                                         CENSY[0] = Double.valueOf(tokens[0]);\r
555                                         CENSY[1] = Double.valueOf(tokens[0]);\r
556                                         validInput = true;\r
557                                     }\r
558                                     catch (Exception e) {};\r
559                                 }\r
560                                 validInput = false;\r
561                                 while (!validInput) {\r
562                                         System.out.print("How many arcs are there on the fundamental boundary section?");\r
563                                         NARCS = input.nextInt();\r
564                                         if (NARCS <= MNARC - 1) {\r
565                                                 validInput = true;\r
566                                         } else {\r
567                                                 System.out.print("NARCS must be <= " + (MNARC - 1));\r
568                                         }\r
569                                 } // while (!validInput)\r
570                         } // if (SYMTY)\r
571                         else { // !SYMTY\r
572                                 validInput = false;\r
573                                 while (!validInput) {\r
574                                         System.out.println("How many arcs are there on the boundary?");\r
575                                         NARCS = input.nextInt();\r
576                                         if (NARCS <= MNARC - 1) {\r
577                                                 validInput = true;\r
578                                         } else {\r
579                                                 System.out.println("NARCS must be <= " + (MNARC - 1));\r
580                                         }\r
581                                 } // while (!validInput)\r
582                         } // else !SYMTY\r
583 \r
584                         GMCO = 0;\r
585                         TXCO = 0;\r
586 \r
587                         for (IA = 1; IA <= NARCS; IA++) {\r
588                                 NUMDER[IA - 1] = false;\r
589                                 validInput = false;\r
590                                 while (!validInput) {\r
591                                         System.out.println("ENTER THE TYPE OF ARC(1-4) for ARC NUMBER " + IA);\r
592                                         TYPE = input.nextInt();\r
593                                         if ((TYPE >= 1) && (TYPE <= 4)) {\r
594                                                 validInput = true;\r
595                                         } else {\r
596                                                 System.out.println("TYPE MUST BE BETWEEN 1 and 4");\r
597                                         }\r
598                                 } // while (!validInput)\r
599                                 if (TYPE == 1) {\r
600                                         ARCTY[IA - 1] = TYPE;\r
601                                         validInput = false;\r
602                                         while (!validInput) {\r
603                                                 try {\r
604                                                         System.out.println("What are the coordinates of the initial point on the line (X,Y)?");\r
605                                                 line = input.next();\r
606                                                 tokens = line.split(",");\r
607                                                 STAPT[IA - 1][0] = Double.valueOf(tokens[0]);\r
608                                                 STAPT[IA - 1][1] = Double.valueOf(tokens[1]);\r
609                                                 validInput = true;\r
610                                                 }\r
611                                                 catch (Exception e) {};\r
612                                         }\r
613                                 } // if (TYPE == 1)\r
614                                 else if (TYPE == 2) {\r
615                                         ARCTY[IA - 1] = TYPE;\r
616                                         validInput = false;\r
617                                         while (!validInput) {\r
618                                                 try {\r
619                                                         System.out.println("What are the coordinates of the initial point on the circle (X,Y)?");\r
620                                                 line = input.next();\r
621                                                 tokens = line.split(",");\r
622                                                 STAPT[IA - 1][0] = Double.valueOf(tokens[0]);\r
623                                                 STAPT[IA - 1][1] = Double.valueOf(tokens[1]);\r
624                                                 validInput = true;\r
625                                                 }\r
626                                                 catch (Exception e) {};\r
627                                         }\r
628                                         validInput = false;\r
629                                         while (!validInput) {\r
630                                                 try {\r
631                                                         System.out.println("What are the coordinates of the center of the circle (X,Y)?");\r
632                                                 line = input.next();\r
633                                                 tokens = line.split(",");\r
634                                                 X = Double.valueOf(tokens[0]);\r
635                                                 Y = Double.valueOf(tokens[1]);\r
636                                                 validInput = true;\r
637                                                 }\r
638                                                 catch (Exception e) {};\r
639                                         }\r
640                                         System.out.println("What is the signed angle subtended at center (in units of PI)?");\r
641                                         ALPHA = input.nextDouble();\r
642                                         GMCO = GMCO + 1;\r
643                                         PGM[IA - 1] = GMCO;\r
644                                         RGM[GMCO - 1] = X;\r
645                                         GMCO = GMCO + 1;\r
646                                         RGM[GMCO - 1] = Y;\r
647                                         GMCO = GMCO + 1;\r
648                                         RGM[GMCO - 1] = ALPHA * PI;\r
649                                 } // else if (TYPE == 2)\r
650                                 else if ((TYPE == 3) || (TYPE == 4)) {\r
651                                         ARCTY[IA - 1] = TYPE;\r
652                                         validInput = false;\r
653                                         while (!validInput) {\r
654                                                 try {\r
655                                                         System.out.println("What are the coordinates of the initial point on the curve (X,Y)?");\r
656                                                 line = input.next();\r
657                                                 tokens = line.split(",");\r
658                                                 STAPT[IA - 1][0] = Double.valueOf(tokens[0]);\r
659                                                 STAPT[IA - 1][1] = Double.valueOf(tokens[1]);\r
660                                                 validInput = true;\r
661                                                 }\r
662                                                 catch (Exception e) {};\r
663                                         }\r
664                                         validInput = false;\r
665                                         while (!validInput) {\r
666                                                 try {\r
667                                                         if (TYPE == 3) {\r
668                                                                 System.out.println("Enter the initial and final parameter values (X,Y)");\r
669                                                         } else {\r
670                                                                 System.out.println("Enter the initial and final polar values (in angles of PI) (X,Y)");\r
671                                                         }\r
672                                                 line = input.next();\r
673                                                 tokens = line.split(",");\r
674                                                 X = Double.valueOf(tokens[0]);\r
675                                                 Y = Double.valueOf(tokens[1]);\r
676                                                 validInput = true;\r
677                                                 }\r
678                                                 catch (Exception e) {};\r
679                                         }\r
680                                         GMCO = GMCO + 1;\r
681                                         PGM[IA - 1] = GMCO;\r
682                                         if (TYPE == 4) {\r
683                                                 RGM[GMCO - 1] = X * PI;\r
684                                                 GMCO = GMCO + 1;\r
685                                                 RGM[GMCO - 1] = Y * PI;\r
686                                         } else {\r
687                                                 RGM[GMCO - 1] = X;\r
688                                                 GMCO = GMCO + 1;\r
689                                                 RGM[GMCO - 1] = Y;\r
690                                         }\r
691                                         for (J = 1; J <= 2; J++) {\r
692                                                 if (J == 1 && TYPE == 3) {\r
693                                                         System.out.println("ENTER JAVA EXPRESSION WITH NO SPACES ENDING IN // FOR PARFUN");\r
694                                                         System.out.println("PUT REAL PART ui IMAGINARY PART");\r
695                                                 } else if (J == 2 && TYPE == 3) {\r
696                                                         System.out.println("ENTER JAVA EXPRESSION WITH NO SPACES ENDING IN // FOR DPARFN");\r
697                                                         System.out.println("PUT REAL PART ui IMAGINARY PART");\r
698                                                 } else if (J == 1 && TYPE == 4) {\r
699                                                         System.out.println("ENTER JAVA EXPRESSION WITH NO SPACES ENDING IN // FOR RADIUS");\r
700                                                 } else {\r
701                                                         System.out.println("ENTER JAVA EXPRESSION WITH NO SPACES ENDING IN // FOR RADIUS DERIVATIVE");\r
702                                                 }\r
703 \r
704                                                 TXCO = TXCO + 1;\r
705                                                 PTX[IA - 1 + (J - 1) * MNARC] = TXCO;\r
706                                                 I = 1;\r
707 \r
708                                                 TXT = input.next();\r
709                                                 L = -1;\r
710                                                 while (L == -1) {\r
711                                                         L = TXT.indexOf("//");\r
712                                                         if (L == -1) {\r
713                                                                 DEFN[TXCO - 1] = TABC + TXT;\r
714                                                                 I = I + 1;\r
715                                                                 TXCO = TXCO + 1;\r
716                                                         } // if (L == -1)\r
717                                                 } // while (L == -1)\r
718                                                 NTX[IA - 1 + (J - 1) * MNARC] = I;\r
719                                                 if (L == 0) {\r
720                                                         DEFN[TXCO - 1] = TABC;\r
721                                                         NUMDER[IA - 1] = true;\r
722                                                 } else {\r
723                                                         DEFN[TXCO - 1] = TABC + TXT.substring(0, L);\r
724                                                 }\r
725                                                 if ((J == 1) && (TYPE == 4)) {\r
726                                                         System.out.println("(... = ZRAD)");\r
727                                                 }\r
728                                         } // for (J = 1; J <= 2; J++)\r
729                                 } // else if ((TYPE == 3) || (TYPE == 4))\r
730                         } // for (IA = 1; IA <= NARCS; IA++)\r
731 \r
732                         if (SYMTY) {\r
733                                 validInput = false;\r
734                                 while (!validInput) {\r
735                                         try {\r
736                                                 System.out.println("ENTER THE COORDINATES OF FINAL POINT ON THIS LAST ARC (X,Y)");\r
737                                         line = input.next();\r
738                                         tokens = line.split(",");\r
739                                         STAPT[NARCS][0] = Double.valueOf(tokens[0]);\r
740                                         STAPT[NARCS][1] = Double.valueOf(tokens[1]);\r
741                                         validInput = true;\r
742                                         }\r
743                                         catch (Exception e) {};\r
744                                 }\r
745                         } else {\r
746                                 STAPT[NARCS][0] = STAPT[0][0];\r
747                                 STAPT[NARCS][1] = STAPT[0][1];\r
748                         }\r
749 \r
750                         validInput = false;\r
751                         while (!validInput) {\r
752                                 System.out.println("END OF INPUT PHASE; CONTINUE WITH PROCESSING (Y/N)?");\r
753                                 String term = input.next();\r
754                                 String firstTerm = term.substring(0, 1);\r
755                                 if (firstTerm.equalsIgnoreCase("Y")) {\r
756                                         validInput = true;\r
757                                 } else if (firstTerm.equalsIgnoreCase("N")) {\r
758                                         validInput = true;\r
759                                         setCompleted(false);\r
760                                         try {\r
761                                                 raFile.close();\r
762                                         } catch (IOException e) {\r
763 \r
764                                         }\r
765                                         return;\r
766                                 } else {\r
767                                         System.out.println(term + " is not a valid response");\r
768                                 }\r
769                         } // while (!validInput)\r
770                 } // if (traditionalInput)\r
771                 HEADER("PARFUN", REDD, raFile);\r
772                 if (SYMTY) {\r
773                         SYINF1(ORDRG, ORDSG, RTUNI, U2, REFLN, CENSY, STAPT[0], STAPT[NARCS], IER);\r
774                         if (IER[0] > 0) {\r
775                                 WRTAIL(6, 0, IER[0], null);\r
776                                 return;\r
777                         }\r
778                         System.out.println("\nN O T E : THE ORDER OF THE SYMMETRY GROUP IS " + ORDSG);\r
779                         if (REFLN) {\r
780                                 System.out.println("          ISYGP = " + (-ORDSG[0]));\r
781                         } else {\r
782                                 System.out.println("          ISYGP = " + (ORDSG[0]));\r
783                         }\r
784                         WRSYM1(NARCS, ORDRG[0], ORDSG[0], RTUNI, U2, CENSY, REFLN, true, REDD, CHNL, raFile);\r
785                         if (REFLN) {\r
786                                 CH = "TS";\r
787                         } else {\r
788                                 CH = "TT";\r
789                         }\r
790                         WRFUN1(NARCS, STAPT, ARCTY, PGM, RGM, PTX, NTX, DEFN, CHNL, "IB", CH, "ZETA  ", REDD, raFile);\r
791                         WRSYM2(NARCS, ORDRG[0], CENSY, REFLN, CHNL, raFile);\r
792                 } else {\r
793                         WRFUN1(NARCS, STAPT, ARCTY, PGM, RGM, PTX, NTX, DEFN, CHNL, "IA", "TT", "PARFUN", REDD, raFile);\r
794                 }\r
795 \r
796                 try {\r
797                         raFile.writeBytes("//\n");\r
798                         raFile.writeBytes("}\n");\r
799 \r
800                         // **** WRITE THE SOURCE CODE FOR DPARFN\r
801 \r
802                         raFile.writeBytes("//...........................................\n");\r
803                 } catch (IOException e) {\r
804                         MipavUtil.displayError("IOException " + e + " in PARGEN");\r
805                         System.exit(-1);\r
806                 }\r
807                 HEADER("DPARFN", REDD, raFile);\r
808                 if (SYMTY) {\r
809                         WRSYM1(NARCS, ORDRG[0], ORDSG[0], RTUNI, U2, CENSY, REFLN, false, REDD, CHNL, raFile);\r
810                         if (REFLN) {\r
811                                 CH = "TS";\r
812                         } else {\r
813                                 CH = "TT";\r
814                         }\r
815                         WRFUN2(NARCS, MNARC, STAPT, ARCTY, PGM, RGM, PTX, NTX, DEFN, CHNL, "IB", CH, "ZETA  ", NUMDER, REDD,\r
816                                         raFile);\r
817                         WRSYM3(NARCS, ORDRG[0], REFLN, CHNL, raFile);\r
818                 } // if (SYMTY)\r
819                 else {\r
820                         WRFUN2(NARCS, MNARC, STAPT, ARCTY, PGM, RGM, PTX, NTX, DEFN, CHNL, "IA", "TT", "DPARFN", NUMDER, REDD,\r
821                                         raFile);\r
822                 }\r
823 \r
824                 try {\r
825                         raFile.writeBytes("//\n");\r
826                         raFile.writeBytes("}\n");\r
827                         raFile.close();\r
828                 } catch (IOException e) {\r
829                         MipavUtil.displayError("IOException " + e + " in PARGEN");\r
830                         System.exit(-1);\r
831                 }\r
832                 WRTAIL(6, 0, IER[0], null);\r
833 \r
834         } // public void PARGEN\r
835 \r
836         private void HEADER(String TXT, String REDD, RandomAccessFile raFile) {\r
837 \r
838                 String TAB6 = "      ";\r
839 \r
840                 String LINE = TAB6 + "private void " + TXT + "(int IA, double TT[]) {\n";\r
841                 try {\r
842                         raFile.writeBytes(LINE);\r
843 \r
844                         LINE = "//" + TAB6 + "IMPLICIT REAL(A-H,O-S),INTEGER(I-N),COMPLEX(T-Z)\n";\r
845                         raFile.writeBytes(LINE);\r
846 \r
847                         raFile.writeBytes("      double PI = " + Math.PI + ";\n");\r
848                         raFile.writeBytes("      double UI[] = new double[]{0.0,1.0};\n");\r
849                         raFile.writeBytes("//\n");\r
850                 } catch (IOException e) {\r
851                         MipavUtil.displayError("IOException " + e + " in HEADER");\r
852                         System.exit(-1);\r
853                 }\r
854 \r
855         } // private void HEADER\r
856 \r
857         private void SYINF1(int ORDRG[], int ORDSG[], double RTUNI[], double U2[], boolean REFLN, double Z0[], double Z1[],\r
858                         double Z2[], int IER[]) {\r
859                 // COMPLEX RTUNI,U2,Z0,Z1,Z2\r
860 \r
861                 // **** GIVEN Z0,THE CENTRE OF SYMMETRY, Z1 AND Z2, THE INITIAL AND\r
862                 // FINAL\r
863                 // **** POINTS ON THE FUNDAMENTAL BOUNDARY SECTION, REFLN, WHICH IS TRUE\r
864                 // **** IF THE SYMMETRY GROUP HAS IMPROPER ROTATIONAL ELEMENTS\r
865                 // **** (I.E. REFLECTIONAL SYMMETRIES), THIS ROUTINE COMPUTES\r
866                 // **** ORDRG - THE ORDER OF THE SUBGROUP OF PROPER ROTATIONS (THIS IS\r
867                 // THE\r
868                 // **** ORDER OF THE SYMMETRY GROUP IF REFLN=.FALSE.)\r
869                 // **** ORDSG - THE ORDER OF THE FULL SYMMETRY GROUP, EITHER ORDRG OR\r
870                 // **** 2*ORDRG DEPENDING ON WHETHER REFLN IS .FALSE. OR .TRUE.\r
871                 // **** RTUNI - THE ROOT OF UNITY FROM WHICH THE PROPER ROTATIONAL\r
872                 // SUBROUP\r
873                 // **** IS GENERATED\r
874                 // **** U2 - THE ADDITIONAL IN-PLANE ROTATION WHICH, WHEN COMBINED WITH\r
875                 // **** CONJUGATION, DEFINES THE IMPROPER ROTATION FOR THE CASE\r
876                 // **** REFLN=.TRUE.\r
877 \r
878                 // LOCAL VARIABLES\r
879 \r
880                 double ALPHA, PI, SQRTEPS;\r
881                 // COMPLEX CT,U\r
882                 double CT[] = new double[2];\r
883                 double U[] = new double[2];\r
884                 double cr[] = new double[1];\r
885                 double ci[] = new double[1];\r
886 \r
887                 PI = Math.PI;\r
888                 SQRTEPS = Math.sqrt(EPS);\r
889                 CT[0] = Z2[0] - Z0[0];\r
890                 CT[1] = Z2[1] - Z0[1];\r
891                 double ABSCT = zabs(CT[0], CT[1]);\r
892                 if (ABSCT < SQRTEPS) {\r
893                         IER[0] = 56;\r
894                         return;\r
895                 }\r
896                 U[0] = CT[0] / ABSCT;\r
897                 U[1] = CT[1] / ABSCT;\r
898                 zmlt(U[0], U[1], U[0], U[1], cr, ci);\r
899                 U2[0] = cr[0];\r
900                 U2[1] = ci[0];\r
901 \r
902                 zmlt(Z1[0] - Z0[0], Z1[1] - Z0[1], U[0], -U[1], cr, ci);\r
903                 CT[0] = cr[0];\r
904                 CT[1] = ci[0];\r
905                 ABSCT = zabs(CT[0], CT[1]);\r
906                 if (ABSCT < SQRTEPS) {\r
907                         IER[0] = 57;\r
908                         return;\r
909                 }\r
910                 ALPHA = Math.atan2(CT[1], CT[0]);\r
911                 ALPHA = Math.abs(ALPHA);\r
912 \r
913                 if (REFLN) {\r
914                         ORDRG[0] = (int) Math.round(PI / ALPHA);\r
915                         ORDSG[0] = 2 * ORDRG[0];\r
916                 } else {\r
917                         ORDRG[0] = 2 * (int) Math.round(PI / ALPHA);\r
918                         ORDSG[0] = ORDRG[0];\r
919                 }\r
920 \r
921                 ALPHA = 2.0 * PI / (double) (ORDRG[0]);\r
922                 RTUNI[0] = Math.cos(ALPHA);\r
923                 RTUNI[1] = Math.sin(ALPHA);\r
924 \r
925                 // NORMAL EXIT\r
926 \r
927                 IER[0] = 0;\r
928                 return;\r
929         } // private void SYINF1\r
930 \r
931         // COMPLEX FUNCTION PARFUN(I,T)\r
932         // INTEGER I\r
933         // COMPLEX T\r
934         double[] PARFUN(int I, double T[]) {\r
935 \r
936                 // DUMMY FUNCTION TO AID LINK-LOADING OF PARGEN\r
937 \r
938                 double result[] = new double[] { 1.0, 0.0 };\r
939                 return result;\r
940         } // double[] PARFUN\r
941 \r
942         // COMPLEX FUNCTION DPARFN(I,T)\r
943         // INTEGER I\r
944         // COMPLEX T\r
945         double[] DPARFN(int I, double T[]) {\r
946 \r
947                 // DUMMY FUNCTION TO AID LINK-LOADING OF PARGEN\r
948 \r
949                 double result[] = new double[] { 1.0, 0.0 };\r
950                 return result;\r
951         } // double[] DPARFN\r
952 \r
953         private void WRHEAD(int I, int CHNL, RandomAccessFile raFile) {\r
954 \r
955                 // **** WRITE A HEADING FOR THE MAIN CONFPACK MODULES JAPHYC (I=1),\r
956                 // **** GQPHYC (I=2), JACANP (I=3), GQCANP (I=4), CNDPLT (I=5), THE\r
957                 // **** PARAMETRIC FUNCTION GENERATOR PARGEN (I=6),THE PARAMETRIC\r
958                 // FUNCTION\r
959                 // **** TESTER TSTPLT (I=7) AND THE LEVEL CURVE ROUTINE LEVCUR (I=8). IF\r
960                 // **** CHNL=0 THEN WRITE ON THE STANDARD OUTPUT CHANNEL, OTHERWISE\r
961                 // WRITE\r
962                 // **** ON THE CHANNEL SPECIFIED BY CHNL.\r
963                 //\r
964                 // LOCAL VARIABLES\r
965                 //\r
966                 String DOTS = ".................................................";\r
967                 String CPHEAD = ": C O N F P A C K    M O D U L E    ";\r
968                 String MOD[] = new String[] { "J A P H Y C :", "G Q P H Y C :", "J A C A N P :", "G Q C A N P :",\r
969                                 "C N D P L T :", "P A R G E N :", "T S T P L T :", "L E V C U R :" };\r
970                 String TXT = CPHEAD + MOD[I - 1];\r
971 \r
972                 if (CHNL == 0 || CHNL == 6) {\r
973                         System.out.println("\n\n      " + DOTS + "\n      " + TXT + "\n      " + DOTS);\r
974                 } else {\r
975                         try {\r
976                                 raFile.writeBytes("\n\n      //" + DOTS + "\n      //" + TXT + "\n      //" + DOTS + "\n");\r
977                         } catch (IOException e) {\r
978                                 MipavUtil.displayError("IOException " + e + " on raFile.writeBytes in WRHEAD");\r
979                                 System.exit(-1);\r
980                         }\r
981                 }\r
982                 return;\r
983         } // private void WRHEAD\r
984 \r
985         private void WRTAIL(int I, int CHNL, int IER, RandomAccessFile raFile) {\r
986 \r
987                 //\r
988                 // **** WRITE A CLOSING MESSAGE FOR THE MAIN CONFPACK MODULES JAPHYC\r
989                 // (I=1)\r
990                 // **** GQPHYC (I=2), JACANP (I=3), GQCANP (I=4), CNDPLT (I=5), THE\r
991                 // PARA-\r
992                 // **** METRIC FUNCTION GENERATOR PARGEN (I=6), THE PARAMETRIC FUNCTION\r
993                 // **** TESTER TSTPLT (I=7) AND THE LEVEL CURVE ROUTINE LEVCUR (I=8). IF\r
994                 // **** CHNL=0 THEN WRITE ON THE STANDARD OUTPUT CHANNEL, OTHERWISE\r
995                 // WRITE\r
996                 // **** ON THE CHANNEL SPECIFIED BY CHNL. THE TEXT OF THE MESSAGE IS\r
997                 // **** DETERMINED BY THE ERROR NUMBER IER VIA THE SUBROUTINE IERTXT.\r
998 \r
999                 // LOCAL VARIABLES\r
1000 \r
1001                 String MOD[] = new String[] { "J A P H Y C :", "G Q P H Y C :", "J A C A N P :", "G Q C A N P :",\r
1002                                 "C N D P L T :", "P A R G E N :", "T S T P L T :", "L E V C U R :" };\r
1003                 String GOOD = "  NORMAL EXIT";\r
1004                 String BAD = "  ABNORMAL EXIT";\r
1005                 String LINE = "__________________________________________________________________&";\r
1006 \r
1007                 String TXT, TXT2;\r
1008                 if (IER == 0) {\r
1009                         TXT = MOD[I - 1] + GOOD;\r
1010                 } else {\r
1011                         TXT = MOD[I - 1] + BAD;\r
1012                 }\r
1013                 TXT2 = IERTXT(IER);\r
1014 \r
1015                 if ((CHNL == 0) || (CHNL == 6)) {\r
1016                         System.out.println("\n\n      " + TXT);\r
1017                         System.out.println("      " + TXT2);\r
1018                         System.out.println(LINE);\r
1019                 } else {\r
1020                         try {\r
1021                                 raFile.writeBytes("\n\n      //" + TXT + "\n");\r
1022                                 raFile.writeBytes("      //" + TXT2 + "\n");\r
1023                                 raFile.writeBytes("//" + LINE + "\n");\r
1024                         } catch (IOException e) {\r
1025                                 MipavUtil.displayError("IOException " + e + " in WRTAIL");\r
1026                                 System.exit(-1);\r
1027                         }\r
1028                 }\r
1029                 return;\r
1030         }\r
1031 \r
1032         private String IERTXT(int IER) {\r
1033 \r
1034                 // **** SUPPLY ERROR MESSAGE TEXT FOR ERROR NUMBER IER\r
1035                 String result = null;\r
1036                 if (IER == 0) {\r
1037                         result = " ";\r
1038                 } else if (IER == 1) {\r
1039                         result = "PARAMETER IBNDS[0] IS TOO SMALL AT START OF JAPHYC";\r
1040                 } else if (IER == 2) {\r
1041                         result = "PARAMETER IBNDS[1] IS TOO SMALL AT START OF JAPHYC";\r
1042                 } else if (IER == 3) {\r
1043                         result = "NQPTS < 1 AT START OF JAPHYC";\r
1044                 } else if (IER == 4) {\r
1045                         result = "FAILURE TO CONVERGE IN EIGSYS; CAN''T SET UP BASIC QUADRATURE RULES";\r
1046                 } else if (IER == 5) {\r
1047                         result = "PARAMETER MNQPT IN IGNLVL MUST BE INCREASED TO AT LEAST NQPTS";\r
1048                 } else if (IER == 6) {\r
1049                         result = "FAILURE TO CONVERGE IN IMTQLH; CAN''T SET UP IGNORE LEVELS";\r
1050                 } else if (IER == 7) {\r
1051                         result = "FAILURE TO CONVERGE IN IMTQLH; CAN''T SET UP COLLOCATION POINTS";\r
1052                 } else if (IER == 8) {\r
1053                         result = "ARGUMENT MNEQN IS TOO SMALL AT START OF JAPHYC";\r
1054                 } else if (IER == 9) {\r
1055                         result = "PARAMETER IBNDS[3] IS TOO SMALL AT START OF JAPHYC";\r
1056                 } else if (IER == 10) {\r
1057                         result = "PARAMETER NMAX IN SUBIN7 MUST BE INCREASED TO AT LEAST 2*NQPTS";\r
1058                 } else if (IER == 11) {\r
1059                         result = "PARAMETER IBNDS[2] IS TOO SMALL AT START OF JAPHYC";\r
1060                 } else if (IER == 12) {\r
1061                         result = "PARAMETER NC IN DEJAC7 AND DELEG7 MUST BE INCREASED";\r
1062                 } else if (IER == 13) {\r
1063                         result = "PARAMETER NR IN DEJAC7 AND DELEG7 MUST BE >= (NQPTS -1)";\r
1064                 } else if (IER == 14) {\r
1065                         result = "A CORNER ANGLE IS TOO SMALL; MAY CAUSE OVERFLOW IN GAMMA FUNCTION";\r
1066                 } else if (IER == 15) {\r
1067                         result = "SINGULAR COLLOCATION MATRIX";\r
1068                 } else if (IER == 16) {\r
1069                         result = "COLLOCATION MATRIX IS EFFECTIVELY SINGULAR";\r
1070                 } else if (IER == 17) {\r
1071                         result = "NUMBER OF SUBARCS EXCEEDS IBNDS[0] DURING REFINEMENT";\r
1072                 } else if (IER == 18) {\r
1073                         result = "NUMBER OF EQUATIONS EXCEEDS MNEQN DURING REFINEMENT";\r
1074                 } else if (IER == 19) {\r
1075                         result = "TOTAL NUMBER OF QUADRATURE PTS EXCEEDS IBNDS[3] DURING REFINEMENT";\r
1076                 } else if (IER == 20) {\r
1077                         result = "NUMBER OF QUADRATURE PANELS EXCEEDS IBNDS[2] DURING REFINEMENT";\r
1078                 } else if (IER == 21) {\r
1079                         result = "FAILURE TO CONVERGE IN IMTQLH; CAN''T SET UP TEST POINTS";\r
1080                 } else if (IER == 22) {\r
1081                         result = "ARGUMENT MQUPH OF GQPHYC MUST BE INCREASED";\r
1082                 } else if (IER == 23) {\r
1083                         result = "PARAMETER MNCOF IN POPQF1 MUST BE >= NQPTS";\r
1084                 } else if (IER == 24) {\r
1085                         result = "NUMBER OF QUADRATURE PANELS EXCEEDS MQIN1 IN GQPHYC";\r
1086                 } else if (IER == 25) {\r
1087                         result = "PARAMETER MNXI IN DEPPJ8 AND DEPPL8 MUST BE INCREASED";\r
1088                 } else if (IER == 26) {\r
1089                         result = "PARAMETER MAXNZ IN DEPPJ9 AND DEPPL9 MUST BE INCREASED";\r
1090                 } else if (IER == 27) {\r
1091                         result = "PARAMETER MXNQD IN PHTCA1 MUST BE INCREASED";\r
1092                 } else if (IER == 28) {\r
1093                         result = "PARAMETER MXCOF IN PHTCA1 MUST BE INCREASED";\r
1094                 } else if (IER == 29) {\r
1095                         result = "PARAMETER MQIN1 IN PHTCA1 MUST BE INCREASED";\r
1096                 } else if (IER == 30) {\r
1097                         result = "PARAMETER MNDG IN JCFIM5 MUST BE INCREASED";\r
1098                 } else if (IER == 31) {\r
1099                         result = "PARAMETER MNQD IN JCFIM5 MUST BE INCREASED";\r
1100                 } else if (IER == 32) {\r
1101                         result = "ARGUMENT IBNDS[1] SUPPLIED TO JACANP MUST BE INCREASED";\r
1102                 } else if (IER == 33) {\r
1103                         result = "ARGUMENT IBNDS[0] SUPPLIED TO JACANP MUST BE INCREASED";\r
1104                 } else if (IER == 34) {\r
1105                         result = "FN HAS SAME SIGN AT INTERVAL ENDS IN BISNEW; CAN''T SOLVE BCF EQN";\r
1106                 } else if (IER == 35) {\r
1107                         result = "DERIVATIVE OF BCF IS ZERO IN BISNEW; CAN''T SOLVE BCF EQN";\r
1108                 } else if (IER == 36) {\r
1109                         result = "ELEMENT OF ARGUMENT ARRAY SVAL IN RHOFN IS +-1; CAN''T CONTINUE";\r
1110                 } else if (IER == 37) {\r
1111                         result = "PARAMETER MXNQD IN CINRAD MUST BE INCREASED";\r
1112                 } else if (IER == 38) {\r
1113                         result = "PARAMETER MXCOF IN CINRAD MUST BE INCREASED";\r
1114                 } else if (IER == 39) {\r
1115                         result = "CENTRE POINT IS PATHOLOGICALLY CLOSE TO BOUNDARY;CAN''T CONTINUE";\r
1116                 } else if (IER == 40) {\r
1117                         result = "PARAMETER MQIN1 IN CINRAD MUST BE INCREASED";\r
1118                 } else if (IER == 41) {\r
1119                         result = "ARGUMENT MQUCA OF GQCANP MUST BE INCREASED";\r
1120                 } else if (IER == 42) {\r
1121                         result = "PARAMETER MNCOF IN POPQG1 MUST BE >= NQPTS";\r
1122                 } else if (IER == 43) {\r
1123                         result = "NUMBER OF QUADRATURE PANELS EXCEEDS MQIN1 IN GQCANP";\r
1124                 } else if (IER == 44) {\r
1125                         result = "PARAMETER MNCOF IN BMPHC1 MUST BE >= NQPTS";\r
1126                 } else if (IER == 45) {\r
1127                         result = "ARGUMENTS IARC, PHYPT OF BMPHYC DON''T DEFINE A BOUNDARY POINT";\r
1128                 } else if (IER == 46) {\r
1129                         result = "PARAMETER MNCOF IN BMCAP1 MUST BE >= NQPTS";\r
1130                 } else if (IER == 47) {\r
1131                         result = "PARAMETER MXNQD IN CATPH4 MUST BE INCREASED";\r
1132                 } else if (IER == 48) {\r
1133                         result = "PARAMETER MNCOF IN CATPH4 MUST BE >= NQPTS";\r
1134                 } else if (IER == 49) {\r
1135                         result = "PARAMETER MQIN1 IN CATPH4 MUST BE INCREASED";\r
1136                 } else if (IER == 50) {\r
1137                         result = "PARAMETER MXCOF IN DIAGN3 MUST BE >= NQPTS";\r
1138                 } else if (IER == 51) {\r
1139                         result = "NON-ANALYTIC ARC DETECTED IN DIAGN3";\r
1140                 } else if (IER == 52) {\r
1141                         result = "PARAMETER MAXSA IN CNDPLT MUST BE INCREASED";\r
1142                 } else if (IER == 53) {\r
1143                         result = "OVERFLOW EXPECTED IN IGNLVL; A CORNER ANGLE IS TOO SMALL";\r
1144                 } else if (IER == 54) {\r
1145                         result = "PARAMETER MXCO IN AXION1 MUST BE INCREASED";\r
1146                 } else if (IER == 55) {\r
1147                         result = "NARCS ISN''T AN INTEGER MULTIPLE OF THE ORDER OF THE SYMMETRY GROUP";\r
1148                 } else if (IER == 56) {\r
1149                         result = "CENTRE OF SYMMETRY IS PATHOLOGICALLY CLOSE TO LAST POINT ON FBS";\r
1150                 } else if (IER == 57) {\r
1151                         result = "CENTRE OF SYMMETRY IS PATHOLOGICALLY CLOSE TO FIRST POINT ON FBS";\r
1152                 } else if (IER == 58) {\r
1153                         result = "NUMBER OF ARCS IS TOO BIG; INCREASE PARAMETER MNARC IN PARGEN";\r
1154                 } else if (IER == 59) {\r
1155                         result = "NUMBER OF ARCS IS TOO BIG; INCREASE PARAMETER MNARC IN TSTPLT";\r
1156                 } else if (IER == 60) {\r
1157                         result = "NON-ANALYTIC ARC (DPARFN=(0.,0.)) DETECTED IN TSTPLT";\r
1158                 } else {\r
1159                         result = "UNRECOGNIZED ERROR NUMBER IN IERTXT ROUTINE !!";\r
1160                 }\r
1161                 return result;\r
1162         }\r
1163 \r
1164         private void WRSYM1(int NARCS, int ORDRG, int ORDSG, double[] RTUNI, double[] U2, double[] CENSY, boolean REFLN,\r
1165                         boolean PARFUN, String REDD, int CHNL, RandomAccessFile raFile) {\r
1166                 // COMPLEX RTUNI,U2,CENSY\r
1167 \r
1168                 // **** TO WRITE THE DIMENSION AND PARAMETER STATEMENTS AND THE CODE TO\r
1169                 // **** TO REDUCE A GIVEN ARC NUMBER TO ITS SYMMETRIC COUNTERPART ON THE\r
1170                 // **** FUNDAMENTAL BOUNDARY SECTION.\r
1171 \r
1172                 // .......................................................................\r
1173                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1174                 // LAST UPDATE: 4 AUG 1990\r
1175                 // .......................................................................C\r
1176                 // LOCAL VARIABLES\r
1177 \r
1178                 int I;\r
1179                 double R, A;\r
1180                 // COMPLEX ZT\r
1181                 double ZT[] = new double[2];\r
1182                 boolean NEEDC;\r
1183                 // String FMT;\r
1184                 double cr[] = new double[1];\r
1185                 double ci[] = new double[1];\r
1186 \r
1187                 // FMT="(A12,"+REDD+",A1,"+REDD+",A2)";\r
1188 \r
1189                 try {\r
1190                         if (PARFUN) {\r
1191                                 NEEDC = ((CENSY[0] != 0.0) || (CENSY[1] != 0.0));\r
1192                                 if (NEEDC || REFLN) {\r
1193                                         raFile.writeBytes("      PARAMETER (\n");\r
1194                                         if (NEEDC && REFLN) {\r
1195                                                 R = U2[0];\r
1196                                                 A = U2[1];\r
1197                                                 raFile.writeBytes("U2[0] = " + R + ";\n");\r
1198                                                 raFile.writeBytes("U2[1] = " + A + ";\n");\r
1199                                                 R = CENSY[0];\r
1200                                                 A = CENSY[1];\r
1201                                                 raFile.writeBytes("ZCEN[0] = " + R + ";\n");\r
1202                                                 raFile.writeBytes("ZCEN[1] = " + A + ";)\n");\r
1203                                         } // if (NEEDC && REFLN)\r
1204                                         else if (NEEDC && (!REFLN)) {\r
1205                                                 R = CENSY[0];\r
1206                                                 A = CENSY[1];\r
1207                                                 raFile.writeBytes("ZCEN[0] = " + R + ";\n");\r
1208                                                 raFile.writeBytes("ZCEN[1] = " + A + ";)\n");\r
1209                                         } // else if (NEEDC && (!REFLN))\r
1210                                         else {\r
1211                                                 R = U2[0];\r
1212                                                 A = U2[1];\r
1213                                                 raFile.writeBytes("U2[0] = " + R + ";\n");\r
1214                                                 raFile.writeBytes("U2[1] = " + A + ";)\n");\r
1215                                         } // else\r
1216                                         raFile.writeBytes("//\n");\r
1217                                 } // if (NEEDC || REFLN)\r
1218                         } else if (REFLN) {\r
1219                                 R = U2[0];\r
1220                                 A = U2[1];\r
1221                                 raFile.writeBytes("      PARAMETER (\n");\r
1222                                 raFile.writeBytes("U2[0] = " + R + ";\n");\r
1223                                 raFile.writeBytes("U2[1] = " + A + ";)\n");\r
1224                                 raFile.writeBytes("//\n");\r
1225                         }\r
1226 \r
1227                         // FMT="(A7,"+REDD+",A1,"+REDD+",A2)";\r
1228 \r
1229                         if (ORDRG >= 2) {\r
1230                                 raFile.writeBytes("double WW[] = new double[" + (ORDRG - 1) + "];\n");\r
1231                                 ZT[0] = 1.0;\r
1232                                 ZT[1] = 0.0;\r
1233                                 for (I = 0; I < ORDRG - 2; I++) {\r
1234                                         zmlt(ZT[0], ZT[1], RTUNI[0], RTUNI[1], cr, ci);\r
1235                                         ZT[0] = cr[0];\r
1236                                         ZT[1] = ci[0];\r
1237                                         raFile.writeBytes("WW[" + I + "][0] = " + ZT[0] + ";\n");\r
1238                                         raFile.writeBytes("WW[" + I + "][1] = " + ZT[1] + ";\n");\r
1239                                 }\r
1240                                 zmlt(ZT[0], ZT[1], RTUNI[0], RTUNI[1], cr, ci);\r
1241                                 ZT[0] = cr[0];\r
1242                                 ZT[1] = ci[0];\r
1243                                 raFile.writeBytes("WW[" + I + "][0] = " + ZT[0] + ";\n");\r
1244                                 raFile.writeBytes("WW[" + I + "][1] = " + ZT[1] + ";)\n");\r
1245                                 raFile.writeBytes("//\n");\r
1246                         } // if (ORDRG >= 2)\r
1247 \r
1248                         if (ORDRG > 19) {\r
1249                                 System.out.println("\n");\r
1250                                 System.out.println("             ****WARNING****");\r
1251                                 System.out.println("MORE THAN 19 CONTINUTATION LINES HAVE BEEN WRITTEN");\r
1252                         }\r
1253 \r
1254                         if (REFLN) {\r
1255                                 if (ORDRG > 1) {\r
1256                                         if (NARCS > 1) {\r
1257                                                 I = 2 * NARCS;\r
1258                                                 raFile.writeBytes("IB = IA%" + I + ";\n");\r
1259                                                 raFile.writeBytes("if (IB == 0) IB = " + I + ";\n");\r
1260                                                 I = I + 1;\r
1261                                                 raFile.writeBytes("if (IB > " + NARCS + ") {\n");\r
1262                                                 raFile.writeBytes("    IB = " + I + " - IB;\n");\r
1263                                                 raFile.writeBytes("    TS[0] = -TT[0];\n");\r
1264                                                 raFile.writeBytes("    TS[1] = TT[1]);\n");\r
1265                                                 raFile.writeBytes("}\n");\r
1266                                                 raFile.writeBytes("else {\n");\r
1267                                                 raFile.writeBytes("    TS[0] = TT[0];\n");\r
1268                                                 raFile.writeBytes("    TS[1] = TT[1];\n");\r
1269                                                 raFile.writeBytes("}\n");\r
1270                                         } // if (NARCS > 1)\r
1271                                         else {\r
1272                                                 raFile.writeBytes("if ((IA%2) == 0) {\n");\r
1273                                                 raFile.writeBytes("    TS[0] = -TT[0];\n");\r
1274                                                 raFile.writeBytes("    TS[1] = TT[1];\n");\r
1275                                                 raFile.writeBytes("}\n");\r
1276                                                 raFile.writeBytes("else {\n");\r
1277                                                 raFile.writeBytes("    TS[0] = TT[0];\n");\r
1278                                                 raFile.writeBytes("    TS[1] = TT[1];\n");\r
1279                                                 raFile.writeBytes("}\n");\r
1280                                         } // else\r
1281                                 } // if (ORDRG > 1)\r
1282                                 else {\r
1283                                         if (NARCS > 1) {\r
1284                                                 I = 2 * NARCS + 1;\r
1285                                                 raFile.writeBytes("if (IA > " + NARCS + "){\n");\r
1286                                                 raFile.writeBytes("    IB = " + I + " -IA;\n");\r
1287                                                 raFile.writeBytes("    TS[0] = -TT[0];\n");\r
1288                                                 raFile.writeBytes("    TS[1] = TT[1]);\n");\r
1289                                                 raFile.writeBytes("}\n");\r
1290                                                 raFile.writeBytes("else {\n");\r
1291                                                 raFile.writeBytes("    IB = IA;\n");\r
1292                                                 raFile.writeBytes("    TS[0] = TT[0];\n");\r
1293                                                 raFile.writeBytes("    TS[1] = TT[1];\n");\r
1294                                                 raFile.writeBytes("}\n");\r
1295                                         } // if (NARCS)\r
1296                                         else {\r
1297                                                 raFile.writeBytes("if (IA == 2) {\n");\r
1298                                                 raFile.writeBytes("    TS[0] = -TT[0];\n");\r
1299                                                 raFile.writeBytes("    TS[1] = TT[1]);\n");\r
1300                                                 raFile.writeBytes("}\n");\r
1301                                                 raFile.writeBytes("else {\n");\r
1302                                                 raFile.writeBytes("    TS[0] = TT[0];\n");\r
1303                                                 raFile.writeBytes("    TS[1] = TT[1];\n");\r
1304                                                 raFile.writeBytes("}\n");\r
1305                                         } // else\r
1306                                 } // else\r
1307                         } // if (REFLN)\r
1308                         else if (NARCS > 1) {\r
1309                                 raFile.writeBytes("IB = IA%" + NARCS + ";\n");\r
1310                                 raFile.writeBytes("if (IB == 0) IB = " + NARCS + ";\n");\r
1311                         } // else if (NARCS > 1)\r
1312 \r
1313                         raFile.writeBytes("//\n");\r
1314                 } // try\r
1315                 catch (IOException e) {\r
1316                         MipavUtil.displayError("IOException " + e + " in WRSYM1");\r
1317                         System.exit(-1);\r
1318                 }\r
1319 \r
1320         }\r
1321 \r
1322         private void WRFUN1(int NARCS, double STAPT[][], int ARCTY[], int PGM[], double RGM[], int PTX[], int NTX[],\r
1323                         String DEFN[], int CHNL, String CHIA, String CHTT, String VAR, String REDD, RandomAccessFile raFile) {\r
1324                 // COMPLEX STAPT(*)\r
1325                 // CHARACTER DEFN(*)*72,CHIA*2,CHTT*2,VAR*6,REDD*6\r
1326 \r
1327                 // **** TO WRITE THE SOURCE CODE FOR PARFUN IN THE CASE WHERE NO\r
1328                 // **** SYMMETRY IS INVOLVED.\r
1329 \r
1330                 // .......................................................................\r
1331                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1332                 // LAST UPDATE: 4 AUG 1990\r
1333                 // .......................................................................\r
1334 \r
1335                 // LOCAL VARIABLES\r
1336 \r
1337                 int IA, I, J, K;\r
1338                 // CHARACTER TX1*16,TX2*21,FMT1*11,FMT2*11\r
1339                 String TX1, TX2;\r
1340                 // String FMT1,FMT2;\r
1341                 // EXTERNAL PTFUN1\r
1342                 TX1 = "     if(" + CHIA + " == ";\r
1343                 TX2 = "     else if(" + CHIA + " == ";\r
1344                 // FMT1="(A16,I3,A6)";\r
1345                 // FMT2="(A21,I3,A6)";\r
1346                 double STAPT2[][];\r
1347                 double RGM2[];\r
1348                 String DEFN2[];\r
1349 \r
1350                 try {\r
1351                         for (IA = 1; IA <= NARCS; IA++) {\r
1352                                 I = PGM[IA - 1];\r
1353                                 J = PTX[IA - 1];\r
1354                                 STAPT2 = new double[STAPT.length - IA + 1][2];\r
1355                                 for (K = IA; K <= STAPT.length; K++) {\r
1356                                         STAPT2[K - IA] = STAPT[K - 1];\r
1357                                 }\r
1358                                 if (ARCTY[IA-1] != 1) {\r
1359                                         RGM2 = new double[RGM.length - I + 1];\r
1360                                         for (K = I; K <= RGM.length; K++) {\r
1361                                                 RGM2[K - I] = RGM[K - 1];\r
1362                                         }\r
1363                                 }\r
1364                                 else {\r
1365                                         RGM2 = null;\r
1366                                 }\r
1367                                 if ((ARCTY[IA-1] == 3) || (ARCTY[IA-1] == 4)) {\r
1368                                         DEFN2 = new String[DEFN.length - J + 1];\r
1369                                         for (K = J; K <= DEFN.length; K++) {\r
1370                                                 DEFN2[K - J] = DEFN[K - 1];\r
1371                                         }\r
1372                                 }\r
1373                                 else {\r
1374                                         // DEFN2 goes to TXT in PTFUN1 which is not used for TYPES 1 and 2\r
1375                                         DEFN2 = null;\r
1376                                 }\r
1377                                 if (NARCS == 1) {\r
1378                                         PTFUN1(ARCTY[IA - 1], STAPT2, RGM2, NTX[IA - 1], DEFN2, CHNL, CHTT, VAR, REDD, raFile);\r
1379                                 } else {\r
1380                                         if (IA == 1) {\r
1381                                                 raFile.writeBytes(TX1 + IA + ") {\n");\r
1382                                         } else if (IA == NARCS) {\r
1383                                                 raFile.writeBytes("      else {\n");\r
1384                                         } else {\r
1385                                                 raFile.writeBytes(TX2 + IA + ") {\n");\r
1386                                         }\r
1387                                         PTFUN1(ARCTY[IA - 1], STAPT2, RGM2, NTX[IA - 1], DEFN2, CHNL, CHTT, VAR, REDD, raFile);\r
1388                                         if (IA == NARCS)\r
1389                                                 raFile.writeBytes("      }\n");\r
1390                                 } // else\r
1391                         } // for (IA=1; IA <= NARCS; IA++)\r
1392                 } // try\r
1393                 catch (IOException e) {\r
1394                         MipavUtil.displayError("IOException " + e + " in WRFUN1");\r
1395                         System.exit(-1);\r
1396                 }\r
1397 \r
1398         } // private void WRFUN1\r
1399 \r
1400         private void PTFUN1(int TYPE, double STAPT[][], double RGM[], int NTX, String TXT[], int CHNL, String CHTT,\r
1401                         String VAR, String REDD, RandomAccessFile raFile) {\r
1402 \r
1403                 // COMPLEX STAPT(*)\r
1404                 // CHARACTER TXT(*)*72,CHTT*2,VAR*6,REDD*6\r
1405 \r
1406                 // .......................................................................\r
1407                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1408                 // LAST UPDATE: 8 AUG 1990\r
1409                 // .......................................................................C\r
1410                 // **** LOCAL VARIABLES\r
1411 \r
1412                 int I;\r
1413                 double HA, MD, RAD;\r
1414                 double C1[] = new double[2];\r
1415                 double C2[] = new double[2];\r
1416                 double CENTR[] = new double[2];\r
1417                 // COMPLEX C1,C2,CENTR\r
1418                 // String TX1, TX1B, TX2;\r
1419                 // String TX2B, CTX1B, FMT1, FMT2, FMT3, FMT4, FMT5;\r
1420                 // CHARACTER TX1*4,TX1B*5,TX2*13,TX2B*14,CTX1B*10,\r
1421                 // +FMT1*25,FMT2*25,FMT3*14,FMT4*25,FMT5*24\r
1422 \r
1423                 // TX1 = "+"+CHTT+"*";\r
1424                 // TX1B = TX1 + "(";\r
1425                 // CTX1B=" " + TX1B;\r
1426                 // TX2=" "+VAR+" = ";\r
1427                 // TX2B=TX2+"(";\r
1428 \r
1429                 // FMT1='(A14,'//REDD//',A1,'//REDD//',A2)'\r
1430                 // FMT2='(A10,'//REDD//',A1,'//REDD//',A1)'\r
1431                 // FMT3='(A6,'//REDD//',A1)'\r
1432                 // FMT4='(A14,'//REDD//',A5,'//REDD//',A3)'\r
1433                 // FMT5='(A8,'//REDD//',A5,'//REDD//',A1)'\r
1434 \r
1435                 try {\r
1436                         if (TYPE == 1) {\r
1437                                 C1[0] = 0.5 * (STAPT[1][0] + STAPT[0][0]);\r
1438                                 C1[1] = 0.5 * (STAPT[1][1] + STAPT[0][1]);\r
1439                                 C2[0] = 0.5 * (STAPT[1][0] - STAPT[0][0]);\r
1440                                 C2[1] = 0.5 * (STAPT[1][1] - STAPT[0][1]);\r
1441                                 raFile.writeBytes("//\n");\r
1442                                 raFile.writeBytes(\r
1443                                                 VAR + "[0] = " + C1[0] + "+" + CHTT + "[0]*" + C2[0] + " - " + CHTT + "[1]*" + C2[1] + ";\n");\r
1444                                 raFile.writeBytes(\r
1445                                                 VAR + "[1] = " + C1[1] + "+" + CHTT + "[0]*" + C2[1] + " + " + CHTT + "[1]*" + C2[0] + ";\n");\r
1446                                 raFile.writeBytes("//\n");\r
1447                         } // if (TYPE == 1)\r
1448                         else if (TYPE == 2) {\r
1449                                 CENTR[0] = RGM[0];\r
1450                                 CENTR[1] = RGM[1];\r
1451                                 C1[0] = STAPT[0][0] - CENTR[0];\r
1452                                 C1[1] = STAPT[0][1] - CENTR[1];\r
1453                                 HA = 0.5 * RGM[2];\r
1454                                 MD = Math.atan2(C1[1], C1[0]) + HA;\r
1455                                 RAD = zabs(C1[0], C1[1]);\r
1456                                 raFile.writeBytes("//\n");\r
1457                                 raFile.writeBytes(VAR + "[0] = " + CENTR[0] + "+" + RAD + " * " + "Math.exp(-" + CHTT + "[1]*" + HA\r
1458                                                 + ")*" + "Math.cos(" + MD + CHTT + "[0]*" + HA + ");\n");\r
1459                                 raFile.writeBytes(VAR + "[1] = " + CENTR[1] + "+" + RAD + " * " + "Math.exp(-" + CHTT + "[1]*" + HA\r
1460                                                 + ")*" + "Math.sin(" + MD + CHTT + "[0]*" + HA + ");\n");\r
1461                                 raFile.writeBytes("//\n");\r
1462                         } // else if (TYPE == 2)\r
1463                         else if (TYPE == 3) {\r
1464                                 MD = 0.5 * (RGM[1] + RGM[0]);\r
1465                                 HA = 0.5 * (RGM[1] - RGM[0]);\r
1466                                 raFile.writeBytes("//\n");\r
1467                                 raFile.writeBytes("T[0] = " + MD + "+" + CHTT + "[0] * " + HA + ";\n");\r
1468                                 raFile.writeBytes("T[1] = " + MD + "+" + CHTT + "[1] * " + HA + ";\n");\r
1469 \r
1470                                 raFile.writeBytes(VAR + "[0] = ");\r
1471                                 // NTX = 1 if statements are entered without newlines for\r
1472                                 // multiple lines\r
1473                                 for (I = 1; I <= NTX; I++) {\r
1474                                         int index = TXT[I - 1].indexOf("ui");\r
1475                                         String realString = null;\r
1476                                         if (index == -1) {\r
1477                                                 realString = TXT[I - 1];\r
1478                                         } else {\r
1479                                                 realString = TXT[I - 1].substring(0, index);\r
1480                                         }\r
1481                                         if ((index == -1) || (index > 0)) {\r
1482                                                 raFile.writeBytes(realString);\r
1483                                         }\r
1484                                         if (I == NTX) {\r
1485                                                 raFile.writeBytes(";\n");\r
1486                                         }\r
1487                                 }\r
1488                                 raFile.writeBytes(VAR + "[1] = ");\r
1489                                 // NTX = 1 if statements are entered without newlines for\r
1490                                 // multiple lines\r
1491                                 for (I = 1; I <= NTX; I++) {\r
1492                                         int index = TXT[I - 1].indexOf("ui");\r
1493                                         String imagString = null;\r
1494                                         if ((index >= 0) && (index + 2 < TXT[I - 1].length())) {\r
1495                                                 imagString = TXT[I - 1].substring(index + 2);\r
1496                                                 raFile.writeBytes(imagString);\r
1497                                         }\r
1498                                         if (I == NTX) {\r
1499                                                 raFile.writeBytes(";\n");\r
1500                                         }\r
1501                                 }\r
1502                                 raFile.writeBytes("//\n");\r
1503                         } // else if (TYPE == 3)\r
1504                         else {\r
1505                                 MD = 0.5 * (RGM[1] + RGM[0]);\r
1506                                 HA = 0.5 * (RGM[1] - RGM[0]);\r
1507                                 raFile.writeBytes("//\n");\r
1508                                 raFile.writeBytes("T[0] = " + MD + "+" + CHTT + "[0] * " + HA + ";\n");\r
1509                                 raFile.writeBytes("T[1] = " + MD + "+" + CHTT + "[1] * " + HA + ";\n");\r
1510                                 raFile.writeBytes("     ZRAD = ");\r
1511                                 // NTX = 1 if statements are entered without newlines for\r
1512                                 // multiple lines\r
1513                                 for (I = 1; I <= NTX; I++) {\r
1514                                         raFile.writeBytes(TXT[I - 1]);\r
1515                                 }\r
1516                                 raFile.writeBytes(VAR + "[0] = Math.exp(-T[1])*(ZRAD * Math.cos(T[0]));\n");\r
1517                                 raFile.writeBytes(VAR + "[1] = Math.exp(-T[1])*(ZRAD * Math.sin(T[0]));\n");\r
1518                                 raFile.writeBytes("//\n");\r
1519                         }\r
1520                 } // try\r
1521                 catch (IOException e) {\r
1522                         MipavUtil.displayError("IOException " + e + " in PTFUN1");\r
1523                         System.exit(-1);\r
1524                 }\r
1525 \r
1526         } // private void PTFUN1\r
1527 \r
1528         private void WRSYM2(int NARCS, int ORDRG, double CENSY[], boolean REFLN, int CHNL, RandomAccessFile raFile) {\r
1529 \r
1530                 // COMPLEX CENSY\r
1531 \r
1532                 // **** TO WRITE THE CODE TO RECOVER THE BOUNDARY POINT FROM ITS\r
1533                 // SYMMETRIC\r
1534                 // **** COUNTERPART ON THE FUNDAMENTAL BOUNDARY SECTION.\r
1535 \r
1536                 // .......................................................................\r
1537                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1538                 // LAST UPDATE: 4 AUG 1990\r
1539                 // .......................................................................C\r
1540                 // LOCAL VARIABLES\r
1541 \r
1542                 int I;\r
1543                 boolean NEEDC;\r
1544 \r
1545                 NEEDC = ((CENSY[0] != 0.0) || (CENSY[1] != 0.0));\r
1546                 try {\r
1547                         raFile.writeBytes("//\n");\r
1548 \r
1549                         if (REFLN) {\r
1550                                 if (ORDRG > 1) {\r
1551                                         I = 2 * NARCS;\r
1552                                         if (NARCS > 1) {\r
1553                                                 raFile.writeBytes("IS = (IA-IB)%" + I + ";\n");\r
1554                                                 raFile.writeBytes("IR = (IA-IB-IS)/" + I + ";\n");\r
1555                                         } // if (NARCS > 1)\r
1556                                         else {\r
1557                                                 raFile.writeBytes("IS = (IA-1)%2;\n");\r
1558                                                 raFile.writeBytes("IR = (IA-1-IS)/2;\n");\r
1559                                         }\r
1560                                         raFile.writeBytes("if ((IR == 0) && (IS == 0)) {\n");\r
1561                                         raFile.writeBytes("    PARFUN[0] = ZETA[0];\n");\r
1562                                         raFile.writeBytes("    PARFUN[1] = ZETA[1];\n");\r
1563                                         raFile.writeBytes("}\n");\r
1564                                         raFile.writeBytes("else if ((IR > 0) && (IS == 0)) {\n");\r
1565                                         if (NEEDC) {\r
1566                                                 raFile.writeBytes("    PARFUN[0] = ZCEN[0] + WW[IR-1][0]*(ZETA[0] - ZCEN[0]) - "\r
1567                                                                 + "WW[IR-1][1]*(ZETA[1] - ZCEN[1]);\n");\r
1568                                                 raFile.writeBytes("    PARFUN[1] = ZCEN[1] + WW[IR-1][0]*(ZETA[1] - ZCEN[1]) + "\r
1569                                                                 + "WW[IR-1][1]*(ZETA[0] - ZCEN[0]);\n");\r
1570                                         } // if (NEEDC)\r
1571                                         else {\r
1572                                                 raFile.writeBytes("PARFUN[0] = WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\n");\r
1573                                         }\r
1574                                         raFile.writeBytes("else if ((IR == 0) && (IS > 0)) {\n");\r
1575                                         if (NEEDC) {\r
1576                                                 raFile.writeBytes(\r
1577                                                                 "    PARFUN[0] = ZCEN[0] + U2[0]*(ZETA[0]-ZCEN[0]) + " + "U2[1]*(ZETA[1]-ZCEN[1]);\n");\r
1578                                                 raFile.writeBytes(\r
1579                                                                 "    PARFUN[1] = ZCEN[1] - U2[0]*(ZETA[1]-ZCEN[1]) + " + "U2[1]*(ZETA[0]-ZCEN[0]);\n");\r
1580                                         } else {\r
1581                                                 raFile.writeBytes("    PARFUN[0] = U2[0]*ZETA[0] + U2[1]*ZETA[1];\n");\r
1582                                                 raFile.writeBytes("    PARFUN[1] = -U2[0]*ZETA[1] + U2[1]*ZETA[0];\n");\r
1583                                         }\r
1584                                         raFile.writeBytes("}\n");\r
1585                                         raFile.writeBytes("else {\n");\r
1586                                         if (NEEDC) {\r
1587                                                 raFile.writeBytes("double realPart = U2[0]*WW[IR-1][0] - U2[1]*WW[IR-1][1];\n");\r
1588                                                 raFile.writeBytes("double imagPart = U2[0]*WW[IR-1][1] + U2[1]*WW[IR-1][0];\n");\r
1589                                                 raFile.writeBytes("PARFUN[0] = ZCEN[0] + realPart*(ZETA[0]-ZCEN[0]) + "\r
1590                                                                 + "imagPart*(ZETA[1]-ZCEN[1]);\n");\r
1591                                                 raFile.writeBytes("PARFUN[1] = ZCEN[1] - realPart*(ZETA[1]-ZCEN[1]) + "\r
1592                                                                 + "imagPart*(ZETA[0]-ZCEN[0]);\n");\r
1593                                         } else {\r
1594                                                 raFile.writeBytes("double realPart = U2[0]*WW[IR-1][0] - U2[1]*WW[IR-1][1];\n");\r
1595                                                 raFile.writeBytes("double imagPart = U2[0]*WW[IR-1][1] + U2[1]*WW[IR-1][0];\n");\r
1596                                                 raFile.writeBytes("PARFUN[0] = realPart * ZETA[0] + imagPart * ZETA[1];\n");\r
1597                                                 raFile.writeBytes("PARFUN[1] = -realPart * ZETA[1] + imagPart * ZETA[0];\n");\r
1598                                         }\r
1599                                         raFile.writeBytes("}\n");\r
1600                                 } // if (ORDRG > 1)\r
1601                                 else { // ORDRG <= 1\r
1602                                         if (NARCS > 1) {\r
1603                                                 raFile.writeBytes("IS = IA - IB;\n");\r
1604                                         } else {\r
1605                                                 raFile.writeBytes("IS = IA - 1;\n");\r
1606                                         }\r
1607                                         raFile.writeBytes("if (IS == 0) {\n");\r
1608                                         raFile.writeBytes("    PARFUN[0] = ZETA[0];\n");\r
1609                                         raFile.writeBytes("    PARFUN[1] = ZETA[1];\n");\r
1610                                         raFile.writeBytes("}\n");\r
1611                                         raFile.writeBytes("else {\n");\r
1612                                         if (NEEDC) {\r
1613                                                 raFile.writeBytes(\r
1614                                                                 "    PARFUN[0] = ZCEN[0] + U2[0]*(ZETA[0]-ZCEN[0]) + " + "U2[1]*(ZETA[1]-ZCEN[1]);\n");\r
1615                                                 raFile.writeBytes(\r
1616                                                                 "    PARFUN[1] = ZCEN[1] - U2[0]*(ZETA[1]-ZCEN[1]) + " + "U2[1]*(ZETA[0]-ZCEN[0]);\n");\r
1617                                         } else {\r
1618                                                 raFile.writeBytes("    PARFUN[0] = U2[0]*(ZETA[0]-ZCEN[0]) + " + "U2[1]*(ZETA[1]-ZCEN[1]);\n");\r
1619                                                 raFile.writeBytes(\r
1620                                                                 "    PARFUN[1] =  -U2[0]*(ZETA[1]-ZCEN[1]) + " + "U2[1]*(ZETA[0]-ZCEN[0]);\n");\r
1621                                         }\r
1622                                         raFile.writeBytes("}\n");\r
1623                                 } // else ORDRG <= 1\r
1624                         } // if (REFLN)\r
1625                         else { // !REFLN\r
1626                                 if (NARCS > 1) {\r
1627                                         raFile.writeBytes("IR = (IA - IB)/" + NARCS + ";\n");\r
1628                                 } else {\r
1629                                         raFile.writeBytes("IR = IA - 1;\n");\r
1630                                 }\r
1631                                 raFile.writeBytes("if (IR == 0) {\n");\r
1632                                 raFile.writeBytes("PARFUN[0] = ZETA[0]);\n");\r
1633                                 raFile.writeBytes("PARFUN[1] = ZETA[1]);\n");\r
1634                                 raFile.writeBytes("}\n");\r
1635                                 raFile.writeBytes("else {\n");\r
1636                                 if (NEEDC) {\r
1637                                         raFile.writeBytes("    PARFUN[0] = ZCEN[0] + WW[IR-1][0]*(ZETA[0] - ZCEN[0]) - "\r
1638                                                         + "WW[IR-1][1]*(ZETA[1] - ZCEN[1]);\n");\r
1639                                         raFile.writeBytes("    PARFUN[1] = ZCEN[1] + WW[IR-1][0]*(ZETA[1] - ZCEN[1]) + "\r
1640                                                         + "WW[IR-1][1]*(ZETA[0] - ZCEN[0]);\n");\r
1641                                 } else {\r
1642                                         raFile.writeBytes("    PARFUN[0] = WW[IR-1][0]*(ZETA[0] - ZCEN[0]) - "\r
1643                                                         + "WW[IR-1][1]*(ZETA[1] - ZCEN[1]);\n");\r
1644                                         raFile.writeBytes("    PARFUN[1] = WW[IR-1][0]*(ZETA[1] - ZCEN[1]) + "\r
1645                                                         + "WW[IR-1][1]*(ZETA[0] - ZCEN[0]);\n");\r
1646                                 }\r
1647                                 raFile.writeBytes("}\n");\r
1648                         } // else !REFLN\r
1649                 } // try\r
1650                 catch (IOException e) {\r
1651                         MipavUtil.displayError("IOException " + e + " in WRSYM2");\r
1652                         System.exit(-1);\r
1653                 }\r
1654 \r
1655         } // private void WRSYM2\r
1656 \r
1657         private void WRFUN2(int NARCS, int MNARC, double STAPT[][], int ARCTY[], int PGM[], double RGM[], int PTX[],\r
1658                         int NTX[], String DEFN[], int CHNL, String CHIA, String CHTT, String VAR, boolean NUMDER[], String REDD,\r
1659                         RandomAccessFile raFile) {\r
1660 \r
1661                 // COMPLEX STAPT(*)\r
1662                 // CHARACTER DEFN(*)*72,CHIA*2,CHTT*2,VAR*6,REDD*6\r
1663 \r
1664                 // **** TO WRITE THE SOURCE CODE FOR DPARFN IN THE CASE WHERE NO\r
1665                 // **** SYMMETRY IS INVOLVED.\r
1666 \r
1667                 // .......................................................................\r
1668                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1669                 // LAST UPDATE: 4 AUG 1990\r
1670                 // .......................................................................\r
1671 \r
1672                 // LOCAL VARIABLES\r
1673 \r
1674                 int IA, I, J1, J2, N1, N2;\r
1675                 String TX1, TX2;\r
1676                 // String FMT1, FMT2;\r
1677                 // CHARACTER TX1*16,TX2*21,FMT1*11,FMT2*11\r
1678                 // EXTERNAL PTFUN2\r
1679                 double STAPT2[][];\r
1680                 double RGM2[];\r
1681                 String DEFN2[];\r
1682                 String DEFN3[];\r
1683                 int K;\r
1684 \r
1685                 TX1 = "      if (" + CHIA + " == ";\r
1686                 TX2 = "      else if (" + CHIA + " == ";\r
1687                 // FMT1="(A16,I3,A6)";\r
1688                 // FMT2="(A21,I3,A6)";\r
1689 \r
1690                 try {\r
1691                         for (IA = 1; IA <= NARCS; IA++) {\r
1692                                 I = PGM[IA - 1];\r
1693                                 J1 = PTX[IA - 1];\r
1694                                 J2 = PTX[IA + MNARC - 1];\r
1695                                 N1 = NTX[IA - 1];\r
1696                                 N2 = NTX[IA + MNARC - 1];\r
1697                                 STAPT2 = new double[STAPT.length - IA + 1][2];\r
1698                                 for (K = IA; K <= STAPT.length; K++) {\r
1699                                         STAPT2[K - IA][0] = STAPT[K - 1][0];\r
1700                                         STAPT2[K - IA][1] = STAPT[K - 1][1];\r
1701                                 }\r
1702                                 if ((ARCTY[IA-1] == 2) || ((!NUMDER[IA-1]) && ((ARCTY[IA-1] == 3) || (ARCTY[IA-1] == 4)))) {\r
1703                                         RGM2 = new double[RGM.length - I + 1];\r
1704                                         for (K = I; K <= RGM.length; K++) {\r
1705                                                 RGM2[K - I] = RGM[K - 1];\r
1706                                         }\r
1707                                 }\r
1708                                 else {\r
1709                                         RGM2 = null;\r
1710                                 }\r
1711                                 if ((!NUMDER[IA-1]) && (ARCTY[IA-1] == 4)) {\r
1712                                         DEFN2 = new String[DEFN.length - J1 + 1];\r
1713                                         for (K = J1; K <= DEFN.length; K++) {\r
1714                                                 DEFN2[K - J1] = DEFN[K - 1];\r
1715                                         }\r
1716                                 }\r
1717                                 else {\r
1718                                         DEFN2 = null;\r
1719                                 }\r
1720                                 if ((!NUMDER[IA-1]) && ((ARCTY[IA-1] == 3) || (ARCTY[IA-1] == 4))) {\r
1721                                         DEFN3 = new String[DEFN.length - J2 + 1];\r
1722                                         for (K = J2; K <= DEFN.length; K++) {\r
1723                                                 DEFN3[K - J2] = DEFN[K - 1];\r
1724                                         }\r
1725                                 }\r
1726                                 else {\r
1727                                         DEFN3 = null;\r
1728                                 }\r
1729                                 if (NARCS == 1) {\r
1730                                         PTFUN2(ARCTY[IA - 1], STAPT2, RGM2, N1, DEFN2, N2, DEFN3, CHNL, CHTT, VAR, " 1", NUMDER[IA - 1],\r
1731                                                         REDD, raFile);\r
1732                                 } else {\r
1733                                         if (IA == 1) {\r
1734                                                 raFile.writeBytes(TX1 + IA + ") {\n");\r
1735                                         } else if (IA == NARCS) {\r
1736                                                 raFile.writeBytes("      else {\n");\r
1737                                         } else {\r
1738                                                 raFile.writeBytes(TX2 + IA + ") {\n");\r
1739                                         }\r
1740                                         PTFUN2(ARCTY[IA - 1], STAPT2, RGM2, N1, DEFN2, N2, DEFN3, CHNL, CHTT, VAR, CHIA, NUMDER[IA - 1],\r
1741                                                         REDD, raFile);\r
1742                                         if (IA == NARCS) {\r
1743                                                 raFile.writeBytes("      }\n");\r
1744                                         }\r
1745                                 } // else\r
1746                         } // for (IA=1; IA <= NARCS; IA++)\r
1747                 } // try\r
1748                 catch (IOException e) {\r
1749                         MipavUtil.displayError("IOException " + e + " in WRFUN2");\r
1750                         System.exit(-1);\r
1751                 }\r
1752 \r
1753         } // private void WRFUN2\r
1754 \r
1755         private void PTFUN2(int TYPE, double STAPT[][], double RGM[], int NTX1, String TXT1[], int NTX2, String TXT2[],\r
1756                         int CHNL, String CHTT, String VAR, String CHIA, boolean NUMDER, String REDD, RandomAccessFile raFile) {\r
1757                 // COMPLEX STAPT(*)\r
1758                 // CHARACTER TXT1(*)*72,TXT2(*)*72,CHTT*2,VAR*6,CHIA*2,REDD*6\r
1759 \r
1760                 // .......................................................................\r
1761                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1762                 // LAST UPDATE: 8 AUG 1990\r
1763                 // .......................................................................C\r
1764                 // **** LOCAL VARIABLES\r
1765 \r
1766                 int I;\r
1767                 double HA, MD, RAD, A, R;\r
1768                 double C1[] = new double[2];\r
1769                 double CENTR[] = new double[2];\r
1770                 // COMPLEX C1,CENTR\r
1771                 // String TX1, TX1B, TX2, TX2B, TX3, FMT1, FMT2, FMT3, FMT4, FMT5;\r
1772                 // CHARACTER TX1*4,TX1B*5,TX2*13,TX2B*14,TX3*39,\r
1773                 // FMT1*25,FMT2*15,FMT3*15,FMT4*25,FMT5*24\r
1774                 // String TX3R, TX3I;\r
1775 \r
1776                 // TX1 = "+"+CHTT+"*";\r
1777                 // TX1B = TX1 + "(";\r
1778                 // TX2 = " "+VAR+" = ";\r
1779                 // TX2B = TX2 + "(";\r
1780                 // TX3R = " "+VAR+"[0] = (ZDER * Math.cos(T) - ZRAD * Math.sin(T))*(";\r
1781                 // TX3I = " "+VAR+"[1] = (ZDER * Math.sin(T) + ZRAD * Math.cos(T))*(";\r
1782                 // TX3=TX2//'(ZDER+UI*ZRAD)*EXP(UI*T)*('\r
1783 \r
1784                 // FMT1='(A14,'//REDD//',A1,'//REDD//',A2)'\r
1785                 // FMT2='(A39,'//REDD//',A1)'\r
1786                 // FMT3='(A13,'//REDD//',A2)'\r
1787                 // FMT4='(A14,'//REDD//',A5,'//REDD//',A3)'\r
1788                 // FMT5='(A8,'//REDD//',A5,'//REDD//',A1)'\r
1789 \r
1790                 try {\r
1791                         if (TYPE == 1) {\r
1792                                 C1[0] = 0.5 * (STAPT[1][0] - STAPT[0][0]);\r
1793                                 C1[1] = 0.5 * (STAPT[1][1] - STAPT[0][1]);\r
1794                                 raFile.writeBytes("//\n");\r
1795                                 R = C1[0];\r
1796                                 A = C1[1];\r
1797                                 raFile.writeBytes("      " + VAR + "[0] = " + R + ";\n");\r
1798                                 raFile.writeBytes("      " + VAR + "[1] = " + A + ";\n");\r
1799                                 raFile.writeBytes("//\n");\r
1800                         } // if (TYPE == 1)\r
1801                         else if (TYPE == 2) {\r
1802                                 CENTR[0] = RGM[0];\r
1803                                 CENTR[1] = RGM[1];\r
1804                                 C1[0] = STAPT[0][0] - CENTR[0];\r
1805                                 C1[1] = STAPT[0][1] - CENTR[1];\r
1806                                 HA = 0.5 * RGM[2];\r
1807                                 MD = Math.atan2(C1[1], C1[0]) + HA;\r
1808                                 RAD = zabs(C1[0], C1[1]);\r
1809                                 raFile.writeBytes("//\n");\r
1810                                 raFile.writeBytes(\r
1811                                                 VAR + "[0] = -" + RAD + "*" + HA + "*Math.sin(" + MD + "+" + CHTT + "*" + HA + ");\n");\r
1812                                 raFile.writeBytes(VAR + "[1] = " + RAD + "*" + HA + "*Math.cos(" + MD + "+" + CHTT + "*" + HA + ");\n");\r
1813                                 raFile.writeBytes("//\n");\r
1814                         } // else if (TYPE == 2)\r
1815                         else if (NUMDER) {\r
1816                                 raFile.writeBytes("//\n");\r
1817                                 raFile.writeBytes("      " + VAR + " = ZDPARF(" + CHIA + "," + CHTT + ");\n");\r
1818                                 raFile.writeBytes("//\n");\r
1819                         } // else if (NUMDER)\r
1820                         else if (TYPE == 3) {\r
1821                                 MD = 0.5 * (RGM[1] + RGM[0]);\r
1822                                 HA = 0.5 * (RGM[1] - RGM[0]);\r
1823                                 raFile.writeBytes("//\n");\r
1824                                 raFile.writeBytes("      T = " + MD + "+" + CHTT + "*" + "(" + HA + ");\n");\r
1825                                 raFile.writeBytes("      " + VAR + " = " + HA + "*(");\r
1826                                 for (I = 1; I <= NTX2; I++) {\r
1827                                         raFile.writeBytes(TXT2[I - 1]);\r
1828                                 }\r
1829                                 raFile.writeBytes(");\n");\r
1830                                 raFile.writeBytes("//\n");\r
1831                         } // else if (TYPE == 3)\r
1832                         else {\r
1833                                 MD = 0.5 * (RGM[1] + RGM[0]);\r
1834                                 HA = 0.5 * (RGM[1] - RGM[0]);\r
1835                                 raFile.writeBytes("//\n");\r
1836                                 raFile.writeBytes("      T = " + MD + "+" + CHTT + "*" + "(" + HA + ");\n");\r
1837                                 raFile.writeBytes("      ZRAD = ");\r
1838                                 for (I = 1; I <= NTX1; I++) {\r
1839                                         raFile.writeBytes(TXT1[I - 1]);\r
1840                                 }\r
1841                                 raFile.writeBytes(";\n");\r
1842                                 raFile.writeBytes("      ZDER = ");\r
1843                                 for (I = 1; I <= NTX2; I++) {\r
1844                                         raFile.writeBytes(TXT2[I - 1]);\r
1845                                 }\r
1846                                 raFile.writeBytes(";\n");\r
1847                                 raFile.writeBytes("      " + VAR + "[0] = (ZDER * Math.cos(T) - ZRAD * Math.sin(T))*(" + HA + ");\n");\r
1848                                 raFile.writeBytes("      " + VAR + "[1] = (ZDER * Math.sin(T) + ZRAD * Math.cos(T))*(" + HA + ");\n");\r
1849                                 raFile.writeBytes("//\n");\r
1850                         } // else\r
1851                 } // try\r
1852                 catch (IOException e) {\r
1853                         MipavUtil.displayError("IOException " + e + " in PTFUN2");\r
1854                         System.exit(-1);\r
1855                 }\r
1856 \r
1857         }\r
1858 \r
1859         private void WRSYM3(int NARCS, int ORDRG, boolean REFLN, int CHN, RandomAccessFile raFile) {\r
1860 \r
1861                 // **** TO WRITE THE CODE TO RECOVER THE DERIVATIVE FROM ITS SYMMETRIC\r
1862                 // **** COUNTERPART ON THE FUNDAMENTAL BOUNDARY SECTION.\r
1863 \r
1864                 // .......................................................................\r
1865                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
1866                 // LAST UPDATE: 4 AUG 1990\r
1867                 // .......................................................................C\r
1868                 // LOCAL VARIABLES\r
1869 \r
1870                 int I;\r
1871 \r
1872                 try {\r
1873                         raFile.writeBytes("//\n");\r
1874 \r
1875                         if (REFLN) {\r
1876                                 if (ORDRG > 1) {\r
1877                                         I = 2 * NARCS;\r
1878                                         if (NARCS > 1) {\r
1879                                                 raFile.writeBytes("      IS = (IA-IB)%" + I + ";\n");\r
1880                                                 raFile.writeBytes("      IR=(IA-IB-IS)/" + I + ";\n");\r
1881                                         } else {\r
1882                                                 raFile.writeBytes("      IS=(IA-1)%2;\n");\r
1883                                                 raFile.writeBytes("      IR=(IA-1-IS)/2;\n");\r
1884                                         }\r
1885                                         raFile.writeBytes("      if ((IR == 0) && (IS == 0)) {\n");\r
1886                                         raFile.writeBytes("          DPARFN[0] =ZETA[0];\n");\r
1887                                         raFile.writeBytes("          DPARFN[1] =ZETA[1];\n");\r
1888                                         raFile.writeBytes("      }\n");\r
1889                                         raFile.writeBytes("      else if ((IR > 0) && (IS == 0)) {\n");\r
1890                                         raFile.writeBytes("          DPARFN[0] = WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\n");\r
1891                                         raFile.writeBytes("          DPARFN[1] = WW[IR-1][0]*ZETA[1] + WW[IR-1][1]*ZETA[0];\n");\r
1892                                         raFile.writeBytes("      }\n");\r
1893                                         raFile.writeBytes("       else if ((IR == 0 && (IS > 0)) {\n");\r
1894                                         raFile.writeBytes("           DPARFN[0] = -U2[0]*ZETA[0] -U2[1]*ZETA[1];\n");\r
1895                                         raFile.writeBytes("           DPARFN[1] = U2[0]*ZETA[1] -U2[1]*ZETA[0];\n");\r
1896                                         raFile.writeBytes("       }\n");\r
1897                                         raFile.writeBytes("       else {\n");\r
1898                                         raFile.writeBytes("           double realPart = -U2[0]*WW[IR-1][0] + U2[1]*WW[IR-1][1]);\n");\r
1899                                         raFile.writeBytes("           double imagPart = -U2[0]*WW[IR-1][1] + U2[1]*WW[IR-1][0]);\n");\r
1900                                         raFile.writeBytes("           DPARFN[0] = realPart*ZETA[0] + imagPart*ZETA[1]);\n");\r
1901                                         raFile.writeBytes("           DPARFN[1] = imagPart*ZETA[0] - realPart*ZETA[1]);\n");\r
1902                                         raFile.writeBytes("       }\n");\r
1903                                 } // if (ORDRG > 1)\r
1904                                 else {\r
1905                                         if (NARCS > 1) {\r
1906                                                 raFile.writeBytes("      IS=IA-IB;\n");\r
1907                                         } else {\r
1908                                                 raFile.writeBytes("      IS=IA-1;\n");\r
1909                                         }\r
1910                                         raFile.writeBytes("      if (IS == 0) {\n");\r
1911                                         raFile.writeBytes("          DPARFN[0] = ZETA[0]\n");\r
1912                                         raFile.writeBytes("          DPARFN[0] = ZETA[0]\n");\r
1913                                         raFile.writeBytes("      }\n");\r
1914                                         raFile.writeBytes("         else {\n");\r
1915                                         raFile.writeBytes("             DPARFN[0] = -U2[0]*ZETA[0] + U2[1]*ZETA[1]);\n");\r
1916                                         raFile.writeBytes("             DPARFN[1] = U2[0]*ZETA[1] - U2[1]*ZETA[0]);\n");\r
1917                                         raFile.writeBytes("         }\n");\r
1918                                 } // else\r
1919                         } // if (REFLN)\r
1920                         else {\r
1921                                 if (NARCS > 1) {\r
1922                                         raFile.writeBytes("      IR=(IA-IB)/" + NARCS + ";\n");\r
1923                                 } else {\r
1924                                         raFile.writeBytes("      IR=IA-1;\n");\r
1925                                 }\r
1926                                 raFile.writeBytes("      if (IR == 0) {\n");\r
1927                                 raFile.writeBytes("          DPARFN[0] = ZETA[0];\n");\r
1928                                 raFile.writeBytes("          DPARFN[1] = ZETA[1];\n");\r
1929                                 raFile.writeBytes("      }\n");\r
1930                                 raFile.writeBytes("      else {\n");\r
1931                                 raFile.writeBytes("          DPARFN[0]= WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\n");\r
1932                                 raFile.writeBytes("          DPARFN[1]= WW[IR-1][0]*ZETA[1] + WW[IR-1][1]*ZETA[0];\n");\r
1933                                 raFile.writeBytes("      }\n");\r
1934                         }\r
1935                 } // try\r
1936                 catch (IOException e) {\r
1937                         MipavUtil.displayError("IOException " + e + " in WRSYM3");\r
1938                         System.exit(-1);\r
1939                 }\r
1940 \r
1941         } // private void WRSYM3\r
1942 \r
1943         private void TSTPLT(String JBNM, double MXMIS, double MXDIF, int NARCS, double PSD, double MINPD, double MAXPD,\r
1944                         int CHNL, int IER[]) {\r
1945                 // CHARACTER*4 JBNM\r
1946 \r
1947                 // ......................................................................\r
1948 \r
1949                 // 1. TSTPLT\r
1950                 // TESTS THE PARAMETRIC FUNCTION ROUTINES PARFUN AND DPARFN\r
1951                 // FOR CONSISTENCY AND OUTPUTS BOUNDARY POINTS FOR PLOTTING.\r
1952 \r
1953                 // 2. PURPOSE\r
1954                 // THE ROUTINE FIRST CHECKS THAT THE PARAMETRIC FUNCTION\r
1955                 // ROUTINE PARFUN IS CONSISTENT WITH RESPECT TO ITS DEFINITION\r
1956                 // OF ANY CORNERS ON THE BOUNDARY. THIS IS DONE BY CHECKING\r
1957                 // THAT THE COMPUTED POINT AT THE END OF EACH ARC MATCHES THE\r
1958                 // COMPUTED POINT AT THE START OF THE NEXT ONE. IF ALL THE\r
1959                 // RELATIVE MISFIT ERRORS AT CORNERS ARE LESS THAN\r
1960                 // 10*(UNIT ROUNDOFF) THEN ALL CORNERS ARE CONSIDERED TO FIT\r
1961                 // SATISFACTORILY, OTHERWISE THE MAXIMUM RELATIVE MISFIT\r
1962                 // ERROR IS REPORTED.\r
1963 \r
1964                 // THE SECOND PURPOSE OF THE ROUTINE IS TO OUTPUT TO A DATA\r
1965                 // FILE THE COORDINATES OF A NUMBER OF POINTS ON THE BOUNDARY.\r
1966                 // THE BOUNDARY POINTS ARE SELECTED ADAPTIVELY TO MEET THE\r
1967                 // PLOTTING RESOLUTION SPECIFICATIONS DEFINED BY THE ARGUMENTS\r
1968                 // PSD,MINPD,MAXPD (SEE BELOW). THE HOPE IS THAT THE USER MAY\r
1969                 // EASILY FEED THESE DATA POINTS TO HIS LOCAL GRAPH PLOTTING\r
1970                 // ROUTINES SO AS TO CONSTRUCT A PLOT OF THE BOUNDARY. THIS\r
1971                 // WILL PROVIDE AN ESSENTIAL VISUAL CHECK ON THE VALIDITY OF\r
1972                 // THE ROUTINE PARFUN. THE OUTPUT DATA FILE IS AUTOMATICALLY\r
1973                 // NAMED <JBNM>zz.\r
1974 \r
1975                 // THE THIRD PURPOSE OF THE ROUTINE IS TO CHECK PARFUN AND\r
1976                 // DPARFN FOR MUTUAL CONSISTECY. THIS IS DONE BY COMPUTING\r
1977                 // TWO POINT FINITE DIFFERENCE APPROXIMATIONS TO DPARFN.\r
1978                 // THESE DIFFERENCE APPROXIMATIONS ARE COMPUTED AT EACH BOUND-\r
1979                 // ARY POINT THAT IS OUTPUT FOR PLOTTING AND ALSO AT NEARBY\r
1980                 // POINTS WHICH LIE JUST O F F THE BOUNDARY. THIS LATTER\r
1981                 // COMPARISON ALSO TESTS PARFUN AND DPARFN FOR CORRECTNESS IN\r
1982                 // ACCEPTING COMPLEX PARAMETER VALUES. A RELATIVE ERROR IN\r
1983                 // THE FINITE DIFFERENCE APPROXIMATION GREATER THAN 0.1 IS\r
1984                 // REPORTED AS A POSSIBLE LOGICAL INCONSISTENCY BETWEEN PARFUN\r
1985                 // AND DPARFN. (THE CRITICAL RELATIVE ERROR VALUE OF 0.1 CAN\r
1986                 // BE ALTERED BY ADJUSTING THE LOCAL PARAMETER DTOL).\r
1987 \r
1988                 // 3. CALLING SEQUENCE\r
1989                 // CALL TSTPLT(JBNM,MXMIS,MXDIF,NARCS,PSD,MINPD,MAXPD,CHNL,IER)\r
1990 \r
1991                 // PARAMETERS\r
1992                 // ON ENTRY\r
1993                 // JBNM - CHARACTER*4\r
1994                 // THE JOB NAME. THIS IS USED TO CREATE THE OUTPUT\r
1995                 // FILE WITH FILENAME\r
1996 \r
1997                 // <JBNM>zz ,\r
1998 \r
1999                 // WHERE <JBNM> DENOTES THE VALUE OF VARIABLE JBNM\r
2000                 // WITH ANY TRAILING SPACES DELETED.\r
2001 \r
2002                 // NARCS - INTEGER\r
2003                 // THE NUMBER OF ANALYTIC ARCS THAT MAKE UP THE\r
2004                 // W H O L E BOUNDARY OF THE PHYSICAL DOMAIN.\r
2005 \r
2006                 // PSD - REAL\r
2007                 // THE PLOTTING SIZE FOR THE DOMAIN IN ANY APPROPR-\r
2008                 // IATE UNITS. IF PSD .LE. 0.0 THEN IT IS ASSIGNED\r
2009                 // THE DEFAULT VALUE OF 160.0 (A REASONBLE WIDTH IN\r
2010                 // MM FOR PLOTTING ON A4 PAPER).\r
2011 \r
2012                 // MINPD - REAL\r
2013                 // THE MINIMUM SIGNIFICANT PLOTTING DISTANCE, IN THE\r
2014                 // SAME UNITS AS PSD. IF PSD .LE. 0.0 THEN IT IS\r
2015                 // ASSIGNED THE DEFAULT VALUE OF 2.0.\r
2016 \r
2017                 // MAXPD - REAL\r
2018                 // THE MAXIMUM ALLOWED PLOTTING DISTANCE, IN THE\r
2019                 // SAME UNITS AS PSD. IF PSD .LE. 0.0 THEN IT IS\r
2020                 // ASSIGNED THE DEFAULT VALUE OF 5.0. THE LARGER\r
2021                 // MAXPD, THE COARSER WILL BE THE RESOLUTION OF THE\r
2022                 // BOUNDARY POINTS OUTPUT TO <JBNM>zz, BUT THE\r
2023                 // QUICKER THEY WILL BE COMPUTED.\r
2024 \r
2025                 // CHNL - INTEGER\r
2026                 // DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
2027                 // WRITING THE FILE <JBNM>zz.\r
2028                 // ON EXIT\r
2029                 // MXMIS - REAL\r
2030                 // THE MAXIMUM RELATIVE MISFIT ERROR OVER ALL\r
2031                 // CORNER POINTS\r
2032 \r
2033                 // MXDIF - REAL\r
2034                 // THE MAXIMUM RELATIVE ERROR IN FINITE DIFFERENCE\r
2035                 // APPROXIMATIONS TO DPARFN OVER ALL BOUNDARY\r
2036                 // POINTS OUTPUT TO <JBNM>zz AND NEARBY POINTS OFF\r
2037                 // THE BOUNDARY.\r
2038 \r
2039                 // PSD - REAL\r
2040                 // IF PSD .LE. 0.0 ON ENTRY THEN IT WILL HAVE\r
2041                 // THE DEFAULT VALUE OF 160.0 ON EXIT.\r
2042 \r
2043                 // MINPD - REAL\r
2044                 // IF PSD .LE. 0.0 ON ENTRY THEN MINPD WILL HAVE\r
2045                 // THE DEFAULT VALUE OF 2.0 ON EXIT\r
2046 \r
2047                 // MAXPD - REAL\r
2048                 // IF PSD .LE. 0.0 ON ENTRY THEN MAXPD WILL HAVE\r
2049                 // THE DEFAULT VALUE OF 5.0 ON EXIT\r
2050 \r
2051                 // IER - INTEGER\r
2052                 // IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
2053                 // A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
2054                 // WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
2055                 // IER=0 - NORMAL EXIT.\r
2056                 // IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
2057                 // BE SELF EXPLANATORY.\r
2058 \r
2059                 // 4. SUBROUTINES OR FUNCTIONS NEEDED\r
2060                 // - THE CONFPACK LIBRARY.\r
2061                 // - THE REAL FUNCTION R1MACH.\r
2062                 // - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN.\r
2063 \r
2064                 // 5. FURTHER COMMENTS\r
2065                 // - A SUMMARY LISTING IS AUTOMATICALLY WRITTEN ON THE\r
2066                 // STANDARD OUTPUT CHANNEL.\r
2067                 // - THE OUTPUT FILE <JBNM>zz CONTAINS COORDINATE PAIRS\r
2068 \r
2069                 // X Y\r
2070 \r
2071                 // FOR POINTS ON THE PHYSICAL BOUNDARY, WITH ONE PAIR\r
2072                 // PER LINE.\r
2073 \r
2074                 // ......................................................................\r
2075                 // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
2076                 // LAST UPDATE: 6 JULY 1990\r
2077                 // ......................................................................C\r
2078                 // LOCAL VARIABLES\r
2079 \r
2080                 int I, IA;\r
2081                 int IMX = 0;\r
2082                 double TINC = 0.0;\r
2083                 double TMX = 0.0;\r
2084                 double A1, DIFF, ERR, HH, MINC, RMAX, RMEAN, RMIN, T, TOL1, TSD;\r
2085                 double TT[] = new double[2];\r
2086                 // REAL TT(2)\r
2087                 double C1[] = new double[2];\r
2088                 double C2[] = new double[2];\r
2089                 double CENTR[] = new double[2];\r
2090                 double ZZ0[] = new double[2];\r
2091                 double DZZ[] = new double[2];\r
2092                 double NDZZ[] = new double[2];\r
2093                 // COMPLEX C1,C2,CENTR,ZZ0,DZZ,NDZZ;\r
2094                 double ZZ[][] = new double[2][2];\r
2095                 // COMPLEX ZZ(2)\r
2096                 // CHARACTER OFL*6\r
2097                 final int MNARC = 200;\r
2098                 final double DTOL = 1.0E-1;\r
2099                 final int NH = 4;\r
2100                 boolean ATEND, FIRST, WARND;\r
2101                 boolean LNSEG[] = new boolean[MNARC];\r
2102                 double PIN[] = new double[2];\r
2103                 double PAROUT[];\r
2104                 int zzindex;\r
2105                 double cr[] = new double[1];\r
2106                 double ci[] = new double[1];\r
2107                 // EXTERNAL DPARFN,LINSEG,PARFUN,R1MACH,WRHEAD,WRTAIL,ZDPARF\r
2108                 // COMPLEX PARFUN,DPARFN,ZDPARF\r
2109 \r
2110                 // **** WRITE CONFPACK HEADING\r
2111 \r
2112                 WRHEAD(7, 0, null);\r
2113 \r
2114                 if (NARCS > MNARC) {\r
2115                         IER[0] = 59;\r
2116                         WRTAIL(7, 0, IER[0], null);\r
2117                 }\r
2118 \r
2119                 // 1 FORMAT(A45)\r
2120                 // 2 FORMAT(A45,I4)\r
2121                 // 3 FORMAT(A45,E10.3)\r
2122                 // 4 FORMAT(//,T17,A)\r
2123 \r
2124                 TOL1 = 10.0 * EPS;\r
2125 \r
2126                 // **** CHECK THAT ALL ARCS MEET AT CORNER POINTS\r
2127 \r
2128                 IER[0] = 0;\r
2129                 CENTR[0] = 0.0;\r
2130                 CENTR[1] = 0.0;\r
2131                 MXMIS = 0.0;\r
2132                 for (IA = 1; IA <= NARCS; IA++) {\r
2133                         if (IA == 1) {\r
2134                                 I = NARCS;\r
2135                         } else {\r
2136                                 I = IA - 1;\r
2137                         }\r
2138                         PIN[0] = -1.0;\r
2139                         PIN[1] = 0.0;\r
2140                         C1 = PARFUN(IA, PIN);\r
2141                         CENTR[0] = CENTR[0] + C1[0];\r
2142                         CENTR[1] = CENTR[1] + C1[1];\r
2143                         A1 = zabs(C1[0], C1[1]);\r
2144                         PIN[0] = 1.0;\r
2145                         PIN[1] = 0.0;\r
2146                         C2 = PARFUN(I, PIN);\r
2147                         ERR = zabs(C1[0] - C2[0], C1[1] - C2[1]);\r
2148                         if (A1 >= 1.0) {\r
2149                                 ERR = ERR / A1;\r
2150                         }\r
2151                         if (ERR > MXMIS) {\r
2152                                 IMX = IA;\r
2153                                 MXMIS = ERR;\r
2154                         }\r
2155                 } // for (IA=1; IA <= NARCS; IA++)\r
2156                 if (MXMIS >= TOL1) {\r
2157                         System.out.println("MAXIMUM CORNER MISFIT: " + MXMIS);\r
2158                         System.out.println("OCCURS AT CORNER: " + IMX);\r
2159                 } else {\r
2160                         System.out.println("ALL ARCS FIT AT CORNERS:");\r
2161                 }\r
2162 \r
2163                 // **** ESTIMATE THE DIAMETER (TSD) OF THE PHYSICAL DOMAIN\r
2164 \r
2165                 CENTR[0] = CENTR[0] / NARCS;\r
2166                 CENTR[1] = CENTR[1] / NARCS;\r
2167                 TSD = 0.0;\r
2168                 HH = 2.0 / (double) (NH);\r
2169                 for (IA = 1; IA <= NARCS; IA++) {\r
2170                         T = -1.0;\r
2171                         for (I = 1; I <= NH; I++) {\r
2172                                 T = T + HH;\r
2173                                 PIN[0] = T;\r
2174                                 PIN[1] = 0;\r
2175                                 PAROUT = PARFUN(IA, PIN);\r
2176                                 C1[0] = PAROUT[0] - CENTR[0];\r
2177                                 C1[1] = PAROUT[1] - CENTR[1];\r
2178                                 A1 = zabs(C1[0], C1[1]);\r
2179                                 TSD = Math.max(TSD, A1);\r
2180                         } // for (I=1; I <= NH; I++)\r
2181                 } // for (IA=1; IA <= NARCS; IA++)\r
2182                 TSD = 2.0 * TSD;\r
2183 \r
2184                 // **** DETERMINE WHICH ARCS (IF ANY) ARE LINE SEGMENTS\r
2185 \r
2186                 LINSEG(LNSEG, NARCS);\r
2187 \r
2188                 // **** OPEN FILE TO RECEIVE BOUNDARY DATA POINTS FOR PLOTTING\r
2189 \r
2190                 // L=INDEX(JBNM,' ')-1\r
2191                 // IF (L.EQ.-1) L=4\r
2192                 // OFL=JBNM(1:L)//'zz'\r
2193                 // OPEN(CHNL,FILE=OFL)\r
2194                 // Use global zzset instead\r
2195 \r
2196                 // **** SET DEFAULT PLOTTING DISTANCES, IF NECESSARY\r
2197 \r
2198  &n