Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 8 Jan 2018 23:15:06 +0000 (23:15 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Mon, 8 Jan 2018 23:15:06 +0000 (23:15 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15336 ba61647d-9d00-f842-95cd-605cb4296b96

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

index bc95cf6..69d4b93 100644 (file)
@@ -2679,6 +2679,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                     MCHEP,MXERM,PI,RQTOL,SSUPE,\r
                     TGTER,TOLNR;\r
                double SOLUNB[][] = new double[MNEQN][1];\r
+               double A[][] = new double[MNEQN][MNEQN];\r
                \r
                double ZMXER[] = new double[2];\r
                //COMPLEX ZMXER\r
@@ -2937,7 +2938,6 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    //**** SOLVE LINEAR SYSTEM BY GAUSSIAN ELIMINATION USING LINPACK\r
                    // LAPACK equivalent of LINPACK SGECO are DLANGE, DGETRF, and DGECON\r
                    // LU factorization and condition estimation of a general matrix.\r
-                   double A[][] = new double[MNEQN][MNEQN];\r
                    for (I = 0; I < MNEQN; I++) {\r
                        for (J = 0; J < MNEQN; J++) {\r
                                A[I][J] = MATRX[I][J][1];\r
@@ -3265,36 +3265,59 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                      } // else !GACPT \r
            } // while (true)\r
                \r
-               /*      IF (OULVL.EQ.2 .OR. OULVL.EQ.3) THEN\r
-                       CALL RSLT71(QIERC,RCOND,RSNPH(SOLUN),NEQNS,ISNPH(LOSUB),\r
-                    +  IWORK(HISUB),RWORK(COLSC),NQPTS,ISNPH(JATYP),IGEOM(PARNT),TNSUA,\r
-                    +  INTER,MQERR,MCQER,RWORK(WORK2),IWORK(AXION),IWORK(NEWDG),NJIND,\r
-                    +  RSNPH(JACIN),IWORK(NQUAD),RWORK(TOLOU),LGTOL,SOLCO,OCH)\r
-                     ENDIF\r
-               C\r
-               C**** ESTIMATE MAXIMUM ERROR IN MODULUS.\r
-               C\r
-                       WRITE(*,*)\r
-                       WRITE(*,1) 'ERRORS IN MODULUS STARTED:'\r
-               C\r
-                     CALL TSJAC3(IWORK(LOTES),IWORK(HITES),RWORK(COLPR),ZWORK(ZCOLL),\r
-                    +NQPTS,NTEST,ORDSG,TNSUA,TSTNG,ISNPH(DGPOL),ISNPH(JATYP),\r
-                    +IGEOM(PARNT),RSNPH(AICOF),RSNPH(BICOF),RWORK(DIAG),RGEOM(HALEN),\r
-                    +RSNPH(JACIN),RGEOM(MIDPT),RWORK(SDIAG),IER)\r
-                     IF (IER .GT. 0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-               C\r
-                     CALL TESMD9(RWORK(WORK2),MATRX(1,1,2),RSNPH(SOLUN),MNEQN,NCOLL,\r
-                    +NTEST,NQPTS,TNSUA,ISNPH(JATYP),IGEOM(PARNT),ISNPH(DGPOL),\r
-                    +ISNPH(LOSUB),IWORK(HISUB),IWORK(LOTES),IWORK(HITES),IWORK(NQUAD),\r
-                    +IWORK(LOQSB),TOLNR,RGEOM(MIDPT),RGEOM(HALEN),RSNPH(H0VAL),\r
-                    +RWORK(COLSC),RSNPH(ACOEF),RSNPH(BCOEF),RWORK(COLPR),RWORK(QCOMX),\r
-                    +RWORK(QCOMW),CENTR,ZWORK(ZCOLL),INTER,LWORK(LNSEG),RWORK(WORK),\r
-                    +QIERR,MQERR,RSNPH(JACIN),RWORK(A1COF),RWORK(B1COF),AQTOL,RQTOL,\r
-                    +RWORK(AQCOF),RWORK(BQCOF),RWORK(CQCOF),MXERM,IMXER,ZMXER,\r
-                    +RSNPH(ERARC),ORDSG,REFLN)\r
-               C\r
+               if (OULVL == 2 || OULVL == 3) {\r
+                   RSLT71(QIERC,RCOND[0],SOLUN,NEQNS,LOSUB,\r
+                    HISUB,COLSC,NQPTS,JATYP,PARNT,TNSUA[0],\r
+                    INTER,MQERR[0],MCQER[0],WORK2,AXION,NEWDG,NJIND,\r
+                    JACIN,NQUAD,TOLOU,LGTOL,SOLCO);\r
+               } // if (OULVL == 2 || OULVL == 3)\r
+               \r
+               // ESTIMATE MAXIMUM ERROR IN MODULUS.\r
+               \r
+               System.out.println();\r
+               System.out.println("ERRORS IN MODULUS STARTED:");\r
+               \r
+               /*TSJAC3(LOTES,HITES,COLPR,ZCOLL,\r
+                 NQPTS,NTEST,ORDSG,TNSUA[0],TSTNG,DGPOL,JATYP,\r
+                 PARNT,AICOF,BICOF,DIAG,HALEN,\r
+                 JACIN,MIDPT,SDIAG,IER);\r
+                 if (IER[0] > 0) {\r
+                         if (SOLCO >= 1) {\r
+                                       \r
+                                   // ****   COMPUTE THE BOUNDARY CORRESPONDENCE COEFFICIENTS BCFSN AND THE\r
+                                   // ****   ARGUMENTS OF ALL SUBARC END POINTS ON THE UNIT DISC,\r
+                                   // ****   AS REQUIRED BY SUBSEQUENT PROCESSING ROUTINES.\r
+                               \r
+                                   BCFVTF(BCFSN,VTARG,DGPOL,JATYP,\r
+                                       LOSUB,PARNT,RFARC,TNSUA[0],H0VAL,JACIN,\r
+                                       RFARG[0],SOLUN);\r
+                               \r
+                                   // ****   OUTPUT DATA REQUIRED FOR POST-PROCESSING.\r
+                               \r
+                                   IGEOM[2]=TNSUA[0];\r
+                                   ISNPH[2]=TNSUA[0];\r
+                                   ISNPH[3]=NEQNS;\r
+                       } // if (SOLCO >= 1)\r
+                                 \r
+                       WRTAIL(1,0,IER[0],null);\r
+                           return;     \r
+                 } // if (IER[0] > 0)\r
+               \r
+                 for (I = 0; I < MNEQN; I++) {\r
+                         for (J = 0; J < MNEQN; J++) {\r
+                                 A[I][J] = MATRX[I][J][1];\r
+                         }\r
+                 }\r
+                 TESMD9(WORK2,A,SOLUN,MNEQN,NCOLL,\r
+                    NTEST,NQPTS,TNSUA[0],JATYP,PARNT,DGPOL,\r
+                    LOSUB,HISUB,LOTES,HITES,NQUAD,\r
+                    LOQSB,TOLNR,MIDPT,HALEN,H0VAL,\r
+                    COLSC,ACOEF,BCOEF,COLPR,QCOMX,\r
+                    QCOMW,CENTR,ZCOLL,INTER,LNSEG,WORK,\r
+                    QIERR,MQERR,JACIN,A1COF,B1COF,AQTOL,RQTOL,\r
+                    AQCOF,BQCOF,CQCOF,MXERM,IMXER,ZMXER,\r
+                    ERARC,ORDSG,REFLN);\r
+               \r
                      WRITE(*,1) 'ERRORS IN MODULUS DONE:'\r
                      WRITE(*,3) 'MAXIMUM ERROR AT TEST POINTS:',MXERM\r
                      DO 60 I=0,6\r
@@ -3344,8 +3367,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                        CLOSE(OCH)\r
                      ENDIF\r
                  \r
-                     CALL WRTAIL(1,0,IER)\r
-               */\r
+                     CALL WRTAIL(1,0,IER)*/\r
+               \r
     } // private void JAPHYC\r
 \r
     private void ANGLE7(double BE[], int NA, boolean IN) {\r
@@ -5302,25 +5325,6 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
     \r
        //   INTEGRAL  [(1+X)**BETA*P(X,I)*LOG(ZZ(J)-X)*dX], I=0,1,...,MAXDG\r
        // TAU1<=X<=TAU2                                     J=1,2\r
-       \r
-       //     WHERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I\r
-       //     ASSOCIATED WITH THE WEIGHT (1+X)**BETA AND -1<TAU1<TAU2<=1.\r
-       //     THE REMAINDER CORRESPONDING TO P(.,I) AND ZZ(J) IS STORED IN \r
-       //     REMND(I+J+MAXDG*(J-1)), I=0,1,...,MAXDG, J=1,2.  THIS ROUTINE USES\r
-       //     THE SIMPLEST POSSIBLE ESTIMATES; I.E. THE LEADING TERM ONLY IN\r
-       //     THE ASYMPTOTIC EXPANSION AND THE WATSON-DOETSCH ESTIMATE FOR ANY\r
-       //     INTEGRALS.\r
-       \r
-       //     THE PURPOSE OF THIS ROUTINE IS TO DETERMINE A VALUE FOR EITHER\r
-       //     TAU2 (TAU1 REMAINS FIXED IF T1FXD IS "TRUE") OR TAU1 (TAU2 REMAINS\r
-       //     FIXED IF T1FXD IS "FALSE") SUCH THAT\r
-       \r
-       //         ABS( REAL(REMND(I)) )*CSCAL(I) < TOL , I=1,2*MAXDG+2\r
-       \r
-       //     AND THAT, IF POSSIBLE, \r
-       \r
-       //        0.5*TOL <= ABS( REAL(REMND(I)) )*CSCAL(I) < TOL\r
-       \r
        //     FOR AT LEAST ONE VALUE OF I.\r
        //     IER=0 - NORMAL EXIT\r
        //     IER=12- LOCAL PARAMETER NC NEEDS INCREASING TO AT LEAST NZZ\r
@@ -5882,7 +5886,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    ABER0 = qmod.getAbserr();\r
                    NVAL = qmod.getNeval();\r
                    IER[0] = qmod.getErrorStatus();\r
-                   QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                   QIERR[IER[0]]=QIERR[IER[0]]+1;\r
                    TYPE=2;\r
                    qmod = new qaws(TC,1.0,routine,0.0,0.0,2,AQTOL,RQTOL,limit);\r
                    qmod.driver();\r
@@ -5890,7 +5894,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    ABER1 = qmod.getAbserr();\r
                    NVAL = qmod.getNeval();\r
                    IER[0] = qmod.getErrorStatus();\r
-                   QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                   QIERR[IER[0]]=QIERR[IER[0]]+1;\r
                    WORK[0]=U0+U1;\r
                    MQERR[0]=Math.max(MQERR[0],ABER0+ABER1);\r
                    if (DG > 0) {\r
@@ -5902,7 +5906,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                            ABER0 = qmod.getAbserr();\r
                            NVAL = qmod.getNeval();\r
                            IER[0] = qmod.getErrorStatus();\r
-                           QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                           QIERR[IER[0]]=QIERR[IER[0]]+1;\r
                        TYPE=4;\r
                        qmod = new qaws(TC,1.0,routine,0.0,0.0,2,AQTOL,RQTOL,limit);\r
                        qmod.driver();\r
@@ -5910,7 +5914,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                            ABER1 = qmod.getAbserr();\r
                            NVAL = qmod.getNeval();\r
                            IER[0] = qmod.getErrorStatus();\r
-                           QIERR[IER[0]-1]=QIERR[IER[0]-1]+1;\r
+                           QIERR[IER[0]]=QIERR[IER[0]]+1;\r
                        WORK[1]=U0+U1;\r
                        MQERR[0]=Math.max(MQERR[0],ABER0+ABER1);\r
        \r
@@ -7073,7 +7077,8 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
        //     LOCAL VARIABLES\r
        \r
-       int I,I0,I1,I2,J,K,JI,JI0,JI1,JI2,P0,P1,P2,HI,LO,QINTS,NQ,DIFF,TNCQP;\r
+       int QINTS[] = new int[1];\r
+       int I,I0,I1,I2,J,K,JI,JI0,JI1,JI2,P0,P1,P2,HI,LO,NQ,DIFF,TNCQP;\r
        double DST[] = new double[2];\r
        double BETA,H1,M1,T0,T2,SUM1,RR,RRB,MEAN,RXI,IXI;\r
        final double ONE[] = new double[]{1.0,0.0};\r
@@ -7210,98 +7215,689 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
        TNCQP=LOQSB[NJIND-1]+NQUAD[NJIND-1]-1;\r
        HI=0;\r
-       /*for (JI=1; JI <= NJIND; JI++) {\r
+       for (JI=1; JI <= NJIND; JI++) {\r
            NQ=NQUAD[JI-1];\r
            if (NEWQU[JI-1]) {\r
                LO=(JI-1)*NQPTS+1;\r
                BETA=JACIN[JI-1];\r
                I2=2*JI-1;\r
-               SUBIN7(XIVAL(I2),2,BETA,MDGPO,NQPTS,ACOEF(LO),BCOEF(LO),\r
-                     H0VAL(JI),COLSC(LO),AQTOL,TOLOU(JI),XENPT,QINTS,MQIN1,\r
+               double ZZin[][] = new double[2][2];\r
+               for (K = 0; K < 2; K++) {\r
+                       ZZin[K][0] = XIVAL[I2-1+K][0];\r
+                       ZZin[K][1] = XIVAL[I2-1+K][1];\r
+               }\r
+               double AJAC[] = new double[MDGPO];\r
+               double BJAC[] = new double[MDGPO];\r
+               for (K = 0; K < MDGPO; K++) {\r
+                       AJAC[K] = ACOEF[LO-1+K];\r
+                       BJAC[K] = BCOEF[LO-1+K];\r
+               }\r
+               double CSCAL[]= new double[MDGPO+1];\r
+               for (K = 0; K < MDGPO+1; K++) {\r
+                       CSCAL[K] = COLSC[LO-1+K];\r
+               }\r
+               double TOLIO[] = new double[1];\r
+               TOLIO[0] = TOLOU[JI-1];\r
+               SUBIN7(ZZin,2,BETA,MDGPO,NQPTS,AJAC,BJAC,\r
+                     H0VAL[JI-1],CSCAL,AQTOL,TOLIO,XENPT,QINTS,MQIN1,\r
                      IER);\r
-                 IF (IER .GT. 0) THEN\r
-                   IF (IER .EQ. 11) THEN\r
-                     IER=20\r
-                   ENDIF\r
-                   RETURN\r
-                 ENDIF\r
+               TOLOU[JI-1] = TOLIO[0];\r
+               if (IER[0] > 0) {\r
+                   if (IER[0] == 11) {\r
+                       IER[0]=20;\r
+                   }\r
+                   return;\r
+               } // if (IER[0] > 0)\r
+       \r
+               DIFF=QINTS[0]*NQPTS-NQ;\r
+               if (TNCQP+DIFF > MNQUA) {\r
+                   IER[0]=19;\r
+                   return;\r
+               }\r
+               I1=HI+NQ+1;\r
+        \r
+               // IF DIFF IS NON-ZERO WE MUST MAKE SPACE IN ARRAYS QCOMX AND\r
+                // QCOMW TO RECEIVE THE NEW DATA.\r
+       \r
+               if (DIFF > 0) {\r
+                   for (I=TNCQP; I >= I1; I--) {\r
+                       J=I+DIFF;\r
+                       QCOMX[J-1]=QCOMX[I-1];\r
+                       QCOMW[J-1]=QCOMW[I-1];\r
+                   } // for (I=TNCQP; I >= I1; I--)\r
+               } // if (DIFF > 0)\r
+               else if (DIFF < 0) {\r
+                   for (I=I1; I <= TNCQP; I++) {\r
+                       J=I+DIFF;\r
+                       QCOMX[J-1]=QCOMX[I-1];\r
+                       QCOMW[J-1]=QCOMW[I-1];\r
+                   } // for (I=I1; I <= TNCQP; I++)\r
+               } // else if (DIFF < 0)\r
+       \r
+               // NOW SET UP THE NEW RULE AND STORE DATA IN QCOMX, QCOMW\r
+       \r
+               TNCQP=TNCQP+DIFF;\r
+               NQUAD[JI-1]=NQ+DIFF;\r
+               LOQSB[JI-1]=HI+1;\r
+               SUM1=BETA+1.0;\r
+               K=HI;\r
+               for (I=1; I <= QINTS[0]; I++) {\r
+                   RR=(XENPT[I]-XENPT[I-1])*0.5;\r
+                   MEAN=(XENPT[I]+XENPT[I-1])*0.5; \r
+                   if (I == 1) {\r
+                       RRB=Math.pow(RR,SUM1);\r
+                       LO=LO-1;\r
+                       for (J=1; J <= NQPTS; J++) {\r
+                           K=K+1;\r
+                           QCOMX[K-1]=MEAN+RR*QUPTS[LO+J-1];\r
+                           QCOMW[K-1]=RRB*QUWTS[LO+J-1];\r
+                       } // for (J=1; J <= NQPTS; J++)\r
+                   } // if (I == 1)\r
+                   else {\r
+                       LO=NARCS*NQPTS;\r
+                       for (J=1; J <= NQPTS; J++) {\r
+                           K=K+1;\r
+                           QCOMX[K-1]=MEAN+RR*QUPTS[LO+J-1];\r
+                           QCOMW[K-1]=RR*QUWTS[LO+J-1]*Math.pow((1.0+QCOMX[K-1]),BETA);\r
+                       } // for (J=1; J <= NQPTS; J++)\r
+                   } // else\r
+               } // for (I=1; I <= QINTS[0]; I++)\r
+               HI=HI+NQUAD[JI-1];\r
+           } // if (NEWQU[JI-1])\r
+           else {\r
+       \r
+               // HERE WE DO NOTHING OTHER THAN UPDATE SOME SUBSCRIPTS.\r
+       \r
+               LOQSB[JI-1]=HI+1;\r
+               HI=HI+NQ;\r
+           }\r
+       \r
+       } // for (JI=1; JI <= NJIND; JI++)\r
+       \r
+       MCQER[0]=0.0;\r
+       for (I=1; I <= NJIND; I++) {\r
+           MCQER[0]=Math.max(MCQER[0],TOLOU[I-1]);\r
+       }\r
+       \r
+       NUQTL[0]=false;\r
+       \r
+       //     NORMAL TERMINATION\r
+       \r
+       IER[0]=0;\r
+       \r
+    } // private void UPOCOQ1\r
+\r
+\r
+    private void TSJAC3(int LOTES[], int HITES[], double TESPR[], double ZTEST[][], int NQPTS,\r
+        int NTEST, int ORDSG, int TNSUA, int TSTNG, int DGPOL[], int JATYP[], int PARNT[],\r
+        double AICOF[], double BICOF[], double DIAG[], double HALEN[], double JACIN[],\r
+        double MIDPT[], double SDIAG[], int IER[]) {\r
+       \r
+        // INTEGER IER,NQPTS,NTEST,ORDSG,TNSUA,TSTNG\r
+       // INTEGER DGPOL(*),HITES(*),JATYP(*),LOTES(*),PARNT(*)\r
+       // REAL AICOF(*),BICOF(*),DIAG(*),HALEN(*),JACIN(*),MIDPT(*),\r
+       // +SDIAG(*),TESPR(*)\r
+        // COMPLEX ZTEST(*)\r
+       \r
+       // TO ASSIGN THE TEST PARAMETERS (STORED IN TESPR),THE TEST POINTS ON\r
+       // THE PHYSICAL BOUNDARY (STORED IN ZTEST) AND THE ARRAYS LOTES AND \r
+       // HITES NEEDED TO ACCESS THIS DATA CORRECTLY. \r
+       \r
+       // IER=0 - NORMAL EXIT\r
+       // IER=21 - FAILURE IN IMTQLH\r
+       \r
+       // LOCAL VARIABLES\r
+       \r
+       int IFAIL[] = new int[1];\r
+       int D,FIRST,I,J,K,K1,P,PREV,TSFBS;\r
+       double S,TT;\r
+       double PIN[] = new double[2];\r
+       double POUT[];\r
+       // COMPLEX PARFUN\r
+       // EXTERNAL ASONJ7,IMTQLH,PARFUN\r
+       \r
+       TSFBS=TNSUA/ORDSG;\r
+        if (TSTNG != 1) {\r
+           for (I=1; I <= TSFBS; I++) {\r
+               TESPR[I-1]=-1.0;\r
+               P=PARNT[I-1];\r
+               TT=MIDPT[I-1]-HALEN[I-1];\r
+               PIN[0] = TT;\r
+               PIN[1] = 0.0;\r
+               POUT = PARFUN(P,PIN);\r
+               ZTEST[I-1][0] = POUT[0];\r
+               ZTEST[I-1][1] = POUT[1];\r
+               LOTES[I-1]=I;\r
+               HITES[I-1]=I;\r
+           } // for (I=1; I <= TSFBS; I++) \r
+        } // if (TSTNG != 1)\r
+        else { // TSTNG == 1           \r
+           LOTES[0]=1;\r
+           HITES[0]=1+DGPOL[0];\r
+           for (I=2; I <= TSFBS; I++) {\r
+               LOTES[I-1]=HITES[I-2]+1;\r
+               HITES[I-1]=LOTES[I-1]+DGPOL[I-1];\r
+           } // for (I=2; I <= TSFBS; I++)\r
+       \r
+           for (I=1; I <= TSFBS; I++) {\r
+               D=DGPOL[I-1];\r
+               P=PARNT[I-1];\r
+               FIRST=LOTES[I-1];\r
+               TESPR[FIRST-1]=-1.0;\r
+               TT=MIDPT[I-1]-HALEN[I-1];\r
+               PIN[0] = TT;\r
+               PIN[1] = 0.0;\r
+               POUT = PARFUN(P,PIN);\r
+               ZTEST[FIRST-1][0] = POUT[0];\r
+               ZTEST[FIRST-1][1] = POUT[1];\r
+               if (D > 0) {\r
+                   J=JATYP[I-1];\r
+                   if (J > 0) {\r
+                       S=1.0;\r
+                   }\r
+                   else {\r
+                       S=-1.0;\r
+                       J=-J;\r
+                   }\r
+                   PREV=(J-1)*NQPTS;\r
+                   for (K=1; K <= D; K++) {\r
+                       K1=PREV+K;\r
+                       DIAG[K-1]=BICOF[K1-1];\r
+                       if (K == 1) {\r
+                           SDIAG[K-1]=0.0;\r
+                       }\r
+                       else {\r
+                           SDIAG[K-1]=AICOF[K1-2];\r
+                       }\r
+                   } // for (K=1; K <= D; K++)\r
+                   IMTQLH(D,DIAG,SDIAG,IFAIL);\r
+                   if (IFAIL[0] > 0) {\r
+                       IER[0]=21;\r
+                       return;\r
+                   }\r
+                   for (K=1; K <= D; K++) {\r
+                       TT=S*DIAG[K-1];\r
+                       K1=FIRST+K;\r
+                       TESPR[K1-1]=TT;\r
+                       TT=MIDPT[I-1]+HALEN[I-1]*TT;\r
+                       PIN[0] = TT;\r
+                       PIN[1] = 0.0;\r
+                       POUT = PARFUN(P,PIN);\r
+                       ZTEST[K1-1][0] = POUT[0];\r
+                       ZTEST[K1-1][1] = POUT[1];\r
+                   } // for (K=1; K <= D; K++)\r
+               } // if (D > 0)\r
+           } // for (I=1; I <= TSFBS; I++) \r
+        } // else TSTNG == 1\r
+       NTEST=HITES[TSFBS-1];\r
+       if (ORDSG > 1) {\r
+           NTEST=NTEST+1;\r
+           LOTES[TSFBS]=NTEST;\r
+           HITES[TSFBS]=NTEST;\r
+           TESPR[NTEST-1]=-1.0;\r
+           P=PARNT[TSFBS];\r
+           PIN[0] = -1.0;\r
+           PIN[1] = 0.0;\r
+           POUT = PARFUN(P,PIN);\r
+           ZTEST[NTEST-1][0] = POUT[0];\r
+           ZTEST[NTEST-1][1] = POUT[1];\r
+       }\r
+       \r
+       IER[0]=0;\r
+       \r
+    } // private void TSJAC3 \r
+    \r
+    private void TESMD9(double ERMOD[], double MATRX[][], double SOLUN[], int MNEQN, int NCOLL,\r
+        int NTEST, int NQPTS, int TNSUA, int JATYP[], int PARNT[], int DGPOL[], int LOSUB[],\r
+        int HISUB[], int LOTES[], int HITES[], int NQUAD[], int LOQSB[], double TOLNR,\r
+        double MIDPT[], double HALEN[], double H0VAL[], double COLSC[], double ACOEF[], \r
+        double BCOEF[], double TESPR[], double QCOMX[], double QCOMW[], double CENTR[], \r
+        double ZTEST[][], boolean INTER, boolean LNSEG[], double WORK[], int QIERR[],\r
+        double MQERR[], double JACIN[], double A1COF[], double B1COF[], double AQTOL,\r
+        double RQTOL, double AQCOF[], double BQCOF[], double CQCOF[], double MXERM,\r
+        int IMXER, double ZMXER[], double ERARC[], int ORDSG, boolean REFLN) {\r
+       \r
+       //INTEGER MNEQN,NCOLL,NTEST,NQPTS,TNSUA,JATYP(*),PARNT(*),DGPOL(*),\r
+       //+LOSUB(*),HISUB(*),LOTES(*),HITES(*),NQUAD(*),LOQSB(*),QIERR(0:6),\r
+       //+IMXER,ORDSG\r
+       \r
+       //REAL MATRX(MNEQN,*),SOLUN(*),TOLNR,MIDPT(*),HALEN(*),H0VAL(*),\r
+       //+COLSC(*),ACOEF(*),BCOEF(*),TESPR(*),QCOMX(*),QCOMW(*),MQERR,\r
+       //+WORK(*),A1COF(*),B1COF(*),AQTOL,JACIN(*),ERARC(*),\r
+       //+RQTOL,AQCOF(*),BQCOF(*),CQCOF(*),ERMOD(*),MXERM\r
+       \r
+       //COMPLEX CENTR,ZMXER,ZTEST(*)\r
+       \r
+       //LOGICAL INTER,LNSEG(*),REFLN\r
+       \r
+       // TO COMPUTE THE ERROR IN MODULUS AT THE VECTOR OF TEST POINTS\r
+       // ZTEST (PARAMETER VALUES IN TESPR) STORING RESULTS IN ERMOD.\r
+       \r
+        // LOCAL VARIABLES\r
+       \r
+       int IER[] = new int[1];\r
+       int I,J,I1,J1,IA,IQ,JI,AJI,PT,DG,LOCLM,HICLM,LOQ,LOD,NQ,IT,\r
+           NVAL,FIRST,LAST,NEQNS,ROW,TSFBS,ORDRG,IB;\r
+       double MD,HL,RH,TQ,WQ,DD,LDD,SG,SJI,TT,C0,U0,U1,CURR,PREV,NEXT,ABER0,ABER1,RLIM,LLIM,SUM;\r
+       // REAL FNVAL\r
+       double GTQ[] = new double[2];\r
+       double ZQ[] = new double[2];\r
+       double ZT[] = new double[2];\r
+       // COMPLEX GTQ,ZQ,ZT,PARFUN,DPARFN\r
+       // COMMON /FNDEF/BETA,A1,B1,P0VAL,SCALE,TYPE\r
+       // EXTERNAL DPARFN,FNVAL,JAPAR7,PARFUN,QAWS,R1MACH\r
+       qaws qmod;\r
+       int routine = Integration2.DQAWSE;\r
+       /**\r
+         * Gives an upper bound on the number of subintervals in the partition of lower, upper.\r
+         */\r
+        int limit = 100;\r
+       \r
+       TSFBS=TNSUA/ORDSG;\r
+       NEQNS=NCOLL+1;\r
+       RLIM=1.0-5.0*EPS;\r
+       LLIM=-RLIM;\r
+       for (J=1; J <= NCOLL; J++) {\r
+           for (I=1; I <= NTEST; I++) {\r
+               MATRX[I-1][J-1]=0.0;\r
+           }\r
+       } // for (J=1; J <= NCOLL; J++)\r
+       for (I=0; I <= 6; I++) {\r
+           QIERR[I]=0;\r
+       }\r
+       MQERR[0]=0.0;\r
+       \r
+       // NOW SELECT THE INTEGRATION ARC\r
+       \r
+       /*for (IA=1; IA <= TNSUA; IA++) {\r
+\r
+       \r
+           // INITIALISE DATA FOR THIS ARC\r
+       \r
+           PT=PARNT[IA-1];\r
+           MD=MIDPT[IA-1];\r
+           HL=HALEN[IA-1];\r
+           DG=DGPOL[IA-1];\r
+           LOCLM=LOSUB[IA-1];\r
+           HICLM=HISUB[IA-1];\r
+           if (ORDSG == 1 || (ORDSG > 1 && IA <= (TSFBS+1))) {\r
+               FIRST=LOTES[IA-1];\r
+               LAST=HITES[IA-1]+1;\r
+           }\r
+           else if (IA == TNSUA && ORDSG > 1) {\r
+               FIRST=NTEST+1;\r
+               LAST=FIRST;\r
+           }\r
+           else {\r
+               FIRST=NTEST+1;\r
+               LAST=NTEST;\r
+           }\r
+           JI=JATYP[IA-1];\r
+           if (JI < 0) {\r
+               SJI=-1.0;\r
+           }\r
+           else {\r
+               SJI=1.0;\r
+           }\r
+           AJI=Math.abs(JI);\r
+           LOD=(AJI-1)*NQPTS+1;\r
+           RH=Math.sqrt(H0VAL[AJI-1]);\r
+           BETA=JACIN[AJI-1];\r
+           A1=A1COF[AJI-1];\r
+           B1=B1COF[AJI-1];\r
+           P0VAL=1.0/RH;\r
+       \r
+           // SET UP THE DIAGONAL BLOCK OF SCALED PRINCIPAL SINGULAR INTEGRALS\r
+       \r
+           for (IT=FIRST; IT <= LAST; IT++) {\r
+               if (ORDSG > 1 && IA == (TSFBS+1) && IT == LAST) continue;\r
+               if (IT == LAST) {\r
+                   TT=SJI;\r
+               }\r
+               else {\r
+                   TT=SJI*TESPR[IT-1];\r
+               }\r
+       \r
+               if (IT > NTEST) {\r
+                   ROW=1;\r
+               }\r
+               else {\r
+                   ROW=IT;\r
+               }\r
+       \r
+               SCALE=COLSC[LOD-1];\r
+               if (TT < LLIM) {\r
+                   U0=0.0;\r
+                   ABER0=0.0;\r
+                   TYPE=1;\r
+                   // Use dqawse which is the same as (quadpack/dqaws) but provides more\r
+                   // information and control\r
+                   qmod = new qaws(-1.0,1.0,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                   qmod.driver();\r
+                   U1 = qmod.getIntegral();\r
+                   ABER1 = qmod.getAbserr();\r
+                   NVAL = qmod.getNeval();\r
+                   IER[0] = qmod.getErrorStatus();\r
+                   QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+               }\r
+               else if (TT >= RLIM) {\r
+                   TYPE=1;\r
+                   qmod = new qaws(-1.0,1.0,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                   qmod.driver();\r
+                   U0 = qmod.getIntegral();\r
+                   ABER0 = qmod.getAbserr();\r
+                   NVAL = qmod.getNeval();\r
+                   IER[0] = qmod.getErrorStatus();\r
+                   QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+                   U1=0.0;\r
+                   ABER1=0.0;\r
+               }\r
+               else {\r
+                   TYPE=1;\r
+                   qmod = new qaws(-1.0,TT,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                   qmod.driver();\r
+                   U0 = qmod.getIntegral();\r
+                   ABER0 = qmod.getAbserr();\r
+                   NVAL = qmod.getNeval();\r
+                   IER[0] = qmod.getErrorStatus();\r
+                   QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+                   TYPE=2;\r
+                   qmod = new qaws(TT,1.0,routine,0.0,0.0,2,AQTOL,RQTOL,limit);\r
+                   qmod.driver();\r
+                   U1 = qmod.getIntegral();\r
+                   ABER1 = qmod.getAbserr();\r
+                   NVAL = qmod.getNeval();\r
+                   IER[0] = qmod.getErrorStatus();\r
+                   QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+               }\r
+               WORK[0]=U0+U1;\r
+               MQERR[0]=Math.max(MQERR[0],ABER0+ABER1);\r
+       \r
+               if (DG > 0) {\r
+                   SCALE=COLSC[LOD];\r
+                   if (TT < LLIM) {\r
+                       U0=0.0;\r
+                       ABER0=0.0;\r
+                       TYPE=3;\r
+                       qmod = new qaws(-1.0,1.0,routine,BETA,0.0,2,AQTOL,RQTOL,limit);\r
+                           qmod.driver();\r
+                           U1 = qmod.getIntegral();\r
+                           ABER1 = qmod.getAbserr();\r
+                           NVAL = qmod.getNeval();\r
+                           IER[0] = qmod.getErrorStatus();\r
+                           QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+                   }\r
+                   else if (TT > RLIM) {\r
+                       TYPE=3;\r
+                       qmod = new qaws(-1.0,1.0,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                           qmod.driver();\r
+                           U0 = qmod.getIntegral();\r
+                           ABER0 = qmod.getAbserr();\r
+                           NVAL = qmod.getNeval();\r
+                           IER[0] = qmod.getErrorStatus();\r
+                           QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+                       U1=0.0;\r
+                       ABER1=0.0;\r
+                   }\r
+                   else {\r
+                       TYPE=3;\r
+                       qmod = new qaws(-1.0,TT,routine,BETA,0.0,3,AQTOL,RQTOL,limit);\r
+                           qmod.driver();\r
+                           U0 = qmod.getIntegral();\r
+                           ABER0 = qmod.getAbserr();\r
+                           NVAL = qmod.getNeval();\r
+                           IER[0] = qmod.getErrorStatus();\r
+                           QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+                       TYPE=4;\r
+                       qmod = new qaws(TT,1.0,routine,0.0,0.0,2,AQTOL,RQTOL,limit);\r
+                           qmod.driver();\r
+                           U1 = qmod.getIntegral();\r
+                           ABER1 = qmod.getAbserr();\r
+                           NVAL = qmod.getNeval();\r
+                           IER[0] = qmod.getErrorStatus();\r
+                           QIERR[IER[0]]=QIERR[IER[0]]+1;\r
+                   }\r
+                     WORK(2)=U0+U1\r
+                     MQERR=MAX(MQERR,ABER0+ABER1)\r
+       C\r
+       C             NOW USE THE (WEAKLY) STABLE FORWARD RECURRENCE SCHEME FOR \r
+       C             WORK(I),I=3,NQPTS (WITH SCALE FACTOR FOR WORK(2))\r
+       C\r
+                     CURR=WORK(2)\r
+                     PREV=SCALE\r
+                     DO 10 I=1,DG-1\r
+                         J=LOD+I-1\r
+                         NEXT=(AQCOF(J)*TT-BQCOF(J))*CURR-CQCOF(J)*PREV\r
+                         WORK(I+2)=NEXT\r
+                         PREV=CURR\r
+                         CURR=NEXT\r
+       10            CONTINUE\r
+       C\r
+       C             ASSIGN CORRECT SCALE FACTORS.\r
        C\r
-                 DIFF=QINTS*NQPTS-NQ\r
-                 IF (TNCQP+DIFF .GT. MNQUA) THEN\r
-                   IER=19\r
-                   RETURN\r
+                     DO 20 I=3,DG+1\r
+                         J=LOD+I-1\r
+                         WORK(I)=WORK(I)*COLSC(J)/SCALE\r
+       20            CONTINUE\r
+               } // if (DG > 0)\r
+       C\r
+                 SG=1E+0\r
+                 DO 30 J=LOCLM,HICLM\r
+                   MATRX(ROW,J)=MATRX(ROW,J)+SG*WORK(J-LOCLM+1)\r
+                   SG=SG*SJI\r
+       30        CONTINUE\r
+       C\r
+           } // for (IT=FIRST; IT <= LAST; IT++)\r
+       C\r
+       C       INITIALISE SOME DATA FOR THE NON-SINGULAR INTEGRALS\r
+       C\r
+               WORK(1)=1E+0/RH\r
+               NQ=NQUAD(AJI)\r
+               LOQ=LOQSB(AJI)\r
+               IF (IA.EQ.TNSUA) THEN\r
+                   I1=2\r
+               ELSE\r
+                   I1=1\r
+               ENDIF\r
+       C\r
+               DO 100 IQ=1,NQ\r
+                 I=LOQ+IQ-1\r
+                 TQ=QCOMX(I)\r
+                 WQ=QCOMW(I)\r
+                 CALL JAPAR7(WORK,TQ,ACOEF(LOD),BCOEF(LOD),DG)\r
+                 IF (JI .LT. 0) THEN\r
+                   TQ=-TQ\r
                  ENDIF\r
-                 I1=HI+NQ+1\r
-       C \r
-       C         IF DIFF IS NON-ZERO WE MUST MAKE SPACE IN ARRAYS QCOMX AND\r
-       C         QCOMW TO RECEIVE THE NEW DATA.\r
+                 GTQ=CMPLX(MD+HL*TQ)\r
+                 ZQ=PARFUN(PT,GTQ)\r
+       C\r
+       C         ACCUMULATE THE ELEMENTS ABOVE THE DIAGONAL BLOCK\r
        C\r
-                 IF (DIFF .GT. 0) THEN\r
-                   DO 40 I=TNCQP,I1,-1\r
-                     J=I+DIFF\r
-                     QCOMX(J)=QCOMX(I)\r
-                     QCOMW(J)=QCOMW(I)\r
+       C\r
+                 DO 50 IT=I1,FIRST-1\r
+                   ZT=ZTEST(IT)\r
+                   DD=ABS(ZT-ZQ)\r
+                   LDD=LOG(DD)*WQ\r
+                   SG=1E+0\r
+                   DO 40 J1=LOCLM,HICLM\r
+                     J=J1-LOCLM+1\r
+                     I=J1-LOCLM+LOD\r
+                     MATRX(IT,J1)=MATRX(IT,J1)+SG*WORK(J)*LDD*COLSC(I)\r
+                     SG=SG*SJI\r
        40          CONTINUE\r
-                 ELSE IF (DIFF .LT. 0) THEN\r
-                   DO 50 I=I1,TNCQP\r
-                     J=I+DIFF\r
-                     QCOMX(J)=QCOMX(I)\r
-                     QCOMW(J)=QCOMW(I)\r
-       50          CONTINUE\r
+       50        CONTINUE\r
+       C\r
+       C         ACCUMULATE THE ELEMENTS BELOW THE DIAGONAL BLOCK\r
+       C\r
+                 DO 70 IT=LAST+1,NTEST\r
+                   ZT=ZTEST(IT)\r
+                   DD=ABS(ZT-ZQ)\r
+                   LDD=LOG(DD)*WQ\r
+                   SG=1E+0\r
+                   DO 60 J1=LOCLM,HICLM\r
+                     J=J1-LOCLM+1\r
+                     I=J1-LOCLM+LOD\r
+                     MATRX(IT,J1)=MATRX(IT,J1)+SG*WORK(J)*LDD*COLSC(I)\r
+                     SG=SG*SJI\r
+       60          CONTINUE\r
+       70        CONTINUE\r
+       C\r
+       C         ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
+       C         THE DIAGONAL BLOCK FOR THE NON-LINE-SEGMENT CASE.\r
+       C\r
+                 IF (.NOT.LNSEG(IA)) THEN\r
+                   DO 90 IT=FIRST,LAST\r
+                     IF (ORDSG.GT.1 .AND. IA.EQ.(TSFBS+1) .AND. IT.EQ.LAST) \r
+            +        GOTO 90\r
+                     IF (IT .EQ. LAST) THEN\r
+                         TT=1E+0\r
+                     ELSE\r
+                         TT=TESPR(IT)\r
+                     ENDIF\r
+       C\r
+                     IF (IT .GT. NTEST) THEN\r
+                         ROW=1\r
+                     ELSE\r
+                         ROW=IT\r
+                     ENDIF\r
+       C\r
+                     DD=ABS(TT-TQ)\r
+                     IF (DD .LE. TOLNR) THEN\r
+                         DD=ABS(DPARFN(PT,GTQ))*HL\r
+                     ELSE\r
+                         IF (IT .GT. NTEST) THEN\r
+                             ZT=ZTEST(1)\r
+                         ELSE\r
+                             ZT=ZTEST(IT)\r
+                         ENDIF\r
+                         DD=ABS(ZT-ZQ)/DD\r
+                         IF (DD.LT.TOLNR) DD=ABS(DPARFN(PT,GTQ))*HL\r
+                     ENDIF\r
+                     LDD=LOG(DD)*WQ\r
+                     SG=1E+0\r
+                     DO 80 J1=LOCLM,HICLM\r
+                       J=J1-LOCLM+1\r
+                       I=J1-LOCLM+LOD\r
+                       MATRX(ROW,J1)=MATRX(ROW,J1)+SG*WORK(J)*LDD*COLSC(I)\r
+                       SG=SG*SJI\r
+       80            CONTINUE\r
+       90          CONTINUE\r
                  ENDIF\r
+       100     CONTINUE\r
        C\r
-       C         NOW SET UP THE NEW RULE AND STORE DATA IN QCOMX, QCOMW\r
+       C       ACCUMULATE THE RESIDUAL NON-SINGULAR CONTRIBUTIONS INTO\r
+       C       THE DIAGONAL BLOCK FOR THE LINE-SEGMENT CASE.\r
        C\r
-                 TNCQP=TNCQP+DIFF\r
-                 NQUAD(JI)=NQ+DIFF\r
-                 LOQSB(JI)=HI+1\r
-                 SUM1=BETA+1E+0\r
-                 K=HI\r
-                 DO 80 I=1,QINTS\r
-                   RR=(XENPT(I+1)-XENPT(I))*5E-1\r
-                   MEAN=(XENPT(I+1)+XENPT(I))*5E-1 \r
-                   IF (I .EQ. 1) THEN\r
-                     RRB=RR**SUM1\r
-                     LO=LO-1\r
-                     DO 60 J=1,NQPTS\r
-                       K=K+1\r
-                       QCOMX(K)=MEAN+RR*QUPTS(LO+J)\r
-                       QCOMW(K)=RRB*QUWTS(LO+J)\r
-       60            CONTINUE\r
+               IF (LNSEG(IA)) THEN\r
+                 ZT=DPARFN(PT,CMPLX(MD))*HL\r
+                 C0=ABS(ZT)\r
+                 C0=RH*LOG(C0)*COLSC(LOD)\r
+                 DO 110 IT=FIRST,LAST\r
+                   IF (ORDSG.GT.1 .AND. IA.EQ.(TSFBS+1) .AND. IT.EQ.LAST) \r
+            +      GOTO 110\r
+                   IF (IT .GT. NTEST) THEN\r
+                       ROW=1\r
                    ELSE\r
-                     LO=NARCS*NQPTS\r
-                     DO 70 J=1,NQPTS\r
-                       K=K+1\r
-                       QCOMX(K)=MEAN+RR*QUPTS(LO+J)\r
-                       QCOMW(K)=RR*QUWTS(LO+J)*(1E+0+QCOMX(K))**BETA\r
-       70            CONTINUE\r
+                       ROW=IT\r
                    ENDIF\r
-       80        CONTINUE\r
-                 HI=HI+NQUAD(JI)\r
-           } // if (NEWQU[JI-1])\r
-           else {\r
+                   MATRX(ROW,LOCLM)=MATRX(ROW,LOCLM)+C0\r
+       110       CONTINUE\r
+               ENDIF\r
+       } // for (IA=1; IA <= TNSUA; IA++)\r
        C\r
-       C         HERE WE DO NOTHING OTHER THAN UPDATE SOME SUBSCRIPTS.\r
+       C     SET UP THE LAST COLUMN\r
        C\r
-                 LOQSB(JI)=HI+1\r
-                 HI=HI+NQ\r
-           }\r
+             DO 130 I=1,NTEST\r
+               MATRX(I,NEQNS)=1E+0\r
+       130   CONTINUE\r
        C\r
-       } // for (JI=1; JI <= NJIND; JI++)\r
+       C     COMPUTE MATRIX-VECTOR PRODUCT\r
        C\r
-             MCQER=0E+0\r
-             DO 100 I=1,NJIND\r
-               MCQER=MAX(MCQER,TOLOU(I))\r
-       100   CONTINUE\r
+             DO 150 I=1,NTEST\r
+               SUM=0E+0\r
+               DO 140 J=1,NEQNS\r
+                 SUM=SUM+MATRX(I,J)*SOLUN(J)\r
+       140     CONTINUE\r
+               ERMOD(I)=SUM\r
+       150   CONTINUE\r
        C\r
-             NUQTL=.FALSE.\r
+       C     FORM THE ERROR IN MODULUS\r
        C\r
-       C     NORMAL TERMINATION\r
+             IF (INTER) THEN\r
+                 DO 160 I=1,NTEST\r
+                   SUM=EXP(ERMOD(I))\r
+                   ERMOD(I)=ABS(1E+0-ABS(ZTEST(I)-CENTR)/SUM)\r
+       160       CONTINUE\r
+             ELSE\r
+                 DO 170 I=1,NTEST\r
+                   SUM=EXP(ERMOD(I))\r
+                   ERMOD(I)=ABS(1E+0-SUM)\r
+       170       CONTINUE\r
+             ENDIF\r
        C\r
-             IER=0\r
-       C */\r
-    } // private void UPOCOQ1\r
-\r
+       C     FIND MAXIMUM ERROR IN MODULUS AND THE POINT AND THE ARC AT WHICH \r
+       C     IT OCCURS\r
+       C\r
+             MXERM=0E+0\r
+             DO 190 IA=1,TSFBS\r
+               FIRST=LOTES(IA)\r
+               LAST=HITES(IA)\r
+               MD=0E+0\r
+               DO 180 IT=FIRST,LAST\r
+                 IF (ERMOD(IT) .GT. MD) THEN\r
+                     MD=ERMOD(IT)\r
+                 ENDIF\r
+                 IF (MD .GT. MXERM) THEN\r
+                     MXERM=MD\r
+                     IMXER=IA\r
+                     ZMXER=ZTEST(IT)\r
+                 ENDIF\r
+       180     CONTINUE\r
+               IF (IA.EQ.TSFBS .AND. ORDSG.EQ.1) THEN\r
+                 IT=1\r
+               ELSE\r
+                 IT=LAST+1\r
+               ENDIF\r
+               IF (ERMOD(IT) .GT. MD) THEN\r
+                   MD=ERMOD(IT)\r
+               ENDIF\r
+               IF (MD .GT. MXERM) THEN\r
+                   MXERM=MD\r
+                   IMXER=IA\r
+                   ZMXER=ZTEST(IT)\r
+               ENDIF\r
+               ERARC(IA)=MD\r
+       190   CONTINUE\r
+       C\r
+       C     IF REGION IS SYMMETRIC, FILL UP THE WHOLE ERARC VECTOR USING\r
+       C     SYMMETRY\r
+       C\r
+             IF (ORDSG.GT.1) THEN\r
+             IF (REFLN) THEN\r
+               DO 200 IA=1,TSFBS\r
+                 IB=2*TSFBS+1-IA\r
+                 ERARC(IB)=ERARC(IA)\r
+       200     CONTINUE\r
+               ORDRG=ORDSG/2\r
+               DO 220 I=2,ORDRG\r
+                 I1=(I-1)*TSFBS*2\r
+                 DO 210 IA=1,TSFBS*2\r
+                   IB=I1+IA\r
+                   ERARC(IB)=ERARC(IA)\r
+       210       CONTINUE\r
+       220     CONTINUE\r
+             ELSE\r
+               DO 240 I=2,ORDSG\r
+                 I1=(I-1)*TSFBS\r
+                 DO 230 IA=1,TSFBS\r
+                   IB=I1+IA\r
+                   ERARC(IB)=ERARC(IA)\r
+       230       CONTINUE\r
+       240     CONTINUE\r
+             ENDIF\r
+             ENDIF\r
+       C*/\r
+    } // private void TESMD9\r
 \r
 \r
       /**\r