testStripmap1() plot okay. Plots are identical to Figure 4.6 except for 2 extra...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 14 Nov 2017 22:17:02 +0000 (22:17 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Tue, 14 Nov 2017 22:17:02 +0000 (22:17 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15252 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 7b38158..e69e49b 100644 (file)
@@ -367,6 +367,7 @@ public class SchwarzChristoffelMapping2 extends AlgorithmBase {
        private void testStripmap1() {\r
                // NLConstrainedEngine and NLConstrainedEP do not work.  Nl2sol does not work for f1, but it \r
                // works for f2.  NESolve works for both f1 and f2.\r
+               // Plots are identical to Figure 4.6 except for 2 extra lines from plotpoly\r
                 double x[] = new double[]{Double.POSITIVE_INFINITY, -1, -2.5, Double.POSITIVE_INFINITY, 2.4,\r
                                                   Double.POSITIVE_INFINITY, 2.4, Double.POSITIVE_INFINITY, -1, -2.5};\r
                 double y[] = new double[]{0, -1, -1, 0, -3.3, 0, -1.3, 0, 1, 1};\r
@@ -376,6 +377,11 @@ public class SchwarzChristoffelMapping2 extends AlgorithmBase {
                 scmap f1 = stripmap(p, endidx);\r
                 int endidx2[] = new int[]{3,7};\r
                 scmap f2 = stripmap(p, endidx2);\r
+                double re[] = new double[]{0};\r
+                double im[] = new double[]{8};\r
+                double axis[] = new double[]{-4.3, 5.7, -6.15, 3.85};\r
+                stripplot(f1, re, im, Integer.MIN_VALUE, false, axis);\r
+                stripplot(f2, re, im, Integer.MIN_VALUE, false, axis);\r
        }\r
        \r
        private void testStripmap2() {\r
@@ -2127,6 +2133,427 @@ public class SchwarzChristoffelMapping2 extends AlgorithmBase {
                return fprime;\r
        }\r
        \r
+       public void stripplot(scmap M, double re[], double im[], int yInvert, boolean closed, double axis[]) {\r
+               // Visualize a Schwarz-Christoffel strip map.\r
+               // stripplot plots the polygon associated with the Schwarz-Christoffel\r
+               // half-plane map M and the images of vertical lines whose real\r
+               // parts are given in re and horizontal lines whose imaginary\r
+               // parts are given in im under the S-C transformation.\r
+               // From 1998 MATLAB plot routine copyright by Toby Driscoll.\r
+               int i;\r
+               polygon p = M.poly;\r
+               double w[][] = p.vertex;\r
+               double beta[] = new double[p.angle.length];\r
+               for (i = 0; i < p.angle.length; i++) {\r
+                       beta[i] = p.angle[i] - 1;\r
+               }\r
+               double z[][] = M.prevertex;\r
+               double c[] = M.constant;\r
+               stplot(w, beta, z, c, re, im, yInvert, closed, axis);\r
+       }\r
+       \r
+       public void stplot(double w[][], double beta[], double z[][], double c[],\r
+                       double re[], double im[], int yInvert, boolean closed, double axis[]) {\r
+           // Image of cartesian grid under Schwarz-Christoffel strip map.\r
+               // stplot will adaptively plot the images under the Schwarz-Christoffel\r
+               // exterior map of vertical lines whose real parts are given in re and\r
+               // horizontal lines whose imaginary parts are given in im in the upper half-plane.\r
+               //The abscissae of the vertical lines will bracket the finite extremes of real(Z).\r
+               // From 1998 MATLAB stplot routine copyright by Toby Driscoll.\r
+               int i, j, m;\r
+               double tol;\r
+               boolean haveInfiniteZ = false;\r
+               double betan[];\r
+               int n = w.length;\r
+               double qdat[][] = null;\r
+               double axlim[] = new double[4];\r
+               // Number of quadrature points per integration.\r
+               // Approximately equals -log10(error).  Increase if plot\r
+               // has false little zigzags in curves. \r
+               int nqpts = 5; \r
+               // Minimum line segment length, as a proportion of the axes box\r
+        double minlen = 0.005;\r
+        // Maximum line segment length, as a proportion of the axes box\r
+        double maxlen = 0.02;\r
+        // Max allowed number of adaptive refinements made to meet other requirements \r
+        int maxrefn = 16;\r
+        Vector<Double>zpReal = new Vector<Double>();\r
+               Vector<Double>zpImag = new Vector<Double>();\r
+               Vector<Double>wpReal = new Vector<Double>();\r
+               Vector<Double>wpImag = new Vector<Double>();\r
+               Vector<Boolean>newlog = new Vector<Boolean>();\r
+               Vector<Double> linhx[][] = null;\r
+               Vector<Double> linhy[][] = null;\r
+               Vector<Double>x1Vector = new Vector<Double>();\r
+               Vector<Double>y1Vector = new Vector<Double>();\r
+               Vector<Double>x2Vector = new Vector<Double>();\r
+               Vector<Double>y2Vector = new Vector<Double>();\r
+        \r
+               // Zero arguments default to 10\r
+               if (((re == null) || (re.length == 0)) && ((im == null) || (im.length == 0))) {\r
+                   re = new double[]{10.0};\r
+                   im = new double[]{10.0};\r
+               }\r
+               \r
+               // Integer arguments must be converted to specific values\r
+               double minre = Double.MAX_VALUE;\r
+               double maxre = -Double.MAX_VALUE;\r
+               for (i = 0; i < z.length; i++) {\r
+                       if (!Double.isInfinite(z[i][0])) {\r
+                               if (z[i][0] < minre) {\r
+                                       minre = z[i][0];\r
+                               }\r
+                               if (z[i][0] > maxre) {\r
+                                       maxre = z[i][0];\r
+                               }\r
+                       }\r
+               }\r
+               if ((re.length == 1) && (re[0] == Math.round(re[0]))) {\r
+                       if (re[0] < 1) {\r
+                           re = null;\r
+                       }\r
+                       else if (re[0] < 2) {\r
+                               // Real parts are given in re.\r
+                               re[0] = (minre + maxre)/2.0;\r
+                       }\r
+                       else {\r
+                           m = (int)re[0];     \r
+                           double dre = (maxre - minre)/(m - 1.0);\r
+                           double spacing = (maxre - minre + 2*dre)/(m - 1.0);\r
+                           re = new double[m];\r
+                           for (i = 0; i < m; i++) {\r
+                               re[i] = minre - dre + i*spacing;\r
+                           }\r
+                       } // else \r
+               } // if ((re.length == 1) && (re[0] == Math.round(re[0])))\r
+               if ((im.length == 1) & (im[0] == Math.round(im[0]))) {\r
+                       if (im[0] < 1) {\r
+                               im = null;\r
+                       }\r
+                       else {\r
+                               m = (int)im[0];\r
+                               double spacing = 1.0/(m + 1.0);\r
+                               im = new double[m];\r
+                               for (i = 1; i < m+1; i++) {\r
+                                       im[i-1] = i*spacing;\r
+                               }\r
+                       }\r
+               } // if ((im.length == 1) & (im[0] == Math.round(im[0])))\r
+               boolean atinf[] = new boolean[w.length];\r
+               int numinfinite = 0;\r
+               int numfinite = 0;\r
+               for (i = 0; i < w.length; i++) {\r
+                       if ((!Double.isInfinite(w[i][0])) && (!Double.isInfinite(w[i][1]))) {\r
+                               numfinite++;\r
+                       }\r
+                       else {\r
+                               numinfinite++;\r
+                               atinf[i] = true;\r
+                       }\r
+               }\r
+               int first = -1;\r
+               for (i = 0; i < n-1 && (first == -1); i++) {\r
+                       if ((!atinf[i]) && (!atinf[i+1])) {\r
+                           first = i;  \r
+                       }\r
+               }\r
+               if ((first == -1) && (!atinf[n-1]) && (!atinf[0])) {\r
+                       first = n-1;\r
+               }\r
+               if (first == -1) {\r
+                       MipavUtil.displayError("There must be two consecutive finite vertices.");\r
+                       return;\r
+               }\r
+               int renum[] = new int[n];\r
+               for (i = 0; i < n - first; i++) {\r
+                       renum[i] = first +i;\r
+               }\r
+               for (i = 0; i < first; i++) {\r
+                       renum[n-first+i] = i;\r
+               }\r
+               double wrenum[][] = new double[n][2];\r
+               for (i = 0; i < n; i++) {\r
+                       wrenum[i][0] = w[renum[i]][0];\r
+                       wrenum[i][1] = w[renum[i]][1];\r
+               }\r
+               float xPointArray[];\r
+               float yPointArray[];\r
+               \r
+               // Put 2 finite points at every infinity in plotpoly\r
+               // Unless last one in wrenum is infinite in which case add 3\r
+               numinfinite = 0;\r
+               for (i = 0; i < w.length; i++) {\r
+                   if (Double.isInfinite(wrenum[i][0]) || Double.isInfinite(wrenum[i][1]))     {\r
+                       if (i == w.length-1) {\r
+                               numinfinite = numinfinite + 2;\r
+                       }\r
+                       else {\r
+                           numinfinite++;      \r
+                       }\r
+                   }\r
+               }\r
+               if (closed) {\r
+                   xPointArray = new float[n+1+numinfinite];\r
+                   yPointArray = new float[n+1+numinfinite];\r
+               }\r
+               else {\r
+                       xPointArray = new float[n+numinfinite];\r
+                       yPointArray = new float[n+numinfinite];\r
+               }\r
+               ViewJFrameGraph pointGraph = scm.plotpoly(xPointArray, yPointArray, w, beta, false, axlim, yInvert, closed, axis);\r
+               ViewJComponentGraph graph = pointGraph.getGraph();\r
+               Rectangle graphBounds = graph.getGraphBounds();\r
+               Graphics g = graph.getGraphics();\r
+               double xScale = graphBounds.width / (axlim[1] - axlim[0]);\r
+        double yScale = graphBounds.height / (axlim[3] - axlim[2]);\r
+        double len = Math.max(axlim[1] - axlim[0], axlim[3] - axlim[2]);\r
+               minlen = len * minlen;\r
+               maxlen = len * maxlen;\r
+               tol = Math.pow(10.0, -nqpts);\r
+               for (i = 0; i < z.length; i++) {\r
+                       if (Double.isInfinite(z[i][0]) || Double.isInfinite(z[i][1])) {\r
+                               haveInfiniteZ = true;\r
+                       }\r
+               }\r
+               qdat = new double[nqpts][2*beta.length+2];\r
+               scm.scqdata(qdat, beta, nqpts); \r
+               \r
+               \r
+               // Plot vertical lines\r
+               if (re != null) {\r
+               linhx = new Vector[re.length][2];\r
+               linhy = new Vector[re.length][2];\r
+               for (i = 0; i < re.length; i++) {\r
+                       for (j = 0; j < 2; j++) {\r
+                               linhx[i][j] = new Vector<Double>();\r
+                               linhy[i][j] = new Vector<Double>();\r
+                       }\r
+               }\r
+               \r
+               for (j = 0; j < re.length; j++) {\r
+                   zpReal.clear();\r
+                   zpImag.clear();\r
+                   double spacing = 1.0/14.0;\r
+                   for (i = 0; i < 15; i++) {\r
+                       zpReal.add(re[j]);\r
+                       zpImag.add(i*spacing);\r
+                   }\r
+                   newlog.clear();\r
+                   for (i = 0; i < 15; i++) {\r
+                       newlog.add(true);\r
+                   }\r
+                   wpReal.clear();\r
+                   wpImag.clear();\r
+                   for (i = 0; i < 15; i++) {\r
+                       wpReal.add(Double.NaN);\r
+                   \r
+                       wpImag.add(0.0);\r
+                   }\r
+                   \r
+                       // The individual points will be shown as they are found.\r
+                               \r
+                       // Adaptive refinement to make curve smooth\r
+                       int iter = 0;\r
+                       while (iter < maxrefn) {\r
+                               int numnewlog = 0;\r
+                               for (i = 0; i < newlog.size(); i++) {\r
+                                   if (newlog.get(i)) {\r
+                                       numnewlog++;\r
+                                   }\r
+                               } // for (i = 0; i < newlog.size(); i++)\r
+                               if (numnewlog == 0) {\r
+                                       break;\r
+                               }\r
+                               double zpnew[][] = new double[numnewlog][2];\r
+                               for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                   if (newlog.get(i)) {\r
+                                       zpnew[m][0] = zpReal.get(i);\r
+                                       zpnew[m++][1] = zpImag.get(i);\r
+                                   }\r
+                               } // for (i = 0, m = 0; i < newlog.length; i++)\r
+                               double neww[][] = scm.stmap(zpnew, w, beta, z, c, qdat);\r
+                               for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                       if (newlog.get(i)) {\r
+                                           wpReal.set(i, neww[m][0]);\r
+                                           wpImag.set(i, neww[m++][1]);\r
+                                       } \r
+                               } // for (i = 0, m = 0; i < newlog.size(); i++)\r
+                               iter = iter + 1;\r
+                               \r
+                               linhx[j][0].clear();\r
+                               linhy[j][0].clear();\r
+                               linhx[j][1].clear();\r
+                               linhy[j][1].clear();\r
+                               // Update the points to show progress\r
+                               for (i = 0; i < wpReal.size(); i++) {\r
+                                   linhx[j][0].add(wpReal.get(i));\r
+                                   linhy[j][0].add(wpImag.get(i));\r
+                                   linhx[j][1].add(zpReal.get(i));\r
+                                   linhy[j][1].add(zpImag.get(i));\r
+                               }\r
+                               \r
+                               // Add points to zp where necessary\r
+                               scm.scpadapt(zpReal, zpImag, wpReal, wpImag, newlog, minlen, maxlen, axlim);\r
+                   } // while (iter < maxrefn)\r
+               } // for (j = 0; j < re.length; j++)\r
+               \r
+               for (i = 0; i < re.length; i++) {\r
+                   if ((linhx[i][0] != null)&& (linhy[i][0] != null)) {\r
+                       for (j = 0; j < linhx[i][0].size()-1; j++) {\r
+                               double posx1 = linhx[i][0].get(j);\r
+                               double posy1 = linhy[i][0].get(j);\r
+                               double posx2 = linhx[i][0].get(j+1);\r
+                               double posy2 = linhy[i][0].get(j+1);\r
+                               x1Vector.add(posx1);\r
+                               y1Vector.add(posy1);\r
+                               x2Vector.add(posx2);\r
+                               y2Vector.add(posy2);\r
+                           int x1 =  (int)Math.round(graphBounds.x + xScale*(posx1 - axlim[0]));\r
+                           int y1 =  (int)Math.round(graphBounds.y + yScale*(posy1 - axlim[2]));\r
+                           y1 = -y1 + 2*graphBounds.y + graphBounds.height;\r
+                           int x2 =  (int)Math.round(graphBounds.x + xScale*(posx2 - axlim[0]));\r
+                           int y2 =  (int)Math.round(graphBounds.y + yScale*(posy2 - axlim[2]));\r
+                           y2 = -y2 + 2*graphBounds.y + graphBounds.height;\r
+                           graph.drawLine(g, x1, y1, x2, y2);\r
+                       }\r
+                   }\r
+               }\r
+               } // if (re != null)\r
+               \r
+               // Plot horizontal lines\r
+               if (im != null) {\r
+               linhx = new Vector[im.length][2];\r
+               linhy = new Vector[im.length][2];\r
+               for (i = 0; i < im.length; i++) {\r
+                       for (j = 0; j < 2; j++) {\r
+                               linhx[i][j] = new Vector<Double>();\r
+                               linhy[i][j] = new Vector<Double>();\r
+                       }\r
+               }\r
+               \r
+               double x1h = Math.min(-5, minre);\r
+               double x2h = Math.max(5, maxre);\r
+               double wl[] = new double[2]; // image of left end\r
+               double wr[] = new double[2]; // image of right end\r
+               for (i = 0; i < z.length; i++) {\r
+                       if (Double.isInfinite(z[i][0])) {\r
+                               if (z[i][0] < 0) {\r
+                                       wl[0] = w[i][0];\r
+                                       wl[1] = w[i][1];\r
+                               }\r
+                               else if (z[i][0] > 0) {\r
+                                       wr[0] = w[i][0];\r
+                                       wr[1] = w[i][1];\r
+                               }\r
+                       }\r
+               }\r
+               for (j = 0; j < im.length; j++) {\r
+                   // Start evenly spaced\r
+                       zpReal.clear();\r
+                       zpImag.clear();\r
+                       //zpReal.add(Double.NEGATIVE_INFINITY);\r
+                       //zpImag.add(im[j]);\r
+                       for (i = 0; i < 15; i++) {\r
+                               zpReal.add(x1h + i*(x2h-x1h)/14.0);\r
+                               zpImag.add(im[j]);\r
+                       }\r
+                       //zpReal.add(Double.POSITIVE_INFINITY);\r
+                       //zpImag.add(im[j]);\r
+                       newlog.clear();\r
+                       //newlog.add(false);\r
+                   for (i = 0; i < 15; i++) {\r
+                       newlog.add(true);\r
+                   }\r
+                   //newlog.add(false);\r
+                   wpReal.clear();\r
+                   wpImag.clear();\r
+                   //wpReal.add(wl[0]);\r
+                   //wpImag.add(wl[1]);\r
+                   for (i = 0; i < 15; i++) {\r
+                       wpReal.add(Double.NaN);\r
+                       wpImag.add(0.0);\r
+                   }\r
+                   //wpReal.add(wr[0]);\r
+                   //wpImag.add(wr[1]);\r
+                   \r
+                   // The individual points will be shown as they are found.\r
+                       \r
+                       // Adaptive refinement to make curve smooth\r
+                       int iter = 0;\r
+                       while (iter < maxrefn) {\r
+                               int numnewlog = 0;\r
+                               for (i = 0; i < newlog.size(); i++) {\r
+                                   if (newlog.get(i)) {\r
+                                       numnewlog++;\r
+                                   }\r
+                               } // for (i = 0; i < newlog.size(); i++)\r
+                               if (numnewlog == 0) {\r
+                                       break;\r
+                               }\r
+                               double zpnew[][] = new double[numnewlog][2];\r
+                               for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                   if (newlog.get(i)) {\r
+                                       zpnew[m][0] = zpReal.get(i);\r
+                                       zpnew[m++][1] = zpImag.get(i);\r
+                                   }\r
+                               } // for (i = 0, m = 0; i < newlog.length; i++)\r
+                               double neww[][] = scm.stmap(zpnew, w, beta, z, c, qdat);\r
+                               for (i = 0, m = 0; i < newlog.size(); i++) {\r
+                                       if (newlog.get(i)) {\r
+                                           wpReal.set(i, neww[m][0]);\r
+                                           wpImag.set(i, neww[m++][1]);\r
+                                       } \r
+                               } // for (i = 0, m = 0; i < newlog.size(); i++)\r
+                               iter = iter + 1;\r
+                               \r
+                               linhx[j][0].clear();\r
+                               linhy[j][0].clear();\r
+                               linhx[j][1].clear();\r
+                               linhy[j][1].clear();\r
+                               // Update the points to show progress\r
+                               for (i = 0; i < wpReal.size(); i++) {\r
+                                   linhx[j][0].add(wpReal.get(i));\r
+                                   linhy[j][0].add(wpImag.get(i));\r
+                                   linhx[j][1].add(zpReal.get(i));\r
+                                   linhy[j][1].add(zpImag.get(i));\r
+                               }\r
+                               \r
+                               // Add points to zp where necessary\r
+                               scm.scpadapt(zpReal, zpImag, wpReal, wpImag, newlog, minlen, maxlen, axlim);\r
+                   } // while (iter < maxrefn)\r
+               } // for (j = 0; j < im.length; j++)\r
+               \r
+               for (i = 0; i < im.length; i++) {\r
+                   if ((linhx[i][0] != null)&& (linhy[i][0] != null)) {\r
+                       for (j = 0; j < linhx[i][0].size()-1; j++) {\r
+                               double posx1 = linhx[i][0].get(j);\r
+                               double posy1 = linhy[i][0].get(j);\r
+                               double posx2 = linhx[i][0].get(j+1);\r
+                               double posy2 = linhy[i][0].get(j+1);\r
+                               x1Vector.add(posx1);\r
+                               y1Vector.add(posy1);\r
+                               x2Vector.add(posx2);\r
+                               y2Vector.add(posy2);\r
+                           int x1 =  (int)Math.round(graphBounds.x + xScale*(posx1 - axlim[0]));\r
+                           int y1 =  (int)Math.round(graphBounds.y + yScale*(posy1 - axlim[2]));\r
+                           y1 = -y1 + 2*graphBounds.y + graphBounds.height;\r
+                           int x2 =  (int)Math.round(graphBounds.x + xScale*(posx2 - axlim[0]));\r
+                           int y2 =  (int)Math.round(graphBounds.y + yScale*(posy2 - axlim[2]));\r
+                           y2 = -y2 + 2*graphBounds.y + graphBounds.height;\r
+                           graph.drawLine(g, x1, y1, x2, y2);\r
+                       }\r
+                   }\r
+               }\r
+               } // if (im != null)\r
+               \r
+               graph.setX1Vector(x1Vector);\r
+               graph.setY1Vector(y1Vector);\r
+               graph.setX2Vector(x2Vector);\r
+               graph.setY2Vector(y2Vector);\r
+               graph.setAddSchwarzChristoffelLines(true);\r
+               graph.paintComponent(g);\r
+       }\r
+       \r
        public void hplot(scmap M, double re[], double im[], int yInvert, boolean closed, double axis[]) {\r
                // Visualize a Schwarz-Christoffel half-plane map.\r
                // hplot plots the polygon associated with the Schwarz-Christoffel\r
@@ -2150,7 +2577,7 @@ public class SchwarzChristoffelMapping2 extends AlgorithmBase {
                        double re[], double im[], int yInvert, boolean closed, double axis[]) {\r
            // Image of cartesian grid under Schwarz-Christoffel half-plane map.\r
                // hpplot will adaptively plot the images under the Schwarz-Christoffel\r
-               // exterior map of vertical lines whose real parrts are given in re and\r
+               // exterior map of vertical lines whose real parts are given in re and\r
                // horizontal lines whose imaginary parts are given in im.\r
                // From 1998 MATLAB hpplot routine copyright by Toby Driscoll.\r
                int i, j, m;\r