Porting code.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 16 Nov 2017 22:22:57 +0000 (22:22 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 16 Nov 2017 22:22:57 +0000 (22:22 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15262 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 63181c2..98bd92e 100644 (file)
@@ -1,5 +1,7 @@
 package gov.nih.mipav.model.algorithms;\r
 \r
+import gov.nih.mipav.view.MipavUtil;\r
+\r
 public class DoublyConnectedSC extends AlgorithmBase {\r
        // This is a port of the FORTRAN ACM TOMS Algorithm 785, collected\r
        // algorithms form ACM.  The original ACM work was published in Transactions\r
@@ -7,8 +9,750 @@ public class DoublyConnectedSC extends AlgorithmBase {
        // The original FORTRAN code was written by Dr. Chenglie Hu as described in\r
        // his article "A Software Package for Computing Schwarz-Christoffel Conformal\r
        // Transformation for Doubly Connected Polyhgonal Regions."\r
+       \r
+       // geometries of the polygon region in the 7 test routines\r
+       private final int SQUARE_SYMMETRIC_REGION = 1;\r
+       private final int MILDLY_CROWDED_INFINITE_REGION = 2;\r
+       private final int HEAVILY_CROWDED_REGION = 3;\r
+       private final int CHINESE_CHARACTER_STRUCTURED_REGION = 4;\r
+       private final int FOUR_DIRECTION_INFINITE_REGION = 5;\r
+       private final int EXAMPLE_GIVEN_FOR_CHECKING_INVERSE_MAP = 6;\r
+       private final int UPPER_HALF_PLANE_WITH_HORIZONTAL_SLIT = 7;\r
+       \r
+       // Below are the 4 variables in the common block PARAM4\r
+       private double DLAM;\r
+       private int IU;\r
+       private double UARY[] = new double[8];\r
+       private double VARY[] = new double[3];\r
+       \r
+       private SchwarzChristoffelMapping scm = new SchwarzChristoffelMapping();\r
+       \r
+       // Specifies the geometry of the polygon region for 7 test examples\r
+       private int IPOLY;\r
+       // The number of Gauss-Jacobi points\r
+       // Recommended values for NPTQ are 2-8.\r
+       private int NPTQ;\r
+       private boolean testRoutine = false;\r
+       private double MACHEP = 2.2204460E-16;\r
+       \r
+       public DoublyConnectedSC() {\r
+               \r
+       }\r
+       \r
+       public DoublyConnectedSC(int IPOLY, int NPTQ) {\r
+               this.IPOLY = IPOLY;\r
+               this.NPTQ = NPTQ;\r
+               testRoutine = true;\r
+       }\r
+       \r
        public void runAlgorithm() {\r
+               int M[] = new int[1];\r
+               int N[] = new int[1];\r
+               double Z0[][] = new double[30][2];\r
+               double Z1[][] = new double[30][2];\r
+               double ALFA0[] = new double[30];\r
+               double ALFA1[] = new double[30];\r
+               double QWORK[] = new double[1660];\r
+               int ISHAPE;\r
+               \r
+double neweps;\r
+               \r
+               // MACHEP is a machine dependent parameter specifying the relative precision of\r
+               // floating point arithmetic.\r
+                // epsilon = D1MACH(4)\r
+        // Machine epsilon is the smallest positive epsilon such that\r
+        // (1.0 + epsilon) != 1.0.\r
+        // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
+        // epsilon = 2.2204460e-16\r
+        // epsilon is called the largest relative spacing\r
+        MACHEP = 1.0;\r
+        neweps = 1.0;\r
+\r
+        while (true) {\r
+\r
+            if (1.0 == (1.0 + neweps)) {\r
+                break;\r
+            } else {\r
+                MACHEP = neweps;\r
+                neweps = neweps / 2.0;\r
+            }\r
+        } // while(true)\r
+               \r
+               if (testRoutine) {\r
+                       DSCDATA(IPOLY, M, N, Z0, Z1, ALFA0, ALFA1);\r
+                       ANGLES(N[0], Z1, ALFA1, 1);\r
+                       if ((IPOLY == 1) || (IPOLY == 3) || (IPOLY == 4) || (IPOLY == 6)) {\r
+                               ANGLES(M[0], Z0, ALFA0, 0);\r
+                       }\r
+                       \r
+                       // Generate the Gauss-Jacobi weights & nodes and check the input\r
+                       QINIT(M[0],N[0],ALFA0, ALFA1, NPTQ, QWORK);\r
+                       ISHAPE = 0;\r
+                       if ((IPOLY == 2) || (IPOLY == 5) || (IPOLY == 7)) {\r
+                               ISHAPE = 1;\r
+                       }\r
+                       CHECK(ALFA0, ALFA1, M[0], N[0], ISHAPE);\r
+               } // if (testRoutine)\r
+               \r
+               \r
+       }\r
+       \r
+       private void DSCDATA(int IPOLY, int M[], int N[], double Z0[][], double Z1[][], double ALFA0[], double ALFA1[]) {\r
+               // DSCDATA generates data.\r
+               double Q;\r
+               double S;\r
+               if (IPOLY == 1) {\r
+                       M[0] = 4;\r
+                       N[0] = 4;\r
+                       Q = Math.sqrt(2.0);\r
+                       Z0[0][0] = 1.0 + Q;\r
+                       Z0[0][1] = 1.0 + Q;\r
+                       Z0[1][0] = -1.0 - Q;\r
+                       Z0[1][1] = 1.0 + Q;\r
+                       Z0[2][0] = -1.0 - Q;\r
+                       Z0[2][1] = -1.0 - Q;\r
+                       Z0[3][0] = 1.0 + Q;\r
+                       Z0[3][1] = -1.0 - Q;\r
+                       Z1[0][0] = Q;\r
+                       Z1[0][1] = 0.0;\r
+                       Z1[1][0] = 0.0;\r
+                       Z1[1][1] = Q;\r
+                       Z1[2][0] = -Q;\r
+                       Z1[2][1] = 0.0;\r
+                       Z1[3][0] = 0.0;\r
+                       Z1[3][1] = -Q;\r
+               } // if (IPOLY == 1)\r
+               else if (IPOLY == 2) {\r
+                       M[0] = 12;\r
+                       N[0] = 6;\r
+                       Z0[0][0] = 0.6875;\r
+                       Z0[0][1] = 0.875;\r
+                       Z0[1][0] = 0.0;\r
+                       Z0[1][1] = 0.0;\r
+                       Z0[2][0] = 1.0;\r
+                       Z0[2][1] = 0.0;\r
+                       Z0[3][0] = 0.875;\r
+                       Z0[3][1] = 0.875;\r
+                       Z0[4][0] = 1.125;\r
+                       Z0[4][1] = 1.375;\r
+                       Z0[5][0] = 2.0;\r
+                       Z0[5][1] = 1.375;\r
+                       Z0[6][0] = 1.25;\r
+                       Z0[6][1] = 2.0;\r
+                       Z0[7][0] = 2.25;\r
+                       Z0[7][1] = 2.0;\r
+                       Z0[8][0] = 2.375;\r
+                       Z0[8][1] = 2.75;\r
+                       Z0[9][0] = 1.625;\r
+                       Z0[9][1] = 2.25;\r
+                       Z0[10][0] = 1.125;\r
+                   Z0[10][1] = 2.625;\r
+                   Z0[11][0] = -0.5;\r
+                   Z0[11][1] = 2.75;\r
+                   Z1[0][0] = 0.375;\r
+                   Z1[0][1] = 1.875;\r
+                   Z1[1][0] = 0.5;\r
+                   Z1[1][1] = 2.0;\r
+                   Z1[2][0] = 1.0;\r
+                   Z1[2][1] = 1.5;\r
+                   Z1[3][0] = 0.5;\r
+                   Z1[3][1] = 2.1875;\r
+                   Z1[4][0] = 0.5;\r
+                   Z1[4][1] = 2.5;\r
+                   Z1[5][0] = Z1[3][0];\r
+                   Z1[5][1] = Z1[3][1];\r
+                   ALFA0[0] = 1.39169261159339475;\r
+            ALFA0[1] = 0.28801540784794967;\r
+            ALFA0[2] = 0.454832764699133488;\r
+            ALFA0[3] = 1.19275085295129979;\r
+            ALFA0[4] = 1.35241638234956651;\r
+            ALFA0[5] = 0.0;\r
+            ALFA0[6] = 2.0;\r
+            ALFA0[7] = 0.552568456711253445;\r
+            ALFA0[8] = 0.260264501477747753;\r
+            ALFA0[9] = 1.39199980651013222;\r
+            ALFA0[10] = 0.819604487273064009;\r
+            ALFA0[11] = 0.295854728586457991;\r
+               } // else if (IPOLY == 2)\r
+               else if (IPOLY == 3) {\r
+                       M[0] = 11;\r
+                       N[0] = 6;\r
+                       Z0[0][0] = 0.5;\r
+                       Z0[0][1] = 2.5;\r
+                       Z0[1][0] = 0.5;\r
+                       Z0[1][1] = 0.5;\r
+                       Z0[2][0] = 1.0;\r
+                       Z0[2][1] = 0.5;\r
+                       Z0[3][0] = 1.0;\r
+                       Z0[3][1] = 1.0;\r
+                       Z0[4][0] = 1.0;\r
+                       Z0[4][1] = 0.5;\r
+                       Z0[5][0] = 0.5;\r
+                       Z0[5][1] = 0.5;\r
+                       Z0[6][0] = 0.5;\r
+                       Z0[6][1] = 2.5;\r
+                       Z0[7][0] = 0.0;\r
+                       Z0[7][1] = 2.5;\r
+                       Z0[8][0] = 0.0;\r
+                       Z0[8][1] = 0.0;\r
+                       Z0[9][0] = 2.0;\r
+                       Z0[9][1] = 0.0;\r
+                       Z0[10][0] = 2.0;\r
+                       Z0[10][1] = 2.5;\r
+                       Z1[0][0] = 1.0;\r
+                       Z1[0][1] = 2.0;\r
+                       Z1[1][0] = 1.0;\r
+                       Z1[1][1] = 1.5;\r
+                       Z1[2][0] = 1.0;\r
+                       Z1[2][1] = 2.0;\r
+                       Z1[3][0] = 1.5;\r
+                       Z1[3][1] = 2.0;\r
+                       Z1[4][0] = 1.5;\r
+                       Z1[4][1] = 0.5;\r
+                       Z1[5][0] = 1.5;\r
+                       Z1[5][0] = 2.0;\r
+                       ALFA0[0] = 1.0/2.0;\r
+            ALFA0[1] = 1.0/2.0;\r
+            ALFA0[2] = 1.0/2.0;\r
+            ALFA0[3] = 2.0;\r
+            ALFA0[4] = 3.0/2.0;\r
+            ALFA0[5] = 3.0/2.0;\r
+            ALFA0[6] = 1.0/2.0;\r
+            ALFA0[7] = 1.0/2.0;\r
+            ALFA0[8] = 1.0/2.0;\r
+            ALFA0[9] = 1.0/2.0;\r
+            ALFA0[10] = 1.0/2.0;\r
+            ALFA1[0] = 3.0/2.0;\r
+            ALFA1[1] = 2.0;\r
+            ALFA1[2] = 1.0/2.0;\r
+            ALFA1[3] = 1.0/2.0;\r
+            ALFA1[4] = 2.0;\r
+            ALFA1[5] = 3.0/2.0;\r
+               } // else if (IPOLY == 3)\r
+               else if (IPOLY == 4) {\r
+                   M[0] = 4;\r
+                   N[0] = 17;\r
+                   Z0[0][0] = -1.0;\r
+                   Z0[0][1] = -1.0;\r
+                   Z0[1][0] = 1.0;\r
+                   Z0[1][1] = -1.0;\r
+                   Z0[2][0] = 1.0;\r
+                   Z0[2][1] = 1.0;\r
+                   Z0[3][0] = -1.0;\r
+                   Z0[3][1] = 1.0;\r
+                   Z1[0][0] = 0.0;\r
+                   Z1[0][1] = 0.5;\r
+                   Z1[1][0] = 0.0;\r
+                   Z1[1][1] = 0.0;\r
+                   Z1[2][0] = -0.5;\r
+                   Z1[2][1] = 0.0;\r
+                   Z1[3][0] = Z1[1][0];\r
+                   Z1[3][1] = Z1[1][1];\r
+                   Z1[4][0] = 0.0;\r
+                   Z1[4][1] = -0.5;\r
+                   Z1[5][0] = -0.5;\r
+                   Z1[5][1] = -0.5;\r
+                   Z1[6][0] = 0.5;\r
+                   Z1[6][1] = -0.5;\r
+                   Z1[7][0] = 0.25;\r
+                   Z1[7][1] = -0.5;\r
+                   Z1[8][0] = 0.25;\r
+                   Z1[8][1] = -0.25;\r
+                   Z1[9][0] = Z1[7][0];\r
+                   Z1[9][1] = Z1[7][1];\r
+                   Z1[10][0] = Z1[4][0];\r
+                   Z1[10][1] = Z1[4][1];\r
+                   Z1[11][0] = Z1[1][0];\r
+                   Z1[11][1] = Z1[1][1];\r
+                   Z1[12][0] = 0.5;\r
+                   Z1[12][1] = 0.0;\r
+                   Z1[13][0] = Z1[1][0];\r
+                   Z1[13][1] = Z1[1][1];\r
+                   Z1[14][0] = Z1[0][0];\r
+                   Z1[14][1] = Z1[0][1];\r
+                   Z1[15][0] = 0.5;\r
+                   Z1[15][1] = 0.5;\r
+                   Z1[16][0] = -0.5;\r
+                   Z1[16][1] = 0.5;\r
+               } // else if (IPOLY == 4)\r
+               else if (IPOLY == 5) {\r
+                       M[0] = 12;\r
+                       N[0] = 4;\r
+                       S = 1.0/6.0;\r
+                       Z0[0][0] = 0.0;\r
+                       Z0[0][1] = -8.0*S;\r
+                       Z0[1][0] = 0.0;\r
+                       Z0[1][1] = 0.0;\r
+                       Z0[2][0] = 1.0;\r
+                       Z0[2][1] = -1.5;\r
+                       Z0[3][0] = 1.0;\r
+                       Z0[3][1] = -0.5;\r
+                       Z0[4][0] = 0.0;\r
+                       Z0[4][1] = 0.0;\r
+                       Z0[5][0] = 0.0;\r
+                       Z0[5][1] = 2.0*S;\r
+                       Z0[6][0] = 0.0;\r
+                       Z0[6][1] = 0.0;\r
+                       Z0[7][0] = -1.0;\r
+                       Z0[7][1] = 0.0;\r
+                       Z0[8][0] = 0.0;\r
+                       Z0[8][1] = 0.0;\r
+                       Z0[9][0] = -7.0*S;\r
+                       Z0[9][1] = -1.0;\r
+                       Z0[10][0] = 0.0;\r
+                       Z0[10][1] = 0.0;\r
+                       Z0[11][0] = -2.0*S;\r
+                       Z0[11][1] = -10.0*S;\r
+                       ALFA0[0] = 1.75;\r
+            ALFA0[1] = 0.0;\r
+            ALFA0[2] = 1.0;\r
+            ALFA0[3] = 1.0;\r
+            ALFA0[4] = (Math.PI-3.5)/Math.PI;\r
+            ALFA0[5] = 3.5/Math.PI;\r
+            ALFA0[6] = 1.5;\r
+            ALFA0[7] = 1.0;\r
+            ALFA0[8] = 0.0;\r
+            ALFA0[9] = 1.75;\r
+            ALFA0[10] = 0.0;\r
+            ALFA0[11] = 1.0;\r
+            Z1[0][0] = 2.0*S;\r
+            Z1[0][1] = -0.5;\r
+            Z1[1][0] = -S;\r
+            Z1[1][1] = -2.0*S;\r
+            Z1[2][0] = -4.0*S;\r
+            Z1[2][1] = -5.0*S;\r
+            Z1[3][0] = 0.0;\r
+            Z1[3][1] = -1.0;\r
+               } // else if (IPOLY == 5)\r
+               else if (IPOLY == 6) {\r
+                   M[0] = 7;\r
+                   N[0] = 2;\r
+                   Z0[0][0] = -2.0;\r
+                   Z0[0][1] = -1.0;\r
+                   Z0[1][0] = 2.0;\r
+                   Z0[1][1] = -1.0;\r
+                   Z0[2][0] = 2.0;\r
+                   Z0[2][1] = 2.0;\r
+                   Z0[3][0] = -0.8;\r
+                   Z0[3][1] = 2.0;\r
+                   Z0[4][0] = 1.0;\r
+                   Z0[4][1] = 0.5;\r
+                   Z0[5][0] = -1.0;\r
+                   Z0[5][1] = 2.0;\r
+                   Z0[6][0] = -2.0;\r
+                   Z0[6][1] = 2.0;\r
+                   Z1[0][0] = 0.0;\r
+                   Z1[0][1] = 0.0;\r
+                   Z1[1][0] = -1.0;\r
+                   Z1[1][1] = 0.0;\r
+               } // else if (IPOLY == 6)\r
+               else {\r
+                       M[0] = 3;\r
+                       N[0] = 2;\r
+                       Z0[0][0] = 1.01;\r
+                       Z0[0][1] = 0.0;\r
+                       Z0[1][0] = 100.0;\r
+                       Z0[1][1] = 100.0;\r
+                       Z0[2][0] = -1.01;\r
+                       Z0[2][1] = 0.0;\r
+                       Z1[0][0] = 0.0;\r
+                       Z1[0][1] = 2.0;\r
+                       Z1[1][0] = 0.0;\r
+                       Z1[1][1] = 1.0;\r
+                       ALFA0[0] = 1.0;\r
+                       ALFA0[1] = -1.0;\r
+                       ALFA0[2] = 1.0;\r
+               } // else               \r
+       }\r
+       \r
+       private void ANGLES(int MN, double Z01[][], double ALFA01[], int I01) {\r
+           // Computes the interior angles of a doubly\r
+               // connected and bounded polygonal region.\r
+               // I01 := 0(outer polygon) := 1(inner polygon)\r
+               // MN: The number of vertices\r
+               // Z01: Array containing vertices\r
+               // ALFA01: Array containing interior turning angles on return\r
+               \r
+               int K, KL, KR;\r
+               double cr[] = new double[1];\r
+               double ci[] = new double[1];\r
+               \r
+               for (K = 1; K <= MN; K++) {\r
+                       KL = (K+MN-2)%MN + 1;\r
+                       KR = K%MN + 1;\r
+                       scm.zdiv(Z01[KL-1][0] - Z01[K-1][0], Z01[KL-1][1] - Z01[K-1][1], \r
+                                        Z01[KR-1][0] - Z01[K-1][0], Z01[KR-1][1] - Z01[K-1][1], cr, ci);\r
+                       ALFA01[K-1] = Math.atan2(ci[0], cr[0])/Math.PI;\r
+                       if (ALFA01[K-1] <= 0.0) {\r
+                               ALFA01[K-1] = ALFA01[K-1] + 2.0;\r
+                       }\r
+               } // for (K = 1; K <= MN; K++)\r
+               if (I01 == 1) {\r
+                       for (K = 1; K <= MN; K++) {\r
+                           if (ALFA01[K-1] != 2.0) {\r
+                               ALFA01[K-1] = 2.0 - ALFA01[K-1];\r
+                           }\r
+                       } // for (K = 1; K <= MN; K++)\r
+               } // if (I01 == 1)\r
+               return;\r
+       }\r
+       \r
+       private void CHECK(double ALFA0[], double ALFA1[], int M, int N, int ISHAPE) {\r
+               //   CHECKS IF THE INPUT DATA AND PARAMETERS ARE CORRECT.\r
+               //   NOTE1:    ANGLE-CHECKING MAKES SENSE ONLY IF ANGLES\r
+               //             ARE USER-PROVIDED.\r
+               //   NOTE2:    USERS ARE RESPONSIBLE FOR CHECKING THE FOLLOWING:\r
+               //   1. EACH COMPONENT OF THE OUTER POLYGON CONTAINS AT LEAST ONE\r
+               //      FINITE VERTEX.\r
+               //   2. Z1(N) IS THE CLOSEST (OR ALMOST CLOSEST)\r
+               //      INNER VERTEX TO Z1(M).\r
+               //   3. COUNTERCLOCKWISE ORIENTATION OF THE VERTICES.\r
+               \r
+               // ALFA0(M), ALFA1(N)\r
+               \r
+               if (!((M >= 3) && (M <= 30) && (N <= 30) && (M+N <= 40))) {\r
+                       System.err.println("M must be no less than 3");\r
+                       System.err.println("M, N must be no greater than 30");\r
+                       System.err.println("M+N must be no greater than 40");\r
+                       System.exit(-1);\r
+               }\r
+       \r
+       }\r
+       \r
+       private void QINIT(int M, int N, double ALFA0[], double ALFA1[], int NPTQ, double QWORK[]) {\r
+           // Computes the Gauss-Jacobi nodes & weights for Gauss-Jacobi quadrature.  Work array\r
+               // QWORK must be dimensioned at least NPTQ*(2(M+N) + 3).  It is divided up into\r
+               // 2(M+N)+3 vectors of length NPTQ:  The first M+N+1 contain quadrature nodes on output,\r
+               // the next M+N+1 contain the corresponding weights on output, and the last one is a scratch\r
+               // vector used by GAUSSJ.  NPTQ is the number of G-J nodes (same as weights) used.  See\r
+               // comment on routine WPROD for the rest of the calling sequence.\r
+               \r
+               // For each finite vertex, compute nodes & weights\r
+               // for one-sided Gauss-Jacobi quadrature:\r
+               \r
+               // ALFA0(M), ALFA1(N), QWORK(1660)\r
+               double ALPHA;\r
+               int INODES, ISCR, IWTS, J, K;\r
+               int i;\r
+               // B is used for input scratch array, T  for nodes output and W for weights output\r
+               double B[] = new double[NPTQ];\r
+               double T[] = new double[NPTQ];\r
+               double W[] = new double[NPTQ];\r
+               \r
+               ISCR = NPTQ * (2 * (M+N) + 2) + 1;\r
+               for (K = 1; K <= M + N; K++) {\r
+                   INODES = NPTQ * (K-1) + 1;\r
+                   IWTS = NPTQ * (M+N+K) + 1;\r
+                   if (K <= M) {\r
+                       ALPHA = ALFA0[K-1] - 1.0;\r
+                       if (ALFA0[K-1] > 0.0) {\r
+                           GAUSSJ(NPTQ, 0.0, ALPHA, B, T, W);\r
+                           for (i = 0; i < NPTQ; i++) {\r
+                               QWORK[INODES-1+i] = T[i];\r
+                               QWORK[IWTS-1+i] = W[i];\r
+                           }\r
+                       } // if (ALFA0[K-1] > 0.0)\r
+                       else {\r
+                               for (J = 0; J < NPTQ; J++) {\r
+                                       QWORK[IWTS+J-1] = 0.0;\r
+                                       QWORK[INODES+J-1] = 0.0;\r
+                               }\r
+                       } // else\r
+                   } // if (K <= M)\r
+                   else {\r
+                       ALPHA = ALFA1[K-M-1] - 1.0;\r
+                       GAUSSJ(NPTQ, 0.0, ALPHA, B, T, W);\r
+                   for (i = 0; i < NPTQ; i++) {\r
+                       QWORK[INODES-1+i] = T[i];\r
+                       QWORK[IWTS-1+i] = W[i];\r
+                   }\r
+                   } // else\r
+                   \r
+                   // Take singularities into account in advance for the purpose of saving\r
+                   // certain amount of calculation in WQSUM:\r
+                   for (J = 0; J < NPTQ; J++) {\r
+                       QWORK[IWTS+J-1] = QWORK[IWTS+J-1] * Math.pow((1.0 + QWORK[INODES+J-1]), (-ALPHA));\r
+                   }\r
+               } // for (K = 1; K <= M + N; K++)\r
+               \r
+               // Compute weights and nodes for pure Gaussian quadrature:\r
+               INODES = NPTQ *(M+N) + 1;\r
+               IWTS = NPTQ * (2* (M+N) + 1) + 1;\r
+               GAUSSJ(NPTQ, 0.0, 0.0, B, T, W);\r
+               for (i = 0; i < NPTQ; i++) {\r
+               QWORK[INODES-1+i] = T[i];\r
+               QWORK[IWTS-1+i] = W[i];\r
+           }\r
+               return;\r
+       }\r
+       \r
+       private void GAUSSJ(int N, double ALPHA, double BETA, double B[], double T[], double W[]) {\r
+           // B(N), T(N), W(N)\r
+           //        THIS ROUTINE COMPUTES THE NODES T(J) AND WEIGHTS\r
+               //        W(J) FOR GAUSS-JACOBI QUADRATURE FORMULAS.\r
+               //        THESE ARE USED WHEN ONE WISHES TO APPROXIMATE\r
+               //\r
+               //                 INTEGRAL (FROM A TO B)  F(X) W(X) DX\r
+               //\r
+               //                              N\r
+               //        BY                   SUM W  F(T )\r
+               //                             J=1  J    J\r
+               //\r
+               //        (W(X) AND W(J) HAVE NO CONNECTION WITH EACH OTHER)\r
+               //        WHERE W(X) IS THE WEIGHT FUNCTION\r
+               //\r
+               //                   W(X) = (1-X)**ALPHA * (1+X)**BETA\r
+               //\r
+               //        ON (-1, 1), ALPHA, BETA .GT. -1.\r
+               //\r
+               //     INPUT:\r
+               //\r
+               //        N        THE NUMBER OF POINTS USED FOR THE QUADRATURE RULE\r
+               //        ALPHA    SEE ABOVE\r
+               //        BETA     SEE ABOVE\r
+               //        B        REAL SCRATCH ARRAY OF LENGTH N\r
+               //\r
+               //     OUTPUT PARAMETERS (BOTH DOUBLE PRECISION ARRAYS OF LENGTH N)\r
+               //\r
+               //        T        WILL CONTAIN THE DESIRED NODES.\r
+               //        W        WILL CONTAIN THE DESIRED WEIGHTS W(J).\r
+               //\r
+               //     SUBROUTINES REQUIRED: CLASS, IMTQL2\r
+               //\r
+               //     REFERENCE:\r
+               //\r
+               //        THE ROUTINE HAS BEEN ADAPTED FROM THE MORE GENERAL\r
+               //        ROUTINE GAUSSQ BY GOLUB AND WELSCH.  SEE\r
+               //        GOLUB, G. H., AND WELSCH, J. H., "CALCULATION OF GAUSSIAN\r
+               //        QUADRATURE RULES," MATHEMATICS OF COMPUTATION 23 (APRIL,\r
+               //        1969), PP. 221-230.\r
+               \r
+               double MUZERO[] = new double[1];\r
+               int I;\r
+               int IERR[] = new int[1];\r
                \r
+               CLASS(N, ALPHA, BETA, B, T, MUZERO);\r
+        W[0] = 1.0;\r
+        for (I = 1; I < N; I++) {\r
+               W[I] = 0.0;\r
+        }\r
+        IMTQL2(N, T, B, W, IERR);\r
+        for (I = 0; I < N; I++) {\r
+               W[I] = MUZERO[0]*W[I]*W[I];\r
+        }\r
+        return;\r
        }\r
        \r
+       private void CLASS(int N, double ALPHA, double BETA, double B[], double A[], double MUZERO[]) {\r
+       //\r
+       //           THIS PROCEDURE SUPPLIES THE COEFFICIENTS A(J), B(J) OF THE\r
+       //        RECURRENCE RELATION\r
+       //\r
+       //             B P (X) = (X - A ) P   (X) - B   P   (X)\r
+       //              J J            J   J-1       J-1 J-2\r
+       //\r
+       //        FOR THE VARIOUS CLASSICAL (NORMALIZED) ORTHOGONAL POLYNOMIALS,\r
+       //        AND THE ZERO-TH MOMENT\r
+       //\r
+       //             MUZERO = INTEGRAL W(X) DX\r
+       //\r
+       //        OF THE GIVEN POLYNOMIAL'S WEIGHT FUNCTION W(X).  SINCE THE\r
+       //        POLYNOMIALS ARE ORTHONORMALIZED, THE TRIDIAGONAL MATRIX IS\r
+       //        GUARANTEED TO BE SYMMETRIC.\r
+       \r
+       //      DOUBLE PRECISION A(N),B(N)\r
+       \r
+             double A2B2,AB,ABI;\r
+             int I,NM1;\r
+       //\r
+       //     .. External Functions ..\r
+       //     double DGAMMA\r
+       //      EXTERNAL DGAMMA\r
+       \r
+             NM1 = N - 1;\r
+       //\r
+             AB = ALPHA + BETA;\r
+             ABI = 2.0 + AB;\r
+             MUZERO[0] = Math.pow(2.0, (AB+1.0))*DGAMMA(ALPHA+1.0)*\r
+                     DGAMMA(BETA+1.0)/DGAMMA(ABI);\r
+             A[0] = (BETA-ALPHA)/ABI;\r
+             B[0] = Math.sqrt(4.0* (1.0+ALPHA)* (1.0+BETA)/\r
+                   ((ABI+1.0)*ABI*ABI));\r
+             A2B2 = BETA*BETA - ALPHA*ALPHA;\r
+             for (I = 2; I <= NM1; I++) {\r
+                 ABI = 2.0*I + AB;\r
+                 A[I-1] = A2B2/ ((ABI-2.0)*ABI);\r
+                 B[I-1] = Math.sqrt(4.0*I* (I+ALPHA)* (I+BETA)* (I+AB)/\r
+                       ((ABI*ABI-1)*ABI*ABI));\r
+             } // for (I = 2; I <= NM1; I++)\r
+             ABI = 2.0*N + AB;\r
+             A[N-1] = A2B2/ ((ABI-2.0)*ABI);\r
+             return;\r
+       }\r
+       \r
+       private void IMTQL2(int N, double D[], double E[], double Z[], int IERR[]) {\r
+               // This is a modified version of the EISPACK routine IMTQL2.\r
+               // It finds the eigenvalues and first components of the eigenvectors\r
+               // of a symmetric tridiagonal matrix by the implicit QL method.\r
+               \r
+               // D(N), E(N), Z(N)\r
+               double B, C, F, G, P, R, S;\r
+               int I, II, K, L, M, MML;\r
+               int J = 0;\r
+               boolean zeroJ = true;\r
+               IERR[0] = 0;\r
+               if (N == 1) {\r
+                       return;\r
+               }\r
+               \r
+               E[N-1] = 0.0;\r
+               for (L = 1; L <= N; L++) {\r
+                       if (zeroJ) {\r
+                       J = 0;\r
+                       }\r
+                       else {\r
+                               zeroJ = true;\r
+                       }\r
+                   // Look for small sub-diagonal element\r
+                   for (M = L; M <= N; M++) {\r
+                       if (M == N) {\r
+                               break;\r
+                       }\r
+                       if (Math.abs(E[M-1]) <= MACHEP * (Math.abs(D[M-1] + Math.abs(D[M])))) {\r
+                               break;\r
+                       }\r
+                   } // for (M = L; M <= N; M++)\r
+                   \r
+                   P = D[L-1];\r
+                   if (M == L) {\r
+                       continue;\r
+                   }\r
+                   if (J == 30) {\r
+                       // Set error -- no convergence to an eigenvalue after 30 iterations.\r
+                       IERR[0] = L;\r
+                       return;\r
+                   }\r
+                   J = J+1;\r
+                   // Form shift\r
+                   G = (D[L] - P)/(2.0*E[L-1]);\r
+                   R = Math.sqrt(G*G + 1.0);\r
+                   if (G >= 0) {\r
+                       G = D[M-1] - P + E[L-1]/(G + Math.abs(R));\r
+                   }\r
+                   else {\r
+                       G = D[M-1] - P + E[L-1]/(G - Math.abs(R));      \r
+                   }\r
+                   S = 1.0;\r
+                   C = 1.0;\r
+                   P = 0.0;\r
+                   MML = M - L;\r
+                   // For I=M-1 step -1 until L do\r
+                   for (II = 1; II <= MML; II++) {\r
+                       I = M - II;\r
+                       F = S * E[I-1];\r
+                       B = C * E[I-1];\r
+                       if (Math.abs(F) >= Math.abs(G)) {\r
+                           C = G/F;\r
+                           R = Math.sqrt(C*C + 1.0);\r
+                           E[I] = F*R;\r
+                           S = 1.0/R;\r
+                           C = C*S;\r
+                       }\r
+                       else {\r
+                           S= F/G;\r
+                           R = Math.sqrt(S*S + 1.0);\r
+                           E[I] = G*R;\r
+                           C = 1.0/R;\r
+                           S = S * C;\r
+                       }\r
+                       G = D[I] - P;\r
+                       R = (D[I-1] - G)*S + 2.0*C*B;\r
+                       P = S*R;\r
+                       D[I] = G + P;\r
+                       G = C*R - B;\r
+                       // Form first component of vector\r
+                       F = Z[I];\r
+                       Z[I] = S * Z[I-1] + C * F;\r
+                       Z[I-1] = C * Z[I-1] - S * F;\r
+                   } // for (II = 1; II <= MML; II++)\r
+                   \r
+                   D[L-1] = D[L-1] - P;\r
+                   E[L-1] = G;\r
+                   E[M-1] = 0.0;\r
+                   zeroJ = false;\r
+               } // for (L = 1; L <= N; L++)\r
+               \r
+               // Order eigenvalues and eigenvectors\r
+               for (II = 2; II <= N; II++) {\r
+                   I = II - 1;\r
+                   K = I;\r
+                   P = D[I-1];\r
+                   \r
+                   for (J = II; J <= N; J++) {\r
+                       if (D[J-1] >= P) {\r
+                               continue;\r
+                       }\r
+                       K = J;\r
+                       P = D[J-1];\r
+                   } // for (J = II; J <= N; J++)\r
+                   \r
+                   if (K == I) {\r
+                       continue;\r
+                   }\r
+                   D[K-1] = D[I-1];\r
+                   D[I-1] = P;\r
+                   P = Z[I-1];\r
+                   Z[I-1] = Z[K-1];\r
+                   Z[K-1] = P;\r
+               } // for (II = 2; II <= N; II++)\r
+               return;\r
+       }\r
+\r
+       private double DGAMMA(double X) {\r
+               //\r
+               // COMPUTES THE GAMMA FUNCTION VIA A SERIES EXPANSION FROM\r
+               // ABRAMOWITZ & STEGUN\r
+               //\r
+               // L. N. TREFETHEN, 1/13/84\r
+               //\r
+           double FAC, G, Y;\r
+           int I, II;\r
+           double result;\r
+           \r
+           double C[] = new double[]{1.0000000000000000,.5772156649015329,\r
+                         -.6558780715202538,-.0420026350340952,\r
+                         .1665386113822915,-.0421977345555443,\r
+                         -.0096219715278770,.0072189432466630,\r
+                         -.0011651675918591,-.0002152416741149,\r
+                         .0001280502823882,-.0000201348547807,\r
+                         -.0000012504934821,.0000011330272320,\r
+                         -.0000002056338417,.0000000061160950,.0000000050020075,\r
+                         -.0000000011812746,.0000000001043427,.0000000000077823,\r
+                         -.0000000000036968,.0000000000005100,\r
+                         -.0000000000000206,-.0000000000000054,\r
+                         .0000000000000014,.0000000000000001};\r
+           \r
+           // ARGUMENT REDUCTION:\r
+               FAC = 1.0;\r
+               Y = X;\r
+               while (Y > 1.5) {\r
+                   Y = Y - 1.0;\r
+                   FAC = FAC*Y;\r
+               }\r
+\r
+            while (Y < 0.5) {\r
+               FAC = FAC/Y;\r
+               Y = Y + 1.0;\r
+            }\r
+         \r
+         // SERIES:\r
+            G = C[25];\r
+            for (I = 1; I <= 25; I++) {\r
+                   II = 26 - I;\r
+                   G = Y*G + C[II-1];\r
+            }\r
+         G = Y*G;\r
+         result = FAC/G;\r
+         return result;\r
+       }\r
 }
\ No newline at end of file