Added PARFUN and DPARFN for examples 1 and 3. Corrected PTFUN2.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 7 Feb 2018 23:14:07 +0000 (23:14 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 7 Feb 2018 23:14:07 +0000 (23:14 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15364 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 2fde9c7..c730a36 100644 (file)
@@ -296,6 +296,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        Vector<Double>Ray[]; // x y pairs\r
        // From TSTPLT\r
        Vector<Double>Boundary; // x y pairs\r
+       private int example = 1;\r
 \r
        public SymmsIntegralMapping() {\r
 \r
@@ -858,6 +859,18 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                WRTAIL(6, 0, IER[0], null);\r
 \r
        } // public void PARGEN\r
+       \r
+       public void DRIVE0() {\r
+        int IER[] = new int[1];\r
+        double MXMIS[] = new double[1];\r
+        double MXDIF[] = new double[1];\r
+        double PSD[] = new double[1];\r
+        double MINPD = 0.0;\r
+        double MAXPD = 0.0;\r
+\r
+       TSTPLT(MXMIS,MXDIF,NARCS,PSD,MINPD,MAXPD,IER);\r
+       }\r
+\r
 \r
        private void HEADER(String TXT, String REDD, RandomAccessFile raFile) {\r
 \r
@@ -871,15 +884,11 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                raFile.writeBytes(LINE);\r
                LINE = TAB6 + "double ZRAD[] = new double[2];\n";\r
                raFile.writeBytes(LINE);\r
-               LINE = TAB6 + "int IA;\n";\r
-               raFile.writeBytes(LINE);\r
                LINE = TAB6 + "double T[] = new double[2];\n";\r
                raFile.writeBytes(LINE);\r
-               LINE = TAB6 + "double TT[] = new double[2];\n";\r
+               LINE = TAB6 + "//double TT[] = new double[2];\n";\r
                raFile.writeBytes(LINE);\r
 \r
-                       raFile.writeBytes("      double PI = " + Math.PI + ";\n");\r
-                       raFile.writeBytes("//\n");\r
                } catch (IOException e) {\r
                        MipavUtil.displayError("IOException " + e + " in HEADER");\r
                        System.exit(-1);\r
@@ -961,27 +970,156 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                return;\r
        } // private void SYINF1\r
 \r
-       // COMPLEX FUNCTION PARFUN(I,T)\r
+       // COMPLEX FUNCTION PARFUN(I,TT)\r
        // INTEGER I\r
-       // COMPLEX T\r
+       // COMPLEX TT\r
        // \r
-       double[] PARFUN(int I, double T[]) {\r
-\r
-               // DUMMY FUNCTION TO AID LINK-LOADING OF PARGEN\r
-\r
-               double result[] = new double[] { 1.0, 0.0 };\r
-               return result;\r
-       } // double[] PARFUN\r
-\r
-       // COMPLEX FUNCTION DPARFN(I,T)\r
-       // INTEGER I\r
-       // COMPLEX T\r
-       double[] DPARFN(int I, double T[]) {\r
-\r
-               // DUMMY FUNCTION TO AID LINK-LOADING OF PARGEN\r
-\r
-               double result[] = new double[] { 1.0, 0.0 };\r
-               return result;\r
+       private double[] PARFUN(int IA, double TT[]) {\r
+        double PARFUNANS[] = new double[2];\r
+           double ZRAD[] = new double[2];\r
+           double T[] = new double[2];\r
+           //double TT[] = new double[2];\r
+       if (example == 1) {\r
+          if(IA == 1) {\r
+              PARFUNANS[0] = 2.0+2.0 * Math.exp(-TT[1]*(-0.7853981633974483))*Math.cos(-2.356194490192345+TT[0]*(-0.7853981633974483));\r
+                  PARFUNANS[1] = 2.0+2.0 * Math.exp(-TT[1]*(-0.7853981633974483))*Math.sin(-2.356194490192345+TT[0]*(-0.7853981633974483));\r
+          }\r
+          else if(IA == 2) {\r
+                  T[0] = 2.3561945+TT[0] * (0.7853981999999999);\r
+                  T[1] = TT[1] * (0.7853981999999999);\r
+                  PARFUNANS[0] = Math.cos(T[0])*Math.cosh(T[1])-2.0*Math.cos(T[0])*Math.sinh(T[1]);\r
+                  PARFUNANS[1] = -Math.sin(T[0])*Math.sinh(T[1])+2.0*Math.sin(T[0])*Math.cosh(T[1]);\r
+          }\r
+          else if(IA == 3) {\r
+                  PARFUNANS[0] = -1.5+TT[0]*(-0.5) - TT[1]*(0.0);\r
+                  PARFUNANS[1] = 0.0+TT[0]*(0.0) + TT[1]*(-0.5);\r
+          }\r
+          else if(IA == 4) {\r
+                  PARFUNANS[0] = 0.0+2.0 * Math.exp(-TT[1]*(0.7853981633974483))*Math.cos(3.9269908169872414+TT[0]*(0.7853981633974483));\r
+                  PARFUNANS[1] = 0.0+2.0 * Math.exp(-TT[1]*(0.7853981633974483))*Math.sin(3.9269908169872414+TT[0]*(0.7853981633974483));\r
+          }\r
+          else {\r
+                  T[0] = 5.497787143782138+TT[0] * (0.7853981633974483);\r
+                  T[1] = TT[1] * (0.7853981633974483);\r
+                  ZRAD[0] = 2.0-0.5*Math.sin(4.0*T[0]-6.0*Math.PI)*Math.cosh(4.0*T[1]);\r
+                  ZRAD[1] = -0.5*Math.cos(4.0*T[0]-6.0*Math.PI)*Math.sinh(4.0*T[1]);\r
+                  PARFUNANS[0] = ZRAD[0]*Math.exp(-T[1])*Math.cos(T[0]) - ZRAD[1]*Math.exp(-T[1])*Math.sin(T[0]);\r
+                  PARFUNANS[1] = ZRAD[0]*Math.exp(-T[1])*Math.sin(T[0]) + ZRAD[1]*Math.exp(-T[1])*Math.cos(T[0]);\r
+           }   \r
+       } // if (example == 1)\r
+       else if (example == 3) {\r
+          int IB, IR;\r
+          double ZETA[] = new double[2];\r
+           double WW[][] = new double[3][2];\r
+          WW[0][0] = 6.123233995736766E-17;\r
+          WW[0][1] = 1.0;\r
+          WW[1][0] = -1.0;\r
+          WW[1][1] = 1.2246467991473532E-16;\r
+          WW[2][0] = -1.8369701987210297E-16;\r
+          WW[2][1] = -1.0;\r
+          IB = IA%2;\r
+          if (IB == 0) IB = 2;\r
+          if(IB == 1) {\r
+              ZETA[0] = 1.5+TT[0]*(0.5) - TT[1]*(0.0);\r
+              ZETA[1] = 0.0+TT[0]*(0.0) + TT[1]*(0.5);\r
+          }\r
+          else {\r
+              T[0] = 0.7853981633974483+TT[0] * (0.7853981633974483);\r
+              T[1] = TT[1] * (0.7853981633974483);\r
+              ZRAD[0] = 2.0-2.0*T[0]/Math.PI;\r
+              ZRAD[1] = -2.0*T[1]/Math.PI;\r
+              ZETA[0] = ZRAD[0]*Math.exp(-T[1])*Math.cos(T[0]) - ZRAD[1]*Math.exp(-T[1])*Math.sin(T[0]);\r
+              ZETA[1] = ZRAD[0]*Math.exp(-T[1])*Math.sin(T[0]) + ZRAD[1]*Math.exp(-T[1])*Math.cos(T[0]);\r
+          }\r
+          IR = (IA - IB)/2;\r
+          if (IR == 0) {\r
+              PARFUNANS[0] = ZETA[0];\r
+              PARFUNANS[1] = ZETA[1];\r
+          }\r
+          else {\r
+              PARFUNANS[0] = WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\r
+              PARFUNANS[1] = WW[IR-1][0]*ZETA[1] + WW[IR-1][1]*ZETA[0];\r
+          }   \r
+       } // else if (example == 3)\r
+       return PARFUNANS;\r
+       } // private double[] PARFUN\r
+\r
+       // COMPLEX FUNCTION DPARFN(IA,TT)\r
+       // INTEGER IA\r
+       // COMPLEX TT\r
+       private double[] DPARFN(int IA, double TT[]) {\r
+        double DPARFNANS[] = new double[2];\r
+           double ZRAD[] = new double[2];\r
+           double ZDER[] = new double[2];\r
+           double T[] = new double[2];\r
+           //double TT[] = new double[2];      \r
+           if (example == 1) {\r
+               if (IA == 1) {\r
+                   DPARFNANS[0] = (-2.0)*(-0.7853981633974483)*Math.exp((-TT[1])*(-0.7853981633974483))*Math.sin(-2.356194490192345+TT[0]*(-0.7853981633974483));\r
+                       DPARFNANS[1] = (2.0)*(-0.7853981633974483)*Math.exp((-TT[1])*(-0.7853981633974483))*Math.cos(-2.356194490192345+TT[0]*(-0.7853981633974483));\r
+               }\r
+               else if (IA == 2) {\r
+                   T[0] = 2.3561945+TT[0]*(0.7853981999999999);\r
+                       T[1] = TT[1]*(0.7853981999999999);\r
+                       DPARFNANS[0] = 0.7853981999999999*(-Math.sin(T[0])*Math.cosh(T[1])+2.0*Math.sin(T[0])*Math.sinh(T[1]));\r
+                       DPARFNANS[1] = 0.7853981999999999*(2.0*Math.cos(T[0])*Math.cosh(T[1])-Math.cos(T[0])*Math.sinh(T[1]));\r
+               }\r
+               else if (IA == 3) {\r
+                   DPARFNANS[0] = -0.5;\r
+                       DPARFNANS[1] = 0.0;\r
+               }\r
+               else if (IA == 4) {\r
+                   DPARFNANS[0] = (-2.0)*(0.7853981633974483)*Math.exp((-TT[1])*(0.7853981633974483))*Math.sin(3.9269908169872414+TT[0]*(0.7853981633974483));\r
+                       DPARFNANS[1] = (2.0)*(0.7853981633974483)*Math.exp((-TT[1])*(0.7853981633974483))*Math.cos(3.9269908169872414+TT[0]*(0.7853981633974483));\r
+               }\r
+               else {\r
+                       T[0] = 5.497787143782138+TT[0]*(0.7853981633974483);\r
+                       T[1] = TT[1]*(0.7853981633974483);\r
+                       ZRAD[0] = 2.0-0.5*Math.sin(4.0*T[0]-6.0*Math.PI)*Math.cosh(4.0*T[1]);\r
+                       ZRAD[1] = -0.5*Math.cos(4.0*T[0]-6.0*Math.PI)*Math.sinh(4.0*T[1]);\r
+                       ZDER[0] = -2.0*Math.cos(4.0*T[0]-6.0*Math.PI)*Math.cosh(4.0*T[1]);\r
+                       ZDER[1] = 2.0*Math.sin(4.0*T[0]-6.0*Math.PI)*Math.sinh(4.0*T[1]);\r
+                       DPARFNANS[0] = ((ZDER[0] - ZRAD[1])* Math.cos(T[0]) - (ZRAD[0] + ZDER[1])* Math.sin(T[0]))*Math.exp(-T[1])*(0.7853981633974483);\r
+                       DPARFNANS[1] = ((ZDER[0] - ZRAD[1])* Math.sin(T[0]) + (ZRAD[0] + ZDER[1])* Math.cos(T[0]))*Math.exp(-T[1])*(0.7853981633974483);\r
+                }      \r
+           } // if (example == 1)\r
+           else if (example == 3) {\r
+               int IB, IR;\r
+               double ZETA[] = new double[2];\r
+               double WW[][] = new double[3][2];\r
+               WW[0][0] = 6.123233995736766E-17;\r
+               WW[0][1] = 1.0;\r
+               WW[1][0] = -1.0;\r
+               WW[1][1] = 1.2246467991473532E-16;\r
+               WW[2][0] = -1.8369701987210297E-16;\r
+               WW[2][1] = -1.0;\r
+               IB = IA%2;\r
+               if (IB == 0) IB = 2;\r
+               if (IB == 1) {\r
+                   ZETA[0] = 0.5;\r
+                   ZETA[1] = 0.0;\r
+               }\r
+               else {\r
+                   T[0] = 0.7853981633974483+TT[0]*(0.7853981633974483);\r
+                   T[1] = TT[1]*(0.7853981633974483);\r
+                   ZRAD[0] = 2.0-2.0*T[0]/Math.PI;\r
+                   ZRAD[1] = -2.0*T[1]/Math.PI;\r
+                   ZDER[0] = -2.0/Math.PI;\r
+                   ZDER[1] = 0.0;\r
+                   ZETA[0] = ((ZDER[0] - ZRAD[1])* Math.cos(T[0]) - (ZRAD[0] + ZDER[1])* Math.sin(T[0]))*Math.exp(-T[1])*(0.7853981633974483);\r
+                   ZETA[1] = ((ZDER[0] - ZRAD[1])* Math.sin(T[0]) + (ZRAD[0] + ZDER[1])* Math.cos(T[0]))*Math.exp(-T[1])*(0.7853981633974483);\r
+               }\r
+               IR=(IA-IB)/2;\r
+               if (IR == 0) {\r
+                   DPARFNANS[0] = ZETA[0];\r
+                   DPARFNANS[1] = ZETA[1];\r
+               }\r
+               else {\r
+                   DPARFNANS[0]= WW[IR-1][0]*ZETA[0] - WW[IR-1][1]*ZETA[1];\r
+                   DPARFNANS[1]= WW[IR-1][0]*ZETA[1] + WW[IR-1][1]*ZETA[0];\r
+               }\r
+           } // else if (example == 3)\r
+           return DPARFNANS;\r
        } // double[] DPARFN\r
 \r
        private void WRHEAD(int I, int CHNL, RandomAccessFile raFile) {\r
@@ -1420,8 +1558,9 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                                                raFile.writeBytes(TX2 + IA + ") {\n");\r
                                        }\r
                                        PTFUN1(ARCTY[IA - 1], STAPT2, RGM2, NTX[IA - 1], DEFN2, CHNL, CHTT, VAR, REDD, raFile, writeReturn);\r
+                                       raFile.writeBytes("      }\n");\r
                                        if (IA == NARCS)\r
-                                               raFile.writeBytes("      }\n");\r
+                                               raFile.writeBytes("    }\n");\r
                                } // else\r
                        } // for (IA=1; IA <= NARCS; IA++)\r
                } // try\r
@@ -1877,8 +2016,9 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                                        }\r
                                        PTFUN2(ARCTY[IA - 1], STAPT2, RGM2, N1, DEFN2, N2, DEFN3, CHNL, CHTT, VAR, CHIA, NUMDER[IA - 1],\r
                                                        REDD, raFile, writeReturn);\r
+                                       raFile.writeBytes("      }\n");\r
                                        if (IA == NARCS) {\r
-                                               raFile.writeBytes("      }\n");\r
+                                               raFile.writeBytes("    }\n");\r
                                        }\r
                                } // else\r
                        } // for (IA=1; IA <= NARCS; IA++)\r
@@ -1949,7 +2089,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                                raFile.writeBytes("//\n");\r
                                raFile.writeBytes(\r
                                                PAD + VAR + "[0] = (-" + RAD + ")*(" + HA + ")*Math.exp((-"+CHTT + "[1])*(" + HA +"))*Math.sin(" + MD + "+" + CHTT + "[0]*(" + HA + "));\n");\r
-                               raFile.writeBytes(PAD + VAR + "[1] = (" + RAD + ")*(" + HA + ")*Math.exp((-"+CHTT + "[1])*Math.cos(" + MD + "+" + CHTT + "[0]*(" + HA + "));\n");\r
+                               raFile.writeBytes(PAD + VAR + "[1] = (" + RAD + ")*(" + HA + ")*Math.exp((-"+CHTT + "[1])*(" + HA + "))*Math.cos(" + MD + "+" + CHTT + "[0]*(" + HA + "));\n");\r
                                raFile.writeBytes("//\n");\r
                        } // else if (TYPE == 2)\r
                        else if (NUMDER) {\r
@@ -2279,8 +2419,8 @@ public class SymmsIntegralMapping extends AlgorithmBase {
 \r
        } // private void WRSYM3\r
 \r
-       private void TSTPLT(String JBNM, double MXMIS, double MXDIF, int NARCS, double PSD, double MINPD, double MAXPD,\r
-                       int CHNL, int IER[]) {\r
+       private void TSTPLT(double MXMIS[], double MXDIF[], int NARCS, double PSD[], double MINPD, double MAXPD,\r
+                       int IER[]) {\r
                // CHARACTER*4 JBNM\r
 \r
                // ......................................................................\r
@@ -2467,7 +2607,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                IER[0] = 0;\r
                CENTR[0] = 0.0;\r
                CENTR[1] = 0.0;\r
-               MXMIS = 0.0;\r
+               MXMIS[0] = 0.0;\r
                for (IA = 1; IA <= NARCS; IA++) {\r
                        if (IA == 1) {\r
                                I = NARCS;\r
@@ -2487,12 +2627,12 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                        if (A1 >= 1.0) {\r
                                ERR = ERR / A1;\r
                        }\r
-                       if (ERR > MXMIS) {\r
+                       if (ERR > MXMIS[0]) {\r
                                IMX = IA;\r
-                               MXMIS = ERR;\r
+                               MXMIS[0] = ERR;\r
                        }\r
                } // for (IA=1; IA <= NARCS; IA++)\r
-               if (MXMIS >= TOL1) {\r
+               if (MXMIS[0] >= TOL1) {\r
                        System.out.println("MAXIMUM CORNER MISFIT: " + MXMIS);\r
                        System.out.println("OCCURS AT CORNER: " + IMX);\r
                } else {\r
@@ -2534,13 +2674,13 @@ public class SymmsIntegralMapping extends AlgorithmBase {
 \r
                // **** SET DEFAULT PLOTTING DISTANCES, IF NECESSARY\r
 \r
-               if (PSD <= 0.0) {\r
-                       PSD = 1.6E+2;\r
+               if (PSD[0] <= 0.0) {\r
+                       PSD[0] = 1.6E+2;\r
                        MINPD = 2.0;\r
                        MAXPD = 5.0;\r
                }\r
-               RMIN = MINPD * TSD / PSD;\r
-               RMAX = MAXPD * TSD / PSD;\r
+               RMIN = MINPD * TSD / PSD[0];\r
+               RMAX = MAXPD * TSD / PSD[0];\r
                RMEAN = 0.5 * (RMIN + RMAX);\r
                MINC = Math.sqrt(EPS);\r
 \r
@@ -2548,7 +2688,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                // AND\r
                // **** TESTING\r
 \r
-               MXDIF = 0.0;\r
+               MXDIF[0] = 0.0;\r
                zzindex = 0;\r
                for (IA = 1; IA <= NARCS; IA++) {\r
                        TT[0] = -1.0;\r
@@ -2605,11 +2745,11 @@ public class SymmsIntegralMapping extends AlgorithmBase {
 \r
                                        zdiv(1.0 - NDZZ[0], -NDZZ[1], DZZ[0], DZZ[1], cr, ci);\r
                                        ERR = zabs(cr[0], ci[0]);\r
-                                       if (ERR > MXDIF) {\r
-                                               MXDIF = ERR;\r
+                                       if (ERR > MXDIF[0]) {\r
+                                               MXDIF[0] = ERR;\r
                                                IMX = IA;\r
                                                TMX = TT[0];\r
-                                       } // if (ERR > MXDIF)\r
+                                       } // if (ERR > MXDIF[0)\r
                                } // for (I=1; I <= 2; I++)\r
 \r
                                if (!LNSEG[IA - 1]) {\r
@@ -2662,11 +2802,11 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                        zzset[zzindex++][1] = ZZ0[1];\r
                }\r
 \r
-               if (MXDIF > DTOL) {\r
+               if (MXDIF[0] > DTOL) {\r
                        System.out.println();\r
                        System.out.println("POSSIBLE PARFUN/DPARFN INCONSISTECY ON ARC: " + IMX);\r
                        System.out.println("OCCURS AT STANDARDISED PARAMETER VALUE: " + TMX);\r
-                       System.out.println("RELATIVE FINITE DIFF ERROR: " + MXDIF);\r
+                       System.out.println("RELATIVE FINITE DIFF ERROR: " + MXDIF[0]);\r
                } else {\r
                        System.out.println();\r
                        System.out.println("PARFUN AND DPARFN ARE CONSISTENT:");\r