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

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

index d133cc7..b46bbe1 100644 (file)
@@ -86,6 +86,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     String DEFN[]= new String[2*MNARC];\r
     private boolean traditionalInput = false;\r
     Scanner input = new Scanner(System.in);\r
+    private double zzset[][] = new double[400][2];\r
     \r
        public SymmsIntegralMapping() {\r
                \r
@@ -222,11 +223,10 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
              int IER[] = new int[1];\r
              int ORDRG[] = new int[1];\r
              int ORDSG[] = new int[1];\r
-             double ALPHA,PI,R1MACH,X,Y;\r
+             double ALPHA,PI,X,Y;\r
              //COMPLEX CENSY,RTUNI,U2\r
              double RTUNI[] = new double[2];\r
              double U2[] = new double[2];\r
-             boolean CHRIN,CHCK;\r
             // CHARACTER TXT*72,TABC*6,FORTFL*72,CH*2,SIG(10)*2,WID(10)*2,REDD*6,\r
             //+FMT1*8,FMT2*9\r
              String TXT;\r
@@ -234,13 +234,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
              String SIG[] = new String[]{"7","8","9","10","11","12","13","14","15","16"};\r
              String WID[] = new String[]{"15","16","17","18","19","20","21","22","23","24"};\r
              String REDD;\r
-             String FMT1;\r
-             String FMT2;\r
+             //String FMT1;\r
+             //String FMT2;\r
        \r
              //PARAMETER (MNARC=100,TABC='     +',CHNL=20,CHIN=21)\r
              final String TABC = "     +";\r
              final int CHNL = 20;\r
-             final int CHIN = 21;\r
+             //final int CHIN = 21;\r
              \r
              File file;\r
              RandomAccessFile raFile = null;\r
@@ -271,8 +271,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                  //**** POINT OUTPUT\r
                \r
                  REDD="E"+WID[SW-1]+"."+SIG[SW-1]; \r
-                 FMT1="("+REDD+")";\r
-                 FMT2="(2"+REDD+")";\r
+                 //FMT1="("+REDD+")";\r
+                 //FMT2="(2"+REDD+")";\r
                  \r
                  if (traditionalInput) {\r
                      System.out.println("ENTER FILENAME TO RECEIVE OUTPUT JAVA CODE");  \r
@@ -549,7 +549,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                  System.exit(-1);\r
              }\r
              HEADER("DPARFN",REDD,raFile);\r
-             /*if (SYMTY) {\r
+             if (SYMTY) {\r
                WRSYM1(NARCS,ORDRG[0],ORDSG[0],RTUNI,U2,CENSY,REFLN,false,REDD,\r
                       CHNL, raFile);\r
                if (REFLN) {\r
@@ -560,20 +560,24 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                }\r
                WRFUN2(NARCS,MNARC,STAPT,ARCTY,PGM,RGM,PTX,NTX,DEFN,\r
                           CHNL,"IB",CH,"ZETA  ",NUMDER,REDD, raFile);\r
-               CALL WRSYM3(NARCS,ORDRG,REFLN,CHNL)\r
+               WRSYM3(NARCS,ORDRG[0],REFLN,CHNL, raFile);\r
              } // if (SYMTY)\r
              else {\r
-               CALL WRFUN2(NARCS,MNARC,STAPT,ARCTY,PGM,RGM,PTX,NTX,DEFN,\r
-            +              CHNL,'IA','TT','DPARFN',NUMDER,REDD)\r
+                 WRFUN2(NARCS,MNARC,STAPT,ARCTY,PGM,RGM,PTX,NTX,DEFN,\r
+                        CHNL,"IA","TT","DPARFN",NUMDER,REDD, raFile);\r
+             }\r
+       \r
+             try {\r
+                     raFile.writeBytes("//\n");\r
+                     raFile.writeBytes("}\n");\r
+                     raFile.close();\r
+             }\r
+             catch (IOException e) {\r
+                 MipavUtil.displayError("IOException " + e + " in PARGEN");\r
+                 System.exit(-1);\r
              }\r
-       C\r
-             WRITE(CHNL,'(A1)') 'C'\r
-             WRITE(CHNL,'(A9)') '      END'\r
-             CLOSE(CHNL)\r
-       C\r
-       999   CONTINUE\r
-             CALL WRTAIL(6,0,IER)\r
-       C*/\r
+             WRTAIL(6,0,IER[0], null);\r
+       \r
        } // public void PARGEN\r
     private void HEADER(String TXT, String REDD, RandomAccessFile raFile) {\r
 \r
@@ -993,11 +997,11 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                //COMPLEX ZT\r
                double ZT[] = new double[2];\r
                boolean NEEDC;\r
-               String FMT;\r
+               //String FMT;\r
                double cr[] = new double[1];\r
                double ci[] = new double[1];\r
                \r
-               FMT="(A12,"+REDD+",A1,"+REDD+",A2)";\r
+               //FMT="(A12,"+REDD+",A1,"+REDD+",A2)";\r
                \r
                try {\r
                        if (PARFUN) {\r
@@ -1038,7 +1042,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                            raFile.writeBytes("//\n");\r
                        }\r
                        \r
-                       FMT="(A7,"+REDD+",A1,"+REDD+",A2)";\r
+                       //FMT="(A7,"+REDD+",A1,"+REDD+",A2)";\r
                        \r
                    if (ORDRG >= 2) {\r
                        raFile.writeBytes("double WW[] = new double["+(ORDRG-1)+"];\n");\r
@@ -1151,12 +1155,13 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
           int IA,I,J,K;\r
          //CHARACTER TX1*16,TX2*21,FMT1*11,FMT2*11\r
-          String TX1,TX2,FMT1,FMT2;\r
+          String TX1,TX2;\r
+          //String FMT1,FMT2;\r
          // EXTERNAL PTFUN1\r
          TX1 = "     if("+CHIA+ " == ";\r
          TX2 = "     else if("+CHIA+ " == ";\r
-          FMT1="(A16,I3,A6)";\r
-          FMT2="(A21,I3,A6)";\r
+          //FMT1="(A16,I3,A6)";\r
+          //FMT2="(A21,I3,A6)";\r
           double STAPT2[][];\r
           double RGM2[];\r
           String DEFN2[];\r
@@ -1218,20 +1223,21 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           //**** LOCAL VARIABLES\r
 \r
       int I;\r
-      double HA,MD,RAD,A,R;\r
+      double HA,MD,RAD;\r
       double C1[] = new double[2];\r
       double C2[] = new double[2];\r
       double CENTR[] = new double[2];\r
       //COMPLEX C1,C2,CENTR\r
-      String TX1, TX1B, TX2, TX2B, CTX1B, FMT1, FMT2, FMT3, FMT4, FMT5;\r
+      //String TX1, TX1B, TX2;\r
+      //String TX2B, CTX1B, FMT1, FMT2, FMT3, FMT4, FMT5;\r
       //CHARACTER TX1*4,TX1B*5,TX2*13,TX2B*14,CTX1B*10,\r
       //+FMT1*25,FMT2*25,FMT3*14,FMT4*25,FMT5*24\r
 \r
-      TX1 = "+"+CHTT+"*";\r
-      TX1B = TX1 + "(";\r
-      CTX1B="     " + TX1B;\r
-      TX2="      "+VAR+" = ";\r
-      TX2B=TX2+"(";\r
+      //TX1 = "+"+CHTT+"*";\r
+      //TX1B = TX1 + "(";\r
+      //CTX1B="     " + TX1B;\r
+      //TX2="      "+VAR+" = ";\r
+      //TX2B=TX2+"(";\r
 \r
       //FMT1='(A14,'//REDD//',A1,'//REDD//',A2)'\r
       //FMT2='(A10,'//REDD//',A1,'//REDD//',A1)'\r
@@ -1482,7 +1488,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
          //     LOCAL VARIABLES\r
                \r
          int IA,I,J1,J2,N1,N2;\r
-         String TX1, TX2, FMT1, FMT2;\r
+         String TX1, TX2;\r
+         //String FMT1, FMT2;\r
          //CHARACTER TX1*16,TX2*21,FMT1*11,FMT2*11\r
          //EXTERNAL PTFUN2\r
          double STAPT2[][];\r
@@ -1493,8 +1500,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
           TX1="      if ("+CHIA+" == ";\r
          TX2="      else if ("+CHIA+" == ";\r
-         FMT1="(A16,I3,A6)";\r
-         FMT2="(A21,I3,A6)";\r
+         //FMT1="(A16,I3,A6)";\r
+         //FMT2="(A21,I3,A6)";\r
                \r
          try {\r
          for (IA=1; IA <= NARCS; IA++) {\r
@@ -1554,81 +1561,624 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
           RandomAccessFile raFile) {\r
           //COMPLEX STAPT(*)\r
           //CHARACTER TXT1(*)*72,TXT2(*)*72,CHTT*2,VAR*6,CHIA*2,REDD*6\r
-               /*C\r
-               C.......................................................................\r
-               C     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
-               C     LAST UPDATE: 8 AUG 1990\r
-               C.......................................................................C\r
-               C**** LOCAL VARIABLES\r
-               C\r
-                     INTEGER I\r
-                     REAL HA,MD,RAD,A,R\r
-                     COMPLEX C1,CENTR\r
-                     CHARACTER TX1*4,TX1B*5,TX2*13,TX2B*14,TX3*39,\r
-                    +FMT1*25,FMT2*15,FMT3*15,FMT4*25,FMT5*24\r
-               C\r
-                     TX1='+'//CHTT//'*'\r
-                     TX1B=TX1//'('\r
-                     TX2='      '//VAR//'='\r
-                     TX2B=TX2//'('\r
-                     TX3=TX2//'(ZDER+UI*ZRAD)*EXP(UI*T)*('\r
-               C\r
-                     FMT1='(A14,'//REDD//',A1,'//REDD//',A2)'\r
-                     FMT2='(A39,'//REDD//',A1)'\r
-                     FMT3='(A13,'//REDD//',A2)'\r
-                     FMT4='(A14,'//REDD//',A5,'//REDD//',A3)'\r
-                     FMT5='(A8,'//REDD//',A5,'//REDD//',A1)'\r
-               C\r
-                     IF (TYPE.EQ.1) THEN\r
-                       C1=5E-1*(STAPT(2)-STAPT(1))\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                       R=REAL(C1)\r
-                       A=AIMAG(C1)\r
-                       WRITE(CHNL,FMT1) TX2B,R,',',A,') '\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                     ELSE IF (TYPE.EQ.2) THEN\r
-                       CENTR=CMPLX(RGM(1),RGM(2))\r
-                       C1=STAPT(1)-CENTR\r
-                       HA=5E-1*RGM(3)\r
-                       MD=ATAN2(AIMAG(C1),REAL(C1))+HA\r
-                       RAD=ABS(C1)\r
-                       WRITE(CHNL,'(A1)') 'C' \r
-                       WRITE(CHNL,FMT4) TX2,RAD,'*UI*(',HA,') *'\r
-                       WRITE(CHNL,FMT4) '     +EXP(UI*(',MD,TX1B,HA,')))'\r
-                       WRITE(CHNL,'(A1)') 'C'  \r
-                     ELSE IF (NUMDER) THEN\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                       WRITE(CHNL,'(A26)') TX2,'ZDPARF(',CHIA,',',CHTT,')'     \r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                     ELSE IF (TYPE.EQ.3) THEN\r
-                       MD=5E-1*(RGM(2)+RGM(1))\r
-                       HA=5E-1*(RGM(2)-RGM(1))\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                       WRITE(CHNL,FMT5) '      T=',MD,TX1B,HA,')'\r
-                       WRITE(CHNL,FMT3) TX2,HA,'*('\r
-                       DO 10 I=1,NTX2\r
-                         WRITE(CHNL,'(A72)') TXT2(I)\r
-               10      CONTINUE\r
-                       WRITE(CHNL,'(A7)') '     +)'\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                     ELSE\r
-                       MD=5E-1*(RGM(2)+RGM(1))\r
-                       HA=5E-1*(RGM(2)-RGM(1))\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                       WRITE(CHNL,FMT5) '      T=',MD,TX1B,HA,')'\r
-                       WRITE(CHNL,'(A11)') '      ZRAD='\r
-                       DO 20 I=1,NTX1\r
-                         WRITE(CHNL,'(A72)') TXT1(I)\r
-               20      CONTINUE\r
-                       WRITE(CHNL,'(A11)') '      ZDER='\r
-                       DO 30 I=1,NTX2\r
-                         WRITE(CHNL,'(A72)') TXT2(I)\r
-               30      CONTINUE\r
-                       WRITE(CHNL,FMT2) TX3,HA,')'\r
-                       WRITE(CHNL,'(A1)') 'C'\r
-                     ENDIF\r
-               C*/\r
-               }\r
+               \r
+         //.......................................................................\r
+         //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+         //     LAST UPDATE: 8 AUG 1990\r
+         //.......................................................................C\r
+         // **** LOCAL VARIABLES\r
+               \r
+         int I;\r
+         double HA,MD,RAD,A,R;\r
+         double C1[] = new double[2];\r
+         double CENTR[] = new double[2];\r
+         //COMPLEX C1,CENTR\r
+         //String TX1, TX1B, TX2, TX2B, TX3, FMT1, FMT2, FMT3, FMT4, FMT5;\r
+         //CHARACTER TX1*4,TX1B*5,TX2*13,TX2B*14,TX3*39,\r
+          //FMT1*25,FMT2*15,FMT3*15,FMT4*25,FMT5*24\r
+         //String TX3R, TX3I;\r
+               \r
+         //TX1 = "+"+CHTT+"*";\r
+         //TX1B = TX1 + "(";\r
+         //TX2 = "      "+VAR+" = ";\r
+         //TX2B = TX2 + "(";\r
+         //TX3R = "      "+VAR+"[0] = (ZDER * Math.cos(T) - ZRAD * Math.sin(T))*(";\r
+         //TX3I = "      "+VAR+"[1] = (ZDER * Math.sin(T) + ZRAD * Math.cos(T))*(";\r
+          //TX3=TX2//'(ZDER+UI*ZRAD)*EXP(UI*T)*('\r
+        \r
+         //FMT1='(A14,'//REDD//',A1,'//REDD//',A2)'\r
+         //FMT2='(A39,'//REDD//',A1)'\r
+         //FMT3='(A13,'//REDD//',A2)'\r
+         //FMT4='(A14,'//REDD//',A5,'//REDD//',A3)'\r
+         //FMT5='(A8,'//REDD//',A5,'//REDD//',A1)'\r
+               \r
+            try {\r
+           if (TYPE == 1) {\r
+                        C1[0] = 0.5*(STAPT[1][0] - STAPT[0][0]);\r
+                        C1[1] = 0.5*(STAPT[1][1] - STAPT[0][1]);\r
+                        raFile.writeBytes("//\n");\r
+                        R=C1[0];\r
+                        A=C1[1];\r
+                        raFile.writeBytes("      "+VAR+"[0] = " + R + ";\n");\r
+                        raFile.writeBytes("      "+VAR+"[1] = " + A + ";\n");\r
+                        raFile.writeBytes("//\n");\r
+                } // if (TYPE == 1)\r
+                else if (TYPE == 2) {\r
+                        CENTR[0] = RGM[0];\r
+                        CENTR[1] = RGM[1];\r
+                        C1[0] = STAPT[0][0] - CENTR[0];\r
+                        C1[1] = STAPT[0][1] - CENTR[1];\r
+                        HA = 0.5 * RGM[2];\r
+                        MD = Math.atan2(C1[1], C1[0]) + HA;\r
+                        RAD = zabs(C1[0],C1[1]);\r
+                        raFile.writeBytes("//\n");\r
+                        raFile.writeBytes(VAR+"[0] = -" + RAD+"*"+HA+"*Math.sin("+MD+"+"+CHTT+"*"+HA+");\n");\r
+                        raFile.writeBytes(VAR+"[1] = " + RAD+"*"+HA+"*Math.cos("+MD+"+"+CHTT+"*"+HA+");\n");\r
+                        raFile.writeBytes("//\n");\r
+                } // else if (TYPE == 2)\r
+                else if (NUMDER) {\r
+                        raFile.writeBytes("//\n");\r
+                        raFile.writeBytes("      "+VAR+" = ZDPARF("+CHIA+","+CHTT+");\n");\r
+                        raFile.writeBytes("//\n");  \r
+                } // else if (NUMDER)\r
+                else if (TYPE == 3) {\r
+                        MD = 0.5*(RGM[1] + RGM[0]);\r
+                        HA = 0.5*(RGM[1] - RGM[0]);\r
+                        raFile.writeBytes("//\n");  \r
+                        raFile.writeBytes("      T = " + MD+"+"+CHTT+"*"+"("+HA+");\n");\r
+                        raFile.writeBytes("      "+VAR+" = " +HA + "*(");\r
+                        for (I = 1; I <= NTX2; I++) {\r
+                                raFile.writeBytes(TXT2[I-1]);\r
+                        }\r
+                        raFile.writeBytes(");\n");\r
+                        raFile.writeBytes("//\n");      \r
+                } // else if (TYPE == 3)\r
+                else {\r
+                        MD = 0.5*(RGM[1] + RGM[0]);\r
+                        HA = 0.5*(RGM[1] - RGM[0]);\r
+                        raFile.writeBytes("//\n");  \r
+                        raFile.writeBytes("      T = " + MD+"+"+CHTT+"*"+"("+HA+");\n");\r
+                        raFile.writeBytes("      ZRAD = ");\r
+                        for (I = 1; I <= NTX1; I++) {\r
+                                raFile.writeBytes(TXT1[I-1]);\r
+                        }\r
+                        raFile.writeBytes(";\n");\r
+                        raFile.writeBytes("      ZDEr = ");\r
+                        for (I = 1; I <= NTX2; I++) {\r
+                                raFile.writeBytes(TXT2[I-1]);\r
+                        }\r
+                        raFile.writeBytes(";\n");      \r
+                        raFile.writeBytes("      "+VAR+"[0] = (ZDER * Math.cos(T) - ZRAD * Math.sin(T))*("+HA+");\n");  \r
+                        raFile.writeBytes("      "+VAR+"[1] = (ZDER * Math.sin(T) + ZRAD * Math.cos(T))*("+HA+");\n");\r
+                        raFile.writeBytes("//\n");      \r
+                } // else\r
+            } // try\r
+            catch(IOException e) {\r
+                MipavUtil.displayError("IOException " + e + " in PTFUN2");\r
+                System.exit(-1);\r
+            }\r
+               \r
+    }\r
+      \r
+   private void WRSYM3(int NARCS,int ORDRG, boolean REFLN, int CHN, RandomAccessFile raFile) {\r
+\r
+    //**** TO WRITE THE CODE TO RECOVER THE DERIVATIVE FROM ITS SYMMETRIC\r
+    //**** COUNTERPART ON THE FUNDAMENTAL BOUNDARY SECTION.\r
+\r
+    //.......................................................................\r
+    // AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+    // LAST UPDATE: 4 AUG 1990\r
+    //.......................................................................C\r
+    //     LOCAL VARIABLES\r
+\r
+    int I;\r
+\r
+\r
+    try {\r
+           raFile.writeBytes("//\n"); \r
+       \r
+           if (REFLN) {\r
+               if (ORDRG > 1) {\r
+                   I=2*NARCS;\r
+                   if (NARCS > 1) {\r
+                       raFile.writeBytes("      IS = (IA-IB)%"+I+";\n");\r
+                       raFile.writeBytes("      IR=(IA-IB-IS)/"+I+";\n");\r
+                   }\r
+                   else {\r
+                       raFile.writeBytes("      IS=(IA-1)%2;\n");\r
+                       raFile.writeBytes("      IR=(IA-1-IS)/2;\n");\r
+                   }\r
+                   raFile.writeBytes("      if ((IR == 0) && (IS == 0)) {\n");\r
+                   raFile.writeBytes("          DPARFN[0] =ZETA[0];\n");\r
+                   raFile.writeBytes("          DPARFN[1] =ZETA[1];\n");\r
+                   raFile.writeBytes("      }\n");\r
+                   raFile.writeBytes("      else if ((IR > 0) && (IS == 0)) {\n");\r
+                   raFile.writeBytes("          DPARFN[0] = WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\n");\r
+                   raFile.writeBytes("          DPARFN[1] = WW[IR-1][0]*ZETA[1] + WW[IR-1][1]*ZETA[0];\n");\r
+                   raFile.writeBytes("      }\n");\r
+                   raFile.writeBytes("       else if ((IR == 0 && (IS > 0)) {\n");\r
+                   raFile.writeBytes("           DPARFN[0] = -U2[0]*ZETA[0] -U2[1]*ZETA[1];\n");\r
+                   raFile.writeBytes("           DPARFN[1] = U2[0]*ZETA[1] -U2[1]*ZETA[0];\n");\r
+                   raFile.writeBytes("       }\n");\r
+                   raFile.writeBytes("       else {\n");\r
+                   raFile.writeBytes("           double realPart = -U2[0]*WW[IR-1][0] + U2[1]*WW[IR-1][1]);\n");\r
+                   raFile.writeBytes("           double imagPart = -U2[0]*WW[IR-1][1] + U2[1]*WW[IR-1][0]);\n");\r
+                   raFile.writeBytes("           DPARFN[0] = realPart*ZETA[0] + imagPart*ZETA[1]);\n");\r
+                   raFile.writeBytes("           DPARFN[1] = imagPart*ZETA[0] - realPart*ZETA[1]);\n"); \r
+                   raFile.writeBytes("       }\n");\r
+               } // if (ORDRG > 1)\r
+               else {\r
+                   if (NARCS > 1) {\r
+                       raFile.writeBytes("      IS=IA-IB;\n");\r
+                   }\r
+                   else {\r
+                       raFile.writeBytes("      IS=IA-1;\n");\r
+                   }\r
+                   raFile.writeBytes("      if (IS == 0) {\n");\r
+                   raFile.writeBytes("          DPARFN[0] = ZETA[0]\n");\r
+                   raFile.writeBytes("          DPARFN[0] = ZETA[0]\n");\r
+                   raFile.writeBytes("      }\n");\r
+                   raFile.writeBytes("         else {\n");\r
+                   raFile.writeBytes("             DPARFN[0] = -U2[0]*ZETA[0] + U2[1]*ZETA[1]);\n");\r
+                   raFile.writeBytes("             DPARFN[1] = U2[0]*ZETA[1] - U2[1]*ZETA[0]);\n");\r
+                   raFile.writeBytes("         }\n");\r
+               } // else\r
+           } // if (REFLN)\r
+           else {\r
+               if (NARCS > 1) {\r
+                   raFile.writeBytes("      IR=(IA-IB)/"+NARCS+";\n");\r
+               }\r
+               else {\r
+                   raFile.writeBytes("      IR=IA-1;\n");\r
+               }\r
+               raFile.writeBytes("      if (IR == 0) {\n");\r
+               raFile.writeBytes("          DPARFN[0] = ZETA[0];\n");\r
+               raFile.writeBytes("          DPARFN[1] = ZETA[1];\n");\r
+               raFile.writeBytes("      }\n");\r
+               raFile.writeBytes("      else {\n");\r
+               raFile.writeBytes("          DPARFN[0]= WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\n");\r
+               raFile.writeBytes("          DPARFN[1]= WW[IR-1][0]*ZETA[1] + WW[IR-1][1]*ZETA[0];\n");\r
+               raFile.writeBytes("      }\n");\r
+           }\r
+    } // try\r
+    catch(IOException e) {\r
+       MipavUtil.displayError("IOException " + e + " in WRSYM3");\r
+       System.exit(-1);\r
+    }\r
+\r
+    } // private void WRSYM3\r
+   \r
+   private void TSTPLT(String JBNM, double MXMIS, double MXDIF, int NARCS, double PSD,\r
+                  double MINPD, double MAXPD, int CHNL, int IER[]) {\r
+   //CHARACTER*4 JBNM\r
+\r
+    //......................................................................\r
+\r
+    // 1.     TSTPLT\r
+    //           TESTS THE PARAMETRIC FUNCTION ROUTINES PARFUN AND DPARFN\r
+    //           FOR CONSISTENCY AND OUTPUTS BOUNDARY POINTS FOR PLOTTING.\r
+           \r
+\r
+    // 2.     PURPOSE\r
+    //           THE ROUTINE FIRST CHECKS THAT THE PARAMETRIC FUNCTION\r
+    //           ROUTINE PARFUN IS CONSISTENT WITH RESPECT TO ITS DEFINITION\r
+    //           OF ANY CORNERS ON THE BOUNDARY.  THIS IS DONE BY CHECKING \r
+    //           THAT THE COMPUTED POINT AT THE END OF EACH ARC MATCHES THE \r
+    //           COMPUTED POINT AT THE START OF THE NEXT ONE.  IF ALL THE\r
+    //           RELATIVE MISFIT ERRORS AT CORNERS ARE LESS THAN\r
+    //           10*(UNIT ROUNDOFF) THEN ALL CORNERS ARE CONSIDERED TO FIT\r
+    //           SATISFACTORILY, OTHERWISE THE MAXIMUM RELATIVE MISFIT\r
+    //           ERROR IS REPORTED.\r
+\r
+    //           THE SECOND PURPOSE OF THE ROUTINE IS TO OUTPUT TO A DATA \r
+    //           FILE THE COORDINATES OF A NUMBER OF POINTS ON THE BOUNDARY. \r
+    //           THE BOUNDARY POINTS ARE SELECTED ADAPTIVELY TO MEET THE \r
+    //           PLOTTING RESOLUTION SPECIFICATIONS DEFINED BY THE ARGUMENTS \r
+    //           PSD,MINPD,MAXPD (SEE BELOW).  THE HOPE IS THAT THE USER MAY \r
+    //           EASILY FEED THESE DATA POINTS TO HIS LOCAL GRAPH PLOTTING  \r
+    //           ROUTINES SO AS TO CONSTRUCT A PLOT OF THE BOUNDARY.  THIS   \r
+    //           WILL PROVIDE AN ESSENTIAL VISUAL CHECK ON THE VALIDITY OF \r
+    //           THE ROUTINE PARFUN.  THE OUTPUT DATA FILE IS AUTOMATICALLY \r
+    //           NAMED <JBNM>zz.\r
\r
+    //           THE THIRD PURPOSE OF THE ROUTINE IS TO CHECK PARFUN AND \r
+    //           DPARFN FOR MUTUAL CONSISTECY.  THIS IS DONE BY COMPUTING \r
+    //           TWO POINT FINITE DIFFERENCE APPROXIMATIONS TO DPARFN.  \r
+    //           THESE DIFFERENCE APPROXIMATIONS ARE COMPUTED AT EACH BOUND- \r
+    //           ARY POINT THAT IS OUTPUT FOR PLOTTING AND ALSO AT NEARBY \r
+    //           POINTS WHICH LIE JUST O F F THE BOUNDARY.  THIS LATTER \r
+    //           COMPARISON ALSO TESTS PARFUN AND DPARFN FOR CORRECTNESS IN \r
+    //           ACCEPTING COMPLEX PARAMETER VALUES.   A RELATIVE ERROR IN \r
+    //           THE FINITE DIFFERENCE APPROXIMATION GREATER THAN 0.1 IS \r
+    //           REPORTED AS A POSSIBLE LOGICAL INCONSISTENCY BETWEEN PARFUN \r
+    //           AND DPARFN.  (THE CRITICAL RELATIVE ERROR VALUE OF 0.1 CAN \r
+    //           BE ALTERED BY ADJUSTING THE LOCAL PARAMETER DTOL).\r
+\r
+    // 3.     CALLING SEQUENCE\r
+    //           CALL TSTPLT(JBNM,MXMIS,MXDIF,NARCS,PSD,MINPD,MAXPD,CHNL,IER)\r
+\r
+    //        PARAMETERS\r
+    //         ON ENTRY\r
+    //            JBNM   - CHARACTER*4\r
+    //                     THE JOB NAME.  THIS IS USED TO CREATE THE OUTPUT \r
+    //                     FILE WITH FILENAME\r
+\r
+    //                                <JBNM>zz ,\r
+\r
+    //                     WHERE <JBNM> DENOTES THE VALUE OF VARIABLE JBNM\r
+    //                     WITH ANY TRAILING SPACES DELETED.\r
+\r
+    //            NARCS  - INTEGER\r
+    //                     THE NUMBER OF ANALYTIC ARCS THAT MAKE UP THE \r
+    //                     W H O L E BOUNDARY OF THE PHYSICAL DOMAIN.\r
+\r
+    //            PSD    - REAL\r
+    //                     THE PLOTTING SIZE FOR THE DOMAIN IN ANY APPROPR-\r
+    //                     IATE UNITS.  IF PSD .LE. 0.0 THEN IT IS ASSIGNED\r
+    //                     THE DEFAULT VALUE OF 160.0 (A REASONBLE WIDTH IN\r
+    //                     MM FOR PLOTTING ON A4 PAPER).\r
+\r
+    //            MINPD  - REAL\r
+    //                     THE MINIMUM SIGNIFICANT PLOTTING DISTANCE, IN THE\r
+    //                     SAME UNITS AS PSD.  IF PSD .LE. 0.0 THEN IT IS\r
+    //                     ASSIGNED THE DEFAULT VALUE OF 2.0.\r
+\r
+    //            MAXPD  - REAL\r
+    //                     THE MAXIMUM ALLOWED PLOTTING DISTANCE, IN THE\r
+    //                     SAME UNITS AS PSD.  IF PSD .LE. 0.0 THEN IT IS\r
+    //                     ASSIGNED THE DEFAULT VALUE OF 5.0.  THE LARGER\r
+    //                     MAXPD, THE COARSER WILL BE THE RESOLUTION OF THE\r
+    //                     BOUNDARY POINTS OUTPUT TO <JBNM>zz, BUT THE\r
+    //                     QUICKER THEY WILL BE COMPUTED. \r
+\r
+    //            CHNL   - INTEGER\r
+    //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+    //                     WRITING THE FILE <JBNM>zz.\r
+    //         ON EXIT\r
+    //            MXMIS  - REAL\r
+    //                     THE MAXIMUM RELATIVE MISFIT ERROR OVER ALL\r
+    //                     CORNER POINTS\r
+\r
+    //            MXDIF  - REAL\r
+    //                     THE MAXIMUM RELATIVE ERROR IN FINITE DIFFERENCE\r
+    //                     APPROXIMATIONS TO DPARFN OVER ALL BOUNDARY \r
+    //                     POINTS OUTPUT TO <JBNM>zz AND NEARBY POINTS OFF\r
+    //                     THE BOUNDARY.\r
+\r
+    //            PSD    - REAL\r
+    //                     IF PSD .LE. 0.0 ON ENTRY THEN IT WILL HAVE\r
+    //                     THE DEFAULT VALUE OF 160.0 ON  EXIT.\r
+\r
+    //            MINPD  - REAL\r
+    //                     IF PSD .LE. 0.0 ON ENTRY THEN MINPD WILL HAVE\r
+    //                     THE DEFAULT VALUE OF 2.0 ON EXIT\r
+\r
+    //            MAXPD  - REAL\r
+    //                     IF PSD .LE. 0.0 ON ENTRY THEN MAXPD WILL HAVE\r
+    //                     THE DEFAULT VALUE OF 5.0 ON EXIT\r
+\r
+    //            IER    - INTEGER\r
+    //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+    //                     A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
+    //                     WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+    //                     IER=0 - NORMAL EXIT.\r
+    //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+    //                             BE SELF EXPLANATORY.\r
+\r
+\r
+    // 4.     SUBROUTINES OR FUNCTIONS NEEDED\r
+    //              - THE CONFPACK LIBRARY.\r
+    //              - THE REAL FUNCTION R1MACH.\r
+    //              - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN.\r
+\r
+\r
+    // 5.     FURTHER COMMENTS\r
+    //              - A SUMMARY LISTING IS AUTOMATICALLY WRITTEN ON THE \r
+    //                STANDARD OUTPUT CHANNEL.\r
+    //              - THE OUTPUT FILE <JBNM>zz CONTAINS COORDINATE PAIRS\r
+\r
+    //                                 X Y\r
+\r
+    //                FOR POINTS ON THE PHYSICAL BOUNDARY, WITH ONE PAIR\r
+    //                PER LINE.\r
+\r
+    // ......................................................................\r
+    //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+    //     LAST UPDATE: 6 JULY 1990\r
+    // ......................................................................C\r
+    //     LOCAL VARIABLES\r
+\r
+    int I,IA,L;\r
+    int IMX = 0;\r
+    double A1,DIFF,ERR,HH,MINC,R1MACH,RMAX,RMEAN,RMIN,T,TINC,TOL1,TMX,TSD;\r
+    double TT[] = new double[2];\r
+    //REAL TT(2)\r
+    double C1[] = new double[2];\r
+    double C2[] = new double[2];\r
+    double CENTR[] = new double[2];\r
+    double ZZ0[] = new double[2];\r
+    double DZZ[] = new double[2];\r
+    double NDZZ[] = new double[2];\r
+    //COMPLEX C1,C2,CENTR,ZZ0,DZZ,NDZZ;\r
+    double ZZ[][] = new double[2][2];\r
+    //COMPLEX ZZ(2)\r
+    String OFL;\r
+    //CHARACTER OFL*6\r
+    final int MNARC = 200;\r
+    final double DTOL = 1.0E-1;\r
+    final int NH = 4;\r
+    boolean ATEND,FIRST,WARND;\r
+    boolean LNSEG[] = new boolean[MNARC];\r
+    double PIN[] = new double[2];\r
+    double PAROUT[];\r
+    int zzindex;\r
+    //EXTERNAL DPARFN,LINSEG,PARFUN,R1MACH,WRHEAD,WRTAIL,ZDPARF\r
+    //COMPLEX PARFUN,DPARFN,ZDPARF\r
+\r
+    //**** WRITE CONFPACK HEADING\r
+\r
+    WRHEAD(7,0,null);\r
+\r
+    if (NARCS > MNARC) {\r
+        IER[0]=59;\r
+        WRTAIL(7,0,IER[0], null);\r
+    }\r
+\r
+    //1     FORMAT(A45)\r
+    //2     FORMAT(A45,I4)\r
+    //3     FORMAT(A45,E10.3)\r
+    //4     FORMAT(//,T17,A)\r
+\r
+    TOL1= 10.0*EPS;\r
+\r
+    //**** CHECK THAT ALL ARCS MEET AT CORNER POINTS\r
+\r
+    IER[0]=0;\r
+    CENTR[0] = 0.0;\r
+    CENTR[1] = 0.0;\r
+    MXMIS=0.0;\r
+    for (IA=1; IA <= NARCS; IA++) {\r
+        if (IA == 1) {\r
+            I=NARCS;\r
+        }\r
+        else {\r
+            I=IA-1;\r
+        }\r
+        PIN[0] = -1.0;\r
+        PIN[1] = 0.0;\r
+        C1 = PARFUN(IA, PIN);\r
+        CENTR[0] = CENTR[0] + C1[0];\r
+        CENTR[1] = CENTR[1] + C1[1];\r
+        A1 = zabs(C1[0],C1[1]);\r
+        PIN[0] = 1.0;\r
+        PIN[1] = 0.0;\r
+        C2=PARFUN(I,PIN);\r
+        ERR = zabs(C1[0]-C2[0],C1[1]-C2[1]);\r
+        if (A1 >= 1.0) {\r
+            ERR=ERR/A1;\r
+        }\r
+        if (ERR > MXMIS) {\r
+            IMX=IA;\r
+            MXMIS=ERR;\r
+        }\r
+    } // for (IA=1; IA <= NARCS; IA++)\r
+       if (MXMIS >= TOL1) {\r
+        System.out.println("MAXIMUM CORNER MISFIT: " + MXMIS);\r
+        System.out.println("OCCURS AT CORNER: " + IMX);\r
+    }\r
+    else {\r
+        System.out.println("ALL ARCS FIT AT CORNERS:");\r
+    }  \r
+\r
+    // **** ESTIMATE THE DIAMETER (TSD) OF THE PHYSICAL DOMAIN\r
+\r
+    CENTR[0]=CENTR[0]/NARCS;\r
+    CENTR[1]=CENTR[1]/NARCS;\r
+    TSD=0.0;\r
+    HH=2.0/(double)(NH);\r
+    for (IA=1; IA <= NARCS; IA++) {\r
+        T=-1.0;\r
+        for (I=1; I <= NH; I++) {\r
+            T=T+HH;\r
+            PIN[0] = T;\r
+            PIN[1] = 0;\r
+            PAROUT = PARFUN(IA,PIN);\r
+            C1[0] = PAROUT[0] - CENTR[0];\r
+            C1[1] = PAROUT[1] - CENTR[1];\r
+            A1=zabs(C1[0],C1[1]);\r
+            TSD=Math.max(TSD,A1);\r
+        } // for (I=1; I <= NH; I++)\r
+    } // for (IA=1; IA <= NARCS; IA++)\r
+   TSD=2.0*TSD;\r
+\r
+   //**** DETERMINE WHICH ARCS (IF ANY) ARE LINE SEGMENTS\r
+\r
+   LINSEG(LNSEG,NARCS);\r
+\r
+   //**** OPEN FILE TO RECEIVE BOUNDARY DATA POINTS FOR PLOTTING\r
+\r
+   //L=INDEX(JBNM,' ')-1\r
+   //IF (L.EQ.-1) L=4\r
+   //OFL=JBNM(1:L)//'zz'\r
+   //OPEN(CHNL,FILE=OFL)\r
+   // Use global zzset instead\r
+\r
+   // **** SET DEFAULT PLOTTING DISTANCES, IF NECESSARY\r
+\r
+   if (PSD <= 0.0) {\r
+       PSD=1.6E+2;\r
+       MINPD=2.0;\r
+       MAXPD=5.0;\r
+   }\r
+   RMIN=MINPD*TSD/PSD;\r
+   RMAX=MAXPD*TSD/PSD;\r
+   RMEAN=0.5*(RMIN+RMAX);\r
+   MINC=Math.sqrt(EPS);\r
+   \r
+   // **** START EVALUATING BOUNDARY POINTS AND DERIVATIVES FOR PLOTTING AND\r
+   // **** TESTING\r
+\r
+   MXDIF=0.0;\r
+   zzindex = 0;\r
+   /*for (IA=1; IA <= NARCS; IA++) {\r
+       TT[0]=-1.0;\r
+       PIN[0] = TT[0];\r
+       PIN[1] = 0.0;\r
+       ZZ[0]=PARFUN(IA,PIN);\r
+       zzset[zzindex][0] = ZZ[0][0];\r
+       zzset[zzindex++][1] = ZZ[0][1];\r
+       if (IA == 1) {\r
+          ZZ0[0] = ZZ[0][0];\r
+          ZZ0[1] = ZZ[0][1];\r
+       }\r
+       FIRST = true;\r
+       WARND= false;\r
+20      CONTINUE\r
+C\r
+C****   TEST THE COMPATIBILTY OF PARFUN AND DPARFN BY ESTIMATING DPARFN\r
+C****   NUMERICALLY AT BOTH REAL AND COMPLEX PARAMETER VALUES.\r
+C\r
+     DO 30 I=1,2\r
+       IF (I.EQ.1) THEN\r
+         C1=CMPLX(TT(1))\r
+       ELSE\r
+         C1=CMPLX(TT(1),MINC)\r
+       ENDIF\r
+       DZZ=DPARFN(IA,C1)\r
+       NDZZ=ZDPARF(IA,C1)\r
+       A1=ABS(DZZ)\r
+C\r
+       IF (A1.EQ.0E+0) THEN\r
+         IER=60\r
+         WRITE(*,*)\r
+         WRITE(*,*) '              ***DPARFN=(0.,0.)***'\r
+         WRITE(*,*) '                             ARC:',IA \r
+         WRITE(*,*) '    STANDARDISED PARAMETER VALUE:',TT(1) \r
+         GOTO 999\r
+       ENDIF\r
+C\r
+       IF (A1.LE.TOL1 .AND. .NOT.WARND) THEN\r
+         WRITE(*,4) '*** W A R N I N G  ***'\r
+         WRITE(*,2) 'PATHOLOGICALLY SMALL DERIVATIVE ON ARC',IA\r
+         WARND=.TRUE.\r
+       ENDIF\r
+C\r
+       IF (FIRST) THEN\r
+         TINC=RMEAN/A1\r
+         TINC=MAX(TINC,MINC)\r
+         FIRST=.FALSE.\r
+       ENDIF\r
+C\r
+       ERR=ABS(1E+0-NDZZ/DZZ)\r
+       IF (ERR.GT.MXDIF) THEN\r
+         MXDIF=ERR\r
+         IMX=IA\r
+         TMX=TT(1)\r
+       ENDIF\r
+30      CONTINUE\r
+C\r
+     IF (.NOT.LNSEG(IA)) THEN\r
+C\r
+C****     DETERMINE THE NEXT BOUNDARY POINT TO BE PLOTTED\r
+C\r
+40        CONTINUE\r
+       TT(2)=TT(1)+TINC\r
+       IF (TT(2) .GE. 1E+0) THEN\r
+         TT(2)=1E+0\r
+         ATEND=.TRUE.\r
+       ELSE\r
+         ATEND=.FALSE.\r
+       ENDIF\r
+C\r
+       ZZ(2)=PARFUN(IA,CMPLX(TT(2)))\r
+       DIFF=ABS(ZZ(2)-ZZ(1))\r
+       IF (DIFF.EQ.0E+0 .AND. .NOT.ATEND) THEN\r
+         TINC=MAX(MINC,2*TINC)\r
+         GOTO 40\r
+       ENDIF\r
+C\r
+       IF (DIFF.GT.RMAX .OR. (DIFF.LT.RMIN .AND. .NOT.ATEND)) THEN\r
+         TINC=RMEAN*TINC/DIFF\r
+         TINC=MAX(TINC,MINC)\r
+         GOTO 40\r
+       ENDIF\r
+C\r
+       WRITE(CHNL,'(2E16.7)') ZZ(2)\r
+       IF (.NOT. ATEND) THEN\r
+         ZZ(1)=ZZ(2)\r
+         TT(1)=TT(2)\r
+         GOTO 20\r
+       ENDIF\r
+     ENDIF\r
+\r
+   } // for (IA=1; IA <= NARCS; IA++)\r
+   IF (LNSEG(NARCS)) WRITE(CHNL,'(2E16.7)') ZZ0\r
+   CLOSE(CHNL)   \r
+C\r
+   IF (MXDIF .GT. DTOL) THEN\r
+     WRITE(*,*)\r
+     WRITE(*,2) 'POSSIBLE PARFUN/DPARFN INCONSISTECY ON ARC:',IMX\r
+     WRITE(*,3) 'OCCURS AT STANDARDISED PARAMETER VALUE:',TMX\r
+     WRITE(*,3) 'RELATIVE FINITE DIFF ERROR:',MXDIF \r
+   ELSE\r
+     WRITE(*,*)\r
+     WRITE(*,1) 'PARFUN AND DPARFN ARE CONSISTENT:'\r
+   ENDIF\r
+C\r
+999   CALL WRTAIL(7,0,IER)*/\r
+\r
+   } // private void TSTPLOT\r
+   \r
+   private void LINSEG(boolean LNSEG[], int NARCS) {\r
+\r
+       //**** TO DETERMINE THE ARRAY LNSEG, WHERE LNSEG(I) IS SET TO TRUE IF THE\r
+       //**** I'TH ARC IS A LINE SEGMENT, I=1,...,NARCS.\r
+\r
+       //**** LOCAL VARIABLES\r
+\r
+       final int NPTS = 9;\r
+          int IA,J,NINTS;\r
+       double DIFF,HH,MXDIF,TOL,R1MACH;\r
+       double SUM[] = new double[2];\r
+       double TT[] = new double[2];\r
+       // COMPLEX SUM,TT\r
+       double DF[][] = new double[NPTS][2];\r
+       //COMPLEX DF(NPTS)\r
+       //EXTERNAL DPARFN,R1MACH\r
+       // COMPLEX DPARFN\r
+\r
+       NINTS=NPTS-1;\r
+       HH=2.0/NINTS;\r
+       TOL=10.0*EPS;\r
+\r
+       for (IA=1; IA <= NARCS; IA++) {\r
+          SUM[0] = 0.0;\r
+          SUM[1] = 0.0;\r
+           for (J=1; J <= NPTS; J++) {\r
+               TT[0]=-1.0+(J-1)*HH;\r
+               TT[1] = 0.0;\r
+               DF[J-1]=DPARFN(IA,TT);\r
+               SUM[0]=SUM[0]+DF[J-1][0];\r
+               SUM[1]=SUM[1]+DF[J-1][1];\r
+           } // for (J=1; J <= NPTS; J++)\r
+           SUM[0]=SUM[0]/NPTS;\r
+           SUM[1]=SUM[1]/NPTS;\r
+\r
+           MXDIF=0.0;\r
+           for (J=1; J <= NPTS; J++) {\r
+               DIFF=zabs(SUM[0]-DF[J-1][0],SUM[1]-DF[J-1][1]);\r
+               MXDIF = Math.max(MXDIF,DIFF);\r
+           } // for (J=1; J <= NPTS; J++)\r
+\r
+           if (MXDIF <= TOL) {\r
+               LNSEG[IA-1] = true;\r
+           }\r
+           else {\r
+               LNSEG[IA-1] = false;\r
+           }\r
+       } // for (IA=1; IA <= NARCS; IA++)\r
+\r
+   } // private void LINSEG\r
 \r
 \r
       \r