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

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

index 13560ca..baa7af0 100644 (file)
@@ -13604,13 +13604,15 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        final int MNCOF = 32;\r
        final int MQIN1 = 21;\r
        final int MXNQD = 144;\r
-       final int NBSCT = 3;\r
+       //final int NBSCT = 3;\r
+       int QINTS[] = new int[1];\r
        int AJTC,DGC,I,IA,IP,K,J,J1,JQ,JTC,LODC,LOL,\r
-            LOMC,NQ,NQUAD,PSA,QINTS;\r
+            LOMC,NQ,NQUAD,PSA;\r
        final double PTHTL = 1.0E-3;\r
        final double DELTA = 2.0E-1;\r
+       double TOLOU[] = new double[1];\r
        double AA,ARG,ARGW,AWW,BB,BETAC,DIST,EFPTL,EPS_CATPH4,HA,ILW,IMXI,MD,\r
-            MEAN,PI,REXI,RLW,RR,RRB,S2C,SCO,SJTC,SUM1,TOLOU,TT,\r
+            MEAN,PI,REXI,RLW,RR,RRB,S2C,SCO,SJTC,SUM1,TT,\r
             TUPI,UPPER;\r
        double CT0[] = new double[2];\r
        double CT2[] = new double[2];\r
@@ -13635,12 +13637,27 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
         double ci[] = new double[1];\r
         double cr2[] = new double[1];\r
         double ci2[] = new double[1];\r
-       \r
+        double var;\r
+        double ang;\r
+        double A[];\r
+        double B[];\r
+        double CO[][];\r
+        int N;\r
+        double CIN[] = new double[2];\r
+        double diff[] = new double[2];\r
+       double POUT[];\r
+       double ebase;\r
+       double expr;\r
+       double expi;\r
+       double JOUT[];\r
+       double logr;\r
+       double logi;\r
+        \r
        EPS_CATPH4=10.0 * EPS;\r
        PI=Math.PI;\r
        TUPI=2.0*PI;\r
        LOL=NARCS*NQPTS;\r
-       /*loopIP: for (IP=1; IP <= NPTS; IP++) {\r
+       loopIP: for (IP=1; IP <= NPTS; IP++) {\r
            WW[0]=CANPT[IP-1][0];\r
            WW[1]=CANPT[IP-1][1];\r
            AWW=zabs(WW[0],WW[1]);\r
@@ -13777,114 +13794,182 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                    TERM[1] = SJTC*XI[1];\r
                    if ((TERM[0] == 0.0) && (TERM[1] == 0.0)) {\r
                        PHYPT[IP-1]=PARFUN(PARNT[PSA-1],CT0);\r
-                       continue IPloop;\r
+                       continue loopIP;\r
                    }\r
     \r
-                   IF (TERM.EQ.(2E+0,0E+0)) THEN\r
-                     PHYPT(IP)=PARFUN(PARNT(PSA),CT2)\r
-                     GOTO 300\r
-                   ENDIF\r
-       C\r
-                   ARG=ATAN2(AIMAG(TERM),REAL(TERM))\r
-                   UPPER=(1E+0+S2C*5E-1)*PI\r
-                   IF (ARG.GT.UPPER) THEN\r
-                     ARG=ARG-TUPI\r
-                   ELSE IF (ARG.LE.(UPPER-TUPI)) THEN\r
-                     ARG=ARG+TUPI\r
-                   ENDIF\r
-       C\r
-                   FP=ABS(TERM)**(BETAC+1E+0)*CMPLX(COS((BETAC+1E+0)*ARG),\r
-            +                                       SIN((BETAC+1E+0)*ARG))\r
-                   IF (INTER) THEN\r
-                     FP=-FP\r
-                   ENDIF\r
-                   PHI=CCJACS(SJTC*XI,DGC-1,A1COC(LODC),B1COC(LODC),\r
-            +                 H1VLC(AJTC),BFSNC(LOMC+1))\r
-                   PHI=(BFSNC(LOMC)-(1E+0-SJTC*XI)*PHI)*SJTC*CMPLX(0E+0,1E+0)\r
-                   PHYPT(IP)=CENTR+(PARFUN(PARNT(PSA),CT0)-CENTR)*EXP(FP*PHI)\r
-                   GOTO 300\r
+                   if ((TERM[0] == 2.0) && (TERM[1] == 0.0)) {\r
+                       PHYPT[IP-1]=PARFUN(PARNT[PSA-1],CT2);\r
+                       continue loopIP;\r
+                   }\r
+       \r
+                   ARG=Math.atan2(TERM[1],TERM[0]);\r
+                   UPPER=(1.0+S2C*0.5)*PI;\r
+                   if (ARG > UPPER) {\r
+                       ARG=ARG-TUPI;\r
+                   }\r
+                   else if (ARG <= (UPPER-TUPI)) {\r
+                       ARG=ARG+TUPI;\r
+                   }\r
+       \r
+                   var = Math.pow(zabs(TERM[0],TERM[1]),(BETAC + 1.0));\r
+                   ang = (BETAC + 1.0)*ARG;\r
+                   FP[0] = var * Math.cos(ang);\r
+                   FP[1] = var * Math.sin(ang);\r
+                   if (INTER) {\r
+                       FP[0]=-FP[0];\r
+                       FP[1]=-FP[1];\r
+                   }\r
+                   CIN[0] = SJTC*XI[0];\r
+                   CIN[1] = SJTC*XI[1];\r
+                   A = new double[DGC-1];\r
+                   B = new double[DGC-1];\r
+                   for (N = 0; N < DGC-1; N++) {\r
+                       A[N] = A1COC[LODC-1+N];\r
+                       B[N] = B1COC[LODC-1+N];\r
+                   }\r
+                   CO = new double[DGC][2];\r
+                   for (N = 0; N < DGC; N++) {\r
+                       CO[N][0] = BFSNC[LOMC+N][0];\r
+                       CO[N][1] = BFSNC[LOMC+N][1];\r
+                   }\r
+                   PHI=CCJACS(CIN,DGC-1,A,B,H1VLC[AJTC-1],CO);\r
+                   zmlt(1.0 - SJTC*XI[0],-SJTC*XI[1], PHI[0],PHI[1],cr,ci);\r
+                   diff[0] = BFSNC[LOMC-1][0] - cr[0];\r
+                   diff[1] = BFSNC[LOMC-1][1] - ci[0];\r
+                   zmlt(diff[0],diff[1],0.0,SJTC,cr,ci);\r
+                   PHI[0] = cr[0];\r
+                   PHI[1] = ci[0];\r
+                   POUT = PARFUN(PARNT[PSA-1],CT0);\r
+                   diff[0] = POUT[0]- CENTR[0];\r
+                   diff[1] = POUT[1] - CENTR[1];\r
+                   zmlt(FP[0],FP[1],PHI[0],PHI[1],cr,ci);\r
+                   ebase = Math.exp(cr[0]);\r
+                   expr = ebase * Math.cos(ci[0]);\r
+                   expi = ebase * Math.sin(ci[0]);\r
+                   zmlt(diff[0],diff[1],expr,expi,cr,ci);\r
+                   PHYPT[IP-1][0] = CENTR[0] + cr[0];\r
+                   PHYPT[IP-1][1] = CENTR[1] + ci[0];\r
+                   continue loopIP;\r
                } // else if (DIST < EFPTL)\r
                else {\r
-       C\r
-       C           SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS\r
-       C           PARTICULAR POINT WW.\r
-       C\r
-       C           INITIALISE SOME DATA\r
-       C\r
-                   DGC=DGPOC(IA)\r
-                   IF (DGC+1 .GT. MNCOF) THEN\r
-                     IER=48\r
-                     RETURN\r
-                   ENDIF\r
-                   JTC=JTYPC(IA)\r
-                   AJTC=ABS(JTC)\r
-                   BETAC=JAINC(AJTC)\r
-                   LOMC=LSUBC(IA)\r
-                   LODC=(AJTC-1)*NQPTS+1\r
-                   SJTC=SIGN(1E+0,REAL(JTC))\r
-                   SCO=SJTC\r
-       C\r
-                   DO 100 J=1,DGC+1\r
-                     J1=LOMC+J-1\r
-                     SCO=SCO*SJTC\r
-                     JCOFC(J)=SOLNC(J1)*SCO\r
-       100         CONTINUE\r
-       C\r
-                   XI=SJTC*CMPLX(REXI,IMXI)            \r
-                   CALL PPSBI1(XI,BETAC,NQPTS,DGC,ACOFC(LODC),BCOFC(LODC),\r
-            +                  H0VLC(AJTC),JCOFC,LGTOL,TOLOU,XENPT,QINTS,\r
-            +                  MQIN1,IER)\r
-                   IF (IER .GT. 0) THEN\r
-                     IF (IER .EQ. 29) THEN\r
-                       IER=49\r
-                     ENDIF\r
-                     RETURN\r
-                   ENDIF\r
-                   NQUAD=QINTS*NQPTS\r
-                   IF (NQUAD .GT. MXNQD) THEN\r
-                     IER=47\r
-                     RETURN\r
-                   ENDIF\r
-                   K=0\r
-                   SUM1=BETAC+1E+0\r
-                   DO 130 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
-                       DO 110 J=1,NQPTS\r
-                         J1=LODC+J-1\r
-                         K=K+1\r
-                         TT=MEAN+RR*QUPTC(J1)\r
-                         WSPEC(K)=RRB*QUWTC(J1)*JACSUC(TT,DGC,ACOFC(LODC),\r
-            +                     BCOFC(LODC),H0VLC(AJTC),JCOFC)\r
-                         TT=MD+TT*SJTC*HA\r
-                         ZSPEC(K)=CMPLX(COS(TT),SIN(TT))\r
-       110             CONTINUE\r
-                     ELSE\r
-                       DO 120 J=1,NQPTS\r
-                         J1=LOL+J\r
-                         K=K+1\r
-                         TT=MEAN+RR*QUPTC(J1)\r
-                         WSPEC(K)=RR*QUWTC(J1)*(1E+0+TT)**BETAC*JACSUC(TT,\r
-            +                     DGC,ACOFC(LODC),BCOFC(LODC),H0VLC(AJTC),\r
-            +                     JCOFC)\r
-                         TT=MD+TT*SJTC*HA\r
-                         ZSPEC(K)=CMPLX(COS(TT),SIN(TT))\r
-       120             CONTINUE\r
-                     ENDIF\r
-       130         CONTINUE\r
-       C\r
-       C           THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS\r
-       C           AND POINTS WSPEC AND ZSPEC.  NOW ESTIMATE THE INTEGRAL.\r
-       C\r
-                   DO 140 K=1,NQUAD\r
-                     CT=WW/ZSPEC(K)\r
-                     IF (.NOT. INTER) THEN\r
-                       CT=1E+0/CT\r
-                     ENDIF\r
-                     CSUM=CSUM+WSPEC(K)*CLOG(1E+0-CT)\r
-       140         CONTINUE\r
+       \r
+                   // SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS\r
+                   // PARTICULAR POINT WW.\r
+       \r
+                   // INITIALISE SOME DATA\r
+       \r
+                   DGC=DGPOC[IA-1];\r
+                   if (DGC+1 > MNCOF) {\r
+                       IER[0]=48;\r
+                       return;\r
+                   }\r
+                   JTC=JTYPC[IA-1];\r
+                   AJTC=Math.abs(JTC);\r
+                   BETAC=JAINC[AJTC-1];\r
+                   LOMC=LSUBC[IA-1];\r
+                   LODC=(AJTC-1)*NQPTS+1;\r
+                   if (JTC >= 0) {\r
+                       SJTC = 1.0;\r
+                   }\r
+                   else {\r
+                       SJTC = -1.0;\r
+                   }\r
+                   SCO=SJTC;\r
+       \r
+                   for (J=1; J <=DGC+1; J++) {\r
+                       J1=LOMC+J-1; \r
+                       SCO=SCO*SJTC;\r
+                       JCOFC[J-1][0]=SOLNC[J1-1][0]*SCO;\r
+                       JCOFC[J-1][1]=SOLNC[J1-1][1]*SCO;\r
+                   }\r
+       \r
+                   XI[0]= SJTC * REXI;\r
+                   XI[1] = SJTC * IMXI; \r
+                   A = new double[DGC];\r
+                   B = new double[DGC];\r
+                   for (N = 0; N < DGC; N++) {\r
+                       A[N] = ACOFC[LODC-1+N];\r
+                       B[N] = BCOFC[LODC-1+N];\r
+                   }\r
+                   PPSBI1(XI,BETAC,NQPTS,DGC,A,B,\r
+                          H0VLC[AJTC-1],JCOFC,LGTOL,TOLOU,XENPT,QINTS,\r
+                          MQIN1,IER);\r
+                   if (IER[0] > 0) {\r
+                       if (IER[0] == 29) {\r
+                           IER[0]=49;\r
+                       }\r
+                       return ;\r
+                   } // if (IER[0] > 0)\r
+                   NQUAD=QINTS[0]*NQPTS;\r
+                   if (NQUAD > MXNQD) {\r
+                         IER[0]=47;\r
+                         return;\r
+                   }\r
+                   K=0;\r
+                   SUM1=BETAC+1.0;\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
+                           for (J=1; J <= NQPTS; J++) {\r
+                               J1=LODC+J-1;\r
+                               K=K+1;\r
+                               TT=MEAN+RR*QUPTC[J1-1];\r
+                               A = new double[DGC];\r
+                               B = new double[DGC];\r
+                               for (N = 0; N < DGC; N++) {\r
+                                       A[N] = ACOFC[LODC-1+N];\r
+                                       B[N] = BCOFC[LODC-1+N];\r
+                               }\r
+                               JOUT = JACSUC(TT,DGC,A,B,H0VLC[AJTC-1],JCOFC);\r
+                               WSPEC[K-1][0]=RRB*QUWTC[J1-1]*JOUT[0];\r
+                               WSPEC[K-1][1]=RRB*QUWTC[J1-1]*JOUT[1];\r
+                               TT=MD+TT*SJTC*HA;\r
+                               ZSPEC[K-1][0]=Math.cos(TT);\r
+                               ZSPEC[K-1][1]=Math.sin(TT);\r
+                           } // for (J=1; J <= NQPTS; J++)\r
+                       } // if (I == 1)\r
+                       else  { // I != 1\r
+                           for (J=1; J <= NQPTS; J++) {\r
+                               J1=LOL+J;\r
+                               K=K+1;\r
+                               TT=MEAN+RR*QUPTC[J1-1];\r
+                               A = new double[DGC];\r
+                               B = new double[DGC];\r
+                               for (N = 0; N < DGC; N++) {\r
+                                       A[N] = ACOFC[LODC-1+N];\r
+                                       B[N] = BCOFC[LODC-1+N];\r
+                               }\r
+                               JOUT = JACSUC(TT,DGC,A,B,H0VLC[AJTC-1],JCOFC);\r
+                               var = RR*QUWTC[J1-1]*Math.pow((1.0+TT), BETAC);\r
+                               WSPEC[K-1][0] = var * JOUT[0];\r
+                               WSPEC[K-1][1] = var * JOUT[1];\r
+                               TT=MD+TT*SJTC*HA;\r
+                               ZSPEC[K-1][0]=Math.cos(TT);\r
+                               ZSPEC[K-1][1]=Math.sin(TT);\r
+                           } // for (J=1; J <= NQPTS; J++)\r
+                       } // else I != 1\r
+                   } // for (I=1; I <= QINTS[0]; I++)\r
+       \r
+                   // THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS\r
+                   // AND POINTS WSPEC AND ZSPEC.  NOW ESTIMATE THE INTEGRAL.\r
+       \r
+                   for (K=1; K <= NQUAD; K++) {\r
+                       zdiv(WW[0],WW[1],ZSPEC[K-1][0],ZSPEC[K-1][1],cr,ci);\r
+                       CT[0] = cr[0];\r
+                       CT[1] = ci[0];\r
+                       if (!INTER) {\r
+                               zdiv(1.0,0.0,CT[0],CT[1],cr,ci);\r
+                           CT[0] = cr[0];\r
+                           CT[1] = ci[0];\r
+                       }\r
+                       logr = Math.log(zabs(1.0-CT[0],-CT[1]));\r
+                       logi = Math.atan2(-CT[1], 1.0-CT[0]);\r
+                       zmlt(WSPEC[K-1][0],WSPEC[K-1][1],logr,logi,cr,ci);\r
+                       CSUM[0] = CSUM[0] + cr[0];\r
+                       CSUM[1] = CSUM[1] + ci[0];\r
+                   } // for (K=1; K <= NQUAD; K++)\r
        \r
                     // END OF ELSE BLOCK RELATING TO SPECIAL QUADRATURE RULE FOR\r
                     // WW NEAR ARC IA\r
@@ -13904,11 +13989,557 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
        \r
            // END OF MAP CALCULATION FOR FIELD POINT NUMBER IP\r
        \r
-       } // loopIP: for (IP=1; IP <= NPTS; IP++) */\r
+       } // loopIP: for (IP=1; IP <= NPTS; IP++) \r
        \r
        IER[0]=0;\r
        \r
     } // private void CATPH4\r
+    \r
+    private void BMPHYC(int IARC, double PHYPT[], double CANPT[], double DERIV[], boolean WANTD,\r
+        boolean WANTM, int IER[]) {\r
+       \r
+       //INTEGER IARC,IER\r
+       //INTEGER IGEOM(*),ISNPH(*),IQUPH(*)\r
+       //REAL RGEOM(*),RSNPH(*),RQUPH(*)\r
+       //COMPLEX PHYPT,CANPT,DERIV\r
+       //COMPLEX ZQUPH(*)\r
+       //LOGICAL WANTD,INTER,WANTM\r
+       \r
+       // ......................................................................\r
+       \r
+       // 1.     BMPHYC\r
+       //           BOUNDARY MAPPING FOR THE PHYSICAL --> CANONICAL MAP.\r
+       \r
+       // 2.     PURPOSE\r
+       //           GIVEN A POINT ON THE BOUNDARY OF THE PHYSICAL DOMAIN, THIS\r
+       //           ROUTINE COMPUTES THE CORRESPONDING APPROXIMATE IMAGE POINT\r
+       //           ON THE UNIT DISC AND, IF REQUESTED, ALSO COMPUTES THE \r
+       //           DERIVATIVE OF THE MAP : PHYSICAL --> CANONICAL AT THE GIVEN\r
+       //           POINT ON THE PHYSICAL BOUNDARY.\r
+       \r
+       // 3.     CALLING SEQUENCE\r
+       //           CALL BMPHYC(IARC,PHYPT,CANPT,DERIV,WANTD,INTER,IGEOM,RGEOM,\r
+       //                       ISNPH,RSNPH,IQUPH,RQUPH,ZQUPH,WANTM,IER)\r
+       \r
+       //        PARAMETERS\r
+       //         ON ENTRY\r
+       //            IARC   - INTEGER\r
+       //                     ALLOWS TWO MODES OF DEFINING THE POINT ON THE \r
+       //                     BOUNDARY OF THE PHYSICAL DOMAIN.  IF IARC > 0 THEN\r
+       //                     THE PHYSICAL POINT LIES ON ANALYTIC ARC NUMBER \r
+       //                     IARC (AS DEFINED IN THE PARAMETRIC FUNCTION \r
+       //                     PARFUN).  IF IARC <= 0 THEN THE PARTICULAR ARC ON\r
+       //                     WHICH THE PHYSICAL POINT LIES IS CONSIDERED TO BE\r
+       //                     UNKNOWN ON ENTRY.\r
+       \r
+       //            PHYPT  - COMPLEX\r
+       //                     IF IARC > 0 THEN PHYPT IS THE (COMPLEX) PARAMETER\r
+       //                     VALUE WHICH DEFINES THE PHYSICAL POINT ON ANALYTIC\r
+       //                     ARC NUMBER IARC.  IF IARC <= 0 THEN PHYPT IS THE\r
+       //                     GIVEN PHYSICAL POINT.\r
+       \r
+       //            WANTD  - LOGICAL\r
+       //                     IF WANTD IS TRUE THEN DERIV IS COMPUTED OTHERWISE\r
+       //                     DERIV ISN'T COMPUTED.\r
+       \r
+       //            INTER  - LOGICAL\r
+       //                     TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE \r
+       //                     OTHERWISE. (AS PREVIOUSLY USED IN JAPHYC, GQPHYC)\r
+       \r
+       //            IGEOM  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR IGEOM PREVIOUSLY SET UP BY \r
+       //                     JAPHYC.\r
+       \r
+       //            RGEOM  - REAL ARRAY\r
+       //                     THE REAL VECTOR RGEOM PREVIOUSLY SET UP BY JAPHYC.\r
+       \r
+       //            ISNPH  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR ISNPH PREVIOUSLY SET UP BY \r
+       //                     JAPHYC.\r
+       \r
+       //            RSNPH  - REAL ARRAY\r
+       //                     THE REAL VECTOR RSNPH PREVIOUSLY SET UP BY JAPHYC.\r
+       \r
+       //            IQUPH  - INTEGER ARRAY\r
+       //                     THE INTEGER VECTOR IQUPH PREVIOUSLY SET UP BY\r
+       //                     GQPHYC.\r
+       \r
+       //            RQUPH  - REAL ARRAY\r
+       //                     THE REAL ARRAY PREVIOUSLY SET UP BY GQPHYC.\r
+       \r
+       //            ZQUPH  - COMPLEX ARRAY\r
+       //                     THE COMPLEX ARRAY PREVIOUSLY SET UP BY GQPHYC.\r
+       \r
+       //            WANTM  - LOGICAL\r
+       //                     IF WANTM IS TRUE THEN, ON AN ABNORMAL EXIT, AN\r
+       //                     ERROR MESSAGE IS WRITTEN ON THE STANDARD OUTPUT\r
+       //                     CHANNEL.  IF WANTM IS FALSE THEN NO MESSAGE IS\r
+       //                     WRITTEN.\r
+       \r
+       \r
+       //         ON EXIT\r
+       //            CANPT  - COMPLEX\r
+       //                     THE POINT ON THE UNIT DISC CORRESPONDING TO THE\r
+       //                     GIVEN PHYSICAL POINT UNDER THE MAP : PHYSICAL -->\r
+       //                     CANONICAL.\r
+       \r
+       //            DERIV  - COMPLEX\r
+       //                     THE COMPUTED DERIVATIVE OF THE MAP : PHYSICAL --> \r
+       //                     CANONICAL AT THE GIVEN POINT ON THE PHYSICAL BOUN-\r
+       //                     DARY.  THIS IS COMPUTED ONLY IF WANTD IS TRUE.\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
+       //           - NOTE THAT THIS ROUTINE CAN ONLY BE USED  A F T E R  THE  \r
+       //             ROUTINES JAPHYC AND GQPHYC HAVE SUCCESSFULLY EXECUTED, \r
+       //             AND THAT MANY INPUT ARGUMENTS FOR BMPHYC ARE OUTPUT VALUES\r
+       //             FROM JAPHYC AND GQPHYC.\r
+       \r
+       // ......................................................................\r
+       //     AUTHOR: DAVID HOUGH, ETH, ZUERICH\r
+       //     LAST UPDATE: 3 JULY 1990\r
+       // ......................................................................\r
+            \r
+       //     LOCAL VARAIBLES\r
+       \r
+       int MNSUA,TNSUA;\r
+       String IERTXT;\r
+       //CHARACTER*6 IERTXT\r
+       \r
+       // EXTERNAL BMPHC1,IERTXT\r
+       \r
+       NARCS=ISNPH[0];\r
+       NJIND=NARCS+1;\r
+       NQPTS=ISNPH[1];\r
+       TNSUA=ISNPH[2];\r
+       MNSUA=ISNPH[4];\r
+       MNEQN=ISNPH[5];\r
+       MQUPH=IQUPH[3];\r
+       TNGQP=NQPTS*NJIND;\r
+       \r
+        // SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC\r
+       \r
+       // SET UP THE POINTERS TO  ISNPH AND RSNPH, AS IN JAPHYC\r
+       \r
+       // SET UP POINTERS TO IQUPH AND RQUPH, AS IN GQPHYC\r
+       \r
+       // GET REQUIRED MAP AND DERIVATIVE\r
+       \r
+        BMPHC1(IARC,PHYPT,CANPT,DERIV,NQPTS,TNSUA,DGPOL,\r
+            JATYP,LOSUB,LQSBF,NPPQF,PARNT,\r
+            AICOF,ACOEF,BICOF,BCFSN,BCOEF,\r
+            H0VAL,HIVAL,HALEN,JACIN,MIDPT,\r
+            SOLUN,TPPQF,TRRAD,VTARG,ZPPQF,\r
+            INTER,WANTD,IER);\r
+       \r
+       // SEND ERROR MESSAGE TO STANDARD OUTPUT OF NECESSARY\r
+       \r
+       if (IER[0] > 0 && WANTM) System.out.println(IERTXT(IER[0]));\r
+       \r
+    } // private void BMPHYC\r
+    \r
+    private void BMPHC1(int IARC, double PHYPT[], double CANPT[], double DERIV[], int NQPTS,\r
+        int TNSUA, int DGPOL[], int JATYP[], int LOSUB[], int LPQSB[], int NPPQF[], int PARNT[],\r
+        double A1COF[], double ACOEF[], double B1COF[], double BCFSN[], double BCOEF[],\r
+        double H0VAL[], double H1VAL[], double HALEN[], double JACIN[], double MIDPT[], \r
+        double SOLUN[], double TPPQF[], double TRRAD[], double VTARG[], double ZPPQF[][],\r
+        boolean INTER, boolean WANTD, int IER[]) {\r
+       //INTEGER IARC,IER,NQPTS,TNSUA\r
+       //INTEGER DGPOL(*),JATYP(*),LOSUB(*),LPQSB(*),NPPQF(*),PARNT(*)\r
+       //REAL A1COF(*),ACOEF(*),B1COF(*),BCFSN(*),BCOEF(*),H0VAL(*),\r
+       //+H1VAL(*),HALEN(*),JACIN(*),MIDPT(*),SOLUN(*),TPPQF(*),TRRAD(*),\r
+       //+VTARG(*)\r
+       //COMPLEX CANPT,DERIV,PHYPT,ZPPQF(*)\r
+       //LOGICAL INTER,WANTD\r
+       \r
+       //     GIVEN A POINT (DEFINED BY IARC AND PHYPT AS EXPLAINED NEXT) ON\r
+       //     THE BOUNDARY OF THE PHYSICAL DOMAIN, TO COMPUTE ITS IMAGE CANPT\r
+       //     ON THE UNIT DISC.  IN CASE WANTD=.TRUE. THEN ALSO COMPUTE THE \r
+       //     DERIVATIVE DERIV OF THE MAP ONTO THE DISC AT THE GIVEN BOUNDARY \r
+       //     POINT.\r
+       \r
+       //     IF IARC > 0 THEN PHYPT IS THE PARAMETER VALUE OF THE \r
+       //     PHYSICAL POINT ON THE PARENT ANALYTIC ARC NUMBER IARC OTHERWISE\r
+       //     PHYPT DEFINES THE COORDINATES OF A PHYSICAL POINT SOMEWHERE ON \r
+       //     THE BOUNDARY.\r
+       \r
+       //     IER=0  - NORMAL EXIT.\r
+       //     IER=44 - LOCAL PARAMETER MNCOF NEEDS INCREASING.\r
+       //     IER=45 - THE PHYSICAL POINT DEFINED BY PHYPT (WITH IARC<0) HAS NOT\r
+       //              BEEN DETECTED AS LYING ON THE BOUNDARY; ACTUALLY, IT HAS\r
+       //              NOT BEEN DETECTED AS LYING PATHOLOGICALY CLOSE TO THE\r
+       //              BOUNDARY.  CHECK THAT PHYPT IS CORRECT AND IF IT IS THEN \r
+       //              CONSIDER INCREASING THE PATHOLOGICAL TOLERANCE PARAMETER \r
+       //              PTHTL IN THE FIRST LINE BELOW.\r
+       \r
+       //     LOCAL VARIABLES\r
+       \r
+       final int MNCOF = 32;\r
+       int AJT,DEG,I,I1,IA,K,J1,JQ,JT,LOD,LOM,NQ,PT;\r
+       double BETA,DIST,HL,JACSUM,MD,MINDS,NEWTL,PHI,PTHTL,R1MACH,RBCF,\r
+            SJT,SUM,TOLSM,TT,TUPI,TXI;\r
+       double BCF[] = new double[2];\r
+       double CT[] = new double[2];\r
+       double DIFF2[] = new double[2];\r
+       double PSI[] = new double[2];\r
+       double WGHT[] = new double[2];\r
+       double XI[] = new double[2];\r
+       double ZXI[] = new double[2];\r
+       double ZZ[] = new double[2];\r
+       //COMPLEX BCF,CJACSU,CT,DIFF2,DPARFN,PARFUN,PSI,WGHT,XI,ZXI,\r
+       //+ZTOB1,ZZ\r
+       \r
+       double JACOF[] = new double[MNCOF];\r
+       double A[];\r
+       double B[];\r
+       double CO[];\r
+       int N;\r
+       double DOUT[];\r
+       double POUT[];\r
+       double base;\r
+       double cr[] = new double[1];\r
+       double ci[] = new double[1];\r
+       double cr2[] = new double[1];\r
+       double ci2[] = new double[1];\r
+       double ZIN[] = new double[2];\r
+       \r
+       // EXTERNAL CJACSU,DPARFN,JACSUM,PARFUN,R1MACH,ZTOB1\r
+       \r
+       NEWTL=Math.sqrt(EPS);\r
+       PTHTL=NEWTL;\r
+       TOLSM=10.0*EPS;\r
+       TUPI=2.0 * Math.PI;\r
+       \r
+       /*if (IARC > 0) {\r
+       \r
+           // *PHYPT* IS THE PARAMETER VALUE OF THE PHYSICAL POINT ON THE\r
+           // PARENT ANALYTIC ARC NUMBER *IARC*.\r
+        \r
+           // FIRST FIND THE CORRESPONDING SUBARC NUMBER AND LOCAL PARAMETER.\r
+       \r
+           TT=PHYPT[0];\r
+           if (TT < -1.0) {\r
+               TT=-1.0;\r
+           }\r
+           if (TT > 1.0) {\r
+               TT=1.0;\r
+           }\r
+           I=1;\r
+           while (true) {\r
+               PT=PARNT[I-1];\r
+               HL=HALEN[I-1];\r
+               MD=MIDPT[I-1];\r
+               SUM=MD+HL;\r
+       \r
+               if (Math.abs(SUM-1.0) < TOLSM) {\r
+                   SUM=1.0;\r
+               }\r
+       \r
+               if (PT == IARC && TT <= SUM) {\r
+                   IA=I;\r
+                   TT=(TT-MD)/HL;\r
+                   break;\r
+               }\r
+               else {\r
+                   I=I+1;\r
+                   continue;\r
+               }\r
+           } // while (true)\r
+       \r
+           if (TT < -1.0) {\r
+               TT=-1.0;\r
+           }\r
+           if (TT > 1.0) {\r
+               TT=1.0;\r
+           }\r
+       \r
+           JT=JATYP[IA-1];\r
+           if (JT >= 0) {\r
+               SJT = 1.0;\r
+           }\r
+           else {\r
+               SJT = -1.0;\r
+           }\r
+           AJT=Math.abs(JT);\r
+           BETA=JACIN[AJT-1];\r
+           DEG=DGPOL[IA-1];\r
+           if (DEG+1 > MNCOF) {\r
+               IER[0]=44;\r
+               return;\r
+           }\r
+           LOM=LOSUB[IA-1];\r
+           LOD=(AJT-1)*NQPTS+1;\r
+       \r
+           TT=SJT*TT;\r
+           A = new double[DEG-1];\r
+           B = new double[DEG-1];\r
+           for (N = 0; N < DEG-1; N++) {\r
+               A[N] = A1COF[LOD-1+N];\r
+               B[N] = B1COF[LOD-1+N];\r
+           }\r
+           CO = new double[DEG];\r
+           for (N = 0; N < DEG; N++) {\r
+               CO[N] = BCFSN[LOM+N];\r
+           }\r
+           PHI=JACSUM(TT,DEG-1,A,B,H1VAL[AJT-1],CO);   \r
+           PHI=Math.pow((1.0+TT),(BETA+1.0))*(BCFSN[LOM-1]-(1.0-TT)*PHI);\r
+          if (JT > 0) {\r
+              RBCF=VTARG[IA-1];\r
+          }\r
+          else {\r
+              RBCF=VTARG[IA];\r
+          }\r
+          RBCF=RBCF+SJT*PHI;\r
+          CANPT[0] = Math.cos(RBCF);\r
+          CANPT[1] = Math.sin(RBCF);\r
+          if (WANTD) {\r
+              if (BETA < 0.0 && (1E+0+TT) <= 0.0) {\r
+       \r
+                   // WE ARE AT A CORNER WITH INFINITE DERIVATIVE.\r
+                  DERIV[0] = Double.MAX_VALUE;\r
+                  DERIV[1] = Double.MAX_VALUE:\r
+              }\r
+              else if (BETA > 0.0 && (1.0+TT) <= 0E+0) {\r
+       \r
+                   // WE ARE AT A CORNER WITH ZERO DERIVATIVE\r
+       \r
+                  DERIV[0] = 0.0;\r
+                  DERIV[1] = 0.0; \r
+              }\r
+              else {\r
+                  for (I=1; I <= DEG+1; I++) {\r
+                      I1=I+LOM-1;\r
+                      JACOF[I-1]=SOLUN[I1-1];\r
+                  } // for (I=1; I <= DEG+1; I++)\r
+                  for (I=2; I <= DEG+1; I +=2) {\r
+                      JACOF[I-1]=SJT*JACOF[I-1];\r
+                  } // for (I=2; I <= DEG+1; I +=2)\r
+                  A = new double[DEG];\r
+                  B = new double[DEG];\r
+                  for (N = 0; N < DEG; N++) {\r
+                       A[N] = ACOEF[LOD-1+N];\r
+                       B[N] = BCOEF[LOD-1+N];\r
+                   }\r
+                   PHI=JACSUM(TT,DEG,A,B,H0VAL[AJT-1],JACOF);\r
+                   DOUT = DPARFN(IAC,PHYPT);\r
+                   base = TUPI * Math.pow((1.0+TT), BETA)*PHI;\r
+                   zmlt(0.0,base,CANPT[0],CANPT[1],cr,ci);\r
+                   zdiv(cr[0],ci[0],DOUT[0]/HL,DOUT[1]/HL,cr2,ci2);\r
+                   DERIV[0] = cr2[0];\r
+                   DERIV[1] = ci2[0];\r
+              } // else\r
+          } // if (WANTD)\r
+       } // if (IARC > 0)\r
+       else { // IARC <= 0\r
+       \r
+           // *PHYPT* IS A POINT SOMEWHERE ON THE PHYSICAL BOUNDARY.\r
+       \r
+           ZZ[0]=PHYPT[0];\r
+           ZZ[1]=PHYPT[1];\r
+           loopIA: for (IA=1; IA <= TNSUA; IA++) {\r
+               PT=PARNT[IA-1];\r
+               JT=JATYP[IA-1];\r
+               NQ=NPPQF[IA-1];\r
+               K=LPQSB[IA-1]-1;\r
+               HL=HALEN[IA-1];\r
+               MD=MIDPT[IA-1];\r
+               for (JQ=1; JQ <= NQ; JQ++) {\r
+                   K=K+1;\r
+                   DIFF2[0]=ZZ[0]-ZPPQF[K-1][0];\r
+                   DIFF2[1]=ZZ[1]-ZPPQF[K-1][1];\r
+                   DIST=zabs(DIFF2[0],DIFF2[1]);\r
+                   if (DIST < TRRAD[K-1]) {\r
+       \r
+                       // ZZ IS CLOSE TO ARC IA.\r
+       \r
+                       J1=JQ;\r
+                       MINDS=DIST;\r
+                       TXI=TPPQF[K-1];\r
+                       ZXI[0]=ZPPQF[K-1][0];\r
+                       ZXI[1]=ZPPQF[K-1][1];\r
+                       while (true) {\r
+                           J1=J1+1;\r
+                           if (J1 <= NQ) {\r
+                               K=K+1;\r
+                               DIFF2[0]=ZZ[0]-ZPPQF[K-1][0];\r
+                               DIFF2[1]=ZZ[1]-ZPPQF[K-1][1];\r
+                               DIST=zabs(DIFF2[0],DIFF2[1]);\r
+                               if (DIST < MINDS) {\r
+                                   MINDS=DIST;\r
+                                   TXI=TPPQF[K-1];\r
+                                   ZXI[0]=ZPPQF[K-1][0];\r
+                                   ZXI[1]=ZPPQF[K-1][1];\r
+                                   continue;\r
+                               } // if (DIST < MINDS)\r
+                           } // if (J1 <= NQ)\r
+                           break;\r
+                       } // while (true)\r
+       \r
+                       // PRELIMINARIES\r
+                 \r
+                       if (JT >= 0) {\r
+                           SJT = 1.0;  \r
+                       }\r
+                       else {\r
+                               SJT = -1.0;\r
+                       }\r
+                       AJT=Math.abs(JT);\r
+                       BETA=JACIN[AJT-1];\r
+                       DEG=DGPOL[IA-1];\r
+                       if (DEG+1 > MNCOF) {\r
+                           IER[0]=44;\r
+                           return;\r
+                       } // if (DEG+1 > MNCOF)\r
+                       LOM=LOSUB[IA-1];\r
+                       LOD=(AJT-1)*NQPTS+1;\r
+       \r
+                       // NOW USE NEWTON'S METHOD TO ESTIMATE THE PARAMETRIC\r
+                       // PRE-IMAGE XI OF ZZ.\r
+       \r
+                       XI[0]=TXI;\r
+                       XI[1] = 0.0;\r
+                       CT[0]=MD+HL*XI[0];\r
+                       CT[1]=HL*XI[1];\r
+                       DOUT = DPARFN(PT,CT);\r
+                       zdiv(ZXI[0]-ZZ[0],ZXI[1]-ZZ[1],DOUT[0]*HL,DOUT[1]*HL,cr,ci);\r
+                       DIFF2[0] = cr[0];\r
+                       DIFF2[1] = ci[0];\r
+                       XI[0]=XI[0]-DIFF2[0];\r
+                       XI[1]=XI[1]-DIFF2[1];\r
+                       while (true) {\r
+                           if (zabs(DIFF2[0],DIFF2[1]) > NEWTL) {\r
+                               CT[0]=MD+HL*XI[0];\r
+                               CT[1] = HL*XI[1];\r
+                               POUT = PARFUN(PT,CT);\r
+                               DOUT = DPARFN(PT,CT);\r
+                               zdiv(POUT[0]-ZZ[0],POUT[1]-ZZ[1],DOUT[0]*HL,DOUT[1]*HL,cr,ci);\r
+                               DIFF2[0] = cr[0];\r
+                               DIFF2[1] = ci[0];\r
+                               XI[0]=XI[0]-DIFF2[0];\r
+                               XI[1]=XI[1]-DIFF2[1];\r
+                               continue;\r
+                           } // if (zabs(DIFF2[0],DIFF2[1]) > NEWTL)\r
+                           else {\r
+       \r
+                               // LAST ITERATION\r
+                               CT[0]=MD+HL*XI[0];\r
+                               CT[1] = HL*XI[1];\r
+                               POUT = PARFUN(PT,CT);\r
+                               DOUT = DPARFN(PT,CT);\r
+                               zdiv(POUT[0]-ZZ[0],POUT[1]-ZZ[1],DOUT[0]*HL,DOUT[1]*HL,cr,ci);\r
+                               DIFF2[0] = cr[0];\r
+                               DIFF2[1] = ci[0];\r
+                               XI[0]=XI[0]-DIFF2[0];\r
+                               XI[1]=XI[1]-DIFF2[1];\r
+                               break;\r
+                           }\r
+                       } // while (true)\r
+                       XI[0]=SJT*XI[0];\r
+                       XI[1]=SJT*XI[1];\r
+       \r
+                       if (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL) {\r
+       \r
+                           // ZZ IS PATHOLOGICALLY CLOSE TO ARC IA AND WE USE THE\r
+                           // CONTINUATION OF THE BOUNDARY CORRESPONDENCE FUNCTION\r
+                           // TO ESTIMATE CANPT AND THE DERIVATIVE\r
+                           A = new double[DEG-1];\r
+                           B = new double[DEG-1];\r
+                           for (N = 0; N < DEG-1; N++) {\r
+                               A[N] = A1COF[LOD-1+N];\r
+                               B[N] = B1COF[LOD-1+N];\r
+                           }\r
+                           CO = new double[DEG];\r
+                           for (N = 0; N < DEG; N++) {\r
+                               CO[N] = BCFSN[LOM+N];\r
+                           }\r
+                           PSI=CJACSU(XI,DEG-1,A,B,H1VAL[AJT-1],CO);\r
+                           ZIN[0] = XI[0]+1.0;\r
+                           ZIN[1] = XI[1];\r
+                           WGHT=ZTOB1(ZIN,BETA+1.0,JT,INTER);\r
+                           zmlt(BCFSN[LOM-1]-(1.0-XI[0]),XI[1],PSI[0],PSI[1],cr,ci);\r
+                           zmlt(WGHT[0],WGHT[1],cr[0],ci[0],cr2,ci2);\r
+                           PSI[0] = cr2[0];\r
+                           PSI[1] = ci2[0];\r
+                           if (JT > 0) {\r
+                               BCF[0]=VTARG[IA-1];\r
+                               BCF[1]= 0.0;\r
+                           }\r
+                           else {\r
+                               BCF[0]=VTARG[IA];\r
+                               BCF[1] = 0.0;\r
+                           }\r
+                           BCF[0]=BCF[0]+SJT*PSI[0];\r
+                           BCF[1]=BCF[1]+SJT*PSI[1];\r
+                           base = Math.exp(-BCF[1]);\r
+                           CANPT[0] = base * Math.cos(BCF[0]);\r
+                           CANPT[1] = base * Math.sin(BCF[0]);\r
+                           if (WANTD) {\r
+                               if (((1.0+XI[0]) == 0.0) && (XI[1] == 0.0) && (BETA < 0.0)) { \r
+       \r
+                                   // WE ARE AT A CORNER WITH INFINITE DERIVATIVE.\r
+       \r
+                                   DERIV[0]= Double.MAX_VALUE;\r
+                                   DERIV[1] = Double.MAX_VALUE;\r
+                               }\r
+                         ELSE IF ((1E+0+XI).EQ.(0E+0,0E+0) .AND. BETA.GT.0E+0) \r
+            +            THEN\r
+       C\r
+       C                   WE ARE AT A CORNER WITH ZERO DERIVATIVE.\r
+       C\r
+                           DERIV=(0E+0,0E+0)\r
+                         ELSE\r
+                           DO 55 I=1,DEG+1\r
+                             I1=I+LOM-1\r
+                             JACOF(I)=SOLUN(I1)\r
+       55                  CONTINUE\r
+                           DO 60 I=2,DEG+1,2\r
+                             JACOF(I)=SJT*JACOF(I)\r
+       60                  CONTINUE\r
+                           PSI=CJACSU(XI,DEG,ACOEF(LOD),BCOEF(LOD),H0VAL(AJT),\r
+            +                         JACOF)\r
+                           CT=MD+HL*SJT*XI\r
+                           WGHT=WGHT/(1E+0+XI)\r
+                           DERIV=TUPI*WGHT*PSI*(0E+0,1E+0)*CANPT\r
+            +                    /DPARFN(PT,CT)/HL\r
+                         ENDIF\r
+                           } // if (WANTD)\r
+                       // NORMAL EXIT\r
+                       IER[0] = 0;\r
+                       return;\r
+                       \r
+                       } // if (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL)\r
+       \r
+                   // END OF *IF (DIST .LT. TRRAD(K)) THEN* FOLLOWS\r
+       \r
+                   } // if (DIST < TRRAD[K-1])\r
+               } // for (JQ=1; JQ <= NQ; JQ++)\r
+           } // loopIA: for (IA=1; IA <= TNSUA; IA++)\r
+           IER[0]=45;\r
+           return;\r
+       } // else IARC <= 0\r
+       \r
+       300   CONTINUE\r
+       \r
+       // NORMAL EXIT\r
+       \r
+       IER[0]=0; */\r
+    } // private void BMPHC1\r
+\r
+\r
 \r
 \r
       /**\r