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

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

index c1b089c..fd90e54 100644 (file)
@@ -7,6 +7,7 @@ import java.util.Scanner;
 \r
 import gov.nih.mipav.model.structures.ModelImage;\r
 import gov.nih.mipav.view.MipavUtil;\r
+import gov.nih.mipav.view.Preferences;\r
 \r
 public class SymmsIntegralMapping extends AlgorithmBase  {\r
        String fileDir;\r
@@ -2804,44 +2805,90 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                //**** ASSIGN THE LOGICAL LINE SEGMENT TYPE TO EACH ARC.\r
                \r
-               /*      CALL LINSEG(LWORK(LNSEG),NARCS)\r
-               C\r
-               C**** LIST THE INPUT ARGUMENTS AND ASSOCIATED QUANTITIES\r
-               C\r
-                     CALL RSLT80(JBNM,HEAD,GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
-                    +INCST,INDEG,RFARC,RFARG,CENTR,RSNPH(JACIN),LWORK(LNSEG),\r
-                    +TSTNG,OULVL,IBNDS,MNEQN,OCH)\r
-                     PI=4E+0*ATAN(1E+0)\r
-                     RFARG=RFARG*PI\r
-               C\r
-               C**** SET UP THE GAUSS-JACOBI AND GAUSS-LEGENDRE QUADRATURE DATA AND \r
-               C**** STORE IN ARRAYS QUPTS AND QUWTS.  SET UP THREE TERM RECURRENCE\r
-               C**** COEFFICIENTS AND STORE IN ACOEF, BCOEF.  DETERMINE ZEROTH\r
-               C**** MOMENTS OF JACOBI DISTRIBUTIONS AND STORE IN H0VAL. \r
-               C**** ALSO SET UP THREE TERM RECURRENCE COEFFICIENTS AND ZEROTH MOMENTS\r
-               C**** FOR THE INTEGRATED POLYNOMIALS, STORING RESULTS IN AICOF,BICOF\r
-               C**** AND HIVAL.\r
-               C      \r
-                     CALL OPQUD1(NJIND,NQPTS,RSNPH(JACIN),RSNPH(ACOEF),RSNPH(BCOEF),\r
-                    +RSNPH(H0VAL),RSNPH(AICOF),RSNPH(BICOF),RSNPH(HIVAL),RSNPH(QUPTS),\r
-                    +RSNPH(QUWTS),RWORK(WORK),IER)\r
-                     IF (IER .GT. 0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-                     J=1-NQPTS\r
-                     DO 10 I=1,NJIND\r
-                       J=J+NQPTS\r
-                       RWORK(A1COF+I-1)=RSNPH(ACOEF+J-1)\r
-                       RWORK(B1COF+I-1)=RSNPH(BCOEF+J-1)\r
-               10    CONTINUE  \r
-                     WRITE(*,1) 'BASIC GAUSS QUADRATURE DATA DONE:'\r
-               C\r
-               C**** SET UP THE COEFFICIENTS IN THE THREE TERM RECURRENCE FORMULAE\r
-               C**** FOR THE PRINCIPAL SINGULAR INTEGRALS ASSOCIATED WITH THE VARIOUS\r
-               C**** JACOBI WEIGHT FUNCTIONS AND THEIR ORTHONORMAL POLYNOMIALS; STORE\r
-               C**** THESE COEFFICIENTS IN AQCOF, BQCOF AND CQCOF\r
-               C\r
-                     CALL ASQUC7(RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),RSNPH(JACIN),\r
+               boolean linout[] = new boolean[NARCS];\r
+               LINSEG(linout, NARCS);\r
+               for (I = 0; I < NARCS; I++) {\r
+                       LWORK[LNSEG-1+I]= linout[I];\r
+               }\r
+               //LINSEG(LWORK(LNSEG),NARCS);\r
+               \r
+               //**** LIST THE INPUT ARGUMENTS AND ASSOCIATED QUANTITIES\r
+               \r
+               double beta[] = new double[NARCS];\r
+               boolean linear[] = new boolean[NARCS];\r
+               for (I = 0; I < NARCS; I++) {\r
+                       beta[I] = RSNPH[JACIN-1+I];\r
+                       linear[I] = LWORK[LNSEG-1+I];\r
+               }\r
+               RSLT80(JBNM,HEAD,GSUPE,MAXER,GAQTL,INTER,NARCS,ORDSG,NQPTS,\r
+                    INCST,INDEG,RFARC,RFARG[0],CENTR,beta,linear,\r
+                    TSTNG,OULVL,IBNDS,MNEQN,OCH);\r
+               PI=Math.PI;\r
+               RFARG[0]=RFARG[0]*PI;\r
+               \r
+               //**** SET UP THE GAUSS-JACOBI AND GAUSS-LEGENDRE QUADRATURE DATA AND \r
+               //**** STORE IN ARRAYS QUPTS AND QUWTS.  SET UP THREE TERM RECURRENCE\r
+               //**** COEFFICIENTS AND STORE IN ACOEF, BCOEF.  DETERMINE ZEROTH\r
+               //**** MOMENTS OF JACOBI DISTRIBUTIONS AND STORE IN H0VAL. \r
+               //**** ALSO SET UP THREE TERM RECURRENCE COEFFICIENTS AND ZEROTH MOMENTS\r
+               //**** FOR THE INTEGRATED POLYNOMIALS, STORING RESULTS IN AICOF,BICOF\r
+               //**** AND HIVAL.\r
+                     \r
+               double JAC[] = new double[NJIND];\r
+               for (I = 0; I < NJIND; I++) {\r
+                       JAC[I] = RSNPH[JACIN+I-1];\r
+               }\r
+               double ACO[] = new double[NJIND*NQPTS];\r
+               double BCO[] = new double[NJIND*NQPTS];\r
+               double H0V[] = new double[NJIND];\r
+               double AIC[] = new double[NJIND*NQPTS+1];\r
+               double BIC[] = new double[NJIND*NQPTS+1];\r
+               double HIV[] = new double[NJIND];\r
+               double QUP[] = new double[NJIND*NQPTS];\r
+               double QUW[] = new double[NJIND*NQPTS];\r
+               double WOR[] = new double[NQPTS];\r
+               OPQUD1(NJIND,NQPTS,JAC,ACO,BCO,\r
+                    H0V,AIC,BIC,HIV,QUP,\r
+                    QUW,WOR,IER);\r
+               for (I = 0; I < NJIND*NQPTS; I++) {\r
+                       RSNPH[ACOEF+I-1] = ACO[I];\r
+                       RSNPH[BCOEF+I-1] = BCO[I];\r
+               }\r
+               for (I = 0; I < NJIND; I++) {\r
+                       RSNPH[H0VAL+I-1] = H0V[I];\r
+               }\r
+               for (I = 0; I < NJIND*NQPTS+1; I++) {\r
+                       RSNPH[AICOF+I-1] = AIC[I];\r
+                       RSNPH[BICOF+I-1] = BIC[I];\r
+               }\r
+               for (I = 0; I < NJIND; I++) {\r
+                       RSNPH[HIVAL+I-1] = HIV[I];\r
+               }\r
+               for (I = 0; I < NJIND*NQPTS; I++) {\r
+                       RSNPH[QUPTS+I-1] = QUP[I];\r
+                       RSNPH[QUWTS+I-1] = QUW[I];\r
+               }\r
+               for (I = 0; I < NQPTS; I++) {\r
+                       RWORK[WORK+I-1] = WOR[I];\r
+               }\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(1,0,IER[0],null);\r
+                   return;            \r
+               }\r
+               J=1-NQPTS;\r
+               for (I=1; I <= NJIND; I++) {\r
+                   J=J+NQPTS;\r
+                   RWORK[A1COF+I-2]=RSNPH[ACOEF+J-2];\r
+                   RWORK[B1COF+I-2]=RSNPH[BCOEF+J-2];\r
+               } // for (I=1; I <= NJIND; I++)\r
+               System.out.println("BASIC GAUSS QUADRATURE DATA DONE:");\r
+               \r
+               //**** SET UP THE COEFFICIENTS IN THE THREE TERM RECURRENCE FORMULAE\r
+               //**** FOR THE PRINCIPAL SINGULAR INTEGRALS ASSOCIATED WITH THE VARIOUS\r
+               //**** JACOBI WEIGHT FUNCTIONS AND THEIR ORTHONORMAL POLYNOMIALS; STORE\r
+               //**** THESE COEFFICIENTS IN AQCOF, BQCOF AND CQCOF\r
+               \r
+/*                   CALL ASQUC7(RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),RSNPH(JACIN),\r
                     +NJIND,NQPTS)\r
                      WRITE(*,1) 'DATA FOR SINGULAR INTEGRALS DONE:'\r
                C\r
@@ -3167,9 +3214,9 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        CALL OUPTPH(ISNPH,RSNPH,OCH)\r
                        CLOSE(OCH)\r
                      ENDIF\r
-               C  \r
-                     CALL WRTAIL(1,0,IER)\r
-               C */\r
+                 \r
+                     CALL WRTAIL(1,0,IER)*/\r
+               \r
     } // private void JAPHYC\r
 \r
     private void ANGLE7(double BE[], int NA, boolean IN) {\r
@@ -3272,6 +3319,518 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         } // for (K=1; K <= NA; K++)\r
 \r
     } // private void ANGLE7\r
+    \r
+    private void RSLT80(String JBNM, String HEAD, double SUPER, double MAXER, double AQTOL, boolean INTER, int NARCS, int ORDSG,\r
+        int NQPTS, boolean INCST, int INDEG, int RFARC, double RFARG,double CENTR[], double BETA[], boolean LINEAR[],\r
+        int TSTNG, int OULVL, int IBNDS[], int MNEQN,int OCH) {\r
+       \r
+       //INTEGER IBNDS(*)\r
+       //REAL BETA(*)\r
+       //COMPLEX CENTR\r
+       //LOGICAL LINEAR(*)\r
+       //CHARACTER JBNM*4,HEAD*72\r
+       \r
+       //**** WRITE THE MAIN ARGUMENTS OF JAPHYC AND ASSOCIATED QUANTITIES ON  \r
+       //**** THE LISTING FILE.\r
+       \r
+       // LOCAL VARIABLES\r
+       \r
+       int I;\r
+       final String TXT1 = "REQUESTED ACCURACY UNREALISTIC";\r
+\r
+\r
+       Preferences.debug("HEAD\n",Preferences.DEBUG_ALGORITHM);\r
+        Preferences.debug("JOB NAME : " + JBNM + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+       if (INTER) {\r
+           Preferences.debug("INTERIOR DOMAIN WITH " + NARCS + " ARCS\n",Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       else {\r
+           Preferences.debug("EXTERIOR DOMAIN WITH " + NARCS + " ARCS\n",Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       if (ORDSG > 1) {\r
+           Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+           Preferences.debug("ORDER OF SYMMETRY GROUP IS : " + ORDSG + "\n",Preferences.DEBUG_ALGORITHM);\r
+           Preferences.debug("NUMBER OF ARCS ON FBS IS   : " + (NARCS/ORDSG)  + "\n",Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       \r
+       Preferences.debug("ACCURACY REQUESTED            : " + MAXER  + "\n",Preferences.DEBUG_ALGORITHM);\r
+       if (MAXER < SUPER) {\r
+           Preferences.debug("WORKING ACCURACY              : " + SUPER + "  " + TXT1 + "\n",Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       Preferences.debug("ABSOLUTE QUADRATURE TOLERENCE : " + AQTOL + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM NUMBER OF SUBARCS           : " + IBNDS[0] + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM NUMBER OF EQUATIONS         : " + MNEQN + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM NUMBER OF QUADRATURE PANELS : " + (IBNDS[2]-1) + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM TOTAL  OF QUADRATURE POINTS : " + IBNDS[3] + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MINIMUM NUMBER OF QUADRATURE POINTS : " + NQPTS + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("MAXIMUM DEGREE OF POLYNOMIAL        : " + (NQPTS-1) + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("INITIAL DEGREE OF POLYNOMIAL        : " + INDEG + (NQPTS-1) + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("INCREMENTAL STRATEGY                : " + INCST + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("REFERENCE ARC         : " + RFARC + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("REFERENCE ARGUMENT/PI : " + RFARG + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("CENTRE POINT          : " + CENTR[0] + ", " + CENTR[1] + "\n",Preferences.DEBUG_ALGORITHM);\r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("CORNER        ANGLE/PI          JACOBI INDEX      LINEAR\n",Preferences.DEBUG_ALGORITHM);\r
+       for (I = 1; I <= NARCS; I++) {\r
+               Preferences.debug("  " + I + "     " + (1.0/(1.0 + BETA[I-1])) + "    " + BETA[I-1] + "      " + LINEAR[I-1] + "\n",Preferences.DEBUG_ALGORITHM);\r
+       }\r
+       \r
+       \r
+       Preferences.debug("\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("TESTING LEVEL : " + TSTNG + "\n",Preferences.DEBUG_ALGORITHM);\r
+       Preferences.debug("OUTPUT  LEVEL : " + OULVL + "\n",Preferences.DEBUG_ALGORITHM);\r
+    } // private void RSLT80\r
+    \r
+    private void OPQUD1(int NJIND, int NQPTS,double JACIN[], double ACOEF[], double BCOEF[],\r
+        double H0VAL[], double AICOF[], double BICOF[], double HIVAL[], double QUPTS[], \r
+        double QUWTS[], double WORK[], int IER[]) {\r
+       //REAL JACIN(*),ACOEF(*),BCOEF(*),H0VAL(*),QUPTS(*),QUWTS(*),\r
+       //+WORK(*),AICOF(*),BICOF(*),HIVAL(*)\r
+       \r
+       //**** TO SET UP THE THREE TERM RECURRENCE COEFFICIENTS (ACOEF AND BCOEF)\r
+       //**** FOR THE ON JACOBI POLYNOMIALS (UP TO DEGREE NQPTS) ASSOCIATED WITH\r
+       //**** THE JACOBI INDECES GIVEN IN JACIN AND TO STORE THE ZEROTH MOMENTS\r
+       //**** OF THE JACOBI DISTRIBUTIONS IN H0VAL.  ALSO TO REPEAT THESE\r
+       //**** CALCULATIONS FOR THE INCREMENTED JACOBI INDECES ARISING IN THE\r
+       //**** EXPRESSIONS FOR THE BOUNDARY CORRESPONDENCE FUNCTION, STORING\r
+       //**** RESULTS IN AICOF, BICOF AND HIVAL.\r
+       \r
+       //**** ALSO TO SET UP THE NQPT POINT GAUSS JACOBI QUADRATURE RULES,\r
+       //**** STORING THE ABSCISSAE IN QUPTS AND THE WEIGHTS IN QUWTS.\r
+       \r
+       //**** IER=0 - NORMAL EXIT\r
+       //**** IER=4 - FAILURE TO CONVERGE IN EIGSYS; CAN'T SET UP BASIC\r
+       //****         QUADRATURE RULES \r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       int I,J,K,LOSUB,M;\r
+       int IFAIL[] = new int[1];\r
+       double BETA,BETA1,C,PROD,S,T;\r
+       //EXTERNAL ASONJ7,EIGSYS\r
+       \r
+       for (I=1; I <= NJIND; I++) {\r
+           BETA=JACIN[I-1];\r
+           BETA1=BETA+1.0;\r
+           PROD=BETA*BETA;\r
+       \r
+           // CALCULATE THE ZEROTH MOMENT FOR THIS BETA\r
+       \r
+           H0VAL[I-1]=Math.pow(2.0,BETA1)/BETA1;\r
+       \r
+           // START ON THE 3-TERM ORTHONORMAL RECURRENCE COEFFICIENTS ACOEF\r
+           // AND BCOEF FOR THIS BETA\r
+       \r
+           T=2.0+BETA;\r
+           S=T*T;\r
+           C=4.0*BETA1/S/(T+1.0);\r
+           LOSUB=(I-1)*NQPTS+1;\r
+           ACOEF[LOSUB-1]=Math.sqrt(C);\r
+           BCOEF[LOSUB-1]=BETA/T;\r
+       \r
+           for (K=2; K <= NQPTS; K++) {\r
+               J=LOSUB+K-1;\r
+               BCOEF[J-1]=PROD/T/(T+2.0);\r
+               T=2.0*K+BETA;\r
+               S=T*T;\r
+               C=4.0*K*K*(BETA+K)*(BETA+K)/S/(S-1.0);\r
+               ACOEF[J-1]=Math.sqrt(C);\r
+           } // for (K=2; K <= NQPTS; K++)\r
+       \r
+           // START ON THE QUADRATURE POINTS AND WEIGHTS FOR THIS BETA\r
+       \r
+           for (K=1; K <= NQPTS; K++) {\r
+               J=LOSUB+K-1;\r
+               QUPTS[J-1]=BCOEF[J-1];\r
+               QUWTS[J-1]=ACOEF[J-1];\r
+               WORK[K-1]=0.0;\r
+           } // for (K=1; K <= NQPTS; K++)\r
+           WORK[0]=1.0;\r
+       \r
+           // AT THIS POINT THE LOCAL SEGMENTS OF QUPTS AND QUWTS ARE THE \r
+           // DIAGONAL AND SUBDIAGONAL OF A SYMMETRIC TRIDIAGONAL MATRIX \r
+           // WHOSE EIGENVALUES ARE THE QUADRATURE POINTS AND WHOSE \r
+           // EIGENVECTORS GIVE THE QUADRATURE WEIGHTS.\r
+       \r
+           double D[] = new double[NQPTS];\r
+           double E[] = new double[NQPTS];\r
+           for (M = 0; M < NQPTS; M++) {\r
+               D[M] = QUPTS[LOSUB+M-1];\r
+               E[M] = QUWTS[LOSUB+M-1];\r
+           }\r
+           EIGSYS(NQPTS,EPS,D,E,WORK,IFAIL);\r
+           for (M = 0; M < NQPTS; M++) {\r
+               QUPTS[LOSUB+M-1] = D[M];\r
+               QUWTS[LOSUB+M-1] = E[M];\r
+           }\r
+           if (IFAIL[0] > 0) {\r
+               IER[0]=4;\r
+               return;\r
+           }\r
+       \r
+           for (K=1; K <= NQPTS; K++) {\r
+               QUWTS[LOSUB+K-2]=WORK[K-1]*WORK[K-1]*H0VAL[I-1];\r
+           } // for (K=1; K <= NQPTS; K++)\r
+       \r
+           // SET UP THE THREE TERM RECURRENCE COEFFICIENTS AIVAL,BIVAL AND\r
+           // THE ZEROTH MOMENT HIVAL FOR THE INTEGRATED POLYNOMIALS.\r
+       \r
+           double A[] = new double[NQPTS];\r
+           double B[] = new double[NQPTS];\r
+           double H[] = new double[1];\r
+           ASONJ7(1.0,BETA1,A,B,H,NQPTS);\r
+           for (M = 0; M < NQPTS; M++) {\r
+               AICOF[LOSUB+M-1] = A[M];\r
+               BICOF[LOSUB+M-1] = B[M];\r
+           }\r
+           HIVAL[I-1] = H[0];\r
+       \r
+       } // for (I=1; I <= NJIND; I++)\r
+       \r
+       //     NORMAL TERMINATION\r
+       \r
+       IER[0]=0;\r
+       \r
+    } // private void OPQUD1\r
+    \r
+    private void EIGSYS(int N, double EPS, double D[], double E[],\r
+               double Z[],int IER[]) { \r
+        // REAL D(N),E(N),Z(N)      \r
+      \r
+        // THIS IS A MODIFIED VERSION OF THE EISPACK ROUTINE IMTQL2.     \r
+        // IT FINDS THE EIGENVALUES AND FIRST COMPONENTS OF THE\r
+        // EIGENVECTORS OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL\r
+        // METHOD.     \r
+\r
+        // **** LOCAL VARIABLES\r
+\r
+        int L,J,M ,MML,II,K,I;\r
+        double P,G,R,S,C,F,B;  \r
+  \r
+        IER[0] = 0;    \r
+\r
+        if (N == 1) {\r
+               return;\r
+        }\r
+\r
+        E[N-1] = 0.0;  \r
+        loopL: for (L = 1; L <= N; L++) {       \r
+            J = 0;    \r
+\r
+            //****    LOOK FOR SMALL SUB-DIAGONAL ELEMENT     \r
+\r
+            while (true) {\r
+                       for (M = L; M <= N; M++) {    \r
+                           if (M == N) break;   \r
+                           if (Math.abs(E[M-1]) <= EPS * (Math.abs(D[M-1]) + Math.abs(D[M])))  break;\r
+                       } // for (M = L; M <= N; M++) \r
+               \r
+                      P = D[L-1]; \r
+                      if (M == L) continue loopL;      \r
+                      if (J == 30) {\r
+                          // **** SET ERROR -- NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS\r
+                          IER[0] = L;\r
+                          return;\r
+                      }\r
+                      J = J + 1;\r
+               \r
+                      // ****    FORM SHIFT\r
+               \r
+                      G = (D[L] - P) / (2. * E[L-1]);\r
+                      R = Math.sqrt(G*G+1.0);\r
+                      double sig;\r
+                      if (G >= 0) {\r
+                          sig = Math.abs(R);\r
+                      }\r
+                      else {\r
+                          sig = -Math.abs(R);\r
+                      }\r
+                      G = D[M-1] - P + E[L-1] / (G + sig); \r
+                      S = 1.0;  \r
+                      C = 1.0 ;  \r
+                      P = 0.0;  \r
+                      MML = M - L;\r
+                      \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
+                          } // if (Math.abs(F) >= Math.abs(G))\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
+                          } // else\r
+                          G = D[I] - P;  \r
+                          R = (D[I-1] - G) * S + 2. * C * B;     \r
+                          P = S * R;       \r
+                          D[I] = G + P;  \r
+                          G = C * R - B;  \r
+                \r
+                          //****       FORM FIRST COMPONENT OF VECTOR\r
+               \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
+            } // while (true)\r
+        } // loopL: for (L = 1; L <= N; L++)   \r
+       \r
+        // **** ORDER EIGENVALUES AND EIGENVECTORS      \r
+\r
+        for (II = 2; II <= N; II++) {      \r
+            I = II - 1; \r
+            K = I;    \r
+            P = D[I-1]; \r
+            for (J = II; J <= N; J++) {   \r
+                if (D[J-1] >= P) continue;\r
+                K = J; \r
+                P = D[J-1];\r
+            } // for (J = II; J <= N; J++)\r
+       \r
+            if (K == I) continue;      \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
+    } // private void EIGSYS\r
+\r
+    private void ASONJ7(double ALFA,double BETA, double A[], double B[],\r
+               double H[], int N) {\r
+        //REAL A(*),B(*),H,ALFA,BETA\r
+\r
+        // ..TO ASSIGN THE COEFFICIENTS A(K) AND B(K) , K=1(1)N, IN THE\r
+        // ..3-TERM RECURRENCE FORMULA FOR THE ORTHONORMAL JACOBI POLYNOMIALS\r
+        // ..WHERE\r
+        // ..\r
+        // ..     A(K)P (X) = (X - B(K))P   (X) - A(K-1)P   (X) , K=1,2,..,N,\r
+        // ..          K                 K-1             K-2\r
+        // ..\r
+        // ..         P  (X) = 0 , P (X) = 1/SQRT(H)\r
+        // ..          -1           0 \r
+        // .. \r
+        // ..AND H IS THE ZEROTH MOMENT OF THE JACOBI WEIGHT FUNCTION\r
+        // ..(1-X)**ALFA*(1+X)**BETA ON [-1,1].\r
+\r
+        // **** AUTHOR: DAVID HOUGH\r
+        // **** LAST UPDATE: 15.09.89\r
+\r
+        // **** ..LOCAL VARIABLES..\r
+        double SUM,DIFF,PROD,TC,T,SC,S,C;\r
+        int K;\r
+        // EXTERNAL GAMMA\r
+        // double GAMMA;\r
+\r
+        SUM=ALFA+BETA;\r
+        DIFF=BETA-ALFA;\r
+        PROD=SUM*DIFF;\r
+\r
+        // ..CALCULATE H.\r
+        TC=SUM+1.0;\r
+        SC=Math.pow(2.0,TC);\r
+        S=GAMMA(ALFA+1.0);\r
+        T=GAMMA(BETA+1.0);\r
+        C=GAMMA(TC+1.0);\r
+        H[0]=SC*S*T/C;\r
+\r
+        // ..START ON A,B ARRAYS.\r
+        if (N > 0) {\r
+            T=2.0+SUM;\r
+            S=T*T;\r
+            C=4.0*(ALFA+1.0)*(BETA+1.0)/S/(T+1.0);\r
+            A[0]=Math.sqrt(C);\r
+            B[0]=DIFF/T;\r
+\r
+            for (K=2; K<= N; K++) {\r
+                B[K-1]=PROD/T/(T+2.0);\r
+                T=2.0*K+SUM;\r
+                S=T*T;\r
+                C=4.0*K*(ALFA+K)*(BETA+K)*(SUM+K)/S/(S-1.0);\r
+                A[K-1]=Math.sqrt(C);\r
+            } // for (K=2; K<= N; K++) \r
+        } // if (N > 0)\r
+\r
+    } // private void ASONJ7\r
+    \r
+    private double GAMMA( double U) {\r
+\r
+        // TO COMPUTE THE GAMMA FUNCTION FOR REAL ARGUMENT U BY USING\r
+        // THE CHEBYSHEV EXPANSION GIVEN IN TABLE 1.3 OF "MATHEMATICAL\r
+        // FUNCTIONS AND THEIR APPROXIMATION" BY Y.L. LUKE ,ACADEMIC PRESS,\r
+        // NEW YORK, 1975.\r
+        // SINCE GAMMA HAS POLES AT U=0,-1,-2,-3,... DIVISION BY ZERO WILL\r
+        // OCCUR FOR THESE ARGUMENT VALUES.\r
+\r
+        // LOCAL VARIABLES\r
+\r
+        int N;\r
+        double UWORK,FACTOR;\r
+        double B0 = 0.0;\r
+        double X,B1,B2;\r
+        double A[] = new double[]{3.65738772508338243850,\r
+                              1.95754345666126826928,\r
+                              0.33829711382616038916,0.4208951276557549199E-1,\r
+                              0.42876504821290877E-2,0.36521216929461767E-3,\r
+                              0.27400642226422E-4,0.181240233365124E-5,\r
+                              0.10965775865997E-6,0.598718404552E-8,\r
+                              0.30769080535E-9,0.143179303E-10,\r
+                              0.65108773E-12,0.259585E-13,0.110789E-14,\r
+                              0.3547E-16,0.169E-17,0.3E-19};  \r
+        double result;\r
+\r
+        UWORK=U;\r
+        FACTOR=1.0;\r
+\r
+           while (true) {\r
+               if (UWORK > 4.0) {\r
+                   UWORK=UWORK-1.0;\r
+                   FACTOR=FACTOR*UWORK;\r
+                   continue;\r
+               }\r
+               else if (UWORK < 3.0) {\r
+                   FACTOR=FACTOR/UWORK;\r
+                   UWORK=UWORK+1.0;\r
+                   continue;\r
+               }\r
+               else {\r
+                       break;\r
+               }\r
+           } // while (true)\r
+       \r
+           X=UWORK-3.0;      \r
+           X=4.0*X-2.0;\r
+           B2=0.0;\r
+           B1=0.0;\r
+           for (N=17; N >= 0; N--) {\r
+               B0=X*B1-B2+A[N];\r
+               if (N > 0) {\r
+                   B2=B1;\r
+                   B1=B0;\r
+               }\r
+           } // for (N=17; N >= 0; N--)\r
+       \r
+           result=5E-1*(B0+A[0]-B2)*FACTOR;\r
+           return result;\r
+    } // private double GAMMA\r
+\r
+\r
+            /* SUBROUTINE RSLT71(QIERC,RCOND,SOLUN,NEQNS,LOSUB,HISUB,COLSC,\r
+            +NQPTS,JATYP,PARNT,TNSUA,INTER,MQERR,MCQER,CINFN,ACTIN,\r
+            +NEWDG,NJIND,JACIN,NQUAD,TOLOU,LGTOL,SOLCO,OUCH1)\r
+             INTEGER NEQNS,TNSUA,OUCH1,NQPTS,NJIND,NEWDG(*),NQUAD(*),LOSUB(*),\r
+            +HISUB(*),QIERC(0:6),JATYP(*),PARNT(*),ACTIN(*),SOLCO\r
+             REAL SOLUN(*),RCOND,COLSC(*),MQERR,MCQER,LGTOL,\r
+            +CINFN(*),JACIN(*),TOLOU(*)\r
+             LOGICAL INTER\r
+             CHARACTER QTEXT(0:6)*22,LINE*72\r
+             PARAMETER (LINE='_________________________________________________\r
+            +________________')\r
+       C \r
+       C     LOCAL VARIABLES\r
+       C\r
+             INTEGER I,J,JI,K,L,LOD,N,H\r
+             REAL S,CAP\r
+       C\r
+             QTEXT(0)='...........NORMAL EXIT'\r
+             QTEXT(1)='.....MAX. SUBDIVISIONS'\r
+             QTEXT(2)='....ROUNDOFF DETECTION'\r
+             QTEXT(3)='.........BAD INTEGRAND'\r
+             QTEXT(6)='.........INVALID INPUT'\r
+       C\r
+             WRITE(OUCH1,*) LINE\r
+             WRITE(OUCH1,*) '             SOLUTION NUMBER =',SOLCO\r
+             WRITE(OUCH1,*) '                       NEQNS =',NEQNS \r
+             WRITE(OUCH1,*) 'RECIPROCAL COND NO. ESTIMATE =',RCOND\r
+             WRITE(OUCH1,*) '   CONDITION NO. LOWER BOUND =',1E+0/RCOND\r
+       C\r
+             WRITE(OUCH1,*) \r
+             WRITE(OUCH1,997) 'JACOBI INDEX','POINTS','TOLERANCE ACHIEVED'\r
+             DO 10 I=1,NJIND\r
+               WRITE(OUCH1,998) I,NQUAD(I),TOLOU(I)\r
+       10    CONTINUE\r
+       C\r
+             WRITE(OUCH1,*) \r
+             WRITE(OUCH1,*) 'QAWS TERMINATIONS WITH......' \r
+             DO 20 I=0,6\r
+               IF (QIERC(I) .GT. 0) THEN\r
+                 WRITE(OUCH1,1000) QTEXT(I),QIERC(I)\r
+               ENDIF\r
+       20    CONTINUE\r
+       C\r
+             WRITE(OUCH1,*) \r
+             WRITE(OUCH1,999) '              MAXIMUM QAWS ERROR =',MQERR\r
+             WRITE(OUCH1,999) 'MAXIMUM COMPOSITE GAUSSIAN ERROR =',MCQER\r
+             WRITE(OUCH1,*) \r
+             DO 40 I=1,TNSUA\r
+                 WRITE(OUCH1,*)\r
+                 WRITE(OUCH1,*) 'SUB ARC =',I,' ON PARENT ARC',PARNT(I)\r
+                 WRITE(OUCH1,990) 'N','SCALED SOLUN','UNSCALED SOLUN','IGNORE L\r
+            +EVEL'\r
+                 L=LOSUB(I)\r
+                 H=HISUB(I)\r
+                 JI=ABS(JATYP(I))\r
+                 LOD=(JI-1)*NQPTS+1\r
+                 DO 30 J=L,H\r
+                     N=J-L\r
+                     K=LOD+N\r
+                     S=SOLUN(J)\r
+                     WRITE(OUCH1,991) N,S,S*COLSC(K),LGTOL/CINFN(J)\r
+       30        CONTINUE\r
+                 IF (ACTIN(I) .EQ. -1) THEN\r
+                     WRITE(OUCH1,*)'ACTION: REDUCE DEGREE TO ',NEWDG(I),' ***'\r
+                 ELSE IF (ACTIN(I) .EQ. 0) THEN\r
+                     WRITE(OUCH1,*)'ACTION: NONE            ***'\r
+                 ELSE IF (ACTIN(I) .EQ. 1) THEN\r
+                     WRITE(OUCH1,*)'ACTION: INCREASE DEGREE TO ',NEWDG(I)\r
+                 ELSE\r
+                     WRITE(OUCH1,*)'ACTION: SUBDIVIDE THIS ARC'\r
+                 ENDIF\r
+       40    CONTINUE\r
+       C\r
+             WRITE(OUCH1,*) 'KAPPA =',SOLUN(NEQNS)\r
+             IF (.NOT.INTER) THEN\r
+                 CAP=EXP(-SOLUN(NEQNS))\r
+                 WRITE(OUCH1,*) 'CAPACITY = ',CAP\r
+             ENDIF\r
+       C\r
+       990   FORMAT(A,T7,A,T26,A,T44,A)\r
+       991   FORMAT(I3,T6,E15.8,T25,E15.8,T44,E10.3)\r
+       992   FORMAT(E15.8)\r
+       993   FORMAT(I3,T8,E15.8,'  (',E14.7,',',E14.7,')')\r
+       994   FORMAT(A,T8,A,T34,A)\r
+       995   FORMAT(A,T6,A,T23,A,T36,A)\r
+       996   FORMAT(I2,T6,E14.7,T23,F10.5,T36,E14.7)\r
+       997   FORMAT(A,T24,A,T40,A)\r
+       998   FORMAT(T5,I3,T26,I3,T45,E9.2)\r
+       999   FORMAT(A,E10.2)\r
+       1000  FORMAT(A,1X,I5)\r
+       C\r
+    }*/\r
+\r
  \r
       /**\r
        * zabs computes the absolute value or magnitude of a double precision complex variable zr + j*zi.\r