Porting continues.
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 18 Jan 2018 22:57:20 +0000 (22:57 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Thu, 18 Jan 2018 22:57:20 +0000 (22:57 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15346 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 1eba8bc..ba4b1e1 100644 (file)
@@ -10600,9 +10600,11 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
                int I,MNCOF,MNEQN,MNSUA,MNSUC,NEQNS,NJCOG,PT,TNSUA,TNSUC;\r
                double COLSC[] = new double[]{-1.0};\r
-               double INNRAD,LGTOL,PI,SUPER,THET0,NEWTL;\r
+               double INNRAD = 0.0;\r
+               double LGTOL,PI,SUPER,THET0,NEWTL;\r
                double DCAP0[] = new double[2];\r
                double FACTR[] = new double[2];\r
+               double mult;\r
                //COMPLEX DCAP0,FACTR,CINRAD\r
                \r
                // EXTERNAL BCFSNG,CINRAD,IGNLVL,JCFIM5,OPQUD1,R1MACH,OUPTCA,OUPTCL,\r
@@ -10644,6 +10646,7 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                // SET UP POINTERS TO ELEMENTS IN IQUPH AND RQUPH, AS IN GQPHYC\r
                \r
                // SET UP POINTERS TO ELEMENTS IN ISNCA, RSNCA AND ZSNCA\r
+               // ZSNCA called CSNCA elsewhere\r
                \r
            // INITIALISE JACOBI INDECES *JAINC* FOR THE INVERSE MAP\r
                \r
@@ -10728,41 +10731,45 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                THET0=VTARG[0];\r
                FACTR[0] = Math.cos(THET0);\r
                FACTR[1] = Math.sin(THET0);\r
-               /*      DCAP0=CINRAD(NARCS,NQPTS,TNSUA,ISNPH(DGPOL),ISNPH(JATYP),\r
-                    +ISNPH(LOSUB),IQUPH(LQSBF),IQUPH(NPPQF),IGEOM(PARNT),RSNPH(ACOEF),\r
-                    +RSNPH(BCOEF),RSNPH(H0VAL),RGEOM(HALEN),RSNPH(JACIN),LGTOL,\r
-                    +RGEOM(MIDPT),RSNPH(QUPTS),RSNPH(QUWTS),RSNPH(SOLUN),RQUPH(TPPQF),\r
-                    +RQUPH(TRRAD),RQUPH(WPPQF),CENTR,FACTR,ZQUPH(ZPPQF),IER)\r
-                     IF (IER .GT. 0) THEN\r
-                       GOTO 999\r
-                     ENDIF\r
-               C\r
-                     IF (INTER) THEN\r
-                       INNRAD=ABS(DCAP0)\r
-                       WRITE(*,2) 'INNER RADIUS:',INNRAD\r
-                     ELSE\r
-                       DCAP0=EXP(-RSNPH(SOLUN+NEQNS-1))*DCAP0/ABS(DCAP0)\r
-                     ENDIF\r
-                     ZSNCA(1)=DCAP0\r
-               C\r
-               C     GET THE COEFFICIENTS BFSNC FOR THE COMPLEX BOUNDARY CORRESPONDENCE\r
-               C     FUNCTION FOR THE INVERSE MAP.\r
-               C\r
-                     CALL BCFSNG(TNSUC,ISNCA(DGPOC),ISNCA(JTYPC),ISNCA(LSUBC),\r
-                    +RSNCA(H0VLC),RSNCA(JAINC),ZSNCA(BFSNC),ZSNCA(SOLNC))\r
-               C\r
-                     CALL OUPTCL(ISNCA(DGPOC),ISNCA(JTYPC),LGTOL,ISNCA(LSUBC),NQPTS,\r
-                    +CHNL,IGEOM(PARNT),ISNCA(PRNSA),RWORK,ZSNCA(SOLNC),TNSUC,INTER,\r
-                    +INNRAD,IER)\r
-               C\r
-               C     OUTPUT ALL RESULTS REQUIRED FOR SUBSEQUENT PROCESSING\r
-               C\r
-                     CALL OUPTCA(ISNCA,RSNCA,ZSNCA,CHNL)\r
-               C\r
-               999   CONTINUE\r
-               C\r
-                     CALL WRTAIL(3,0,IER)\r
-               C */\r
+               DCAP0=CINRAD(NARCS,NQPTS,TNSUA,DGPOL,JATYP,\r
+                    LOSUB,LQSBF,NPPQF,PARNT,ACOEF,\r
+                    BCOEF,H0VAL,HALEN,JACIN,LGTOL,\r
+                    MIDPT,QUPTS,QUWTS,SOLUN,TPPQF,\r
+                    TRRAD,WPPQF,CENTR,FACTR,ZPPQF,IER);\r
+               if (IER[0] > 0) {\r
+                       WRTAIL(3,0,IER[0],null);\r
+                       return;\r
+               }\r
+               \r
+               if (INTER) {\r
+                   INNRAD=zabs(DCAP0[0],DCAP0[1]);\r
+                   System.out.println("INNER RADIUS: " + INNRAD);\r
+                   Preferences.debug("INNER RADIUS: " + INNRAD + "\n", Preferences.DEBUG_ALGORITHM);\r
+               }\r
+               else {\r
+                   mult=Math.exp(-SOLUN[NEQNS-1])/zabs(DCAP0[0],DCAP0[1]);\r
+                   DCAP0[0] = DCAP0[0] * mult;\r
+                   DCAP0[1] = DCAP0[1] * mult;\r
+               }\r
+               CSNCA[0]=DCAP0[0];\r
+               CSNCA[1]=DCAP0[1];\r
+               \r
+               // GET THE COEFFICIENTS BFSNC FOR THE COMPLEX BOUNDARY CORRESPONDENCE\r
+               // FUNCTION FOR THE INVERSE MAP.\r
+               \r
+               BCFSNG(TNSUC,DGPOC,JTYPC,LSUBC,\r
+                    H0VLC,JAINC,BFSNC,SOLNC);\r
+               \r
+               OUPTCL(DGPOC,JTYPC,LGTOL,LSUBC,NQPTS,\r
+                    PARNT,PRNSA,RWORK,SOLNC,TNSUC,INTER,\r
+                    INNRAD,IER[0]);\r
+               \r
+               // OUTPUT ALL RESULTS REQUIRED FOR SUBSEQUENT PROCESSING\r
+               \r
+               // CALL OUPTCA(ISNCA,RSNCA,ZSNCA,CHNL)\r
+               \r
+               WRTAIL(3,0,IER[0],null);\r
+           return;\r
     } // private void JACANP\r
 \r
 \r
@@ -11487,6 +11494,564 @@ public class SymmsIntegralMapping extends AlgorithmBase  {
                \r
     } // private void RHOFN\r
 \r
+   // COMPLEX FUNCTION CINRAD\r
+   private double[] CINRAD(int NARCS, int NQPTS, int TNSUA, int DGPOL[], int JATYP[],\r
+       int LOSUB[], int LPQSB[], int NPPQF[], int PARNT[], double ACOEF[], double BCOEF[],\r
+       double H0VAL[], double HALEN[], double JACIN[], double LGTOL, double MIDPT[],\r
+       double QUPTS[], double QUWTS[], double SOLUN[], double TPPQF[], double TRRAD[],\r
+       double WPPQF[], double CENTR[], double FACTR[], double ZPPQF[][], int IER[]) {\r
+          // INTEGER IER,NARCS,NQPTS,TNSUA\r
+          // INTEGER DGPOL(*),JATYP(*),LOSUB(*),LPQSB(*),NPPQF(*),PARNT(*)\r
+          // REAL LGTOL\r
+          // REAL ACOEF(*),BCOEF(*),\r
+          // +H0VAL(*),HALEN(*),JACIN(*),MIDPT(*),QUPTS(*),QUWTS(*),SOLUN(*),\r
+          // +TPPQF(*),TRRAD(*),WPPQF(*)\r
+          // COMPLEX CENTR,FACTR\r
+          // COMPLEX ZPPQF(*)\r
+               \r
+          //     TO COMPUTE THE COMPLEX INNER RADIUS (I.E. THE RECIPROCAL OF THE\r
+          //     DERIVATIVE OF THE INTERIOR MAP AT THE CENTRE POINT OF THE PHYSICAL\r
+          //     DOMAIN)\r
+               \r
+          //      IER=0  - NORMAL EXIT\r
+          //     IER=37 - LOCAL PARAMETER MXNQD NEEDS INCREASING\r
+          //     IER=38 - LOCAL PARAMETER MNCOF NEEDS INCREASING\r
+          //     IER=39 - THE CENTRE POINT IS PATHOLOGICALLY CLOSE TO THE\r
+          //              BOUNDARY\r
+          //     IER=40 - LOCAL PARAMETER MQIN1 MUST BE INCREASED    \r
+       \r
+          // LOCAL VARIABLES\r
+               \r
+          final int MNCOF = 32;\r
+          final int MQIN1 = 11;\r
+          final int MXNQD = 80;\r
+          int QINTS[] = new int[1];\r
+          int AJT,DEG,I,IA,K,J,J1,J2,JQ,JT,LIM,LOD,LOL,LOM,\r
+                    NQ,NQUAD,PT;\r
+          final double PTHTL = 1.0E-3;\r
+          final double LIMIT = 2.3562;\r
+          double TOLOU[] = new double[1];\r
+          double THET1 = 0.0;\r
+          double RT1 = 0.0;\r
+          double CURARG = 0.0;\r
+          double STARG = 0.0;\r
+          double STRT1 = 0.0;\r
+          double STTH1 = 0.0;\r
+          double AISUM,ANGLE,ARGBR,ARSUM,BETA,DIST,HL,ISUM,\r
+                    MD,MEAN,MINDS,NEWTL,PI,RR,RRB,RSUM,\r
+                    RT2,SCO,SS,SUM1,THET2,TT,TXI,TUPI,WT;\r
+          double CT[] = new double[2];\r
+          double XI[] = new double[2];\r
+          double DIFF1[] = new double[2];\r
+          double DIFF2[] = new double[2];\r
+          double STDF1[] = new double[2];\r
+          double ZXI[] = new double[2];\r
+          double ZZ[] = new double[2];\r
+          // COMPLEX CT,DPARFN,PARFUN,XI,DIFF1,DIFF2,\r
+          // +STDF1,ZXI,ZZ\r
+          boolean FIRST;\r
+          // EXTERNAL ARGIN1,DPARFN,JACSUM,PARFUN,PPSBI1\r
+          double JACOF[] = new double[MNCOF];\r
+          double TSPEC[] = new double[MXNQD];\r
+          double WSPEC[] = new double[MXNQD];\r
+          double XENPT[] = new double[MQIN1];\r
+          double JCOFC[][] = new double[MNCOF][2];\r
+          double ZSPEC[][] = new double[MXNQD][2];\r
+          // COMPLEX JCOFC(MNCOF),ZSPEC(MXNQD)\r
+          double cr[] = new double[1];\r
+          double ci[] = new double[1];\r
+          double DIFF1IN[] = new double[2];\r
+          double DIFF2IN[] = new double[2];\r
+          double result[] = new double[2];\r
+          double DOUT[];\r
+          double POUT[];\r
+          int I2;\r
+          double PIN[] = new double[2];\r
+          double expCT;\r
+          \r
+          NEWTL=Math.sqrt(EPS);\r
+          PI=Math.PI;\r
+          TUPI=2.0*PI;\r
+          LOL=NARCS*NQPTS;\r
+          ZZ[0]=CENTR[0];\r
+          ZZ[1]=CENTR[1];\r
+          RSUM=0.0;\r
+          ISUM=0.0;\r
+          FIRST=true;\r
+          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
+                  ARSUM=0.0;\r
+                  AISUM=0.0;\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
+                      RT2=MD+HL*TPPQF[K-1];\r
+                      DIST=zabs(DIFF2[0],DIFF2[1]);\r
+                      if (DIST >= TRRAD[K-1]) {\r
+                          WT=WPPQF[K-1];\r
+                          if (WT != 0.0) {\r
+                              ARSUM=ARSUM+WT*Math.log(DIST);\r
+                              if (FIRST) {\r
+                                  CURARG=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                                  THET2=CURARG;\r
+                                  FIRST=false;\r
+                                  STARG=CURARG;\r
+                              } // if (FIRST)\r
+                              else { // !FIRST\r
+                                  // CT=DIFF2/DIFF1\r
+                                  // CT=DIFF2*CONJG(DIFF1)\r
+                                  // ANGLE=ATAN2(AIMAG(CT),REAL(CT))\r
+                                  THET2=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                                  ANGLE=THET2-THET1;\r
+                                  if (ANGLE <= -PI || ANGLE > PI) {\r
+                                      if (ANGLE > PI) {\r
+                                          ANGLE = ANGLE - TUPI;\r
+                                      }\r
+                                      else {\r
+                                          ANGLE = ANGLE + TUPI;\r
+                                      }\r
+                                  } // if (ANGLE <= -PI || ANGLE > PI)\r
+                                  if (Math.abs(ANGLE) >= LIMIT) {\r
+                                          DIFF1IN[0] = -DIFF1[0];\r
+                                          DIFF1IN[1] = -DIFF1[1];\r
+                                          DIFF2IN[0] = -DIFF2[0];\r
+                                          DIFF2IN[1] = -DIFF2[1];\r
+                                      ANGLE=ARGIN1(RT1,RT2,PT,DIFF1IN,DIFF2IN,ZZ,LIMIT);\r
+                                  }\r
+                                  CURARG=CURARG+ANGLE;\r
+                              } // else !FIRST\r
+                              AISUM=CURARG*WT+AISUM;\r
+                              RT1=RT2;\r
+                              DIFF1[0]=DIFF2[0];\r
+                              DIFF1[1]=DIFF2[1];\r
+                              THET1=THET2;\r
+                          } // if (WT != 0.0)\r
+                      } // if (DIST >= TRRAD[K-1])\r
+                      else { // DIST < TRRAD[K-1]\r
+               \r
+                          // ZZ IS TOO CLOSE TO ARC IA TO USE THE STANDARD RULE.\r
+                          // FIND THE QUADRATURE POINT NEAREST TO ZZ.\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
+                              SS=1.0;\r
+                          }\r
+                          else {\r
+                              SS=-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]=38;\r
+                              return result;\r
+                          }\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
+                              break;\r
+                          } // while (true)\r
+                          XI[0]=SS*XI[0];\r
+                          XI[1]=SS*XI[1];\r
+               \r
+                          if (Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL) {\r
+               \r
+                              // THE CENTRE OF THE DOMAIN (I.E. ZZ) IS PATHOLOGICALLY \r
+                              // CLOSE TO ARC IA AND WE DO NOT ALLOW THIS.\r
+               \r
+                              IER[0]=39;\r
+                              return result;\r
+                          }\r
+                          else {\r
+               \r
+                              // SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS\r
+                              // PARTICULAR POINT ZZ.\r
+               \r
+                              SCO=SS;\r
+                              for (J=1; J <= DEG+1; J++) {\r
+                                  J1=LOM+J-1;\r
+                                  SCO=SCO*SS;\r
+                                  JACOF[J-1]=SOLUN[J1-1]*SCO;\r
+                                  JCOFC[J-1][0]=SOLUN[J1-1]*SCO;\r
+                                  JCOFC[J-1][1] = 0.0;\r
+                              } // for (J=1; J <= DEG+1; J++)\r
+                              double ACOEFIN[] = new double[DEG];\r
+                              double BCOEFIN[] = new double[DEG];\r
+                              for (I = 0; I < DEG; I++) {\r
+                                  ACOEFIN[I] = ACOEF[LOD+I-1];\r
+                                  BCOEFIN[I] = BCOEF[LOD+I-1];\r
+                              }\r
+                              PPSBI1(XI,BETA,NQPTS,DEG,ACOEFIN,BCOEFIN,\r
+                              H0VAL[AJT-1],JCOFC,LGTOL,TOLOU,XENPT,QINTS,MQIN1,IER);\r
+                              if (IER[0] > 0) {\r
+                                  if (IER[0] == 29) {\r
+                                      IER[0]=40;\r
+                                  }\r
+                                  return result;\r
+                              } // if (IER[0] > 0)\r
+                              NQUAD=QINTS[0]*NQPTS;\r
+                              if (NQUAD > MXNQD) {\r
+                                  IER[0]=37;\r
+                                  return result;\r
+                              }\r
+                              K=0;\r
+                              SUM1=BETA+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=LOD+J-1;\r
+                                          K=K+1;\r
+                                          TT=(MEAN+RR*QUPTS[J1-1]);\r
+                                          ACOEFIN = new double[DEG];\r
+                                          BCOEFIN = new double[DEG];\r
+                                          for (I2 = 0; I2 < DEG; I2++) {\r
+                                                  ACOEFIN[I2] = ACOEF[LOD+I2-1];\r
+                                                  BCOEFIN[I2] = BCOEF[LOD+I2-1];\r
+                                          }\r
+                                          WSPEC[K-1]=RRB*QUWTS[J1-1]*JACSUM(TT,DEG,ACOEFIN,\r
+                                             BCOEFIN,H0VAL[AJT-1],JACOF);\r
+                                          TT=TT*SS;\r
+                                          TSPEC[K-1]=MD+TT*HL;\r
+                                          CT[0]=TSPEC[K-1];\r
+                                          CT[1] = 0.0;\r
+                                          ZSPEC[K-1]=PARFUN(PT,CT);\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*QUPTS[J1-1]);\r
+                                          ACOEFIN = new double[DEG];\r
+                                          BCOEFIN = new double[DEG];\r
+                                          for (I2 = 0; I2 < DEG; I2++) {\r
+                                                  ACOEFIN[I2] = ACOEF[LOD+I2-1];\r
+                                                  BCOEFIN[I2] = BCOEF[LOD+I2-1];\r
+                                          }\r
+                                          WSPEC[K-1]=RR*QUWTS[J1-1]*Math.pow((1.0+TT),BETA)*JACSUM(TT,\r
+                                            DEG,ACOEFIN,BCOEFIN,H0VAL[AJT-1],JACOF);\r
+                                          TT=TT*SS;\r
+                                          TSPEC[K-1]=MD+TT*HL;\r
+                                          CT[0]=TSPEC[K-1];\r
+                                          CT[1] = 0.0;\r
+                                          ZSPEC[K-1]=PARFUN(PT,CT);\r
+                                      } // for (J=1; J <= NQPTS; J++)\r
+                                  } // else I != 1\r
+                              } // for (I=1; I <= QINTS[0]; I++)\r
+                              if (SS < 0.0) {\r
+                                  LIM=NQUAD;\r
+                                  if ((LIM%2) == 0) {\r
+                                      LIM=LIM/2;\r
+                                  }\r
+                                  else {\r
+                                      LIM=(LIM-1)/2;\r
+                                  }\r
+                                  J1=0;\r
+                                  J2=NQUAD+1;\r
+                                  for (J=1; J <= LIM; J++) {\r
+                                      J1=J1+1;\r
+                                      J2=J2-1;\r
+                                      TT=WSPEC[J1-1];\r
+                                      WSPEC[J1-1]=WSPEC[J2-1];\r
+                                      WSPEC[J2-1]=TT;\r
+                                      TT=TSPEC[J1-1];\r
+                                     TSPEC[J1-1]=TSPEC[J2-1];\r
+                                     TSPEC[J2-1]=TT;\r
+                                     CT[0]=ZSPEC[J1-1][0];\r
+                                     CT[1]=ZSPEC[J1-1][1];\r
+                                     ZSPEC[J1-1][0]=ZSPEC[J2-1][0];\r
+                                     ZSPEC[J1-1][1]=ZSPEC[J2-1][1];\r
+                                     ZSPEC[J2-1][0]=CT[0];\r
+                                     ZSPEC[J2-1][1]=CT[1];\r
+                                  } // for (J=1; J <= LIM; J++)\r
+                              } // if (SS < 0.0)\r
+               \r
+                              // THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS\r
+                              // AND POINTS WSPEC AND ZSPEC.  NOW ESTIMATE THE INTEGRAL.\r
+               \r
+                              ARSUM=0.0;\r
+                              AISUM=0.0;\r
+                              if (IA == 1) {\r
+                                  FIRST=true;\r
+                              }\r
+                              else {\r
+                                  CURARG=STARG;\r
+                                  RT1=STRT1;\r
+                                  DIFF1[0]=STDF1[0];\r
+                                  DIFF1[1]=STDF1[1];\r
+                                  THET1=STTH1;\r
+                              }\r
+                              for (K=1; K <= NQUAD; K++) {\r
+                                  WT=WSPEC[K-1];\r
+                                  DIFF2[0]=ZZ[0]-ZSPEC[K-1][0];\r
+                                  DIFF2[1]=ZZ[1]-ZSPEC[K-1][1];\r
+                                  RT2=TSPEC[K-1];\r
+                                  DIST=zabs(DIFF2[0],DIFF2[1]);\r
+                                  ARSUM=ARSUM+WT*Math.log(DIST);\r
+                                  if (FIRST) {\r
+                                      CURARG=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                                      THET2=CURARG;\r
+                                      FIRST=false;\r
+                                  } // if (FIRST)\r
+                                  else { // !FIRST\r
+                                      // CT=DIFF2/DIFF1\r
+                                      // CT=DIFF2*CONJG(DIFF1)\r
+                                      // ANGLE=ATAN2(AIMAG(CT),REAL(CT))\r
+                                      THET2=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                                      ANGLE=THET2-THET1;\r
+                                      if (ANGLE <= -PI || ANGLE > PI) {\r
+                                          if (ANGLE > PI) {\r
+                                                  ANGLE = ANGLE - TUPI;\r
+                                          }\r
+                                          else {\r
+                                                  ANGLE = ANGLE + TUPI;\r
+                                          }\r
+                                      } // if (ANGLE <= -PI || ANGLE > PI)\r
+                                      if (Math.abs(ANGLE) >= LIMIT) {\r
+                                          DIFF1IN[0] = -DIFF1[0];\r
+                                          DIFF1IN[1] = -DIFF1[1];\r
+                                          DIFF2IN[0] = -DIFF2[0];\r
+                                          DIFF2IN[1] = -DIFF2[1];\r
+                                          ANGLE=ARGIN1(RT1,RT2,PT,DIFF1IN,DIFF2IN,ZZ,LIMIT);\r
+                                      } // if (Math.abs(ANGLE) >= LIMIT)\r
+                                      CURARG=CURARG+ANGLE;\r
+                                  } // else !FIRST\r
+                                  AISUM=CURARG*WT+AISUM;\r
+                                  RT1=RT2;\r
+                                  DIFF1[0]=DIFF2[0];\r
+                                  DIFF1[1]=DIFF2[1];\r
+                                  THET1=THET2;\r
+                              } // for (K=1; K <= NQUAD; K++)\r
+                              break;\r
+                          } // else !(Math.abs(XI[1]) < PTHTL && Math.abs(XI[0]) < 1.0+PTHTL)             \r
+                      } // else DIST < TRRAD[K-1]\r
+               \r
+                      // END OF QUADRATURE SUM LOOP \r
+               \r
+                  } // for (JQ=1; JQ <= NQ; JQ++)\r
+               \r
+                  RSUM=RSUM+ARSUM;\r
+                  ISUM=ISUM+AISUM;\r
+                  if (JT < 0) {\r
+               \r
+                      // BRING THE ARGUMENT FORWARD TO THE CORNER POINT AND REPLACE\r
+                      // THE INCREMENTED CURARG VALUE BY AN INVERSE TANGENT\r
+                      // EVALUATION\r
+                      PIN[0] = 1.0;\r
+                      PIN[1] = 0.0;\r
+                      POUT = PARFUN(PT,PIN);\r
+                      DIFF2[0] = ZZ[0] - POUT[0];\r
+                      DIFF2[1] = ZZ[1] - POUT[1];\r
+                      RT2=1.0;\r
+                      THET2=Math.atan2(DIFF2[1],DIFF2[0]);\r
+                      ANGLE=THET2-THET1;\r
+                      if (ANGLE <= -PI || ANGLE > PI) {\r
+                          if (ANGLE > PI) {\r
+                                  ANGLE = ANGLE - TUPI;\r
+                          }\r
+                          else {\r
+                                  ANGLE = ANGLE + TUPI;\r
+                          }\r
+                      } // if (ANGLE <= -PI || ANGLE > PI)\r
+                      if (Math.abs(ANGLE) >= LIMIT) {\r
+                          DIFF1IN[0] = -DIFF1[0];\r
+                  DIFF1IN[1] = -DIFF1[1];\r
+                  DIFF2IN[0] = -DIFF2[0];\r
+                  DIFF2IN[1] = -DIFF2[1];\r
+                          ANGLE=ARGIN1(RT1,RT2,PT,DIFF1IN,DIFF2IN,ZZ,LIMIT);\r
+                      } // if (Math.abs(ANGLE) >= LIMIT)\r
+                      CURARG=CURARG+ANGLE;\r
+                      ARGBR=(int)Math.round((CURARG-THET2)/TUPI);\r
+                      CURARG=THET2+TUPI*ARGBR;\r
+                      RT1=-1.0;\r
+                      DIFF1[0]=DIFF2[0];\r
+                      DIFF1[1]=DIFF2[1];\r
+                      THET1=THET2;\r
+                  } // if (JT < 0)          \r
+                  STARG=CURARG;\r
+                  STRT1=RT1;\r
+                  STDF1[0]=DIFF1[0];\r
+                  STDF1[1]=DIFF1[1];\r
+                  STTH1=THET1;\r
+               \r
+                  // END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA\r
+               \r
+          } // for (IA=1; IA <= TNSUA; IA++)\r
+          CT[0] = RSUM;\r
+          CT[1] = ISUM;\r
+          expCT = Math.exp(CT[0]);\r
+          CT[0] = expCT * Math.cos(CT[1]);\r
+          CT[1] = expCT * Math.sin(CT[1]);\r
+          zdiv(CT[0],CT[1],FACTR[0],FACTR[1],cr,ci);\r
+          result[0] = cr[0];\r
+          result[1] = ci[0];\r
+               \r
+          IER[0]=0;\r
+          return result;\r
+   } // private double[] CINRAD\r
+\r
+   private void BCFSNG(int TNSUC, int DGPOC[], int JTYPC[], int LSUBC[],\r
+       double H0VLC[], double JAINC[], double BFSNC[][], double SOLNC[][]) {\r
+       // INTEGER TNSUC\r
+       //  INTEGER DGPOC(*),JTYPC(*),LSUBC(*)\r
+       // REAL H0VLC(*),JAINC(*)\r
+       // COMPLEX BFSNC(*),SOLNC(*)\r
+\r
+       // PERFORMS VARIOUS PRELIMINARY TASKS TO PREPARE FOR THE POST-\r
+       // PROCESSING QUADRATURE CALCULATIONS.\r
+\r
+       // SETS UP THE ARRAY OF COEFFICIENTS BFSNC NEEDED FOR THE CALCULATION\r
+       // OF THE COMPLEX BOUNDARY CORRESPONDENCE FUNCTIONS FOR THE MAP\r
+       // CANONICAL --> PHYSICAL; THESE ARE SIMPLY RELATED TO THE SOLUTION \r
+       // COEFFICIENT ARRAY SOLNC.\r
+\r
+       // LOCAL VARIABLES\r
+\r
+       int AJTC,DEG,I,J,J1,JTC,LO;\r
+       double B1,RTH0,TUPI;\r
+\r
+       TUPI=2.0*Math.PI;\r
+\r
+       for (I=1; I <= TNSUC; I++) {\r
+           JTC=JTYPC[I-1];\r
+           AJTC=Math.abs(JTC);\r
+           RTH0=Math.sqrt(H0VLC[AJTC-1]);\r
+           B1=JAINC[AJTC-1]+1.0;\r
+           DEG=DGPOC[I-1];\r
+           LO=LSUBC[I-1];\r
+\r
+           BFSNC[LO-1][0]=TUPI*SOLNC[LO-1][0]/(B1*RTH0);\r
+           BFSNC[LO-1][1]=TUPI*SOLNC[LO-1][1]/(B1*RTH0);\r
+           for (J=1; J <= DEG; J++) {\r
+               J1=LO+J;\r
+               BFSNC[J1-1][0]=TUPI*SOLNC[J1-1][0]/Math.sqrt(J*(J+B1));\r
+               BFSNC[J1-1][1]=TUPI*SOLNC[J1-1][1]/Math.sqrt(J*(J+B1));\r
+           } // for (J=1; J <= DEG; J++)\r
+\r
+           if (JTC < 0) {\r
+               for (J=1; J <= DEG; J += 2) {\r
+                   J1=LO+J;\r
+                   BFSNC[J1-1][0]=-BFSNC[J1-1][0];\r
+                   BFSNC[J1-1][1]=-BFSNC[J1-1][1];\r
+               } // for (J=1; J <= DEG; J += 2)\r
+           } // if (JTC < 0)\r
+\r
+       } // for (I=1; I <= TNSUC; I++)\r
+\r
+    } // private void BCFSNG\r
+   \r
+   private void OUPTCL(int DGPOC[], int JTYPC[], double LGTOL,int LSUBC[], int NQPTS,\r
+       int PARNT[], int PRNSA[], double RIGLL[], double SOLNC[][], int TNSUC,\r
+       boolean INTER, double INNRAD, int IER) {\r
+          // INTEGER NQPTS,OC,TNSUC,IER\r
+          // INTEGER DGPOC(*),JTYPC(*),LSUBC(*),PARNT(*),PRNSA(*)\r
+          // REAL LGTOL,INNRAD\r
+          // REAL RIGLL(*)\r
+          // COMPLEX SOLNC(*)\r
+          // LOGICAL INTER\r
+               \r
+          // LOCAL VARIABLES\r
+               \r
+          int AJT,DG,I,IC,L,LOM,LOD,PSA,PT;\r
+          double MOD;\r
+          double COF[] = new double[2];\r
+          // COMPLEX COF\r
+          // CHARACTER JBNM*4,OFL*6\r
+          // EXTERNAL WRHEAD,WRTAIL\r
+               \r
+          \r
+          WRHEAD(3,0,null);\r
+               \r
+          if (INTER) {\r
+                  Preferences.debug("INNER RADIUS = " + INNRAD + "\n", Preferences.DEBUG_ALGORITHM);\r
+          }\r
+                       \r
+       /*            WRITE(OC,'(/,A,/)') 'JACOBI COEFFICIENTS FOR INVERSE DENSITY FUNCT\r
+                    +IONS'\r
+                     DO 20 IC=1,TNSUC\r
+                       DG=DGPOC(IC)\r
+                       PSA=PRNSA(IC)\r
+                       PT=PARNT(PSA)\r
+                       WRITE(OC,*)\r
+                       WRITE(OC,101) IC,PSA,PT\r
+                       LOM=LSUBC(IC)\r
+                       AJT=ABS(JTYPC(IC))\r
+                       LOD=(AJT-1)*NQPTS+1\r
+                       DO 10 I=0,DG\r
+                         COF=SOLNC(LOM+I)\r
+                         MOD=ABS(COF)\r
+                         WRITE(OC,102) I,COF,MOD,LGTOL/RIGLL(LOD+I)\r
+               10      CONTINUE\r
+               20    CONTINUE\r
+               C\r
+               101   FORMAT('SUB ARC = ',I3,'; PHYSICAL PARENTAL SUB ARC = ',I3,' ON GL\r
+                    +OBAL ARC ',I2/,' N',T6,'REAL PART',T24,'IMAGINARY PART',T42,'MODUL\r
+                    +US',T56,'IGNORE LVL')\r
+               102   FORMAT(I2,T4,E16.8,T22,E16.8,T40,E11.3,T54,E11.3)\r
+               C\r
+                     CALL WRTAIL(3,OC,IER)\r
+                     CLOSE(OC)\r
+               C */\r
+   } // private void OUPTCL\r
 \r
 \r
 \r