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

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

index e4865b7..8c5db88 100644 (file)
@@ -4,6 +4,7 @@ import java.io.File;
 import java.io.IOException;\r
 import java.io.RandomAccessFile;\r
 import java.util.Scanner;\r
+import java.util.Vector;\r
 \r
 import gov.nih.mipav.model.structures.ModelImage;\r
 import gov.nih.mipav.model.structures.jama.GeneralizedEigenvalue;\r
@@ -291,8 +292,8 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        private double MD;\r
        private double HL;\r
        //From LEVCUR:\r
-       String NEWD;\r
-       private double ZZ[][] = new double[2][2];\r
+       Vector<Double> Contour[];\r
+       Vector<Double>Ray[];\r
 \r
        public SymmsIntegralMapping() {\r
 \r
@@ -14444,6 +14445,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        IER[0]=0;\r
     } // private void BMCAP1\r
        \r
+       @SuppressWarnings("unchecked")\r
        private void LEVCUR(int NCONT, double RADII[], int NARGS, double THETA[], double RAD1[],\r
         double RAD2[], double PSD[], double MINPD[], double MAXPD[], String NEWD[], int IER[]) {\r
                \r
@@ -14657,14 +14659,17 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                //     LOCAL VARAIBLES\r
                \r
                final int MXTRY = 5;\r
-               int COARG,CLNO,IC,IR,L,MNSUA,NRAYS,PPNO,PT,TRY,VTARG;\r
-               double DIAPHY,DIFF,DPD,HS,HT,MXRAD,PI,R1MACH,RMAX,RMEAN,RMIN,SC1,\r
+               //int L, MNSUA;\r
+               int CLNO,IC,IR,NRAYS,PPNO,PT,TRY;\r
+               // double PI;\r
+               double DIFF,DPD,HS,HT,MXRAD,RMAX,RMEAN,RMIN,SC1,\r
                     SC2,TH1,TH2,THET0,TMIN,TINC;\r
                double WW[][] = new double[2][2];\r
                double WEND[] = new double[2];\r
+               double ZZ[][] = new double[2][2];\r
                //COMPLEX WW(2),WEND,ZZ(2)\r
                boolean ATEND,NECES,WANTC,WANTR;\r
-               String OFL, JBNM;\r
+               //String OFL, JBNM;\r
                //CHARACTER OFL*6,JBNM*4\r
                double PHYPT[][] = new double[1][2];\r
                double CANPT[][] = new double[1][2];\r
@@ -14678,12 +14683,12 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                \r
                // INITIALISE SOME VARIABLES\r
                \r
-               /*IER[0]=0;\r
+               IER[0]=0;\r
                WANTC=(NCONT >= 1);\r
                WANTR=(RAD2[0] > RAD1[0]);\r
-               PI=Math.PI;\r
+               //PI=Math.PI;\r
                NARCS=ISNCA[0];\r
-               MNSUA=IGEOM[3];\r
+               //MNSUA=IGEOM[3];\r
                NJIND=NARCS+1;\r
                TNGQP=ISNCA[1]*NJIND;\r
                THET0=VTARG[0];\r
@@ -14725,19 +14730,21 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                    MINPD[0]=2.0;\r
                    MAXPD[0]=5.0;\r
                } // if (PSD <= 0.0)\r
-               RMIN=MINPD*DPD/PSD[0];\r
-               RMAX=MAXPD*DPD/PSD[0];\r
+               RMIN=MINPD[0]*DPD/PSD[0];\r
+               RMAX=MAXPD[0]*DPD/PSD[0];\r
                RMEAN=0.5*(RMIN+RMAX);\r
                TMIN=Math.sqrt(EPS);\r
                \r
                // THE DO 50 LOOP DETERMINES THE IMAGE OF A CONTOUR\r
                \r
                if (WANTC) {\r
+                       Contour = new Vector[NCONT];\r
                    for (IC=1; IC <= NCONT; IC++) {\r
+                       Contour[IC-1] = new Vector<Double>();\r
                        PPNO=0;\r
                        CLNO=0;\r
                        SC1=RADII[IC-1];\r
-                       HT=SC1*(MINPD+MAXPD)/PSD[0];\r
+                       HT=SC1*(MINPD[0]+MAXPD[0])/PSD[0];\r
                        //WRITE(CHNL,*) NEWD\r
                        ATEND=false;\r
                        PT=1;\r
@@ -14745,7 +14752,7 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                        WW[0][0] = SC1 * Math.cos(TH1);\r
                        WW[0][1] = SC1 * Math.sin(TH1);\r
                        CLNO=CLNO+1;\r
-                       DMCANP(1,ZZ,WW,INTER,false,IER);\r
+                       DMCANP(1,ZZ,WW,false,IER);\r
                        if (IER[0] > 0) {\r
                                WRTAIL(8,0,IER[0],null);\r
                                return;\r
@@ -14753,6 +14760,8 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                        PPNO=PPNO+1;\r
                        //WRITE(CHNL,20) ZZ(1)\r
                //20      FORMAT(2E16.8)\r
+                       Contour[IC-1].add(ZZ[0][0]);\r
+                       Contour[IC-1].add(ZZ[0][1]);\r
                        TINC=HT;\r
                        TRY=0;\r
                \r
@@ -14794,17 +14803,20 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                             } //if (TRY < MXTRY)\r
                             PPNO=PPNO+1;\r
                             //WRITE(CHNL,20) ZZ(2)\r
+                            Contour[IC-1].add(ZZ[1][0]);\r
+                            Contour[IC-1].add(ZZ[1][1]);\r
                \r
                            ATEND=(PT > NARCS);\r
                            if (!ATEND) {\r
                                ZZ[0][0]=ZZ[1][0];\r
                                ZZ[0][1]=ZZ[1][1];\r
-                               WW[0[0]]=WW[1][0];\r
+                               WW[0][0]=WW[1][0];\r
                                WW[0][1]=WW[1][1];\r
                                TH1=TH2;\r
                                TRY=0;\r
                                continue;\r
                            } // if (!ATEND)\r
+                           break;\r
                        } // while (true)\r
                \r
                        System.out.println("CONTOUR " + IC + " DONE: " + PPNO + " PTS " + CLNO + " TRIES");\r
@@ -14814,72 +14826,95 @@ public class SymmsIntegralMapping extends AlgorithmBase {
                \r
                // THE DO 70 LOOP DETERMINES THE IMAGE OF A RAY\r
                \r
-                     if (WANTR) {\r
-                         HS=(RAD2[0]-RAD1[0])*(MINPD[0]+MAXPD[0])/PSD;\r
-                         if (NARGS >= 1) {\r
-                       NRAYS=NARGS\r
-                     ELSE\r
-                       NRAYS=NARCS\r
-                     ENDIF\r
-                     DO 70 IR=1,NRAYS\r
-                       IF (NARGS.GE.1) THEN\r
-                         WEND=CMPLX(COS(THETA(IR)),SIN(THETA(IR)))\r
-                       ELSE\r
-                         WEND=CMPLX(COS(RSNCA(COARG-1+IR)),SIN(RSNCA(COARG-1+IR)))\r
-                       ENDIF\r
-                       CLNO=0\r
-                       PPNO=0\r
-                       WRITE(CHNL,*) NEWD\r
-                       SC1=RAD1\r
-                       WW(1)=WEND*SC1\r
-                       CLNO=CLNO+1\r
-                       CALL DMCANP(1,ZZ,WW,INTER,CENTR,IGEOM,RGEOM,ISNCA,RSNCA,ZSNCA,\r
-                    +              IQUCA,ZQUCA,.FALSE.,IER)\r
-                       IF (IER.GT.0) GOTO 999\r
-                       PPNO=PPNO+1\r
-                       WRITE(CHNL,20) ZZ(1)\r
-                       TINC=HS\r
-                       TRY=0\r
-               C\r
-               60      CONTINUE\r
-                       SC2=SC1+TINC\r
-                       IF (SC2 .GT. RAD2) THEN\r
-                         SC2=RAD2\r
-                         ATEND=.TRUE.\r
-                       ELSE\r
-                         ATEND=.FALSE.\r
-                       ENDIF\r
-                       WW(2)=WEND*SC2\r
-                       CLNO=CLNO+1\r
-                       TRY=TRY+1\r
-                       CALL DMCANP(1,ZZ(2),WW(2),INTER,CENTR,IGEOM,RGEOM,ISNCA,RSNCA,\r
-                    +              ZSNCA,IQUCA,ZQUCA,.FALSE.,IER)\r
-                       IF (IER.GT.0) GOTO 999\r
-                       IF (TRY.LT.MXTRY) THEN\r
-                         DIFF=ABS(ZZ(2)-ZZ(1))\r
-                         IF (DIFF.GT.RMAX .OR. (DIFF.LT.RMIN .AND. (.NOT.ATEND))) THEN\r
-                           TINC=RMEAN*TINC/DIFF\r
-                           GOTO 60\r
-                         ENDIF\r
-                       ENDIF\r
-                       PPNO=PPNO+1\r
-                       WRITE(CHNL,20) ZZ(2)\r
-               C\r
-                       IF (.NOT. ATEND) THEN\r
-                         ZZ(1)=ZZ(2)\r
-                         WW(1)=WW(2)\r
-                         SC1=SC2\r
-                         TRY=0\r
-                         GOTO 60\r
-                       ENDIF\r
-               C\r
-                       WRITE(*,2) 'RAY ',IR,' DONE:',PPNO,' PTS  ',CLNO,' TRIES'\r
-               70    CONTINUE\r
-                     } // if (WANTR)\r
-               C\r
-               999   CALL WRTAIL(8,0,IER)\r
-                     CLOSE(CHNL)\r
-               C*/\r
+               if (WANTR) {\r
+                   HS=(RAD2[0]-RAD1[0])*(MINPD[0]+MAXPD[0])/PSD[0];\r
+                   if (NARGS >= 1) {\r
+                       NRAYS=NARGS;\r
+                   }\r
+                   else {\r
+                       NRAYS=NARCS;\r
+                   }\r
+                   Ray = new Vector[NRAYS];\r
+                   for (IR=1; IR <= NRAYS; IR++) {\r
+                       Ray[IR-1] = new Vector<Double>();\r
+                       if (NARGS >= 1) {\r
+                           WEND[0] = Math.cos(THETA[IR-1]);\r
+                           WEND[1] = Math.sin(THETA[IR-1]);\r
+                       }\r
+                       else {\r
+                               WEND[0] = Math.cos(COARG[IR-1]);\r
+                               WEND[1] = Math.sin(COARG[IR-1]);\r
+                       }\r
+                       CLNO=0;\r
+                       PPNO=0;\r
+                       //WRITE(CHNL,*) NEWD\r
+                       SC1=RAD1[0];\r
+                       WW[0][0]=WEND[0]*SC1;\r
+                       WW[0][1]=WEND[1]*SC1;\r
+                       CLNO=CLNO+1;\r
+                       DMCANP(1,ZZ,WW,false,IER);\r
+                       if (IER[0] > 0) {\r
+                               WRTAIL(8,0,IER[0],null);\r
+                               return; \r
+                       }\r
+                       PPNO=PPNO+1;\r
+                       //WRITE(CHNL,20) ZZ(1)\r
+                       Ray[IR-1].add(ZZ[0][0]);\r
+                       Ray[IR-1].add(ZZ[0][1]);\r
+                       TINC=HS;\r
+                       TRY=0;\r
+               \r
+                   while (true) {\r
+                           SC2=SC1+TINC;\r
+                           if (SC2 > RAD2[0]) {\r
+                               SC2=RAD2[0];\r
+                               ATEND=true;\r
+                           }\r
+                           else {\r
+                               ATEND=false;\r
+                           }\r
+                          WW[1][0]=WEND[0]*SC2;\r
+                          WW[1][1]=WEND[1]*SC2;\r
+                          CLNO=CLNO+1;\r
+                          TRY=TRY+1;\r
+                          CANPT[0][0] = WW[1][0];\r
+                          CANPT[0][1] = WW[1][1];\r
+                          DMCANP(1,PHYPT,CANPT,false,IER);\r
+                          ZZ[1][0] = PHYPT[0][0];\r
+                          ZZ[1][1] = PHYPT[0][1];\r
+                          if (IER[0] > 0) {\r
+                                  WRTAIL(8,0,IER[0],null);\r
+                                  return;         \r
+                          }\r
+                          if (TRY < MXTRY) {\r
+                              DIFF=zabs(ZZ[1][0]-ZZ[0][0],ZZ[1][1]-ZZ[0][1]);\r
+                              if (DIFF > RMAX || (DIFF < RMIN && (!ATEND))) {\r
+                                  TINC=RMEAN*TINC/DIFF;\r
+                                  continue;\r
+                              }\r
+                          } // if (TRY < MXTRY)\r
+                          PPNO=PPNO+1;\r
+                          //WRITE(CHNL,20) ZZ(2)\r
+                          Ray[IR-1].add(ZZ[1][0]);\r
+                          Ray[IR-1].add(ZZ[1][1]);\r
+                         if (!ATEND) {\r
+                             ZZ[0][0]=ZZ[1][0];\r
+                             ZZ[0][1]=ZZ[1][1];\r
+                             WW[0][0]=WW[1][0];\r
+                             WW[0][1]=WW[1][1];\r
+                             SC1=SC2;\r
+                             TRY=0;\r
+                             continue;\r
+                         } // if (!ATEND)\r
+                         break;\r
+                   } // while (true)\r
+               \r
+                       System.out.println("RAY "+ IR + " DONE: "+ PPNO +" PTS, " + CLNO + " TRIES");\r
+                       Preferences.debug("RAY "+ IR + " DONE: "+ PPNO +" PTS, " + CLNO + " TRIES\n", Preferences.DEBUG_ALGORITHM);\r
+                   } // for (IR=1; IR <= NRAYS; IR++)\r
+               } // if (WANTR)\r
+               \r
+               WRTAIL(8,0,IER[0],null);\r
     } // private void LEVCUR\r
        \r
        private double DIAPHY(int NARCS) {\r
@@ -14939,6 +14974,401 @@ public class SymmsIntegralMapping extends AlgorithmBase {
        } // private double DIAPHY\r
 \r
 \r
+       private void TSTPLT(double MXMIS[], double MXDIF[],double PSD[], double MINPD[],\r
+        double MAXPD[],int IER[]) {\r
+        //INTEGER NARCS,CHNL,IER\r
+        //REAL MXMIS,MXDIF,PSD,MINPD,MAXPD\r
+        //CHARACTER*4 JBNM\r
+\r
+        // ......................................................................\r
+\r
+        // 1.     TSTPLT\r
+        //           TESTS THE PARAMETRIC FUNCTION ROUTINES PARFUN AND DPARFN\r
+        //           FOR CONSISTENCY AND OUTPUTS BOUNDARY POINTS FOR PLOTTING.\r
+           \r
+\r
+        // 2.     PURPOSE\r
+        //           THE ROUTINE FIRST CHECKS THAT THE PARAMETRIC FUNCTION\r
+        //           ROUTINE PARFUN IS CONSISTENT WITH RESPECT TO ITS DEFINITION\r
+        //           OF ANY CORNERS ON THE BOUNDARY.  THIS IS DONE BY CHECKING \r
+        //           THAT THE COMPUTED POINT AT THE END OF EACH ARC MATCHES THE \r
+        //           COMPUTED POINT AT THE START OF THE NEXT ONE.  IF ALL THE\r
+        //           RELATIVE MISFIT ERRORS AT CORNERS ARE LESS THAN\r
+        //           10*(UNIT ROUNDOFF) THEN ALL CORNERS ARE CONSIDERED TO FIT\r
+        //           SATISFACTORILY, OTHERWISE THE MAXIMUM RELATIVE MISFIT\r
+        //           ERROR IS REPORTED.\r
+\r
+        //           THE SECOND PURPOSE OF THE ROUTINE IS TO OUTPUT TO A DATA \r
+        //           FILE THE COORDINATES OF A NUMBER OF POINTS ON THE BOUNDARY. \r
+        //           THE BOUNDARY POINTS ARE SELECTED ADAPTIVELY TO MEET THE \r
+        //           PLOTTING RESOLUTION SPECIFICATIONS DEFINED BY THE ARGUMENTS \r
+        //           PSD,MINPD,MAXPD (SEE BELOW).  THE HOPE IS THAT THE USER MAY \r
+        //           EASILY FEED THESE DATA POINTS TO HIS LOCAL GRAPH PLOTTING  \r
+        //           ROUTINES SO AS TO CONSTRUCT A PLOT OF THE BOUNDARY.  THIS   \r
+        //           WILL PROVIDE AN ESSENTIAL VISUAL CHECK ON THE VALIDITY OF \r
+        //           THE ROUTINE PARFUN.  THE OUTPUT DATA FILE IS AUTOMATICALLY \r
+        //           NAMED <JBNM>zz.\r
\r
+        //           THE THIRD PURPOSE OF THE ROUTINE IS TO CHECK PARFUN AND \r
+        //           DPARFN FOR MUTUAL CONSISTECY.  THIS IS DONE BY COMPUTING \r
+        //           TWO POINT FINITE DIFFERENCE APPROXIMATIONS TO DPARFN.  \r
+        //           THESE DIFFERENCE APPROXIMATIONS ARE COMPUTED AT EACH BOUND- \r
+        //           ARY POINT THAT IS OUTPUT FOR PLOTTING AND ALSO AT NEARBY \r
+        //           POINTS WHICH LIE JUST O F F THE BOUNDARY.  THIS LATTER \r
+        //           COMPARISON ALSO TESTS PARFUN AND DPARFN FOR CORRECTNESS IN \r
+        //           ACCEPTING COMPLEX PARAMETER VALUES.   A RELATIVE ERROR IN \r
+        //           THE FINITE DIFFERENCE APPROXIMATION GREATER THAN 0.1 IS \r
+        //           REPORTED AS A POSSIBLE LOGICAL INCONSISTENCY BETWEEN PARFUN \r
+        //           AND DPARFN.  (THE CRITICAL RELATIVE ERROR VALUE OF 0.1 CAN \r
+        //           BE ALTERED BY ADJUSTING THE LOCAL PARAMETER DTOL).\r
+\r
+        // 3.     CALLING SEQUENCE\r
+        //           CALL TSTPLT(JBNM,MXMIS,MXDIF,NARCS,PSD,MINPD,MAXPD,CHNL,IER)\r
+\r
+        //        PARAMETERS\r
+        //         ON ENTRY\r
+        //            JBNM   - CHARACTER*4\r
+        //                     THE JOB NAME.  THIS IS USED TO CREATE THE OUTPUT \r
+        //                     FILE WITH FILENAME\r
+\r
+        //                                <JBNM>zz ,\r
+\r
+        //                     WHERE <JBNM> DENOTES THE VALUE OF VARIABLE JBNM\r
+        //                     WITH ANY TRAILING SPACES DELETED.\r
+\r
+        //            NARCS  - INTEGER\r
+        //                     THE NUMBER OF ANALYTIC ARCS THAT MAKE UP THE \r
+        //                     W H O L E BOUNDARY OF THE PHYSICAL DOMAIN.\r
+\r
+        //            PSD    - REAL\r
+        //                     THE PLOTTING SIZE FOR THE DOMAIN IN ANY APPROPR-\r
+        //                     IATE UNITS.  IF PSD .LE. 0.0 THEN IT IS ASSIGNED\r
+        //                     THE DEFAULT VALUE OF 160.0 (A REASONBLE WIDTH IN\r
+        //                     MM FOR PLOTTING ON A4 PAPER).\r
+\r
+        //            MINPD  - REAL\r
+        //                     THE MINIMUM SIGNIFICANT PLOTTING DISTANCE, IN THE\r
+        //                     SAME UNITS AS PSD.  IF PSD .LE. 0.0 THEN IT IS\r
+        //                     ASSIGNED THE DEFAULT VALUE OF 2.0.\r
+\r
+        //            MAXPD  - REAL\r
+        //                     THE MAXIMUM ALLOWED PLOTTING DISTANCE, IN THE\r
+        //                     SAME UNITS AS PSD.  IF PSD .LE. 0.0 THEN IT IS\r
+        //                     ASSIGNED THE DEFAULT VALUE OF 5.0.  THE LARGER\r
+        //                     MAXPD, THE COARSER WILL BE THE RESOLUTION OF THE\r
+        //                     BOUNDARY POINTS OUTPUT TO <JBNM>zz, BUT THE\r
+        //                     QUICKER THEY WILL BE COMPUTED. \r
+\r
+        //            CHNL   - INTEGER\r
+        //                     DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR\r
+        //                     WRITING THE FILE <JBNM>zz.\r
+        //         ON EXIT\r
+        //            MXMIS  - REAL\r
+        //                     THE MAXIMUM RELATIVE MISFIT ERROR OVER ALL\r
+        //                     CORNER POINTS\r
+\r
+        //            MXDIF  - REAL\r
+        //                     THE MAXIMUM RELATIVE ERROR IN FINITE DIFFERENCE\r
+        //                     APPROXIMATIONS TO DPARFN OVER ALL BOUNDARY \r
+        //                     POINTS OUTPUT TO <JBNM>zz AND NEARBY POINTS OFF\r
+        //                     THE BOUNDARY.\r
+\r
+        //            PSD    - REAL\r
+        //                     IF PSD .LE. 0.0 ON ENTRY THEN IT WILL HAVE\r
+        //                     THE DEFAULT VALUE OF 160.0 ON  EXIT.\r
+\r
+        //            MINPD  - REAL\r
+        //                     IF PSD .LE. 0.0 ON ENTRY THEN MINPD WILL HAVE\r
+        //                     THE DEFAULT VALUE OF 2.0 ON EXIT\r
+\r
+        //            MAXPD  - REAL\r
+        //                     IF PSD .LE. 0.0 ON ENTRY THEN MAXPD WILL HAVE\r
+        //                     THE DEFAULT VALUE OF 5.0 ON EXIT\r
+\r
+        //            IER    - INTEGER\r
+        //                     IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED;\r
+        //                     A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY\r
+        //                     WRITTEN ON THE STANDARD OUTPUT CHANNEL.\r
+        //                     IER=0 - NORMAL EXIT.\r
+        //                     IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD\r
+        //                             BE SELF EXPLANATORY.\r
+\r
+\r
+        // 4.     SUBROUTINES OR FUNCTIONS NEEDED\r
+        //              - THE CONFPACK LIBRARY.\r
+        //              - THE REAL FUNCTION R1MACH.\r
+        //              - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN.\r
+\r
+\r
+        // 5.     FURTHER COMMENTS\r
+        //              - A SUMMARY LISTING IS AUTOMATICALLY WRITTEN ON THE \r
+        //                STANDARD OUTPUT CHANNEL.\r
+        //              - THE OUTPUT FILE <JBNM>zz CONTAINS COORDINATE PAIRS\r
+\r
+        //                                 X Y\r
+\r
+        //                FOR POINTS ON THE PHYSICAL BOUNDARY, WITH ONE PAIR\r
+        //                PER LINE.\r
+\r
+        // ......................................................................\r
+        //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+        //     LAST UPDATE: 6 JULY 1990\r
+        // ......................................................................C\r
+        //     LOCAL VARIABLES\r
+\r
+        final int MNARC = 200;\r
+        final int NH = 4;\r
+        int IMX = 0;\r
+               int I,IA,L;\r
+               final double DTOL = 0.1;\r
+        double A1,DIFF,ERR,HH,MINC,R1MACH,RMAX,RMEAN,RMIN,T,TINC,\r
+            TOL1,TMX,TSD;\r
+        double TT[] = new double[2];\r
+        //REAL TT(2)\r
+        double C1[] = new double[2];\r
+        double C2[] = new double[2];\r
+        double CENTR[] = new double[2];\r
+        double ZZ0[] = new double[2];\r
+        double DZZ[] = new double[2];\r
+        double NDZZ[] = new double[2];\r
+        // COMPLEX C1,C2,CENTR,ZZ0,PARFUN,DZZ,DPARFN,ZDPARF,NDZZ\r
+        double ZZ[][] = new double[2][2];\r
+        //COMPLEX ZZ(2)\r
+        //CHARACTER OFL*6\r
+        boolean ATEND,FIRST,WARND;\r
+        boolean LNSEG[]  = new boolean[MNARC];\r
+        double PIN[] = new double[2];\r
+        double POUT[];\r
+        //EXTERNAL DPARFN,LINSEG,PARFUN,R1MACH,WRHEAD,WRTAIL,ZDPARF\r
+\r
+        // WRITE CONFPACK HEADING\r
+\r
+        WRHEAD(7,0,null);\r
+\r
+        if (NARCS > MNARC) {\r
+            IER[0]=59;\r
+            WRTAIL(7,0,IER[0],null);\r
+            return;\r
+        }\r
+\r
+        //1     FORMAT(A45)\r
+        //2     FORMAT(A45,I4)\r
+        //3     FORMAT(A45,E10.3)\r
+        //4     FORMAT(//,T17,A)\r
+\r
+        TOL1=10.0*EPS;\r
+\r
+        // CHECK THAT ALL ARCS MEET AT CORNER POINTS\r
+\r
+        IER[0]=0;\r
+        CENTR[0] = 0.0;\r
+        CENTR[1] = 0.0;\r
+        MXMIS[0]=0.0;\r
+        for (IA=1; IA <= NARCS; IA++) {\r
+            if (IA == 1) {\r
+                I=NARCS;\r
+            }\r
+            else {\r
+                I=IA-1;\r
+            }\r
+            PIN[0] = -1.0;\r
+            PIN[1] = 0.0;\r
+            C1=PARFUN(IA,PIN);\r
+            CENTR[0]=CENTR[0]+C1[0];\r
+            CENTR[1]=CENTR[1]+C1[1];\r
+            A1=zabs(C1[0],C1[1]);\r
+            PIN[0] = 1.0;\r
+            PIN[1] = 0.0;\r
+            C2=PARFUN(I,PIN);\r
+            ERR=zabs(C1[0]-C2[0],C1[1]-C2[1]);\r
+            if (A1 >= 1.0) {\r
+                ERR=ERR/A1;\r
+            }\r
+            if (ERR > MXMIS[0]) {\r
+                IMX=IA;\r
+                MXMIS[0]=ERR;\r
+            }\r
+        } // for (IA=1; IA <= NARCS; IA++)\r
+        if (MXMIS[0] >= TOL1) {\r
+            System.out.println("MAXIMUM CORNER MISFIT: " + MXMIS[0]);\r
+            Preferences.debug("MAXIMUM CORNER MISFIT: " + MXMIS[0] + "\n", Preferences.DEBUG_ALGORITHM);\r
+            System.out.println("OCCURS AT CORNER: " + IMX);\r
+            Preferences.debug("OCCURS AT CORNER: " + IMX + "\n", Preferences.DEBUG_ALGORITHM);\r
+        }\r
+        else {\r
+            System.out.println("ALL ARCS FIT AT CORNERS:");\r
+            Preferences.debug("ALL ARCS FIT AT CORNERS:\n",Preferences.DEBUG_ALGORITHM);\r
+        }\r
+\r
+        // ESTIMATE THE DIAMETER (TSD) OF THE PHYSICAL DOMAIN\r
+\r
+        CENTR[0]=CENTR[0]/NARCS;\r
+        CENTR[1]=CENTR[1]/NARCS;\r
+        TSD=0.0;\r
+        HH=2.0/(double)(NH);\r
+        for (IA=1; IA <= NARCS; IA++) {\r
+            T=-1.0;\r
+            for (I=1; I <= NH; I++) {\r
+                T=T+HH;\r
+                PIN[0] = T;\r
+                PIN[1] = 0.0;\r
+                POUT = PARFUN(IA,PIN);\r
+                C1[0]=POUT[0]-CENTR[0];\r
+                C1[1]=POUT[1]-CENTR[1];\r
+                A1=zabs(C1[0],C1[1]);\r
+                TSD=Math.max(TSD,A1);\r
+            } // for (I=1; I <= NH; I++)\r
+        } // for (IA=1; IA <= NARCS; IA++) \r
+        TSD=2.0*TSD;\r
+\r
+        // DETERMINE WHICH ARCS (IF ANY) ARE LINE SEGMENTS\r
+\r
+        LINSEG(LNSEG,NARCS);\r
+\r
+        // OPEN FILE TO RECEIVE BOUNDARY DATA POINTS FOR PLOTTING\r
+\r
+        //L=INDEX(JBNM,' ')-1\r
+        //IF (L.EQ.-1) L=4\r
+        //OFL=JBNM(1:L)//'zz'\r
+        //OPEN(CHNL,FILE=OFL)\r
+\r
+        //SET DEFAULT PLOTTING DISTANCES, IF NECESSARY\r
+\r
+        if (PSD[0] <= 0.0) {\r
+            PSD[0]=160.0;\r
+            MINPD[0]=2.0;\r
+            MAXPD[0]=5.0;\r
+        }\r
+        RMIN=MINPD[0]*TSD/PSD[0];\r
+        RMAX=MAXPD[0]*TSD/PSD[0];\r
+        RMEAN=0.5*(RMIN+RMAX);\r
+        MINC=Math.sqrt(EPS);\r
+\r
+        // START EVALUATING BOUNDARY POINTS AND DERIVATIVES FOR PLOTTING AND\r
+        // TESTING\r
+\r
+        MXDIF[0]=0.0;\r
+        /*for (IA=1; IA <= NARCS; IA++) {\r
+            TT[0]=-1.0;\r
+            PIN[0] = TT[0];\r
+            PIN[1] = 0.0;\r
+            ZZ[0]=PARFUN(IA,PIN);\r
+            //WRITE(CHNL,'(2E16.7)') ZZ(1)\r
+            if (IA==1) {\r
+               ZZ0[0]=ZZ[0][0];\r
+               ZZ0[1]=ZZ[0][1];\r
+            }\r
+            FIRST=true;\r
+            WARND=false;\r
+            while (true) {\r
+\r
+                // TEST THE COMPATIBILTY OF PARFUN AND DPARFN BY ESTIMATING DPARFN\r
+                // NUMERICALLY AT BOTH REAL AND COMPLEX PARAMETER VALUES.\r
+\r
+                for (I=1; I <= 2; I++) {\r
+                    if (I == 1) {\r
+                       C1[0] = TT[0];\r
+                       C1[1] = TT[1];\r
+                    }\r
+                    else {\r
+                       C1[0] = TT[0];\r
+                       C1[1] = MINC;\r
+                    }\r
+                    DZZ=DPARFN(IA,C1);\r
+                    NDZZ=ZDPARF(IA,C1);\r
+                    A1=zabs(DZZ[0],DZZ[1]);\r
+\r
+                    if (A1 == 0.0) {\r
+                        IER[0]=60;\r
+                        System.out.println();\r
+                        Preferences.debug("\n", Preferences.DEBUG_ALGORITHM);\r
+                        System.out.println("              ***DPARFN=(0.,0.)***");\r
+                        Preferences.debug("              ***DPARFN=(0.,0.)***\n", Preferences.DEBUG_ALGORITHM);\r
+                        System.out.println("                             ARC: " + IA); \r
+                        Preferences.debug("                             ARC: " + IA + "\n", Preferences.DEBUG_ALGORITHM);\r
+                        System.out.println("STANDARDISED PARAMETER VALUE: " + TT[0][0] + " , " + TT[0][1]);\r
+                        Preferences.debug("STANDARDISED PARAMETER VALUE: " + TT[0][0] + " , " + TT[0][1] + "\n", Preferences.DEBUG_ALGORITHM);\r
+                        WRTAIL(7,0,IER[0],null);\r
+                        return;\r
+                    } // if (A1 == 0.0)\r
+\r
+                    if (A1 <= TOL1 && !WARND) {\r
+                        System.out.println("*** W A R N I N G  ***");\r
+                        Preferences.debug("*** W A R N I N G  ***\n", Preferences.DEBUG_ALGORITHM);\r
+                        System.out.println("PATHOLOGICALLY SMALL DERIVATIVE ON ARC " + IA);\r
+                        Preferences.debug("PATHOLOGICALLY SMALL DERIVATIVE ON ARC " + IA + "\n", Preferences.DEBUG_ALGORITHM);\r
+                        WARND=true;\r
+                    } // if (A1 <= TOL1 && !WARND)\r
+C\r
+        IF (FIRST) THEN\r
+          TINC=RMEAN/A1\r
+          TINC=MAX(TINC,MINC)\r
+          FIRST=.FALSE.\r
+        ENDIF\r
+C\r
+        ERR=ABS(1E+0-NDZZ/DZZ)\r
+        IF (ERR.GT.MXDIF) THEN\r
+          MXDIF=ERR\r
+          IMX=IA\r
+          TMX=TT(1)\r
+        ENDIF\r
+                } // for (I=1; I <= 2; I++)\r
+C\r
+      IF (.NOT.LNSEG(IA)) THEN\r
+C\r
+C****     DETERMINE THE NEXT BOUNDARY POINT TO BE PLOTTED\r
+C\r
+40        CONTINUE\r
+        TT(2)=TT(1)+TINC\r
+        IF (TT(2) .GE. 1E+0) THEN\r
+          TT(2)=1E+0\r
+          ATEND=.TRUE.\r
+        ELSE\r
+          ATEND=.FALSE.\r
+        ENDIF\r
+C\r
+        ZZ(2)=PARFUN(IA,CMPLX(TT(2)))\r
+        DIFF=ABS(ZZ(2)-ZZ(1))\r
+        IF (DIFF.EQ.0E+0 .AND. .NOT.ATEND) THEN\r
+          TINC=MAX(MINC,2*TINC)\r
+          GOTO 40\r
+        ENDIF\r
+C\r
+        IF (DIFF.GT.RMAX .OR. (DIFF.LT.RMIN .AND. .NOT.ATEND)) THEN\r
+          TINC=RMEAN*TINC/DIFF\r
+          TINC=MAX(TINC,MINC)\r
+          GOTO 40\r
+        ENDIF\r
+C\r
+        WRITE(CHNL,'(2E16.7)') ZZ(2)\r
+        IF (.NOT. ATEND) THEN\r
+          ZZ(1)=ZZ(2)\r
+          TT(1)=TT(2)\r
+          continue;\r
+        ENDIF\r
+      ENDIF\r
+      break;\r
+            } // while(true)\r
+C\r
+        } // for (IA=1; IA <= NARCS; IA++)\r
+    IF (LNSEG(NARCS)) WRITE(CHNL,'(2E16.7)') ZZ0\r
+    CLOSE(CHNL)   \r
+C\r
+    IF (MXDIF .GT. DTOL) THEN\r
+      WRITE(*,*)\r
+      WRITE(*,2) 'POSSIBLE PARFUN/DPARFN INCONSISTECY ON ARC:',IMX\r
+      WRITE(*,3) 'OCCURS AT STANDARDISED PARAMETER VALUE:',TMX\r
+      WRITE(*,3) 'RELATIVE FINITE DIFF ERROR:',MXDIF \r
+    ELSE\r
+      WRITE(*,*)\r
+      WRITE(*,1) 'PARFUN AND DPARFN ARE CONSISTENT:'\r
+    ENDIF\r
+C\r
+999   CALL WRTAIL(7,0,IER)\r
+C */\r
+    } // private void TSTPLT\r
+       \r
+       \r
 \r
        /**\r
         * zabs computes the absolute value or magnitude of a double precision\r