SUBROUTINE ANGLE7(BE,NA,IN) INTEGER NA REAL BE(*) LOGICAL IN C C**** BE=JACIN,NA=NARCS,IN=INTER C C**** TO COMPUTE THE ARRAY OF JACOBI INDECES CORRESPONDING TO THE C**** CORNER ANGLES ON THE BOUNDARY C C**** LOCAL VARIABLES INTEGER K,B0,B1,B2 REAL X,Y,ANG,PI,R1MACH,EPS,XI,APP COMPLEX U,V,DPARFN EXTERNAL DPARFN,R1MACH C PI=4E+0*ATAN(1E+0) EPS=SQRT(R1MACH(4)) DO 1 K=1,NA U=DPARFN(K,(1.0,0.0)) U=-U/ABS(U) IF (K.EQ.NA) THEN V=DPARFN(1,(-1.0,0.0)) ELSE V=DPARFN(K+1,(-1.0,0.0)) ENDIF V=V/ABS(V) V=U/V X=REAL(V) Y=AIMAG(V) ANG=ATAN2(Y,X) IF (ANG.LT.0) THEN ANG=ANG+2E+0*PI ENDIF ANG=ANG/PI IF (.NOT.IN) THEN ANG=2E+0-ANG ENDIF ANG=-1E+0+1E+0/ANG C C**** TRY TO DETECT SIMPLE RATIONAL INDECES AND FORCE BEST REAL C**** APPROXIMATIONS C IF (ABS(ANG) .LT. EPS) THEN ANG=0E+0 ELSE XI=ABS(ANG) B0=INT(XI) XI=XI-REAL(B0) IF (ABS(XI) .LT. EPS) THEN APP=REAL(B0) ELSE XI=1E+0/XI B1=INT(XI) XI=XI-REAL(B1) IF (ABS(XI) .LT. EPS) THEN APP=REAL(1+B0*B1)/REAL(B1) ELSE XI=1E+0/XI B2=INT(XI) APP=REAL(B0*(1+B1*B2)+B2)/REAL(1+B1*B2) ENDIF ENDIF APP=SIGN(1E+0,ANG)*APP IF (ABS(ANG-APP) .LT. EPS) ANG=APP ENDIF C IF (K.EQ.NA) THEN BE(1)=ANG ELSE BE(K+1)=ANG ENDIF 1 CONTINUE END REAL FUNCTION ARGIN1(RT1,RT2,PT,DIFF1,DIFF2,ZZ,LIMIT) INTEGER PT REAL RT1,RT2,LIMIT COMPLEX DIFF1,DIFF2,ZZ C C ZZ IS A GIVEN FIELD POINT AND DIFF1, DIFF2 ARE THE DIFFERENCES C BETWEEN ZZ AND CONSECUTIVE POINTS ON THE BOUNDARY C (DIFF1=PARFUN(PT,RT1)-ZZ, ZET2=PARFUN(PT,RT2)-ZZ). THE C PURPOSE OF THIS ROUTINE IS TO CALCULATE THE INCREASE IN THE C ARGUMENT ARG(ZZ-Z) AS Z MOVES ALONG THE BOUNDARY FROM THE POINT C WITH PARAMETER VALUE RT1 TO THE POINT WITH PARAMETER VALUE RT2. C C LOCAL VARIABLES C INTEGER NANGS,NINTS REAL ANGLE,T1,T2 COMPLEX D1,D2,PARFUN,V EXTERNAL PARFUN C C LIMIT IS CURRENTLY SET TO 3*PI/4, APPROXIMATELY C T1=RT1 T2=(RT1+RT2)*5E-1 D1=DIFF1 D2=PARFUN(PT,CMPLX(T2))-ZZ NANGS=0 NINTS=2 ARGIN1=0E+0 C 10 CONTINUE V=D2*CONJG(D1) ANGLE=ATAN2(AIMAG(V),REAL(V)) IF (ABS(ANGLE) .GE. LIMIT) THEN T2=(T1+T2)*5E-1 D2=PARFUN(PT,CMPLX(T2))-ZZ NINTS=NINTS+1 GOTO 10 ELSE ARGIN1=ARGIN1+ANGLE NANGS=NANGS+1 IF (NANGS .NE. NINTS) THEN T1=T2 T2=RT2 D1=D2 D2=DIFF2 GOTO 10 ENDIF ENDIF C END SUBROUTINE ASONJ7(ALFA,BETA,A,B,H,N) INTEGER N REAL A(*),B(*),H,ALFA,BETA C ..TO ASSIGN THE COEFFICIENTS A(K) AND B(K) , K=1(1)N, IN THE C ..3-TERM RECURRENCE FORMULA FOR THE ORTHONORMAL JACOBI POLYNOMIALS C ..WHERE C .. C .. A(K)P (X) = (X - B(K))P (X) - A(K-1)P (X) , K=1,2,..,N, C .. K K-1 K-2 C .. C .. P (X) = 0 , P (X) = 1/SQRT(H) C .. -1 0 C .. C ..AND H IS THE ZEROTH MOMENT OF THE JACOBI WEIGHT FUNCTION C ..(1-X)**ALFA*(1+X)**BETA ON [-1,1]. C**** AUTHOR: DAVID HOUGH C**** LAST UPDATE: 15.09.89 C**** ..LOCAL VARIABLES.. REAL SUM,DIFF,PROD,TC,T,SC,S,GAMMA,C INTEGER K EXTERNAL GAMMA SUM=ALFA+BETA DIFF=BETA-ALFA PROD=SUM*DIFF C ..CALCULATE H. TC=SUM+1.0 SC=2.0**TC S=GAMMA(ALFA+1.0) T=GAMMA(BETA+1.0) C=GAMMA(TC+1.0) H=SC*S*T/C C ..START ON A,B ARRAYS. IF (N.GT.0) THEN T=2.0+SUM S=T*T C=4.0*(ALFA+1.0)*(BETA+1.0)/S/(T+1.0) A(1)=SQRT(C) B(1)=DIFF/T DO 10 K=2,N B(K)=PROD/T/(T+2.0) T=2.0*K+SUM S=T*T C=4.0*K*(ALFA+K)*(BETA+K)*(SUM+K)/S/(S-1.0) A(K)=SQRT(C) 10 CONTINUE ENDIF END SUBROUTINE ASQUC7(AQCOF,BQCOF,CQCOF,JACIN,NJIND,NQPTS) INTEGER NJIND,NQPTS REAL AQCOF(*),BQCOF(*),CQCOF(*),JACIN(*) C C ..TO ASSIGN THE COEFFICIENTS A(J), B(J) AND C(J) ,J=1,MN IN THE C ..3-TERM RECURRENCE FORMULA C .. C .. Q (Z) = (A(J)Z - B(J))Q (Z) - C(J)Q (Z) , J=2,...,MN C .. J+1 J J-1 C .. C .. Q (Z) = (A(1)Z - B(1))Q (Z) - C(1) C .. 2 1 C .. C ..WHERE Q (Z):=

, P IS THE ORTHONORMAL JACOBI POLY C .. J J J C ..OF DEGREE J AND THE INNER PRODUCT IS WITH RESPECT TO THE JACOBI C ..DISTRIBUTION OVER (-1,1). HERE A(J,I) STORES "A(J)" FOR THE ITH C ..ARC ON THE BOUNDARY (WITH SIMILAR ROLE FOR ARRAYS B AND C) AND C ..THE JACOBI WEIGHT FUNCTION FOR THE ITH ARC IS C ..(1-X)**AB(I,1)*(1+X)**AB(I,2), J=1,MN, I=1,NA. C C**** AUTHOR: DAVID HOUGH C**** LAST UPDATE: 15.09.89 C C**** LOCAL VARIABLES INTEGER I,J,K,LOSUB REAL BE,H,D,N,N1,N2,F EXTERNAL ASONJ7 DO 10 I=1,NJIND LOSUB=(I-1)*NQPTS+1 BE=JACIN(I) CALL ASONJ7(1E+0,BE+1E+0,AQCOF(LOSUB),BQCOF(LOSUB),H,NQPTS) DO 20 K=1,NQPTS J=LOSUB+K-1 N=REAL(K) D=(N+1E+0)*(N+BE+2E+0) N1=N*(N+BE+1E+0) N2=(N-1E+0)*(N+BE) F=SQRT(N1/D) AQCOF(J)=F/AQCOF(J) BQCOF(J)=BQCOF(J)*AQCOF(J) IF (K.GT.1) THEN CQCOF(J)=AQCOF(J)*N2/AQCOF(J-1)/N1 ELSE CQCOF(J)=-AQCOF(J)*SQRT(H/N1) ENDIF 20 CONTINUE 10 CONTINUE C END SUBROUTINE AXION1(AXION,NEWDG,SOLUN,MDGPO,TNSUA,DGPOL,LOSUB,HISUB, +RIGLL,LGTOL,ACCPT,JACIN,JATYP,NJIND,NEWHL,ESTOL,IER) INTEGER AXION(*),NEWDG(*),MDGPO,TNSUA,DGPOL(*),LOSUB(*),HISUB(*), +JATYP(*),NJIND,IER REAL SOLUN(*),RIGLL(*),LGTOL,JACIN(*),NEWHL(*),ESTOL LOGICAL ACCPT C C TO DETERMINE THE ARRAY "AXION" WHICH SPECIFIES THE ACTIONS THAT C ARE TO TAKEN ON EACH SUBARC, AS FOLLOWS: C AXION(I)=-1 - REDUCE THE DEGREE ON ARC I C AXION(I)=0 - LEAVE THE DEGREE ON ARC I UNCHANGED C AXION(I)=1 - INCREASE THE DEGREE ON ARC I C AXION(I)=2 - SUBDIVIDE ARC I. C C IN CASE ABS(AXION(I))=1, NEWDG(I) RECORDS THE NEW DEGREE TO BE C USED ON ARC I. C C IN CASE AXION(I)<=0 FOR ALL I=1,TNSUA, THE SOLUTION IS DEEMED TO C BE ACCEPTABLE TO THE REQUIRED ACCURACY AND WE SET ACCPT=.TRUE.; C ACCPT=.FALSE. OTHERWISE. C C ALSO TO DETERMINE THE EFFECTIVE STOPPING TOLERANCE ESTOL; IF THE C USER HAD INPUT ESTOL AS THE VALUE FOR THE ARGUMENT MAXER IN C JAPHYC, THEN THE CURRENT SOLUTION WOULD BE ACCEPTED. C C IER=0 - NORMAL EXIT C IER=54 - LOCAL PARAMETER MXCO MUST BE INCREASED C C LOCAL VARIABLES C INTEGER AJT,CI,CRITCO,DG,I,IA,IFL,J,LOD,LOM,MINDG,MXCO,PNDG REAL AA,BETA,CONF,EXA,LAM,POW,SAFEF,TERM,THSLN,XX,VAR REAL COVAR(2,2) LOGICAL CONSV EXTERNAL CRITCO,STATS1 PARAMETER (SAFEF=5E+0,MINDG=0,MXCO=32,CONSV=.FALSE.) REAL POSCO(MXCO) C ESTOL=0E+0 DO 60 IA=1,TNSUA DG=DGPOL(IA) IF (DG+1 .GT. MXCO) THEN IER=54 RETURN ENDIF LOM=LOSUB(IA) AJT=ABS(JATYP(IA)) LOD=(AJT-1)*(MDGPO+1)+1 DO 10 I=0,DG J=LOM+I POSCO(I+1)=ABS(SOLUN(J))*RIGLL(LOD+I)/LGTOL 10 CONTINUE CI=CRITCO(DG+1,POSCO)-1 PNDG=CI+1+MINDG IF (CI .EQ. -1) THEN C C**** ALL COEFFICIENTS ARE ACCEPTABLE. C IF (DG .EQ. MINDG) THEN AXION(IA)=0 NEWDG(IA)=DG ELSE C C**** PROBABLY DECREASE THE DEGREE, BUT ONLY IGNORE THOSE C**** COEFFICIENTS WHICH ARE 'SAFELY' BELOW THE TOLERANCE C**** LIMIT. C DO 20 I=0,DG POSCO(I+1)=POSCO(I+1)*SAFEF 20 CONTINUE PNDG=CRITCO(DG+1,POSCO)+MINDG IF (PNDG .GE. DG) THEN AXION(IA)=0 NEWDG(IA)=DG ELSE IF (PNDG .LE. MINDG) THEN AXION(IA)=-1 NEWDG(IA)=MINDG ELSE AXION(IA)=-1 NEWDG(IA)=PNDG ENDIF ENDIF ELSE IF (PNDG .EQ. DG) THEN AXION(IA)=0 NEWDG(IA)=DG ELSE IF (PNDG .LT. DG) THEN C C**** PROBABLY DECREASE THE DEGREE, BUT ONLY IGNORE THOSE C**** COEFFICIENTS WHICH ARE 'SAFELY' BELOW THE TOLERANCE C**** LIMIT. DO 30 I=0,DG POSCO(I+1)=POSCO(I+1)*SAFEF 30 CONTINUE PNDG=CRITCO(DG+1,POSCO)+MINDG IF (PNDG .GE. DG) THEN AXION(IA)=0 NEWDG(IA)=DG ELSE IF (PNDG .LE. MINDG) THEN AXION(IA)=-1 NEWDG(IA)=MINDG ELSE AXION(IA)=-1 NEWDG(IA)=PNDG ENDIF ELSE IF (DG .EQ. MDGPO) THEN C C**** ARC SUBDIVISION IS REQUIRED AND ASSIGN NEW HALF-LENGTH. C AXION(IA)=2 BETA=JACIN(AJT) LAM=MIN(1E+0,1E+0+BETA) IF (AJT .NE. NJIND) THEN NEWHL(IA)=(1E+0/POSCO(CI+1))**(1E+0/(1E+0+LAM+BETA)) ELSE NEWHL(IA)=5E-1 ENDIF ELSE C C**** MUST DECIDE BETWEEN ARC SUBDIVISION ARC DEGREE INCREASE. C**** FIRST MAKE CONSERVATIVE ESTIMATES FOR THE COEFFICIENTS C**** FOR DEGREES DG+1 TO MDGPO. C CALL STATS1(SOLUN(LOM),DG+1,AA,POW,EXA,COVAR,CONF,IFL) IF (IFL .EQ. 1) THEN C C**** NOT ENOUGH DATA TO MAKE ESTIMATES, WHICH PRESUMES DG C**** IS RATHER SMALL; THEREFORE INCREASE THE DEGREE. C AXION(IA)=1 NEWDG(IA)=MIN(DG+2,MDGPO) ELSE DO 40 I=DG+1,MDGPO IF (CONSV) THEN XX=LOG(1E+0+I) VAR=COVAR(1,1)+XX*XX*COVAR(2,2)+2E+0*XX*COVAR(1,2) THSLN=EXA*(I+1)**POW*EXP(CONF*SQRT(VAR)) ELSE THSLN=EXA*(I+1)**POW ENDIF POSCO(I+1)=ABS(THSLN)*RIGLL(LOD+I)/LGTOL 40 CONTINUE PNDG=CRITCO(MDGPO+1,POSCO)+MINDG IF (PNDG .LE. MDGPO) THEN C C**** INCREASE DEGREE C AXION(IA)=1 NEWDG(IA)=PNDG ELSE C C**** SUBDIVIDE ARC AND ASSIGN NEW HALF-LENGTH. C AXION(IA)=2 BETA=JACIN(AJT) LAM=MIN(1E+0,1E+0+BETA) IF (AJT .NE. NJIND) THEN NEWHL(IA)=(1E+0/POSCO(CI+1))** + (1E+0/(1E+0+LAM+BETA)) ELSE NEWHL(IA)=5E-1 ENDIF ENDIF ENDIF ENDIF C C**** NOW UPDATE THE EFFECTIVE STOPPING TOLERANCE C J=HISUB(IA)-MINDG DO 50 I=J,HISUB(IA) TERM=ABS(SOLUN(I))*RIGLL(LOD+I-LOM) ESTOL=MAX(ESTOL,EXP(TERM)-1E+0) 50 CONTINUE 60 CONTINUE C ACCPT=.TRUE. DO 70 I=1,TNSUA IF (AXION(I) .GT. 0) THEN ACCPT=.FALSE. GOTO 80 ENDIF 70 CONTINUE C 80 CONTINUE C IER=0 C END SUBROUTINE BCFSNG(TNSUC,DGPOC,JTYPC,LSUBC,H0VLC,JAINC,BFSNC,SOLNC) INTEGER TNSUC INTEGER DGPOC(*),JTYPC(*),LSUBC(*) REAL H0VLC(*),JAINC(*) COMPLEX BFSNC(*),SOLNC(*) C C PERFORMS VARIOUS PRELIMINARY TASKS TO PREPARE FOR THE POST- C PROCESSING QUADRATURE CALCULATIONS. C C SETS UP THE ARRAY OF COEFFICIENTS BFSNC NEEDED FOR THE CALCULATION C OF THE COMPLEX BOUNDARY CORRESPONDENCE FUNCTIONS FOR THE MAP C CANONICAL --> PHYSICAL; THESE ARE SIMPLY RELATED TO THE SOLUTION C COEFFICIENT ARRAY SOLNC. C C LOCAL VARIABLES C INTEGER AJTC,DEG,I,J,J1,JTC,LO REAL B1,RTH0,TUPI C TUPI=8E+0*ATAN(1E+0) C DO 30 I=1,TNSUC JTC=JTYPC(I) AJTC=ABS(JTC) RTH0=SQRT(H0VLC(AJTC)) B1=JAINC(AJTC)+1E+0 DEG=DGPOC(I) LO=LSUBC(I) C BFSNC(LO)=TUPI*SOLNC(LO)/(B1*RTH0) DO 10 J=1,DEG J1=LO+J BFSNC(J1)=TUPI*SOLNC(J1)/SQRT(J*(J+B1)) 10 CONTINUE C IF (JTC .LT. 0) THEN DO 20 J=1,DEG,2 J1=LO+J BFSNC(J1)=-BFSNC(J1) 20 CONTINUE ENDIF C 30 CONTINUE C END SUBROUTINE BCFVTF(BCFSN,VTARG,DGPOL,JATYP,LOSUB,PARNT,RFARC,TNSUA, +H0VAL,JACIN,RFARG,SOLUN) INTEGER RFARC,TNSUA INTEGER DGPOL(*),JATYP(*),LOSUB(*),PARNT(*) REAL RFARG REAL BCFSN(*),H0VAL(*),JACIN(*),SOLUN(*),VTARG(*) C C PERFORMS VARIOUS PRELIMINARY TASKS TO PREPARE FOR THE POST- C PROCESSING QUADRATURE CALCULATIONS. C C SETS UP THE ARRAY OF COEFFICIENTS BCFSN NEEDED FOR THE C CALCULATION OF THE BOUNDARY CORRESPONDENCE FUNCTIONS FOR THE MAP C PHYSICAL --> CANONICAL; THESE ARE SIMPLY RELATED TO THE SOLUTION C COEFFICIENT ARRAY SOLUN. C C CALCULATES THE ARGUMENTS VTARG OF THE IMAGES ON THE UNIT CIRCLE C OF THE END POINTS OF ALL SUBARCS, ENFORCING THE USER'S ROTATION C CONDITION THAT PARFUN(RFARC,(-1.0,0.0)) BE MAPPED TO THE POINT C WITH ARGUMENT RFARG. C C LOCAL VARIABLES C INTEGER AJT,DEG,I,J,J1,JT,LO REAL B1,RTH0,THET0,TUPI C TUPI=8E+0*ATAN(1E+0) VTARG(1)=0E+0 C DO 30 I=1,TNSUA JT=JATYP(I) AJT=ABS(JT) RTH0=SQRT(H0VAL(AJT)) B1=JACIN(AJT)+1E+0 DEG=DGPOL(I) LO=LOSUB(I) VTARG(I+1)=VTARG(I)+SOLUN(LO)*TUPI*RTH0 C BCFSN(LO)=TUPI*SOLUN(LO)/(B1*RTH0) DO 10 J=1,DEG J1=LO+J BCFSN(J1)=TUPI*SOLUN(J1)/SQRT(J*(J+B1)) 10 CONTINUE C IF (JT .LT. 0) THEN DO 20 J=1,DEG,2 J1=LO+J BCFSN(J1)=-BCFSN(J1) 20 CONTINUE ENDIF C 30 CONTINUE C I=1 40 IF (PARNT(I) .EQ. RFARC) THEN THET0=RFARG-VTARG(I) ELSE I=I+1 GOTO 40 ENDIF C DO 50 I=1,TNSUA VTARG(I)=VTARG(I)+THET0 50 CONTINUE VTARG(TNSUA+1)=VTARG(1)+TUPI C END SUBROUTINE BISNEW(IER,LL,TT,UU,A1COF,ACOEF,B1COF,BCFSN,BCOEF,BETA, +DEG,H0VAL,H1VAL,JACOF,NEWTL,SJT,RRHS,TOLIW) INTEGER DEG,IER REAL A1COF(*),ACOEF(*),B1COF(*),BCFSN(*),BCOEF(*),BETA,H0VAL, +H1VAL,JACOF,LL,NEWTL,TT,UU,SJT,RRHS,TOLIW C C A MIXTURE OF BISECTION AND NEWTON'S METHOD TO SOLVE THE NON-LINEAR C BOUNDARY CORRESPONDENCE EQUATION C C THETA(TT) = CONST C C FOR REAL PARAMETER TT GIVEN REAL CONST; SEE RB#50 P134. THE C INTERVAL (LL,UU) SHOULD BRACKET TT. C C IER=0 - NORMAL EXIT C IER=34 - FUNCTION HAS SAME SIGN AT LL AND UU C IER=35 - ZERO FUNCTION DERIVATIVE DETECTED C C LOCAL VARIABLES C INTEGER MNITS,NBSCT,NITS,STEPS REAL DFT,EPS,FL,FT,FU,JACSUM,RDIFF,TUPI PARAMETER (MNITS=10,NBSCT=3) EXTERNAL JACSUM C TUPI=8E+0*ATAN(1E+0) C FL=JACSUM(SJT*LL,DEG-1,A1COF,B1COF,H1VAL,BCFSN(2)) FL=BCFSN(1)-(1E+0-SJT*LL)*FL FL=(1E+0+SJT*LL)**(1E+0+BETA)*FL-RRHS C FU=JACSUM(SJT*UU,DEG-1,A1COF,B1COF,H1VAL,BCFSN(2)) FU=BCFSN(1)-(1E+0-SJT*UU)*FU FU=(1E+0+SJT*UU)**(1E+0+BETA)*FU-RRHS C IF (FL*FU .GT. 0E+0) THEN IER=34 RETURN ENDIF C C ENTER NEWTON ITERATION MODE C 10 CONTINUE TT=(UU+LL)*5E-1 NITS=0 20 CONTINUE FT=JACSUM(SJT*TT,DEG-1,A1COF,B1COF,H1VAL,BCFSN(2)) FT=BCFSN(1)-(1E+0-SJT*TT)*FT FT=(1E+0+SJT*TT)**(1E+0+BETA)*FT-RRHS DFT=JACSUM(SJT*TT,DEG,ACOEF,BCOEF,H0VAL,JACOF) DFT=TUPI*(1E+0+SJT*TT)**BETA*DFT*SJT IF (DFT.EQ.0E+0) THEN IER=35 RETURN ENDIF RDIFF=FT/DFT TT=TT-RDIFF NITS=NITS+1 IF (ABS(RDIFF) .LT. NEWTL) THEN C C NEWTON ITERATIONS HAVE CONVERGED FOR TT C IER=0 RETURN ELSE IF (TT.LE.LL .OR. TT.GE.UU .OR. NITS.EQ.MNITS) THEN C C PERFORM NBSCT BISECTION STEPS C STEPS=0 30 CONTINUE EPS=(UU-LL)*5E-1 IF (EPS.LT.TOLIW) THEN C C BISECTION ITERATIONS HAVE CONVERGED FOR TT C IER=0 RETURN ENDIF TT=(UU+LL)*5E-1 FT=JACSUM(SJT*TT,DEG-1,A1COF,B1COF,H1VAL,BCFSN(2)) FT=BCFSN(1)-(1E+0-SJT*TT)*FT FT=(1E+0+SJT*TT)**(1E+0+BETA)*FT-RRHS C IF (FT*FL.LT.0E+0) THEN UU=TT FU=FT ELSE LL=TT FL=FT ENDIF STEPS=STEPS+1 IF (STEPS .EQ. NBSCT) THEN C C RE-START NEWTON MODE C GOTO 10 ELSE C C CONTINUE WITH BISECTION C GOTO 30 ENDIF ELSE C C CONTINUE WITH NEWTON MODE C GOTO 20 ENDIF C END SUBROUTINE BMCANP(THETA,PHYPT,DERIV,WANTD,INTER,CENTR,IGEOM,RGEOM, +ISNCA,RSNCA,ZSNCA,WANTM,IER) C INTEGER IER INTEGER IGEOM(*),ISNCA(*) REAL THETA REAL RGEOM(*),RSNCA(*) COMPLEX PHYPT,DERIV,CENTR COMPLEX ZSNCA(*) LOGICAL WANTD,INTER,WANTM C C ...................................................................... C C 1. BMCANP C BOUNDARY MAPPING FOR THE CANONICAL --> PHYSICAL MAP. C C 2. PURPOSE C GIVEN A POINT ON THE UNIT CIRCLE, THIS ROUTINE COMPUTES THE C CORRESPONDING APPROXIMATE IMAGE POINT ON THE BOUNDARY OF THE C PHYSICAL DOMAIN AND, IF REQUESTED, ALSO COMPUTES THE C DERIVATIVE OF THE MAP : CANONICAL --> PHYSICAL AT THE GIVEN C POINT ON THE UNIT CIRCLE. C C 3. CALLING SEQUENCE C CALL BMCANP(THETA,PHYPT,DERIV,WANTD,INTER,CENTR,IGEOM,RGEOM, C ISNCA,RSNCA,ZSNCA,WANTM,IER) C C PARAMETERS C ON ENTRY C THETA - REAL C THE ARGUMENT OF THE GIVEN POINT ON THE UNIT C CIRCLE. C C WANTD - LOGICAL C IF WANTD IS TRUE THEN DERIV IS COMPUTED OTHERWISE C DERIV ISN'T COMPUTED. C C INTER - LOGICAL C TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE C OTHERWISE. (AS PREVIOUSLY USED IN JAPHYC, GQPHYC) C C IGEOM - INTEGER ARRAY C THE INTEGER VECTOR IGEOM PREVIOUSLY SET UP BY C JAPHYC. C C RGEOM - REAL ARRAY C THE REAL VECTOR RGEOM PREVIOUSLY SET UP BY JAPHYC. C C ISNCA - INTEGER ARRAY C THE INTEGER VECTOR PREVIOUSLY SET UP BY JACANP. C C RSNCA - REAL ARRAY C THE REAL VECTOR PREVIOUSLY SET UP BY JACANP. C C ZSNCA - COMPLEX ARRAY C THE COMPLEX VECTOR PREVIOUSLY SET UP BY JACANP. C C WANTM - LOGICAL C IF WANTM IS TRUE THEN, ON AN ABNORMAL EXIT, AN C ERROR MESSAGE IS WRITTEN ON THE STANDARD OUTPUT C CHANNEL. IF WANTM IS FALSE THEN NO MESSAGE IS C WRITTEN. C C ON EXIT C PHYPT - COMPLEX C THE COMPUTED POINT ON THE PHYSICAL BOUNDARY CORR- C ESPONDING TO THE POINT WITH ARGUMENT THETA ON THE C UNIT CIRCLE. C C DERIV - COMPLEX C THE COMPUTED DERIVATIVE OF THE MAP : CANONICAL --> C PHYSICAL AT THE GIVEN POINT ON THE UNIT CIRCLE. C THIS IS COMPUTED ONLY IF WANTD IS TRUE. C C IER - INTEGER C IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED; C A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY C WRITTEN ON THE STANDARD OUTPUT CHANNEL. C IER=0 - NORMAL EXIT. C IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD C BE SELF EXPLANATORY. C C C 4. SUBROUTINES OR FUNCTIONS NEEDED C - THE CONFPACK LIBRARY. C - THE REAL FUNCTION R1MACH. C - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN. C C C 5. FURTHER COMMENTS C - NOTE THAT THIS ROUTINE CAN ONLY BE USED A F T E R THE C ROUTINE JACANP HAS SUCCESSFULLY EXECUTED, AND THAT SOME C INPUT ARGUMENTS FOR BMCANP ARE OUTPUT VALUES FROM JACANP. C C C ...................................................................... C AUTHOR: DAVID HOUGH, ETH, ZUERICH C LAST UPDATE: 3 JULY 1990 C ...................................................................... C C LOCAL VARAIBLES C INTEGER ACOFC,AICOC,BCOFC,BFSNC,BICOC,DGPOC,H0VLC,HALEN,HIVLC, +JAINC,JTYPC,LSUBC,NARCS,NJIND,NQPTS,MIDPT,MNCOF,MNSUA,MNSUC,PARNT, +PHPAS,PRNSA,SOLNC,TNGQP,TNSUC,VARGC,VTARG CHARACTER IERTXT*6 C EXTERNAL BMCAP1,IERTXT C NARCS=ISNCA(1) NQPTS=ISNCA(2) TNSUC=ISNCA(3) MNSUC=ISNCA(5) MNCOF=ISNCA(6) MNSUA=IGEOM(4) C NJIND=NARCS+1 TNGQP=NQPTS*NJIND C C**** SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC C PARNT=5 HALEN=3 MIDPT=MNSUA+3 VTARG=2*MNSUA+3 C C**** SET UP POINTERS TO ISNCA, RSNCA AND ZSNCA, AS IN JACANP C DGPOC=7 JTYPC=MNSUC+7 LSUBC=2*MNSUC+7 PRNSA=3*MNSUC+7 ACOFC=2 BCOFC=TNGQP+2 AICOC=2*TNGQP+2 BICOC=3*TNGQP+2 H0VLC=6*TNGQP+2 HIVLC=NJIND+6*TNGQP+2 JAINC=2*NJIND+6*TNGQP+2 PHPAS=4*NJIND+6*TNGQP+2 VARGC=MNSUC+4*NJIND+6*TNGQP+2 BFSNC=2 SOLNC=MNCOF+2 C C**** GET REQUIRED MAP AND DERIVATIVE C CALL BMCAP1(PHYPT,DERIV,THETA,NQPTS,TNSUC,ISNCA(DGPOC), +ISNCA(JTYPC),ISNCA(LSUBC),IGEOM(PARNT),ISNCA(PRNSA),RSNCA(AICOC), +RSNCA(ACOFC),RSNCA(BICOC),RSNCA(BCOFC),RSNCA(H0VLC),RSNCA(HIVLC), +RGEOM(HALEN),RSNCA(JAINC),RGEOM(MIDPT),RSNCA(PHPAS),RGEOM(VTARG), +RSNCA(VARGC),ZSNCA(BFSNC),CENTR,ZSNCA(SOLNC),INTER,WANTD,IER) C C**** SEND ERROR MESSAGE TO STANDARD OUTPUT OF NECESSARY C IF (IER.GT.0 .AND. WANTM) WRITE(*,*) IERTXT(IER) C END SUBROUTINE BMCAP1(PHYPT,DERIV,THETA,NQPTS,TNSUC,DGPOC,JTYPC,LSUBC, +PARNT,PRNSA,A1COC,ACOFC,B1COC,BCOFC,H0VLC,H1VLC,HALEN,JAINC,MIDPT, +PHPAS,THET0,VARGC,BFSNC,CENTR,SOLNC,INTER,WANTD,IER) INTEGER IER,NQPTS,TNSUC INTEGER DGPOC(*),JTYPC(*),LSUBC(*),PARNT(*),PRNSA(*) REAL THET0,THETA REAL A1COC(*),ACOFC(*),B1COC(*),BCOFC(*),H0VLC(*),H1VLC(*), +HALEN(*),JAINC(*),MIDPT(*),PHPAS(*),VARGC(*) COMPLEX CENTR,DERIV,PHYPT COMPLEX BFSNC(*),SOLNC(*) LOGICAL INTER,WANTD C C GIVEN A POINT ON THE UNIT CIRCLE DEFINED BY ITS ARGUMENT THETA C TO COMPUTE ITS IMAGE POINT PHYPT ON THE BOUNDARY OF THE C PHYSICAL DOMAIN. IN CASE WANTD=.TRUE. THEN ALSO COMPUTE THE C DERIVATIVE DERIV OF THE MAP ONTO THE PHYSICAL DOMAIN AT THE GIVEN C BOUNDARY POINT ON THE UNIT DISC. C C IER=0 - NORMAL TERMINATION C IER=46 - LOCAL PARAMETER MNCOF NEEDS INCREASING C C LOCAL VARIABLES C INTEGER AJTC,DGC,I,I1,IA,JTC,LODC,LOMC,MNCOF,PSA REAL AA,ARGW,BB,BETAC,HA,MD,R1MACH,SJTC,TT,TUPI COMPLEX CT0,FP,JACSUC,PARFUN,PHI PARAMETER (MNCOF=32) COMPLEX JACOF(MNCOF) EXTERNAL JACSUC,PARFUN,R1MACH C TUPI=8E+0*ATAN(1E+0) ARGW=THETA C 10 CONTINUE IF (ARGW .GT. THET0+TUPI) THEN ARGW=ARGW-TUPI GOTO 10 ENDIF 20 CONTINUE IF (ARGW .LT. THET0) THEN ARGW=ARGW+TUPI GOTO 20 ENDIF C IA=1 30 CONTINUE IF (VARGC(IA).LE.ARGW .AND. ARGW.LE.VARGC(IA+1)) THEN C C WW IS ON ARC IA AND WE USE THE CONTINUATION OF THE BOUNDARY C CORRESPONDENCE FUNCTION TO ESTIMATE PHYPT AND (IF WANTED) THE C DERIVATIVE. C HA=(VARGC(IA+1)-VARGC(IA))*5E-1 MD=(VARGC(IA+1)+VARGC(IA))*5E-1 TT=(ARGW-MD)/HA IF (TT .GT. 1E+0) THEN TT=1E+0 ELSE IF (TT .LT. -1E+0) THEN TT=-1E+0 ENDIF C DGC=DGPOC(IA) IF (DGC+1 .GT. MNCOF) THEN IER=46 RETURN ENDIF JTC=JTYPC(IA) AJTC=ABS(JTC) BETAC=JAINC(AJTC) LOMC=LSUBC(IA) LODC=(AJTC-1)*NQPTS+1 SJTC=SIGN(1E+0,REAL(JTC)) PSA=PRNSA(IA) C IF (PHPAS(IA+1) .LE. PHPAS(IA)) THEN BB=1E+0 ELSE BB=PHPAS(IA+1) ENDIF AA=PHPAS(IA) C IF (SJTC.GT.0) THEN CT0=CMPLX(MIDPT(PSA)+AA*HALEN(PSA)) ELSE CT0=CMPLX(MIDPT(PSA)+BB*HALEN(PSA)) ENDIF C TT=SJTC*TT FP=(1E+0+TT)**(BETAC+1E+0) IF (INTER) THEN FP=-FP ENDIF PHI=JACSUC(TT,DGC-1,A1COC(LODC),B1COC(LODC), + H1VLC(AJTC),BFSNC(LOMC+1)) PHI=(BFSNC(LOMC)-(1E+0-TT)*PHI)*SJTC*CMPLX(0E+0,1E+0) PHYPT=CENTR+(PARFUN(PARNT(PSA),CT0)-CENTR)*EXP(FP*PHI) C IF (WANTD) THEN IF (BETAC.LT.0E+0 .AND. (1E+0+TT).LE.0E+0) THEN C C WE ARE AT A CORNER WITH INFINITE DERIVATIVE. C DERIV=CMPLX(R1MACH(2),R1MACH(2)) ELSE IF (BETAC.GT.0E+0 .AND. (1E+0+TT).LE.0E+0) THEN C C WE ARE AT A CORNER WITH ZERO DERIVATIVE C DERIV=(0E+0,0E+0) ELSE DO 40 I=1,DGC+1 I1=I+LOMC-1 JACOF(I)=SOLNC(I1) 40 CONTINUE DO 50 I=2,DGC+1,2 JACOF(I)=SJTC*JACOF(I) 50 CONTINUE PHI=JACSUC(TT,DGC,ACOFC(LODC),BCOFC(LODC),H0VLC(AJTC), + JACOF) DERIV=TUPI*FP*PHI*(PHYPT-CENTR)* + CMPLX(COS(THETA),-SIN(THETA))/HA/(1E+0+TT) ENDIF ENDIF ELSE IA=IA+1 GOTO 30 ENDIF C C NORMAL EXIT C IER=0 END SUBROUTINE BMPHC1(IARC,PHYPT,CANPT,DERIV,NQPTS,TNSUA,DGPOL,JATYP, +LOSUB,LPQSB,NPPQF,PARNT,A1COF,ACOEF,B1COF,BCFSN,BCOEF,H0VAL,H1VAL, +HALEN,JACIN,MIDPT,SOLUN,TPPQF,TRRAD,VTARG,ZPPQF,INTER,WANTD, +IER) INTEGER IARC,IER,NQPTS,TNSUA INTEGER DGPOL(*),JATYP(*),LOSUB(*),LPQSB(*),NPPQF(*),PARNT(*) REAL A1COF(*),ACOEF(*),B1COF(*),BCFSN(*),BCOEF(*),H0VAL(*), +H1VAL(*),HALEN(*),JACIN(*),MIDPT(*),SOLUN(*),TPPQF(*),TRRAD(*), +VTARG(*) COMPLEX CANPT,DERIV,PHYPT,ZPPQF(*) LOGICAL INTER,WANTD C C GIVEN A POINT (DEFINED BY IARC AND PHYPT AS EXPLAINED NEXT) ON C THE BOUNDARY OF THE PHYSICAL DOMAIN, TO COMPUTE ITS IMAGE CANPT C ON THE UNIT DISC. IN CASE WANTD=.TRUE. THEN ALSO COMPUTE THE C DERIVATIVE DERIV OF THE MAP ONTO THE DISC AT THE GIVEN BOUNDARY C POINT. C C IF IARC > 0 THEN PHYPT IS THE PARAMETER VALUE OF THE C PHYSICAL POINT ON THE PARENT ANALYTIC ARC NUMBER IARC OTHERWISE C PHYPT DEFINES THE COORDINATES OF A PHYSICAL POINT SOMEWHERE ON C THE BOUNDARY. C C IER=0 - NORMAL EXIT. C IER=44 - LOCAL PARAMETER MNCOF NEEDS INCREASING. C IER=45 - THE PHYSICAL POINT DEFINED BY PHYPT (WITH IARC<0) HAS NOT C BEEN DETECTED AS LYING ON THE BOUNDARY; ACTUALLY, IT HAS C NOT BEEN DETECTED AS LYING PATHOLOGICALY CLOSE TO THE C BOUNDARY. CHECK THAT PHYPT IS CORRECT AND IF IT IS THEN C CONSIDER INCREASING THE PATHOLOGICAL TOLERANCE PARAMETER C PTHTL IN THE FIRST LINE BELOW. C C LOCAL VARIABLES C INTEGER AJT,DEG,I,I1,IA,K,J1,JQ,JT,LOD,LOM,MNCOF,NQ,PT REAL BETA,DIST,HL,JACSUM,MD,MINDS,NEWTL,PHI,PTHTL,R1MACH,RBCF, +SJT,SUM,TOLSM,TT,TUPI,TXI COMPLEX BCF,CJACSU,CT,DIFF2,DPARFN,PARFUN,PSI,WGHT,XI,ZXI, +ZTOB1,ZZ C PARAMETER (MNCOF=32) REAL JACOF(MNCOF) C EXTERNAL CJACSU,DPARFN,JACSUM,PARFUN,R1MACH,ZTOB1 C NEWTL=SQRT(R1MACH(4)) PTHTL=NEWTL TOLSM=1E+1*R1MACH(4) TUPI=8E+0*ATAN(1E+0) C IF (IARC .GT. 0) THEN C C *PHYPT* IS THE PARAMETER VALUE OF THE PHYSICAL POINT ON THE C PARENT ANALYTIC ARC NUMBER *IARC*. C C FIRST FIND THE CORRESPONDING SUBARC NUMBER AND LOCAL PARAMETER. C TT=REAL(PHYPT) IF (TT .LT. -1E+0) THEN TT=-1E+0 ENDIF IF (TT .GT. 1E+0) THEN TT=1E+0 ENDIF I=1 10 CONTINUE PT=PARNT(I) HL=HALEN(I) MD=MIDPT(I) SUM=MD+HL C IF (ABS(SUM-1E+0) .LT. TOLSM) THEN SUM=1E+0 ENDIF C IF (PT.EQ.IARC .AND. TT.LE.SUM) THEN IA=I TT=(TT-MD)/HL ELSE I=I+1 GOTO 10 ENDIF C IF (TT .LT. -1E+0) THEN TT=-1E+0 ENDIF IF (TT .GT. 1E+0) THEN TT=1E+0 ENDIF C JT=JATYP(IA) SJT=SIGN(1E+0,REAL(JT)) AJT=ABS(JT) BETA=JACIN(AJT) DEG=DGPOL(IA) IF (DEG+1 .GT. MNCOF) THEN IER=44 RETURN ENDIF LOM=LOSUB(IA) LOD=(AJT-1)*NQPTS+1 C TT=SJT*TT PHI=JACSUM(TT,DEG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) PHI=(1E+0+TT)**(BETA+1E+0)*(BCFSN(LOM)-(1E+0-TT)*PHI) IF (JT .GT. 0) THEN RBCF=VTARG(IA) ELSE RBCF=VTARG(IA+1) ENDIF RBCF=RBCF+SJT*PHI CANPT=CEXP((0E+0,1E+0)*RBCF) IF (WANTD) THEN IF (BETA.LT.0E+0 .AND. (1E+0+TT).LE.0E+0) THEN C C WE ARE AT A CORNER WITH INFINITE DERIVATIVE. C DERIV=CMPLX(R1MACH(2),R1MACH(2)) ELSE IF (BETA.GT.0E+0 .AND. (1E+0+TT).LE.0E+0) THEN C C WE ARE AT A CORNER WITH ZERO DERIVATIVE C DERIV=(0E+0,0E+0) ELSE DO 15 I=1,DEG+1 I1=I+LOM-1 JACOF(I)=SOLUN(I1) 15 CONTINUE DO 20 I=2,DEG+1,2 JACOF(I)=SJT*JACOF(I) 20 CONTINUE PHI=JACSUM(TT,DEG,ACOEF(LOD),BCOEF(LOD),H0VAL(AJT),JACOF) DERIV=TUPI*(1E+0+TT)**BETA*PHI*(0E+0,1E+0)*CANPT/ + DPARFN(IARC,PHYPT)/HL ENDIF ENDIF ELSE C C *PHYPT* IS A POINT SOMEWHERE ON THE PHYSICAL BOUNDARY. C ZZ=PHYPT DO 200 IA=1,TNSUA PT=PARNT(IA) JT=JATYP(IA) NQ=NPPQF(IA) K=LPQSB(IA)-1 HL=HALEN(IA) MD=MIDPT(IA) DO 100 JQ=1,NQ K=K+1 DIFF2=ZZ-ZPPQF(K) DIST=ABS(DIFF2) IF (DIST .LT. TRRAD(K)) THEN C C ZZ IS CLOSE TO ARC IA. C J1=JQ MINDS=DIST TXI=TPPQF(K) ZXI=ZPPQF(K) 40 CONTINUE J1=J1+1 IF (J1 .LE. NQ) THEN K=K+1 DIFF2=ZZ-ZPPQF(K) DIST=ABS(DIFF2) IF (DIST .LT. MINDS) THEN MINDS=DIST TXI=TPPQF(K) ZXI=ZPPQF(K) GOTO 40 ENDIF ENDIF C C PRELIMINARIES C SJT=SIGN(1E+0,REAL(JT)) AJT=ABS(JT) BETA=JACIN(AJT) DEG=DGPOL(IA) IF (DEG+1 .GT. MNCOF) THEN IER=44 RETURN ENDIF LOM=LOSUB(IA) LOD=(AJT-1)*NQPTS+1 C C NOW USE NEWTON'S METHOD TO ESTIMATE THE PARAMETRIC C PRE-IMAGE XI OF ZZ. C XI=CMPLX(TXI) CT=MD+HL*XI DIFF2=(ZXI-ZZ)/(DPARFN(PT,CT)*HL) XI=XI-DIFF2 50 CONTINUE IF (ABS(DIFF2) .GT. NEWTL) THEN CT=MD+HL*XI DIFF2=(PARFUN(PT,CT)-ZZ)/(DPARFN(PT,CT)*HL) XI=XI-DIFF2 GOTO 50 ELSE C C LAST ITERATION C CT=MD+HL*XI DIFF2=(PARFUN(PT,CT)-ZZ)/(DPARFN(PT,CT)*HL) XI=XI-DIFF2 ENDIF XI=SJT*XI C IF (ABS(AIMAG(XI)) .LT. PTHTL .AND. ABS(REAL(XI)) .LT. 1E+ + 0+PTHTL) THEN C C ZZ IS PATHOLOGICALLY CLOSE TO ARC IA AND WE USE THE C CONTINUATION OF THE BOUNDARY CORRESPONDENCE FUNCTION C TO ESTIMATE CANPT AND THE DERIVATIVE C PSI=CJACSU(XI,DEG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) WGHT=ZTOB1(XI+1E+0,BETA+1E+0,JT,INTER) PSI=WGHT*(BCFSN(LOM)-(1E+0-XI)*PSI) IF (JT .GT. 0) THEN BCF=VTARG(IA) ELSE BCF=VTARG(IA+1) ENDIF BCF=BCF+SJT*PSI CANPT=CEXP((0E+0,1E+0)*BCF) IF (WANTD) THEN IF ((1E+0+XI).EQ.(0E+0,0E+0) .AND. BETA.LT.0E+0) THEN C C WE ARE AT A CORNER WITH INFINITE DERIVATIVE. C DERIV=CMPLX(R1MACH(2),R1MACH(2)) ELSE IF ((1E+0+XI).EQ.(0E+0,0E+0) .AND. BETA.GT.0E+0) + THEN C C WE ARE AT A CORNER WITH ZERO DERIVATIVE. C DERIV=(0E+0,0E+0) ELSE DO 55 I=1,DEG+1 I1=I+LOM-1 JACOF(I)=SOLUN(I1) 55 CONTINUE DO 60 I=2,DEG+1,2 JACOF(I)=SJT*JACOF(I) 60 CONTINUE PSI=CJACSU(XI,DEG,ACOEF(LOD),BCOEF(LOD),H0VAL(AJT), + JACOF) CT=MD+HL*SJT*XI WGHT=WGHT/(1E+0+XI) DERIV=TUPI*WGHT*PSI*(0E+0,1E+0)*CANPT + /DPARFN(PT,CT)/HL ENDIF ENDIF GOTO 300 ENDIF C C END OF *IF (DIST .LT. TRRAD(K)) THEN* FOLLOWS C ENDIF 100 CONTINUE 200 CONTINUE IER=45 RETURN ENDIF C 300 CONTINUE C C NORMAL EXIT C IER=0 END SUBROUTINE BMPHYC(IARC,PHYPT,CANPT,DERIV,WANTD,INTER,IGEOM,RGEOM, +ISNPH,RSNPH,IQUPH,RQUPH,ZQUPH,WANTM,IER) C INTEGER IARC,IER INTEGER IGEOM(*),ISNPH(*),IQUPH(*) REAL RGEOM(*),RSNPH(*),RQUPH(*) COMPLEX PHYPT,CANPT,DERIV COMPLEX ZQUPH(*) LOGICAL WANTD,INTER,WANTM C C ...................................................................... C C 1. BMPHYC C BOUNDARY MAPPING FOR THE PHYSICAL --> CANONICAL MAP. C C 2. PURPOSE C GIVEN A POINT ON THE BOUNDARY OF THE PHYSICAL DOMAIN, THIS C ROUTINE COMPUTES THE CORRESPONDING APPROXIMATE IMAGE POINT C ON THE UNIT DISC AND, IF REQUESTED, ALSO COMPUTES THE C DERIVATIVE OF THE MAP : PHYSICAL --> CANONICAL AT THE GIVEN C POINT ON THE PHYSICAL BOUNDARY. C C 3. CALLING SEQUENCE C CALL BMPHYC(IARC,PHYPT,CANPT,DERIV,WANTD,INTER,IGEOM,RGEOM, C ISNPH,RSNPH,IQUPH,RQUPH,ZQUPH,WANTM,IER) C C PARAMETERS C ON ENTRY C IARC - INTEGER C ALLOWS TWO MODES OF DEFINING THE POINT ON THE C BOUNDARY OF THE PHYSICAL DOMAIN. IF IARC > 0 THEN C THE PHYSICAL POINT LIES ON ANALYTIC ARC NUMBER C IARC (AS DEFINED IN THE PARAMETRIC FUNCTION C PARFUN). IF IARC <= 0 THEN THE PARTICULAR ARC ON C WHICH THE PHYSICAL POINT LIES IS CONSIDERED TO BE C UNKNOWN ON ENTRY. C C PHYPT - COMPLEX C IF IARC > 0 THEN PHYPT IS THE (COMPLEX) PARAMETER C VALUE WHICH DEFINES THE PHYSICAL POINT ON ANALYTIC C ARC NUMBER IARC. IF IARC <= 0 THEN PHYPT IS THE C GIVEN PHYSICAL POINT. C C WANTD - LOGICAL C IF WANTD IS TRUE THEN DERIV IS COMPUTED OTHERWISE C DERIV ISN'T COMPUTED. C C INTER - LOGICAL C TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE C OTHERWISE. (AS PREVIOUSLY USED IN JAPHYC, GQPHYC) C C IGEOM - INTEGER ARRAY C THE INTEGER VECTOR IGEOM PREVIOUSLY SET UP BY C JAPHYC. C C RGEOM - REAL ARRAY C THE REAL VECTOR RGEOM PREVIOUSLY SET UP BY JAPHYC. C C ISNPH - INTEGER ARRAY C THE INTEGER VECTOR ISNPH PREVIOUSLY SET UP BY C JAPHYC. C C RSNPH - REAL ARRAY C THE REAL VECTOR RSNPH PREVIOUSLY SET UP BY JAPHYC. C C IQUPH - INTEGER ARRAY C THE INTEGER VECTOR IQUPH PREVIOUSLY SET UP BY C GQPHYC. C C RQUPH - REAL ARRAY C THE REAL ARRAY PREVIOUSLY SET UP BY GQPHYC. C C ZQUPH - COMPLEX ARRAY C THE COMPLEX ARRAY PREVIOUSLY SET UP BY GQPHYC. C C WANTM - LOGICAL C IF WANTM IS TRUE THEN, ON AN ABNORMAL EXIT, AN C ERROR MESSAGE IS WRITTEN ON THE STANDARD OUTPUT C CHANNEL. IF WANTM IS FALSE THEN NO MESSAGE IS C WRITTEN. C C C ON EXIT C CANPT - COMPLEX C THE POINT ON THE UNIT DISC CORRESPONDING TO THE C GIVEN PHYSICAL POINT UNDER THE MAP : PHYSICAL --> C CANONICAL. C C DERIV - COMPLEX C THE COMPUTED DERIVATIVE OF THE MAP : PHYSICAL --> C CANONICAL AT THE GIVEN POINT ON THE PHYSICAL BOUN- C DARY. THIS IS COMPUTED ONLY IF WANTD IS TRUE. C C IER - INTEGER C IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED; C A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY C WRITTEN ON THE STANDARD OUTPUT CHANNEL. C IER=0 - NORMAL EXIT. C IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD C BE SELF EXPLANATORY. C C C 4. SUBROUTINES OR FUNCTIONS NEEDED C - THE CONFPACK LIBRARY. C - THE REAL FUNCTION R1MACH. C - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN. C C C 5. FURTHER COMMENTS C - NOTE THAT THIS ROUTINE CAN ONLY BE USED A F T E R THE C ROUTINES JAPHYC AND GQPHYC HAVE SUCCESSFULLY EXECUTED, C AND THAT MANY INPUT ARGUMENTS FOR BMPHYC ARE OUTPUT VALUES C FROM JAPHYC AND GQPHYC. C C ...................................................................... C AUTHOR: DAVID HOUGH, ETH, ZUERICH C LAST UPDATE: 3 JULY 1990 C ...................................................................... C C LOCAL VARAIBLES C INTEGER ACOEF,AICOF,BCFSN,BCOEF,BICOF,DGPOL,H0VAL,HALEN,HIVAL, +JACIN,JATYP,LOSUB,LQSBF,MIDPT,MNSUA,MNEQN,MQUPH,NARCS,NJIND,NPPQF, +NQPTS,PARNT,SOLUN,TNGQP,TNSUA,TPPQF,TRRAD,VTARG,ZPPQF CHARACTER*6 IERTXT C EXTERNAL BMPHC1,IERTXT C NARCS=ISNPH(1) NJIND=NARCS+1 NQPTS=ISNPH(2) TNSUA=ISNPH(3) MNSUA=ISNPH(5) MNEQN=ISNPH(6) MQUPH=IQUPH(4) TNGQP=NQPTS*NJIND C C**** SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC C PARNT=5 HALEN=3 MIDPT=MNSUA+3 VTARG=2*MNSUA+3 C C**** SET UP THE POINTERS TO ISNPH AND RSNPH, AS IN JAPHYC C DGPOL=7 JATYP=MNSUA+7 LOSUB=2*MNSUA+7 ACOEF=1 BCOEF=TNGQP+1 AICOF=2*TNGQP+1 BICOF=3*TNGQP+1 H0VAL=6*TNGQP+1 HIVAL=NJIND+6*TNGQP+1 JACIN=2*NJIND+6*TNGQP+1 BCFSN=MNSUA+3*NJIND+6*TNGQP+1 SOLUN=MNEQN+MNSUA+3*NJIND+6*TNGQP+1 C C**** SET UP POINTERS TO IQUPH AND RQUPH, AS IN GQPHYC C LQSBF=5 NPPQF=MNSUA+5 TPPQF=2 TRRAD=MQUPH+2 ZPPQF=2 C C**** GET REQUIRED MAP AND DERIVATIVE C CALL BMPHC1(IARC,PHYPT,CANPT,DERIV,NQPTS,TNSUA,ISNPH(DGPOL), +ISNPH(JATYP),ISNPH(LOSUB),IQUPH(LQSBF),IQUPH(NPPQF),IGEOM(PARNT), +RSNPH(AICOF),RSNPH(ACOEF),RSNPH(BICOF),RSNPH(BCFSN),RSNPH(BCOEF), +RSNPH(H0VAL),RSNPH(HIVAL),RGEOM(HALEN),RSNPH(JACIN),RGEOM(MIDPT), +RSNPH(SOLUN),RQUPH(TPPQF),RQUPH(TRRAD),RGEOM(VTARG),ZQUPH(ZPPQF), +INTER,WANTD,IER) C C**** SEND ERROR MESSAGE TO STANDARD OUTPUT OF NECESSARY C IF (IER.GT.0 .AND. WANTM) WRITE(*,*) IERTXT(IER) C END SUBROUTINE CATPH4(NPTS,PHYPT,CANPT,NARCS,NQPTS,TNSUC,DGPOC,JTYPC, +LSUBC,LQSBG,NPPQG,PARNT,PRNSA,A1COC,ACOFC,B1COC,BCOFC,H0VLC, +H1VLC,HALEN,JAINC,LGTOL,MIDPT,PHPAS,QUPTC,QUWTC,THET0, +VARGC,BFSNC,CENTR,FACTR,SOLNC,WPPQG,ZPPQG,INTER,IER) INTEGER IER,NPTS,NARCS,NQPTS,TNSUC INTEGER DGPOC(*),JTYPC(*),LSUBC(*),LQSBG(*),NPPQG(*),PARNT(*), +PRNSA(*) REAL LGTOL,THET0 REAL A1COC(*),ACOFC(*),B1COC(*),BCOFC(*),H0VLC(*),H1VLC(*), +HALEN(*),JAINC(*),MIDPT(*),PHPAS(*),QUPTC(*),QUWTC(*),VARGC(*) COMPLEX CENTR,FACTR COMPLEX BFSNC(*),CANPT(*),PHYPT(*),SOLNC(*),WPPQG(*),ZPPQG(*) LOGICAL INTER C C GIVEN THE ARRAY CANPT OF NPTS POINTS IN THE CANONICAL PLANE, THIS C ROUTINE COMPUTES THE ARRAY PHYPT OF IMAGES IN THE PHYSICAL C PLANE. C C IER=0 - NORMAL EXIT C IER=47 - LOCAL PARAMETER MXNQD NEEDS INCREASING C IER=48 - LOCAL PARAMETER MNCOF NEEDS INCREASING C IER=49 - LOCAL PARAMETER MQIN1 NEEDS INCREASING C C....................................................................... C AUTHOR: DAVID HOUGH, ETH, ZUERICH C LAST UPDATE: 7 JULY 90 C....................................................................... C C LOCAL VARIABLES C INTEGER AJTC,DGC,I,IA,IP,K,J,J1,JQ,JTC,LODC,LOL, +LOMC,MNCOF,MQIN1,MXNQD,NBSCT,NQ,NQUAD,PSA,QINTS REAL AA,ARG,ARGW,AWW,BB,BETAC,DELTA,DIST,EFPTL,EPS,HA,ILW,IMXI,MD, +MEAN,PI,PTHTL,R1MACH,REXI,RLW,RR,RRB,S2C,SCO,SJTC,SUM1,TOLOU,TT, +TUPI,UPPER COMPLEX CT0,CT2,CTA,CTB,FP,CCJACS,CSUM,CT,JACSUC,PARFUN,PHI,TERM, +XI,WW PARAMETER (MNCOF=32,MQIN1=21,MXNQD=144,NBSCT=3, +PTHTL=1E-3,DELTA=2E-1) REAL XENPT(MQIN1) COMPLEX JCOFC(MNCOF),WSPEC(MXNQD),ZSPEC(MXNQD) EXTERNAL CCJACS,JACSUC,PARFUN,PPSBI1,R1MACH C EPS=1E+1*R1MACH(4) PI=4E+0*ATAN(1E+0) TUPI=2E+0*PI LOL=NARCS*NQPTS DO 300 IP=1,NPTS WW=CANPT(IP) AWW=ABS(WW) IF (AWW.LE.EPS) THEN PHYPT(IP)=CENTR GOTO 300 ENDIF RLW=LOG(AWW) ILW=ATAN2(AIMAG(WW),REAL(WW)) 10 CONTINUE IF (ILW .GT. THET0+TUPI) THEN ILW=ILW-TUPI GOTO 10 ENDIF 20 CONTINUE IF (ILW .LT. THET0) THEN ILW=ILW+TUPI GOTO 20 ENDIF CSUM=(0E+0,0E+0) DO 200 IA=1,TNSUC C C PRELIMINARIES FOR ARC IA C HA=(VARGC(IA+1)-VARGC(IA))*5E-1 MD=(VARGC(IA+1)+VARGC(IA))*5E-1 EFPTL=MAX(PTHTL,EPS/HA) IMXI=-RLW/HA C IF (ILW .GT. (MD+PI)) THEN ARGW=ILW-TUPI ELSE IF (ILW .LT. (MD-PI)) THEN ARGW=ILW+TUPI ELSE ARGW=ILW ENDIF C REXI=(ARGW-MD)/HA IF (REXI .GT. 1E+0) THEN DIST=SQRT(IMXI**2+(REXI-1E+0)**2) ELSE IF (REXI .LT. -1E+0) THEN DIST=SQRT(IMXI**2+(REXI+1E+0)**2) ELSE DIST=ABS(IMXI) ENDIF C IF (DIST .GE. DELTA) THEN C C USE THE STANDARD PRECOMPUTED COMPOSITE GAUSSIAN RULE C NQ=NPPQG(IA) K=LQSBG(IA)-1 DO 30 JQ=1,NQ K=K+1 CT=WW/ZPPQG(K) IF (.NOT. INTER) THEN CT=1E+0/CT ENDIF CSUM=CSUM+WPPQG(K)*CLOG(1E+0-CT) 30 CONTINUE C ELSE IF (DIST.LT.EFPTL) THEN C C WW IS PATHOLOGICALLY CLOSE TO ARC IA (OR MAYBE JUST CLOSE TO C ARC IA AND APPARENTLY OUTSIDE THE CANONICAL DOMAIN) AND WE C USE THE CONTINUATION OF THE BOUNDARY CORRESPONDENCE FUNCTION C TO ESTIMATE PHYPT. C C INITIALISE SOME DATA C XI=CMPLX(REXI,IMXI) DGC=DGPOC(IA) IF (DGC+1 .GT. MNCOF) THEN IER=48 RETURN ENDIF JTC=JTYPC(IA) AJTC=ABS(JTC) BETAC=JAINC(AJTC) LOMC=LSUBC(IA) LODC=(AJTC-1)*NQPTS+1 SJTC=SIGN(1E+0,REAL(JTC)) PSA=PRNSA(IA) C IF (PHPAS(IA+1) .LE. PHPAS(IA)) THEN BB=1E+0 ELSE BB=PHPAS(IA+1) ENDIF AA=PHPAS(IA) C IF (INTER) THEN S2C=SJTC ELSE S2C=-SJTC ENDIF C CTA=CMPLX(MIDPT(PSA)+AA*HALEN(PSA)) CTB=CMPLX(MIDPT(PSA)+BB*HALEN(PSA)) IF (SJTC.GT.0) THEN CT0=CTA CT2=CTB ELSE CT0=CTB CT2=CTA ENDIF C TERM=1E+0+SJTC*XI IF (TERM.EQ.(0E+0,0E+0)) THEN PHYPT(IP)=PARFUN(PARNT(PSA),CT0) GOTO 300 ENDIF C IF (TERM.EQ.(2E+0,0E+0)) THEN PHYPT(IP)=PARFUN(PARNT(PSA),CT2) GOTO 300 ENDIF C ARG=ATAN2(AIMAG(TERM),REAL(TERM)) UPPER=(1E+0+S2C*5E-1)*PI IF (ARG.GT.UPPER) THEN ARG=ARG-TUPI ELSE IF (ARG.LE.(UPPER-TUPI)) THEN ARG=ARG+TUPI ENDIF C FP=ABS(TERM)**(BETAC+1E+0)*CMPLX(COS((BETAC+1E+0)*ARG), + SIN((BETAC+1E+0)*ARG)) IF (INTER) THEN FP=-FP ENDIF PHI=CCJACS(SJTC*XI,DGC-1,A1COC(LODC),B1COC(LODC), + H1VLC(AJTC),BFSNC(LOMC+1)) PHI=(BFSNC(LOMC)-(1E+0-SJTC*XI)*PHI)*SJTC*CMPLX(0E+0,1E+0) PHYPT(IP)=CENTR+(PARFUN(PARNT(PSA),CT0)-CENTR)*EXP(FP*PHI) GOTO 300 ELSE C C SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS C PARTICULAR POINT WW. C C INITIALISE SOME DATA C DGC=DGPOC(IA) IF (DGC+1 .GT. MNCOF) THEN IER=48 RETURN ENDIF JTC=JTYPC(IA) AJTC=ABS(JTC) BETAC=JAINC(AJTC) LOMC=LSUBC(IA) LODC=(AJTC-1)*NQPTS+1 SJTC=SIGN(1E+0,REAL(JTC)) SCO=SJTC C DO 100 J=1,DGC+1 J1=LOMC+J-1 SCO=SCO*SJTC JCOFC(J)=SOLNC(J1)*SCO 100 CONTINUE C XI=SJTC*CMPLX(REXI,IMXI) CALL PPSBI1(XI,BETAC,NQPTS,DGC,ACOFC(LODC),BCOFC(LODC), + H0VLC(AJTC),JCOFC,LGTOL,TOLOU,XENPT,QINTS, + MQIN1,IER) IF (IER .GT. 0) THEN IF (IER .EQ. 29) THEN IER=49 ENDIF RETURN ENDIF NQUAD=QINTS*NQPTS IF (NQUAD .GT. MXNQD) THEN IER=47 RETURN ENDIF K=0 SUM1=BETAC+1E+0 DO 130 I=1,QINTS RR=(XENPT(I+1)-XENPT(I))*5E-1 MEAN=(XENPT(I+1)+XENPT(I))*5E-1 IF (I .EQ. 1) THEN RRB=RR**SUM1 DO 110 J=1,NQPTS J1=LODC+J-1 K=K+1 TT=MEAN+RR*QUPTC(J1) WSPEC(K)=RRB*QUWTC(J1)*JACSUC(TT,DGC,ACOFC(LODC), + BCOFC(LODC),H0VLC(AJTC),JCOFC) TT=MD+TT*SJTC*HA ZSPEC(K)=CMPLX(COS(TT),SIN(TT)) 110 CONTINUE ELSE DO 120 J=1,NQPTS J1=LOL+J K=K+1 TT=MEAN+RR*QUPTC(J1) WSPEC(K)=RR*QUWTC(J1)*(1E+0+TT)**BETAC*JACSUC(TT, + DGC,ACOFC(LODC),BCOFC(LODC),H0VLC(AJTC), + JCOFC) TT=MD+TT*SJTC*HA ZSPEC(K)=CMPLX(COS(TT),SIN(TT)) 120 CONTINUE ENDIF 130 CONTINUE C C THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS C AND POINTS WSPEC AND ZSPEC. NOW ESTIMATE THE INTEGRAL. C DO 140 K=1,NQUAD CT=WW/ZSPEC(K) IF (.NOT. INTER) THEN CT=1E+0/CT ENDIF CSUM=CSUM+WSPEC(K)*CLOG(1E+0-CT) 140 CONTINUE C C END OF ELSE BLOCK RELATING TO SPECIAL QUADRATURE RULE FOR C WW NEAR ARC IA C ENDIF C C END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA C 200 CONTINUE PHYPT(IP)=CENTR+FACTR*WW*CEXP(CSUM) C C END OF MAP CALCULATION FOR FIELD POINT NUMBER IP C 300 CONTINUE C IER=0 C END COMPLEX FUNCTION CCJACS(X,N,A,B,H,CO) INTEGER N REAL A(*),B(*),H COMPLEX CO(*),X C ..TO CALCULATE SUMMATION{CO(K+1)*P(K,X)},K=0(1)N, WHERE P(K,X) C ..DENOTES THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE K C ..EVALUATED AT X, ARRAY CO STORES A GIVEN SET OF COEFFICIENTS, C ..ARRAYS A,B STORE THE COEFFICIENTS IN THE THREE-TERM C ..RECURRENCE FORMULA FOR THE JACOBI POLYNOMIALS (SEE ASONJ7) C ..AND H IS THE SQUARED 2-NORM OF UNITY. COMPLEX PREV,CURR,NEXT INTEGER K C IF (N .EQ. 0) THEN CCJACS=CO(1)/SQRT(H) ELSE IF (N .GT. 0) THEN PREV=CO(N+1) CURR=CO(N)+(X-B(N))*PREV/A(N) DO 10 K=N-2,0,-1 NEXT=CO(K+1)+(X-B(K+1))*CURR/A(K+1)-A(K+1)*PREV/A(K+2) PREV=CURR CURR=NEXT 10 CONTINUE CCJACS=CURR/SQRT(H) ELSE CCJACS=(0E+0,0E+0) ENDIF C END LOGICAL FUNCTION CHRIN(A,B) CHARACTER*1 A,B C C**** TEST WHETHER THE EITHER CHARACTER A OR B IS CONTAINED IN THE FIRST C**** 6 CHARACTERS OF THE NEXT RECORD FROM STANDARD INPUT. C C LOCAL VARIABLES C CHARACTER*6 C C READ(*,'(A6)') C CHRIN=(INDEX(C,A).NE.0 .OR. INDEX(C,B).NE.0) C END COMPLEX FUNCTION CINRAD(NARCS,NQPTS,TNSUA,DGPOL,JATYP,LOSUB,LPQSB, +NPPQF,PARNT,ACOEF,BCOEF,H0VAL,HALEN,JACIN,LGTOL,MIDPT,QUPTS,QUWTS, +SOLUN,TPPQF,TRRAD,WPPQF,CENTR,FACTR,ZPPQF,IER) INTEGER IER,NARCS,NQPTS,TNSUA INTEGER DGPOL(*),JATYP(*),LOSUB(*),LPQSB(*),NPPQF(*),PARNT(*) REAL LGTOL REAL ACOEF(*),BCOEF(*), +H0VAL(*),HALEN(*),JACIN(*),MIDPT(*),QUPTS(*),QUWTS(*),SOLUN(*), +TPPQF(*),TRRAD(*),WPPQF(*) COMPLEX CENTR,FACTR COMPLEX ZPPQF(*) C C TO COMPUTE THE COMPLEX INNER RADIUS (I.E. THE RECIPROCAL OF THE C DERIVATIVE OF THE INTERIOR MAP AT THE CENTRE POINT OF THE PHYSICAL C DOMAIN) C C IER=0 - NORMAL EXIT C IER=37 - LOCAL PARAMETER MXNQD NEEDS INCREASING C IER=38 - LOCAL PARAMETER MNCOF NEEDS INCREASING C IER=39 - THE CENTRE POINT IS PATHOLOGICALLY CLOSE TO THE C BOUNDARY C IER=40 - LOCAL PARAMETER MQIN1 MUST BE INCREASED C C LOCAL VARIABLES C INTEGER AJT,DEG,I,IA,K,J,J1,J2,JQ,JT,LIM,LOD,LOL,LOM,MNCOF, +MQIN1,MXNQD,NQ,NQUAD,PT,QINTS REAL AISUM,ANGLE,ARGBR,ARGIN1,ARSUM,BETA,CURARG,DIST,HL,ISUM, +JACSUM,LIMIT,MD,MEAN,MINDS,NEWTL,PI,PTHTL,R1MACH,RR,RRB,RSUM,RT1, +RT2,SCO,SS,STARG,STRT1,STTH1,SUM1,THET1,THET2,TOLOU,TT,TXI,TUPI,WT COMPLEX CT,DPARFN,PARFUN,XI,DIFF1,DIFF2, +STDF1,ZXI,ZZ LOGICAL FIRST EXTERNAL ARGIN1,DPARFN,JACSUM,PARFUN,PPSBI1,R1MACH PARAMETER (MNCOF=32,MQIN1=11,MXNQD=80,PTHTL=1E-3) PARAMETER (LIMIT=2.3562E+0) REAL JACOF(MNCOF),TSPEC(MXNQD),WSPEC(MXNQD),XENPT(MQIN1) COMPLEX JCOFC(MNCOF),ZSPEC(MXNQD) C NEWTL=SQRT(R1MACH(4)) PI=4E+0*ATAN(1E+0) TUPI=2E+0*PI LOL=NARCS*NQPTS ZZ=CENTR RSUM=0E+0 ISUM=0E+0 FIRST=.TRUE. DO 200 IA=1,TNSUA PT=PARNT(IA) JT=JATYP(IA) NQ=NPPQF(IA) K=LPQSB(IA)-1 HL=HALEN(IA) MD=MIDPT(IA) ARSUM=0E+0 AISUM=0E+0 DO 100 JQ=1,NQ K=K+1 DIFF2=ZZ-ZPPQF(K) RT2=MD+HL*TPPQF(K) DIST=ABS(DIFF2) IF (DIST .GE. TRRAD(K)) THEN WT=WPPQF(K) IF (WT .NE. 0E+0) THEN ARSUM=ARSUM+WT*LOG(DIST) IF (FIRST) THEN CURARG=ATAN2(AIMAG(DIFF2),REAL(DIFF2)) THET2=CURARG FIRST=.FALSE. STARG=CURARG ELSE C CT=DIFF2/DIFF1 C CT=DIFF2*CONJG(DIFF1) C ANGLE=ATAN2(AIMAG(CT),REAL(CT)) THET2=ATAN2(AIMAG(DIFF2),REAL(DIFF2)) ANGLE=THET2-THET1 IF (ANGLE .LE. -PI .OR. ANGLE .GT. PI) THEN ANGLE=ANGLE-SIGN(TUPI,ANGLE) ENDIF IF (ABS(ANGLE) .GE. LIMIT) THEN ANGLE=ARGIN1(RT1,RT2,PT,-DIFF1,-DIFF2,ZZ, + LIMIT) ENDIF CURARG=CURARG+ANGLE ENDIF AISUM=CURARG*WT+AISUM RT1=RT2 DIFF1=DIFF2 THET1=THET2 ENDIF ELSE C C ZZ IS TOO CLOSE TO ARC IA TO USE THE STANDARD RULE. C FIND THE QUADRATURE POINT NEAREST TO ZZ. C J1=JQ MINDS=DIST TXI=TPPQF(K) ZXI=ZPPQF(K) 40 CONTINUE J1=J1+1 IF (J1 .LE. NQ) THEN K=K+1 DIFF2=ZZ-ZPPQF(K) DIST=ABS(DIFF2) IF (DIST .LT. MINDS) THEN MINDS=DIST TXI=TPPQF(K) ZXI=ZPPQF(K) GOTO 40 ENDIF ENDIF C C PRELIMINARIES C IF (JT .GT. 0) THEN SS=1E+0 ELSE SS=-1E+0 ENDIF AJT=ABS(JT) BETA=JACIN(AJT) DEG=DGPOL(IA) IF (DEG+1 .GT. MNCOF) THEN IER=38 RETURN ENDIF LOM=LOSUB(IA) LOD=(AJT-1)*NQPTS+1 C C NOW USE NEWTON'S METHOD TO ESTIMATE THE PARAMETRIC C PRE-IMAGE XI OF ZZ. C XI=CMPLX(TXI) CT=MD+HL*XI DIFF2=(ZXI-ZZ)/(DPARFN(PT,CT)*HL) XI=XI-DIFF2 50 CONTINUE IF (ABS(DIFF2) .GT. NEWTL) THEN CT=MD+HL*XI DIFF2=(PARFUN(PT,CT)-ZZ)/(DPARFN(PT,CT)*HL) XI=XI-DIFF2 GOTO 50 ENDIF XI=SS*XI C IF (ABS(AIMAG(XI)) .LT. PTHTL .AND. ABS(REAL(XI)) .LT. 1E+ + 0+PTHTL) THEN C C THE CENTRE OF THE DOMAIN (I.E. ZZ) IS PATHOLOGICALLY C CLOSE TO ARC IA AND WE DO NOT ALLOW THIS. C IER=39 RETURN ELSE C C SET UP A SPECIAL COMPOSITE GAUSSIAN RULE TO HANDLE THIS C PARTICULAR POINT ZZ. C SCO=SS DO 55 J=1,DEG+1 J1=LOM+J-1 SCO=SCO*SS JACOF(J)=SOLUN(J1)*SCO JCOFC(J)=CMPLX(SOLUN(J1)*SCO) 55 CONTINUE CALL PPSBI1(XI,BETA,NQPTS,DEG,ACOEF(LOD),BCOEF(LOD), + H0VAL(AJT),JCOFC,LGTOL,TOLOU,XENPT,QINTS, + MQIN1,IER) IF (IER .GT. 0) THEN IF (IER .EQ. 29) THEN IER=40 ENDIF RETURN ENDIF NQUAD=QINTS*NQPTS IF (NQUAD .GT. MXNQD) THEN IER=37 RETURN ENDIF K=0 SUM1=BETA+1E+0 DO 70 I=1,QINTS RR=(XENPT(I+1)-XENPT(I))*5E-1 MEAN=(XENPT(I+1)+XENPT(I))*5E-1 IF (I .EQ. 1) THEN RRB=RR**SUM1 DO 60 J=1,NQPTS J1=LOD+J-1 K=K+1 TT=(MEAN+RR*QUPTS(J1)) WSPEC(K)=RRB*QUWTS(J1)*JACSUM(TT,DEG,ACOEF(LOD), + BCOEF(LOD),H0VAL(AJT),JACOF) TT=TT*SS TSPEC(K)=MD+TT*HL CT=CMPLX(TSPEC(K)) ZSPEC(K)=PARFUN(PT,CT) 60 CONTINUE ELSE DO 65 J=1,NQPTS J1=LOL+J K=K+1 TT=(MEAN+RR*QUPTS(J1)) WSPEC(K)=RR*QUWTS(J1)*(1E+0+TT)**BETA*JACSUM(TT, + DEG,ACOEF(LOD),BCOEF(LOD),H0VAL(AJT), + JACOF) TT=TT*SS TSPEC(K)=MD+TT*HL CT=CMPLX(TSPEC(K)) ZSPEC(K)=PARFUN(PT,CT) 65 CONTINUE ENDIF 70 CONTINUE IF (SS .LT. 0E+0) THEN LIM=NQUAD IF (MOD(LIM,2) .EQ. 0) THEN LIM=LIM/2 ELSE LIM=(LIM-1)/2 ENDIF J1=0 J2=NQUAD+1 DO 75 J=1,LIM J1=J1+1 J2=J2-1 TT=WSPEC(J1) WSPEC(J1)=WSPEC(J2) WSPEC(J2)=TT TT=TSPEC(J1) TSPEC(J1)=TSPEC(J2) TSPEC(J2)=TT CT=ZSPEC(J1) ZSPEC(J1)=ZSPEC(J2) ZSPEC(J2)=CT 75 CONTINUE ENDIF C C THIS COMPLETES THE SETTING UP OF THE SPECIAL WEIGHTS C AND POINTS WSPEC AND ZSPEC. NOW ESTIMATE THE INTEGRAL. C ARSUM=0E+0 AISUM=0E+0 IF (IA .EQ. 1) THEN FIRST=.TRUE. ELSE CURARG=STARG RT1=STRT1 DIFF1=STDF1 THET1=STTH1 ENDIF DO 80 K=1,NQUAD WT=WSPEC(K) DIFF2=ZZ-ZSPEC(K) RT2=TSPEC(K) DIST=ABS(DIFF2) ARSUM=ARSUM+WT*LOG(DIST) IF (FIRST) THEN CURARG=ATAN2(AIMAG(DIFF2),REAL(DIFF2)) THET2=CURARG FIRST=.FALSE. ELSE C CT=DIFF2/DIFF1 C CT=DIFF2*CONJG(DIFF1) C ANGLE=ATAN2(AIMAG(CT),REAL(CT)) THET2=ATAN2(AIMAG(DIFF2),REAL(DIFF2)) ANGLE=THET2-THET1 IF (ANGLE .LE. -PI .OR. ANGLE .GT. PI) THEN ANGLE=ANGLE-SIGN(TUPI,ANGLE) ENDIF IF (ABS(ANGLE) .GE. LIMIT) THEN ANGLE=ARGIN1(RT1,RT2,PT,-DIFF1,-DIFF2,ZZ, + LIMIT) ENDIF CURARG=CURARG+ANGLE ENDIF AISUM=CURARG*WT+AISUM RT1=RT2 DIFF1=DIFF2 THET1=THET2 80 CONTINUE GOTO 150 ENDIF ENDIF C C END OF QUADRATURE SUM LOOP C 100 CONTINUE C 150 CONTINUE RSUM=RSUM+ARSUM ISUM=ISUM+AISUM IF (JT .LT. 0) THEN C C BRING THE ARGUMENT FORWARD TO THE CORNER POINT AND REPLACE C THE INCREMENTED CURARG VALUE BY AN INVERSE TANGENT C EVALUATION C DIFF2=ZZ-PARFUN(PT,(1E+0,0E+0)) RT2=1E+0 THET2=ATAN2(AIMAG(DIFF2),REAL(DIFF2)) ANGLE=THET2-THET1 IF (ANGLE .LE. -PI .OR. ANGLE .GT. PI) THEN ANGLE=ANGLE-SIGN(TUPI,ANGLE) ENDIF IF (ABS(ANGLE) .GE. LIMIT) THEN ANGLE=ARGIN1(RT1,RT2,PT,-DIFF1,-DIFF2,ZZ, + LIMIT) ENDIF CURARG=CURARG+ANGLE ARGBR=ANINT((CURARG-THET2)/TUPI) CURARG=THET2+TUPI*ARGBR RT1=-1E+0 DIFF1=DIFF2 THET1=THET2 ENDIF STARG=CURARG STRT1=RT1 STDF1=DIFF1 STTH1=THET1 C C END OF LOOP FOR CONTRIBUTIONS FROM ARC NUMBER IA C 200 CONTINUE CT=CMPLX(RSUM,ISUM) CT=CEXP(CT) CINRAD=CT/FACTR C IER=0 C END SUBROUTINE CNDPLT(MAP11,RESMN,UPHYC,UCANP,CRRES,IGEOM,RGEOM,ISNPH, +RSNPH,CH0,CH1,DASH,NEWD,IER) C INTEGER CH0,CH1,IER INTEGER IGEOM(*),ISNPH(*) REAL RESMN,UPHYC,UCANP,CRRES REAL RGEOM(*),RSNPH(*) LOGICAL MAP11 CHARACTER DASH*(*),NEWD*(*) C C ...................................................................... C C 1. CNDPLT C REPORTS ON THE CONDITION OF THE PROBLEM OF EVALUATING THE C MAPPING FUNCTIONS AND ALSO OUTPUTS DATA FOR GRAPH PLOTTING. C C 2. PURPOSE C THE ROUTINE COMPUTES CONDITION NUMBERS FOR THE PROBLEMS OF C EVALUATING THE TWO MAPS PHYSICAL --> CANONICAL AND C CANONICAL --> PHYSICAL AND COMPUTES THE ERROR THAT MAY BE C EXPECTED (IN THE WORST CASE) IN THE RANGE OF EACH APPROX- C IMATE MAP FROM A MACHINE PRECISION LEVEL ROUNDING ERROR IN C THE DOMAIN OF EACH MAP. C THE ROUTINE ALSO COMPUTES THE LEAST RESOLUTION OF THE C COMPUTED MAP : PHYSICAL --> CANONICAL OVER ALL SUB-ARCS ON C THE PHYSICAL BOUNDARY. THE RESOLUTION OF THE MAP FOR ANY C PHYSICAL SUB-ARC IS DEFINED AS THE COMPUTED ANGULAR WIDTH OF C THE IMAGE SUB-ARC ON THE UNIT DISC DIVIDED BY THE ESTIMATED C MAXIMUM ERROR IN THE MODULUS OF THE MAP ON THE GIVEN SUB- C ARC. A LEAST RESOLUTION OF LESS THAN, SAY, 10 INDICATES C THAT THERE ARE REGIONS OF SEVERE CROWDING AND THAT IT WILL C BE PRACTICALLY IMPOSSIBLE TO COMPUTE THE INVERSE MAP C EVERYWHERE ON THE CANONICAL DOMAIN. C THE ROUTINE ALSO SEARCHES (NOT VERY EXHAUSTIVELY) FOR C CHANGES OF SIGN IN THE COMPUTED BOUNDARY CORRESPONDENCE C DERIVATIVE FOR THE MAP : PHYSICAL --> CANONICAL. SUCH C SIGN CHANGES MEAN THAT THE COMPUTED MAP IS NOT ONE-TO-ONE C AND HENCE ONE SHOULD EXPECT DIFFICULTIES IN TRYING TO C COMPUTE THE INVERSE MAP : CANONICAL --> PHYSICAL. C FINALLY THREE OUTPUT FILES C cn, p0, p1 C ARE WRITTEN. THE FIRST OF THESE IS A SUMMARY OF THE ABOVE C RESULTS INTENDED TO BE READ BY THE USER. THE TWO FILES C p0 AND p1 ARE NOT INTENDED TO BE READ BY THE C USER, BUT COULD BE USED TO CREATE PLOTS OF THE BOUNDARY C CORRESPONDENCE FUNCTION AND ITS DERIVATIVE; SEE FURTHER C COMMENTS BELOW. C C 3. CALLING SEQUENCE C CALL CNDPLT(MAP11,RESMN,UPHYC,UCANP,CRRES,IGEOM,RGEOM,ISNPH, C RSNPH,CH0,CH1,DASH,NEWD,IER) C C PARAMETERS C ON ENTRY C CRRES - REAL C THE CRITICAL RESOLUTION. IF THE COMPUTED RESOL- C UTION OF THE PHYSICAL-->CANONICAL MAP ON ANY ARC C FALLS BELOW CRRES THAN A WARNING MESSAGE IS C OUTPUT. C C IGEOM - INTEGER ARRAY C THE INTEGER VECTOR IGEOM PREVIOUSLY SET UP BY C JAPHYC. C C RGEOM - REAL ARRAY C THE REAL VECTOR RGEOM PREVIOUSLY SET UP BY JAPHYC. C C ISNPH - INTEGER ARRAY C THE INTEGER VECTOR ISNPH PREVIOUSLY SET UP BY C JAPHYC. C C RSNPH - REAL ARRAY C THE REAL VECTOR RSNPH PREVIOUSLY SET UP BY JAPHYC. C C CH0 - INTEGER C DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR C WRITING THE FILES cn AND p0. C C CH1 - INTEGER C DEFINES AN OUTPUT CHANNEL THAT MAY BE USED FOR C WRITING THE FILE p1; MUST HAVE CH0.NE.CH1. C C DASH - CHARACTER C A CHARACTER VARIABLE OF USER-DEFINED LENGTH WHICH C DEFINES THE DASH-PATTERN THAT THE USER MAY REQUIRE C FOR GRAPH PLOTTING; SEE FURTHER COMMENTS BELOW C C NEWD - CHARACTER C A CHARACTER VARIABLE OF USER-DEFINED LENGTH WHICH C DENOTES THE START OF A NEW DATA GROUP THAT THE C USER MAY REQUIRE FOR GRAPH PLOTTING; SEE FURTHER C COMMENTS BELOW C ON EXIT C MAP11 - LOGICAL C IF BOUNDARY REVERSALS ARE DETECTED THEN MAP11 IS C SET TO .FALSE. (THE COMPUTED PHYSICAL-->CANONICAL C MAP ISN'T 1-1) OTHERWISE MAP11 IS SET TO .TRUE. C C RESMN - REAL C THE MINIMUM COMPUTED RESOLUTION OF THE PHYSICAL--> C CANONICAL MAP OVER ALL SUBARCS ON THE PHYSICAL C BOUNDARY. IF RESMN IS LESS THAN CRRES THEN C A WARNING MESSAGE IS OUTPUT. C C UPHYC - REAL C ESTIMATED MAXIMUM POSSIBLE ERROR IN THE RANGE OF C THE PHYSICAL-->CANONICAL MAP DUE UNIT ROUNDOFF IN C THE PHYSICAL DOMAIN. C C UCANP - REAL C ESTIMATED MAXIMUM POSSIBLE ERROR IN THE RANGE OF C THE CANONICAL-->PHYSICAL MAP DUE UNIT ROUNDOFF IN C THE CANONICAL DOMAIN. C C IER - INTEGER C IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED; C A MESSAGE TO DESCRIBE THE ERROR IS AUTOMATICALLY C WRITTEN ON THE STANDARD OUTPUT CHANNEL AND THE C LISTING FILE cn. C IER=0 - NORMAL EXIT. C IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD C BE SELF EXPLANATORY. C C C 4. SUBROUTINES OR FUNCTIONS NEEDED C - THE CONFPACK LIBRARY. C - THE REAL FUNCTION R1MACH. C - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN. C C C 5. FURTHER COMMENTS C - NOTE THAT THIS ROUTINE CAN ONLY BE USED A F T E R THE C ROUTINE JAPHYC HAS SUCCESSFULLY EXECUTED, AND THAT MOST C INPUT ARGUMENTS FOR CNDPLT ARE OUTPUT VALUES FROM JAPHYC. C - A DETAILED LISTING OF RESULTS IS WRITTEN ON THE FILE C cn. C - DATA FOR PLOTTING A GRAPH OF THE DIMENSIONLESS BOUNDARY C CORRESPONDENCE FUNCTION AGAINST DIMENSIONLESS ARC LENGTH C ARE WRITTEN ON THE FILE p0. THE CONTENTS OF THIS C FILE ARE AS FOLLOWS: C 1. ABOUT 200 COORDINATE PAIRS X Y, ONE PAIR PER LINE, C WHERE X = DIMENSIONLESS ARC LENGTH (0 <= X <=1) AND C Y = DIMENSIONLESS BOUNDARY CORRESPONDENCE FUNCTION C (0 <= Y <=1); THE NUMBER OF COORDINATE PAIRS IS C CONTROLLED BY THE LOCAL PARAMETER NXINT. C 2. THE SINGLE LINE C C WHERE DENOTES THE VALUE OF THE ARGUMENT DASH; C THIS CAN BE USED TO INDICATE A CHANGE OF DASH PATTERN C TO THE LOCAL GRAPH PLOTTER. C 3. SEVERAL REPETITIONS OF THE FOLLOWING 3-LINE GROUP: C C X 0E+0 C X 1E+0 C HERE DENOTES THE VALUE OF THE ARGUMENT NEWD AND C X (WITH 0 < X < 1) IS THE DIMENSIONLESS ARC C LENGTH OF A CORNER POINT. THE ABOVE GROUP MAY THEN C BE USED TO CONSTRUCT A DASHED LINE FROM (X,0) TO C (X,1). THE NUMBER OF REPETITIONS IS EQUAL TO THE C NUMBER OF CORNERS WITH ARC LENGTH IN THE INTERVAL C 0 < X < 1. C - DATA FOR PLOTTING A GRAPH OF THE DERIVATIVE OF THE DIMEN- C SIONLESS BOUNDARY CORRESPONDENCE FUNCTION WITH RESPECT TO C DIMENSIONLESS ARC LENGTH ARE WRITTEN ON THE FILE p1. C THE CONTENTS OF THIS FILE ARE AS FOLLOWS: C 1. ABOUT 200 COORDINATE PAIRS X Y, ONE PAIR PER LINE, C WHERE X = DIMENSIONLESS ARC LENGTH (0 <= X <=1) AND C Y = DIMENSIONLESS BOUNDARY CORRESPONDENCE DERIVATIVE C (0 <= Y <=1); THE NUMBER OF COORDINATE PAIRS IS C CONTROLLED BY THE LOCAL PARAMETER NXINT. C 2. THE SINGLE LINE C C WHERE DENOTES THE VALUE OF THE ARGUMENT DASH; C THIS CAN BE USED TO INDICATE A CHANGE OF DASH PATTERN C TO THE LOCAL GRAPH PLOTTER. C 3. SEVERAL REPETITIONS OF THE FOLLOWING 3-LINE GROUP: C C X 0E+0 C X 4.4E+0 C HERE DENOTES THE VALUE OF THE ARGUMENT NEWD AND C X (WITH 0 < X < 1) IS THE DIMENSIONLESS ARC C LENGTH OF A RE-ENTRANT CORNER POINT. THE ABOVE GROUP C MAY THEN BE USED TO CONSTRUCT A DASHED LINE FROM (X,0) C TO (X,4.4), TO INDICATE THE PRESENCE OF AN ASYMPTOTE. C SINCE THE AVERAGE VALUE OF THE DIMENSIONLESS BCF DERI- C VATIVE IS 1, 4.4 IS AN ARBITRARILY CHOSEN BUT RELATI- C VELY LARGE HEIGHT AT WHICH TO TERMINATE AN ASYMPTOTE; C THIS HEIGHT IS CONTROLLED BY THE LOCAL PARAMETER BIG. C THE NUMBER OF REPETITIONS OF THIS GROUP IS EQUAL TO C THE NUMBER OF RE-ENTRANT CORNERS WITH DIMENSIONLESS C ARC LENGTH IN THE INTERVAL 0 < X < 1. C - A SUMMARY LISTING OF RESULTS IS AUTOMATICALLY C WRITTEN ON THE STANDARD OUTPUT CHANNEL. C C ...................................................................... C AUTHOR: DAVID HOUGH, ETH, ZUERICH C LAST UPDATE: 17 JULY 1990 C ...................................................................... C C LOCAL VARAIBLES C INTEGER ACOEF,AICOF,BCFSN,BCOEF,BICOF,DGPOL,ERARC,H0VAL,HALEN, +HIVAL,JACIN,JATYP,LOSUB,MIDPT,PARNT,QUPTS,QUWTS,SOLUN,VTARG INTEGER I,IMNLA,J,L,LODP,LODW,MAXSA,MNEQN,MNSUA,NARCS,NASYM,NCRVS, +NINFD,NJIND,NPRVS,NQPTS,NXINT,NZERD,TNGQP,TNSUA REAL ANGSP,BIG,CCAPH,COCAP,COPHC,CPHCA,CR,EXCAP,EXPHC,LA, +OFLOW,PI,R1MACH,TOTLN,MCHEP CHARACTER OFLC*6,OFP0*6,OFP1*6,JBNM*4,CHPC*2,CHCP*2 C PARAMETER(NXINT=200,MAXSA=100,BIG=4.4E+0) C C**** NXINT = GLOBAL NUMBER OF INTERVALS ON [0,1] FOR SAMPLING THE C**** DIMENSIONLESS DERIVATE OF THE BOUNDARY CORRESPONDENCE C**** FUNCTION. C**** IER=52 - LOCAL PARAMETER MAXSA MUST BE INCEASED TO AT LEAST THE C**** VALUE OF ARGUMENT ISNPH(3)=TOTAL NUMBER OF SUBARCS ON C**** PHYSICAL BOUNDARY. C INTEGER ICRVS(MAXSA),IPRVS(MAXSA) REAL ARCLN(MAXSA),ASYMP(MAXSA),BCDMN(MAXSA),CORXX(MAXSA) EXTERNAL DIAGN4,R1MACH,WRHEAD,WRTAIL C C**** WRITE HEADING TO STANDARD OUTPUT CHANNEL C CALL WRHEAD(5,0) C C GET JOBNAME FROM FILE *jbnm* C OPEN(CH0,FILE='jbnm') READ(CH0,'(A4)') JBNM CLOSE(CH0) L=INDEX(JBNM,' ')-1 IF (L.EQ.-1) L=4 C OFLC=JBNM(1:L)//'cn' OFP1=JBNM(1:L)//'p1' OFP0=JBNM(1:L)//'p0' C NARCS=ISNPH(1) NQPTS=ISNPH(2) TNSUA=ISNPH(3) MNSUA=ISNPH(5) MNEQN=ISNPH(6) C NJIND=NARCS+1 TNGQP=NJIND*NQPTS C IF (TNSUA .GT. MAXSA) THEN IER=52 GOTO 999 ENDIF C C**** COPY POINTERS FROM JAPHYC C PARNT=5 HALEN=3 MIDPT=MNSUA+3 VTARG=2*MNSUA+3 DGPOL=7 JATYP=MNSUA+7 LOSUB=2*MNSUA+7 ACOEF=1 BCOEF=TNGQP+1 AICOF=2*TNGQP+1 BICOF=3*TNGQP+1 QUPTS=4*TNGQP+1 QUWTS=5*TNGQP+1 H0VAL=6*TNGQP+1 HIVAL=NJIND+6*TNGQP+1 JACIN=2*NJIND+6*TNGQP+1 ERARC=3*NJIND+6*TNGQP+1 BCFSN=MNSUA+3*NJIND+6*TNGQP+1 SOLUN=MNEQN+MNSUA+3*NJIND+6*TNGQP+1 C OPEN(CH0,FILE=OFP0) OPEN(CH1,FILE=OFP1) WRITE(*,5) 'EVALUATION OF BCF STARTED:' 5 FORMAT(A45) LODP=QUPTS+NARCS*NQPTS LODW=QUWTS+NARCS*NQPTS CALL DIAGN4(CCAPH,COCAP,COPHC,CPHCA,EXCAP,EXPHC,ICRVS,IER, +IPRVS,NASYM,NCRVS,NINFD,NPRVS,NZERD,ARCLN,ASYMP,BCDMN,CORXX, +TOTLN,RGEOM(VTARG),MAP11,ISNPH(DGPOL),ISNPH(JATYP),ISNPH(LOSUB), +NARCS,NQPTS,NXINT,CH0,CH1,IGEOM(PARNT),TNSUA,RSNPH(AICOF), +RSNPH(ACOEF),RSNPH(BICOF),RSNPH(BCFSN),RSNPH(BCOEF),RSNPH(H0VAL), +RSNPH(HIVAL),RGEOM(HALEN),RSNPH(JACIN),RGEOM(MIDPT),RSNPH(SOLUN), +RSNPH(LODP),RSNPH(LODW)) WRITE(*,5) 'EVALUATION OF BCF DONE:' C IF (IER .GT. 0) THEN GOTO 999 ENDIF C IF (NASYM .GT. 0) THEN WRITE(CH1,*) DASH DO 10 I=1,NASYM WRITE(CH1,*) NEWD WRITE(CH1,20) ASYMP(I),0E+0 WRITE(CH1,20) ASYMP(I),BIG 10 CONTINUE 20 FORMAT(2E16.8) ENDIF CLOSE(CH1) C WRITE(CH0,*) DASH DO 30 I=2,NARCS WRITE(CH0,*) NEWD WRITE(CH0,20) CORXX(I),0E+0 WRITE(CH0,20) CORXX(I),1E+0 30 CONTINUE CLOSE(CH0) WRITE(*,5) 'DATA FOR PLOTS DONE:' C OFLOW=R1MACH(2) MCHEP=R1MACH(4) UPHYC=MCHEP*CPHCA UCANP=MCHEP*CCAPH WRITE(*,*) WRITE(*,35) 'PHYSICAL ROUNDOFF MAGNIFIES TO:',UPHYC WRITE(*,35) 'CANONICAL ROUNDOFF MAGNIFIES TO:',UCANP 35 FORMAT(A45,2X,E9.2) C OPEN(CH0,FILE=OFLC) C C**** WRITE CONFPACK HEADING ON LISTING FILE C CALL WRHEAD(5,CH0) C WRITE(CH0,*) WRITE(CH0,40) 40 FORMAT(T4,'MAP',T18,'ESTIMATED EVALUATION',T42,'ESTIMATED MAXIMUM' +,/,T18,'CONDITION NUMBER',T42,'ROUNDOFF ERROR *',/) C IF (NINFD .GT. 0) THEN CHPC='**' ELSE CHPC=' ' ENDIF IF (NZERD .GT. 0) THEN CHCP='**' ELSE CHCP=' ' ENDIF C WRITE(CH0,50) CPHCA,CHPC,UPHYC WRITE(CH0,60) CCAPH,CHCP,UCANP 50 FORMAT('PHY --> CAN',T20,E11.3,A2,T44,E11.3,/) 60 FORMAT('CAN --> PHY',T20,E11.3,A2,T44,E11.3,/) WRITE(CH0,*) '* BASED ON UNIT ROUNDOFF IN DOMAIN OF MAP' IF (NINFD.GT.0 .OR. NZERD.GT.0) THEN WRITE(CH0,*)'** CONDITION NUMBER DEPENDS ON UNIT ROUNDOFF,U:' IF (NINFD .GT. 0) THEN WRITE(CH0,70) COPHC,EXPHC ENDIF IF (NZERD .GT. 0) THEN WRITE(CH0,80) COCAP,EXCAP ENDIF ENDIF 70 FORMAT(T4,'PHY --> CAN : CONDTN NO = ',E11.3,'*U**',E11.3) 80 FORMAT(T4,'CAN --> PHY : CONDTN NO = ',E11.3,'*U**',E11.3) C PI=4E+0*ATAN(1E+0) WRITE(CH0,90) 'END PT.','PARENT','ARGUMENT/PI' 90 FORMAT(//,A,T10,A,T18,A) DO 100 I=1,TNSUA WRITE(CH0,110) I,IGEOM(PARNT+I-1),RGEOM(VTARG+I-1)/PI 100 CONTINUE 110 FORMAT(I3,T10,I3,T18,E16.8) C WRITE(CH0,120) 'SUBARC','% PHYSICAL','% CIRCLE' 120 FORMAT(/,A,T10,A,T29,A) DO 130 I=1,TNSUA ANGSP=RGEOM(VTARG+I)-RGEOM(VTARG+I-1) WRITE(CH0,140) I,ARCLN(I)/TOTLN,ANGSP/2E+0/PI 130 CONTINUE 140 FORMAT(I4,T10,E14.7,T29,E14.7) C WRITE(CH0,150) 'SUB','ACHIEVED','CROWDING','ARC','RESOLUTION', +'FACTOR' 150 FORMAT(/,A,T7,A,T19,A,/,A,T7,A,T19,A) RESMN=OFLOW DO 160 I=1,TNSUA ANGSP=RGEOM(VTARG+I)-RGEOM(VTARG+I-1) IF (ANGSP.EQ.0E+0) THEN CR=OFLOW LA=0E+0 ELSE CR=2E+0*PI*ARCLN(I)/ABS(ANGSP)/TOTLN IF (RSNPH(ERARC+I-1).EQ.0E+0) THEN LA=OFLOW ELSE LA=ABS(ANGSP)/(2E+0*RSNPH(ERARC+I-1)) ENDIF ENDIF IF (LA .LT. RESMN) THEN RESMN=LA IMNLA=I ENDIF WRITE(CH0,170) I,LA,CR 160 CONTINUE 170 FORMAT(I2,T4,2E12.3) C WRITE(CH0,180) RESMN,IMNLA 180 FORMAT(/,'MINIMUM SUBARC RESOLUTION IS ',E11.3,' ON SUBARC ',I2) WRITE(*,*) WRITE(*,35) 'MINIMUM SUBARC RESOLUTION:',RESMN C WRITE(CH0,*) IF (.NOT.MAP11 .OR. RESMN.LT.CRRES) THEN C C**** MESSAGE TO STANDARD OUTPUT C WRITE(*,185) '*** W A R N I N G ***' 185 FORMAT(//,T20,A) IF (RESMN.LT.CRRES) THEN WRITE(*,5) 'THE ABOVE RESOLUTION IS TOO SMALL:' ENDIF IF (.NOT.MAP11) THEN WRITE(*,5) 'BCF DERIVATIVE CHANGES SIGN:' ENDIF C WRITE(CH0,*) ' *** W A R N I N G ***' WRITE(CH0,*) I=0 C IF (RESMN.LT.1) THEN I=I+1 WRITE(CH0,190) I,'. THE ABOVE SUBARC RESOLUTION MEANS THAT IT + WILL BE PRACTICALLY' WRITE(CH0,*) ' IMPOSSIBLE FOR THE INVERSE MAP TO DISCRIMINA +TE CORRECTLY' WRITE(CH0,200) ' BETWEEN NEIGHBOURING POINTS NEAR SUB ARC ' +,IMNLA ELSE IF (RESMN.LT.CRRES) THEN I=I+1 WRITE(CH0,190) I,'. THE ABOVE SUBARC RESOLUTION MEANS THAT THE + INVERSE MAP MAY NOT' WRITE(CH0,*) ' BE ABLE TO RELIABLY DISCRIMINATE CORRECTLY B +ETWEEN' WRITE(CH0,200) ' NEIGHBOURING POINTS NEAR ARC ',IMNLA ENDIF 190 FORMAT(/,I1,A) 200 FORMAT(A,I3) C IF (NCRVS .GT. 0) THEN I=I+1 WRITE(CH0,190) I,'. THERE IS A COMPLETE REVERSAL OF DIRECTION + ON THE FOLLOWING SUB ARCS:' WRITE(CH0,'(T10,I3)') (ICRVS(J),J=1,NCRVS) ENDIF C IF (NPRVS .GT. 0) THEN I=I+1 WRITE(CH0,190) I,'. THERE IS A REVERSAL OF DIRECTION WITHIN T +HE FOLLOWING SUB ARCS:' WRITE(CH0,'(T10,I3)') (IPRVS(J),J=1,NPRVS) WRITE(CH0,*) ' THE CORRESPONDING MINIMUM VALUES OF THE BOUN +DARY CORRESPONDENCE' WRITE(CH0,*) ' DERIVATIVE ARE:' WRITE(CH0,'(T10,E9.2)') (BCDMN(J),J=1,NPRVS) ENDIF ENDIF CLOSE(CH0) 999 CONTINUE C C**** WRITE CLOSING MESSAGE TO STANDARD OUTPUT CHANNEL AND LISTING FILE C CALL WRTAIL(5,0,IER) CALL WRTAIL(5,CH0,IER) C END SUBROUTINE CPJAC3(NARCS,NQPTS,INDEG,DGPOL,JACIN,ACOEF,BCOEF,DIAG, +SDIAG,TNSUA,LOSUB,HISUB,JATYP,PARNT,MIDPT,HALEN,COLPR,ZCOLL,LNSEG, +LOOLD,HIOLD,EPS,IER,INIBT) INTEGER NARCS,NQPTS,INDEG,TNSUA,IER INTEGER DGPOL(*),LOSUB(*),HISUB(*),JATYP(*),PARNT(*) REAL EPS,JACIN(*),ACOEF(*),BCOEF(*),DIAG(*),SDIAG(*), +MIDPT(*),HALEN(*),COLPR(*),LOOLD(*),HIOLD(*) COMPLEX ZCOLL(*) LOGICAL LNSEG(*),INIBT C C**** TO MAKE THE INITIAL ASSIGNMENT OF THE COLLOCATION PARAMETERS C**** (STORED IN COLPR), THE COLLOCATION POINTS ON THE PHYSICAL C**** BOUNDARY (STORED IN ZCOLL) AND THE ARRAYS LOSUB AND HISUB C**** NEEDED TO ACCESS THIS DATA CORRECTLY. ALSO TO SET UP THE C**** ARRAYS C**** JATYP - THE JACOBI INDEX TYPE OF EACH SUBARC C**** PARNT - THE ORIGINAL PARENT ARC OF EACH SUBARC C**** MIDPT - THE GLOBAL PARAMETRIC MIDPOINT OF EACH SUBARC C**** HALEN - THE GLOBAL PARAMETRIC HALF-LENGTH OF EACH SUBARC C**** DGPOL - THE INITIAL POLYNOMIAL DEGREE ON EACH SUBARC C**** LNSEG - THE INITIAL LINE SEGMENT TYPE OF EACH SUBARC. C C**** IER=0 - NORMAL EXIT C**** IER=7 - FAILURE TO CONVERGE IN EIGENVALUE ROUTINE IMTQLH C C LOCAL VARIABLES C INTEGER D,D1,FIRST,I,J,K,K1,K2,P,PREV,IFAIL REAL S,TC,ALFA,BETA,MED,MDNBT COMPLEX PARFUN EXTERNAL PARFUN,IMTQLH,MDNBT C TNSUA=2*NARCS DO 10 I=1,NARCS BETA=JACIN(I) IF (I .EQ. NARCS) THEN ALFA=JACIN(1) ELSE ALFA=JACIN(I+1) ENDIF IF (INIBT) THEN MED=MDNBT(ALFA,BETA) ELSE MED=0E+0 ENDIF J=2*I-1 MIDPT(J)=5E-1*(MED-1E+0) HALEN(J)=5E-1*(MED+1E+0) PARNT(J)=I JATYP(J)=I J=J+1 MIDPT(J)=5E-1*(MED+1E+0) HALEN(J)=5E-1*(1E+0-MED) PARNT(J)=I IF (I .EQ. NARCS) THEN JATYP(J)=-1 ELSE JATYP(J)=-I-1 ENDIF 10 CONTINUE C DO 11 I=NARCS,1,-1 J=2*I LNSEG(J)=LNSEG(I) LNSEG(J-1)=LNSEG(I) 11 CONTINUE C DO 12 I=1,TNSUA DGPOL(I)=INDEG 12 CONTINUE C LOSUB(1)=1 HISUB(1)=1+DGPOL(1) DO 20 I=2,TNSUA LOSUB(I)=HISUB(I-1)+1 HISUB(I)=LOSUB(I)+DGPOL(I) 20 CONTINUE C DO 25 I=1,TNSUA LOOLD(I)=0 HIOLD(I)=-1 25 CONTINUE C DO 50 I=1,TNSUA J=JATYP(I) P=PARNT(I) D=DGPOL(I) D1=D+1 IF (J .GT. 0) THEN S=1E+0 ELSE S=-1E+0 J=-J ENDIF PREV=(J-1)*NQPTS FIRST=LOSUB(I) DO 30 K=1,D1 K1=PREV+K DIAG(K)=BCOEF(K1) IF (K .EQ. 1) THEN SDIAG(K)=0E+0 ELSE SDIAG(K)=ACOEF(K1-1) ENDIF 30 CONTINUE CALL IMTQLH(D1,DIAG,SDIAG,IFAIL) IF (IFAIL .GT. 0) THEN IER=7 RETURN ENDIF DO 40 K=1,D1 TC=S*DIAG(K) K2=FIRST+K-1 COLPR(K2)=TC TC=MIDPT(I)+HALEN(I)*TC ZCOLL(K2)=PARFUN(P,CMPLX(TC)) 40 CONTINUE 50 CONTINUE C C NORMAL EXIT C IER=0 C END COMPLEX FUNCTION CJACSU(X,N,A,B,H,CO) INTEGER N REAL A(*),B(*),H,CO(*) COMPLEX X C ..TO CALCULATE SUMMATION{CO(K+1)*P(K,X)},K=0(1)N, WHERE P(K,X) C ..DENOTES THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE K C ..EVALUATED AT X, ARRAY CO STORES A GIVEN SET OF COEFFICIENTS, C ..ARRAYS A,B STORE THE COEFFICIENTS IN THE THREE-TERM C ..RECURRENCE FORMULA FOR THE JACOBI POLYNOMIALS (SEE ASONJ7) C ..AND H IS THE SQUARED 2-NORM OF UNITY. COMPLEX PREV,CURR,NEXT INTEGER K C IF (N .EQ. 0) THEN CJACSU=CO(1)/SQRT(H) ELSE IF (N .GT. 0) THEN PREV=CO(N+1) CURR=CO(N)+(X-B(N))*PREV/A(N) DO 10 K=N-2,0,-1 NEXT=CO(K+1)+(X-B(K+1))*CURR/A(K+1)-A(K+1)*PREV/A(K+2) PREV=CURR CURR=NEXT 10 CONTINUE CJACSU=CURR/SQRT(H) ELSE CJACSU=(0E+0,0E+0) ENDIF C END INTEGER FUNCTION CRITCO(N,POSCO) INTEGER N REAL POSCO(*) C C GIVEN THE NON-NEGATIVE NUMBERS POSCO(I), I=1,...,N, TO FIND THE C INDEX CRITCO SUCH THAT POSCO(CRITCO) > 1 AND POSCO(I) <=1 FOR C I=CRITCO+1,...,N. IN CASE POSCO(I) <=1 FOR ALL I=1,...,N THEN C POSCO=0. C C**** LOCAL VARIABLE C INTEGER I C I=N 10 CONTINUE IF (I.EQ.0) THEN CRITCO=0 ELSE IF (POSCO(I) .GT. 1E+0) THEN CRITCO=I ELSE I=I-1 GOTO 10 ENDIF C END SUBROUTINE CSCAL3(COLSC,NQPTS,NJIND,ACOEF,BCOEF,H0VAL,QUPTS,QUWTS, +JACIN,MU,TT,QQ) INTEGER NQPTS,NJIND REAL COLSC(*),ACOEF(*),BCOEF(*),H0VAL(*),QUPTS(*),QUWTS(*), +JACIN(*),MU(*),TT(*),QQ(*) C C TO SET UP THE A PRIORI COLUMN SCALE FACTORS USING TRUNCATED C CHEBYSHEV EXAPANSIONS FOR THE LOGARITHMIC KERNEL, GAUSS-JACOBI C QUADRATURE AND GAUSS-JACOBI TEST POINTS. C C LOCAL VARIABLES C INTEGER I,J,J1,JI,K,K1,KMAX,LO,LO1,M REAL BETA,ROOTH,P0SCL,X,MAXVL,SUM1,SUM2 EXTERNAL JAPAR7 C K1=0 MU(1)=-LOG(2E+0) DO 10 I=2,2*NQPTS MU(I)=-2E+0/REAL(I-1) 10 CONTINUE C DO 80 JI=1,NJIND BETA=JACIN(JI) ROOTH=SQRT(H0VAL(JI)) P0SCL=1E+0/ROOTH LO=(JI-1)*NQPTS LO1=LO+1 C DO 30 J=1,NQPTS X=QUPTS(LO+J) QQ(1+(J-1)*NQPTS)=P0SCL CALL JAPAR7(QQ(1+(J-1)*NQPTS),X,ACOEF(LO1),BCOEF(LO1), + NQPTS-1) TT(J)=1E+0 TT(J+NQPTS)=X DO 20 K=3,2*NQPTS TT(J+(K-1)*NQPTS)=2E+0*X*TT(J+(K-2)*NQPTS)- + TT(J+(K-3)*NQPTS) 20 CONTINUE 30 CONTINUE C DO 70 K=1,NQPTS MAXVL=0E+0 DO 60 I=1,NQPTS SUM2=0E+0 KMAX=2*NQPTS+1-K DO 50 M=K,KMAX SUM1=0E+0 DO 40 J=1,NQPTS J1=LO+J SUM1=SUM1+QUWTS(J1)*QQ(K+(J-1)*NQPTS)* + TT(J+(M-1)*NQPTS) 40 CONTINUE SUM2=SUM2+MU(M)*TT(I+(M-1)*NQPTS)*SUM1 50 CONTINUE MAXVL=MAX(MAXVL,ABS(SUM2)) 60 CONTINUE K1=K1+1 COLSC(K1)=1E+0/MAXVL 70 CONTINUE C 80 CONTINUE C END SUBROUTINE DEJAC7(ZZ,NZZ,BETA,TAU,MAXDG,NQUAD,ACOEF,BCOEF,H0VAL, +REMND,CSCAL,TOL,MAXRM,IER) INTEGER MAXDG,NQUAD,NZZ,IER REAL BETA,TAU,H0VAL,TOL,MAXRM REAL ACOEF(*),BCOEF(*),CSCAL(*) COMPLEX ZZ(*),REMND(*) C C WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN C USING AN NQUAD-POINT GAUSS-JACOBI RULE TO ESTIMATE THE INTEGRALS C C INTEGRAL [(1+X)**BETA*P(X,I)*LOG(ZZ(J)-X)*dX], I=0,MAXDG C -1<=X<=TAU J=1,NZZ C C WHERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I C ASSOCIATED WITH THE WEIGHT (1+X)**BETA. THE REMAINDER C CORRESPONDING TO P(.,I) AND ZZ(J) IS STORED IN C REMND(I+J+MAXDG*(J-1)), I=0,MAXDG, J=1,NZZ. THIS ROUTINE USES C THE SIMPLEST POSSIBLE ESTIMATES; I.E. THE LEADING TERM ONLY IN C THE ASYMPTOTIC EXPANSION AND THE WATSON-DOETSCH ESTIMATE FOR ANY C INTEGRALS. C C THE PURPOSE OF THIS ROUTINE IS THEN TO DETERMINE A VALUE FOR TAU C SUCH THAT C C ABS( REAL(REMND(I)) )*CSCAL(I) < TOL , I=1,NZZ*(MAXDG+1) C C AND THAT, IF POSSIBLE, C C 0.5*TOL <= ABS( REAL(REMND(I)) )*CSCAL(I) < TOL C C FOR AT LEAST ONE VALUE OF I. C C IER=0 - NORMAL EXIT C IER=12- LOCAL PARAMETER NC NEEDS INCREASING TO AT LEAST NZZ C (THIS ERROR CAN'T ARISE IN THE PRESENT VERSION, SINCE C NZZ IS FIXED AT 2) C IER=13- LOCAL PARAMETER NR NEEDS INCREASING TO AT LEAST MAXDG C (AT PRESENT MAXDG=NQPTS-1) C IER=14- A JACOBI INDEX MAY BE LARGE ENOUGH TO CAUSE OVERFLOW IN C THE GAMMA FUNCTION (AN ANGLE ON THE PHYSICAL BOUNDARY C MUST BE LESS THAN ABOUT 6 DEGREES) C C LOCAL VARIABLES.. C INTEGER I,J,K,LIM,NC,NR REAL S,KK,GAMMA,LGGAM,SUM1,RI,TURI,RN,OFLOW,P0SCL,TUK,LOWER, +UPPER,R1MACH,TERM,HTOL COMPLEX XI,Z1,XI1,FF,PRE,CUR,NXT PARAMETER (NC=8,NR=30) COMPLEX GG(NC),CONST(NR,NC) EXTERNAL GAMMA,LGGAM,R1MACH C IF (NZZ .GT. NC ) THEN IER=12 RETURN ENDIF C IF (MAXDG .GE. NR) THEN IER=13 RETURN ENDIF C S=BETA+4E+0 IF (S .GT. 2E+1) THEN C C TEST FOR POSSIBLE OVERFLOW IN GAMMA FUNCTION C OFLOW=LOG(R1MACH(2)) KK=LGGAM(S) IF (KK.GT.OFLOW) THEN IER=14 RETURN ELSE KK=EXP(-KK) ENDIF ELSE KK=1E+0/GAMMA(S) ENDIF C C FIRST WE COMPUTE THE FACTORS WHICH ARE INDEPENDENT OF TAU C S=S-1E+0 KK=4E+0**S*KK*GAMMA(BETA+2E+0)/(S-1E+0) SUM1=BETA+1E+0 DO 100 I=2,NQUAD RI=REAL(I) TURI=2E+0*RI KK=KK*16E+0*(RI+BETA)/(TURI+SUM1) KK=KK*RI/(TURI+BETA) KK=KK*(RI+BETA)/(TURI+BETA) KK=KK*RI/(TURI+BETA-1E+0) 100 CONTINUE RN=REAL(NQUAD) TUK=2E+0*RN+SUM1 KK=-KK/TUK/2E+0 C DO 125 I=1,NZZ GG(I)=(1E+0+ZZ(I))**BETA*KK 125 CONTINUE C C NOW GIVE THE JACOBI POLYNOMIALS THE SCALING CORRESPONDING TO C [-1,1] AS STANDARD INTERVAL C P0SCL=1E+0/SQRT(H0VAL) DO 225 J=1,NZZ PRE=CMPLX(P0SCL) CUR=PRE CONST(1,J)=CUR*GG(J)*CSCAL(1) IF (MAXDG .GE. 1) THEN CUR=(ZZ(J)-BCOEF(1))*PRE/ACOEF(1) CONST(2,J)=CUR*GG(J)*CSCAL(2) DO 200 I=2,MAXDG NXT=((ZZ(J)-BCOEF(I))*CUR-ACOEF(I-1)*PRE)/ACOEF(I) PRE=CUR CUR=NXT CONST(I+1,J)=CUR*GG(J)*CSCAL(I+1) 200 CONTINUE ENDIF 225 CONTINUE C C NOW COME THE FACTORS DEPENDENT ON TAU C TAU=1E+0 LOWER=-1E+0 UPPER=1E+0 LIM=NZZ*(MAXDG+1) C 250 CONTINUE C HTOL=5E-1*TOL K=0 DO 325 J=1,NZZ XI=(2E+0*ZZ(J)+1E+0-TAU)/(1E+0+TAU) Z1=SQRT(XI*XI-1E+0) XI1=XI+Z1 IF (ABS(XI1) .LT. 1E+0) THEN XI1=XI-Z1 ENDIF FF=XI1**(-TUK-1E+0)*(XI1*XI1-1E+0)*(1E+0+TAU)*5E-1 DO 300 I=0,MAXDG K=K+1 REMND(K)=CONST(I+1,J)*FF 300 CONTINUE 325 CONTINUE C MAXRM=0E+0 DO 600 I=1,LIM TERM=ABS(REAL(REMND(I))) MAXRM=MAX(MAXRM,TERM) 600 CONTINUE C IF (MAXRM .LT. TOL) THEN C C ACCURACY IS ACHIEVED, BUT MAYBE TAU COULD BE INCREASED. C IF (MAXRM .LT. HTOL) THEN C C TAU NEEDS INCREASING, BUT THIS IS ONLY POSSIBLE IF TAU<1. C IF (TAU .LT. 1E+0) THEN LOWER=TAU TAU=5E-1*(LOWER+UPPER) GOTO 250 ENDIF ENDIF ELSE C C ACCURACY NOT ACHIEVED AND TAU NEEDS DECREASING. C IF (TAU .EQ. 1E+0) THEN TOL=HTOL ENDIF UPPER=TAU TAU=5E-1*(LOWER+UPPER) GOTO 250 ENDIF C C NORMAL TERMINATION C IER=0 C END SUBROUTINE DELEG7(ZZ,NZZ,BETA,TAU1,TAU2,T1FXD,MAXDG,NQUAD,ACOEF, +BCOEF,H0VAL,REMND,CSCAL,TOL,MAXRM,IER) INTEGER MAXDG,NQUAD,IER,NZZ REAL BETA,TAU1,TAU2,H0VAL,TOL,MAXRM REAL ACOEF(*),BCOEF(*),CSCAL(*) LOGICAL T1FXD COMPLEX ZZ(*),REMND(*) C C WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN C USING AN NQUAD-POINT GAUSS-LEGENDRE RULE TO ESTIMATE THE INTEGRALS C C INTEGRAL [(1+X)**BETA*P(X,I)*LOG(ZZ(J)-X)*dX], I=0,1,...,MAXDG C TAU1<=X<=TAU2 J=1,2 C C WHERE P(.,I) IS THE ORTHONORMAL JACOBI POLYNOMIAL OF DEGREE I C ASSOCIATED WITH THE WEIGHT (1+X)**BETA AND -1-1) C IF (T1FXD .AND. TAU2 .LT. 1E+0) THEN LOWER=TAU2 TAU2=5E-1*(LOWER+UPPER) GOTO 250 ELSE IF (.NOT.T1FXD .AND. TAU1 .GT. -1E+0) THEN UPPER=TAU1 TAU1=5E-1*(LOWER+UPPER) GOTO 250 ENDIF ENDIF ELSE C C ACCURACY NOT ACHIEVED AND TAU2 NEEDS DECREASING OR TAU1 NEEDS C INCREASING. C IF (FIRST) THEN TOL=HTOL FIRST=.FALSE. ENDIF IF (T1FXD) THEN UPPER=TAU2 TAU2=5E-1*(LOWER+UPPER) ELSE LOWER=TAU1 TAU1=5E-1*(LOWER+UPPER) ENDIF GOTO 250 ENDIF C C NORMAL TERMINATION C IER=0 C END SUBROUTINE DEPPJ8(BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL, +SOLUN,TOL,MAXRM,NINTS,DELTA,IER) INTEGER NQUAD,IER,DGPOL,NINTS REAL BETA,TAU,TOL,MAXRM,ACOEF(*),BCOEF(*),H0VAL,DELTA COMPLEX SOLUN(*) C C WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN C USING AN NQUAD-POINT GAUSS-JACOBI RULE TO ESTIMATE THE INTEGRALS C C INTEGRAL [(1+X)**BETA*FNPHI(X)*LOG(ZZ-X)*dX] C -1<=X<=TAU C C WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY C CORRESPONDENCE DERIVATIVE / JACOBI WEIGHT QUOTIENT. C C ZZ IS ANY POINT ON A "DELTA-CONTOUR" IN THE UPPER HALF PLANE, C THIS CONTOUR BEING DEFINED BY THE PARAMETE DELTA. TEST VALUES C FOR ZZ ARE ASSIGNED IN THE BODY OF THE ROUTINE AND THE LOCAL ARRAY C XIVAL IS USED FOR STORING THESE TEST VALUES. C C THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO C DEFINE FNPHI AND ARE PASSED TO DEPPJ9 FOR THIS PURPOSE. C C MAXRM RECORDS THE MAXIMUM OF THE ABSOLUTE VALUES OF THE REMAINDER C ESTIMATES ASSOCIATED WITH THE DELTA-CONTOUR. C C ON ENTRY NINTS IS THE NUMBER OF INITIAL TEST INTERVALS TO BE USED C ON THE UPPER HALF DELTA-CONTOUR; ON EXIT NINTS GIVES THE FINAL C NUMBEROF INTERVALS REQUIRED FOR CONVERGENCE TO TAU. C C THE PURPOSE OF THIS ROUTINE IS TO DETERMINE A VALUE FOR TAU C SUCH THAT C C MAXRM < TOL C C AND THAT, IF POSSIBLE, C C 0.5*TOL <= MAXRM < TOL. C C IER=0 - NORMAL EXIT C IER=25- THE LOCAL ARRAY BOUND PARAMETER MNXI NEEDS INCREASING. C C C LOCAL VARIABLES C INTEGER NXI,MNXI REAL SS,PI,SS1,SS2,SS3,DP,LEN,TAUI,SINC,STAR,RXI,DT,HDP COMPLEX XI PARAMETER (MNXI=100) COMPLEX XIVAL(MNXI) EXTERNAL DEPPJ9 C C INITIALISATION C PI=4E+0*ATAN(1E+0) TAUI=1E+0 DP=DELTA*PI DT=DELTA*2E+0 HDP=5E-1*DP LEN=2E+0+DP SINC=LEN/NINTS STAR=SINC C C START OF LOOP FOR INTERVAL HALVING C 10 CONTINUE NXI=0 SS1=-STAR+SINC SS2=LEN DO 30 SS=SS1,SS2,SINC IF (SS .LT. HDP) THEN SS3=SS/DELTA XI=1E+0+DELTA*CMPLX(COS(SS3),SIN(SS3)) ELSE IF (SS .GT. 2E+0+HDP) THEN SS3=(SS-2E+0-HDP)/DELTA XI=-1E+0+DELTA*(0E+0,1E+0)*CMPLX(COS(SS3),SIN(SS3)) ELSE SS3=SS-HDP XI=CMPLX(1E+0-SS3,DELTA) ENDIF RXI=REAL(XI) IF (RXI .LT. (TAUI+DT)) THEN NXI=NXI+1 IF (NXI .GT. MNXI) THEN IER=25 RETURN ENDIF XIVAL(NXI)=XI ENDIF 30 CONTINUE C IF (NXI .EQ. 0) THEN SINC=STAR STAR=5E-1*STAR NINTS=NINTS*2 GOTO 10 ENDIF C TAU=TAUI CALL DEPPJ9(XIVAL,NXI,BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL, + SOLUN(1),TOL,MAXRM,IER) IF (IER .GT. 0) THEN RETURN ENDIF C IF (TAU .NE. TAUI) THEN TAUI=TAU SINC=STAR STAR=5E-1*STAR NINTS=NINTS*2 GOTO 10 ENDIF C IER=0 C END SUBROUTINE DEPPJ9(ZZ,NZZ,BETA,TAU,NQUAD,DGPOL,ACOEF,BCOEF,H0VAL, +SOLUN,TOL,MAXRM,IER) INTEGER NQUAD,NZZ,IER,DGPOL REAL BETA,TAU,TOL,MAXRM,ACOEF(*),BCOEF(*),H0VAL COMPLEX SOLUN(*),ZZ(*) C C WE COMPUTE THE DONALDSON-ELLIOTT ESTIMATES FOR THE REMAINDERS IN C USING AN NQUAD-POINT GAUSS-JACOBI RULE TO ESTIMATE THE INTEGRALS C C INTEGRAL [(1+X)**BETA*FNPHI(X)*LOG(ZZ(I)-X)*dX], I=1,NZZ C -1<=X<=TAU C C WHERE FNPHI IS A POLYNOMIAL APPROXIMATION TO THE BOUNDARY C CORRESPONDENCE DERIVATIVE - JACOBI WEIGHT QUOTIENT AND ZZ IS C A GIVEN ARRAY OF POINTS. C C THE PARAMETERS DGPOL,ACOEF,BCOEF,H0VAL AND SOLUN ARE USED TO C DEFINE FNPHI. C C THE MAXIMUM ABSOLUTE VALUE OF ALL THE REMAINDERS CORRESPONDING TO C ZZ(I) , I=1,NZZ, IS STORED IN MAXRM. C C THIS ROUTINE USES THE SIMPLEST POSSIBLE ESTIMATES; I.E. THE C LEADING TERM ONLY IN THE ASYMPTOTIC EXPANSION AND THE WATSON- C DOETSCH ESTIMATE FOR ANY INTEGRALS. C C THE PURPOSE OF THIS ROUTINE IS THEN TO DETERMINE A VALUE FOR TAU C SUCH THAT C C MAXRM < TOL C C AND THAT, IF POSSIBLE, C C 0.5*TOL <= MAXRM < TOL. C C IER=0 - NORMAL EXIT C IER=26- TOO MANY TEST POINTS ON DELTA CONTOUR; INCREASE C PARAMETER MAXNZ BELOW C IER=14- BETA MAY CAUSE OVERFLOW IN GAMMA FUNCTION; AN ANGLE C ON THE BOUNDARY IS TOO SMALL C C LOCAL VARIABLES C INTEGER I,MAXNZ REAL S,KK,GAMMA,SUM1,RI,TURI,RN,TUK,LOWER,UPPER,TERM,HTOL,TAUI, +OFLOW,LGGAM,R1MACH COMPLEX XI,Z1,XI1,FF,FNPHI,CCJACS,REMND PARAMETER (MAXNZ=200) COMPLEX GG(MAXNZ) EXTERNAL GAMMA,CCJACS,LGGAM,R1MACH C IF (NZZ .GT. MAXNZ ) THEN C C SOME LOCAL ARRAY BOUNDS MUST BE INCREASED C IER=26 RETURN ENDIF C S=BETA+4E+0 IF (S .GT. 2E+1) THEN C C TEST FOR POSSIBLE OVERFLOW IN GAMMA FUNCTION C OFLOW=LOG(R1MACH(2)) KK=LGGAM(S) IF (KK.GT.OFLOW) THEN IER=14 RETURN ELSE KK=EXP(-KK) ENDIF ELSE KK=1E+0/GAMMA(S) ENDIF C C FIRST WE COMPUTE THE FACTORS WHICH ARE INDEPENDENT OF TAU C S=S-1E+0 KK=4E+0**S*KK*GAMMA(BETA+2E+0)/(S-1E+0) SUM1=BETA+1E+0 DO 100 I=2,NQUAD RI=REAL(I) TURI=2E+0*RI KK=KK*16E+0*(RI+BETA)/(TURI+SUM1) KK=KK*RI/(TURI+BETA) KK=KK*(RI+BETA)/(TURI+BETA) KK=KK*RI/(TURI+BETA-1E+0) 100 CONTINUE RN=REAL(NQUAD) TUK=2E+0*RN+SUM1 KK=-KK/TUK/2E+0 C DO 125 I=1,NZZ FNPHI=CCJACS(ZZ(I),DGPOL,ACOEF,BCOEF,H0VAL,SOLUN) GG(I)=(1E+0+ZZ(I))**BETA*KK*FNPHI 125 CONTINUE C C NOW COME THE FACTORS DEPENDENT ON TAU C LOWER=-1E+0 UPPER=TAU TAUI=TAU C 250 CONTINUE C MAXRM=0E+0 HTOL=5E-1*TOL DO 325 I=1,NZZ XI=(2E+0*ZZ(I)+1E+0-TAU)/(1E+0+TAU) Z1=SQRT(XI*XI-1E+0) XI1=XI+Z1 IF (ABS(XI1) .LT. 1E+0) THEN XI1=XI-Z1 ENDIF FF=XI1**(-TUK-1E+0)*(XI1*XI1-1E+0)*(1E+0+TAU)*5E-1 REMND=GG(I)*FF TERM=ABS(REMND) MAXRM=MAX(MAXRM,TERM) 325 CONTINUE C IF (MAXRM .LT. TOL) THEN C C ACCURACY IS ACHIEVED, BUT MAYBE TAU COULD BE INCREASED. C IF (MAXRM .LT. HTOL) THEN C C TAU NEEDS INCREASING, BUT THIS IS ONLY POSSIBLE IF TAUTAU1I) C IF (T1FXD .AND. TAU2 .LT. TAU2I) THEN LOWER=TAU2 TAU2=5E-1*(LOWER+UPPER) GOTO 250 ELSE IF (.NOT.T1FXD .AND. TAU1 .GT. TAU1I) THEN UPPER=TAU1 TAU1=5E-1*(LOWER+UPPER) GOTO 250 ENDIF ENDIF ELSE C C ACCURACY NOT ACHIEVED AND TAU2 NEEDS DECREASING OR TAU1 NEEDS C INCREASING. C IF (FIRST) THEN TOL=HTOL FIRST=.FALSE. ENDIF IF (T1FXD) THEN UPPER=TAU2 TAU2=5E-1*(LOWER+UPPER) ELSE LOWER=TAU1 TAU1=5E-1*(LOWER+UPPER) ENDIF GOTO 250 ENDIF C C NORMAL EXIT C IER=0 C END SUBROUTINE DIAGN4(CCAPH,COCAP,COPHC,CPHCA,EXCAP,EXPHC,ICRVS,IER, +IPRVS,NASYM,NCRVS,NINFD,NPRVS,NZERD,ARCLN,ASYMP,BCDMN,CORXX, +TOTLN,VTARG,MAP11,DGPOL,JATYP,LOSUB,NARCS,NQPTS,NXINT,OUCH0,OUCH1, +PARNT,TNSUA,A1COF,ACOEF,B1COF,BCFSN,BCOEF,H0VAL,H1VAL,HALEN,JACIN, +MIDPT,SOLUN,QUPTS,QUWTS) INTEGER IER,NARCS,NASYM,NCRVS,NINFD,NPRVS,NQPTS,NXINT,NZERD,OUCH0, +OUCH1,TNSUA INTEGER DGPOL(*),ICRVS(*),IPRVS(*),JATYP(*),LOSUB(*),PARNT(*) REAL CCAPH,COCAP,COPHC,CPHCA,EXCAP,EXPHC,TOTLN REAL A1COF(*),ACOEF(*),ARCLN(*),ASYMP(*),B1COF(*), +BCDMN(*),BCFSN(*),BCOEF(*),CORXX(*),JACIN(*),MIDPT(*),H0VAL(*), +H1VAL(*),HALEN(*),SOLUN(*),VTARG(*),QUPTS(*),QUWTS(*) LOGICAL MAP11 C C IER=0 - NORMAL EXIT C IER=50 - LOCAL PARAMETER MXCOF MUST BE >= NQPTS. C IER=51 - NON-ANALYTIC ARC DETECTED C C**** LOCAL VARIABLES C INTEGER AJT,DG,I,I1,IA,JT,K,LOD,LOM,MININ,MXCOF,NINTS,PT,QP REAL AL,ATOL,BT,CC,COF,D,DSDT,H0,HH,HL,JACSUM,MD,MCHEP,MPT, +PHI,R1MACH,RTOL,SEND,SINC,SJT,SS,SUM,TERM,TINC,TT,TUPI,X,XX,YMAX, +YMIN,YY COMPLEX PARFUN,T1,T2 COMMON /DSDTDA/PT,MD,HL PARAMETER (MININ=20,MXCOF=32,QP=4) REAL JACOF(MXCOF) EXTERNAL DSDT,JACSUM,PARFUN,R1MACH C C INITIALISE SOME CONSTANTS C TUPI=8E+0*ATAN(1E+0) MCHEP=R1MACH(4) RTOL=1E+1*MCHEP ATOL=1E+2*MCHEP NCRVS=0 NPRVS=0 CCAPH=0E+0 CPHCA=0E+0 MAP11=.TRUE. YMAX=R1MACH(2) NASYM=0 C C START TO COMPUTE THE ARC LENGTHS OF EACH SUBARC (ARCLN) AND THE C TOTAL LENGTH (TOTLN) OF THE BOUNDARY C TOTLN=0E+0 DO 10 IA=1,TNSUA PT=PARNT(IA) MD=MIDPT(IA) HL=HALEN(IA) T1=CMPLX(MD+HL) T2=CMPLX(MD-HL) C C**** COMPOSITE QP-PANEL GAUSS-LEGENDRE ESTIMATE FOR ARCLN(IA) C HH=1E+0/QP SUM=0E+0 DO 6 K=1,QP MPT=-1E+0+(2E+0*K-1E+0)*HH DO 3 I=1,NQPTS X=MPT+HH*QUPTS(I) SUM=SUM+QUWTS(I)*DSDT(X) 3 CONTINUE 6 CONTINUE ARCLN(IA)=HH*SUM TOTLN=TOTLN+ARCLN(IA) 10 CONTINUE C C TEST FOR COMPLETE REVERSAL OF DIRECTION OF A BOUNDARY SUBARC ON C THE UNIT DISC. C DO 20 IA=1,TNSUA IF (VTARG(IA+1) .LT. VTARG(IA)) THEN NCRVS=NCRVS+1 ICRVS(NCRVS)=IA MAP11=.FALSE. ENDIF 20 CONTINUE C C COMPUTE THE NUMBERS *NINFD* (*NZERD*) OF POINTS WHERE THE C DERIVATIVE OF THE MAP PHYSICAL --> CANONICAL IS RESPECTIVELY C INFINITE (ZERO). C NINFD=0 NZERD=0 DO 25 I=1,NARCS IF (JACIN(I) .LT. 0E+0) THEN NINFD=NINFD+1 ELSE IF (JACIN(I) .GT. 0E+0) THEN NZERD=NZERD+1 ENDIF 25 CONTINUE C C NOW START TO EVALUATE THE DIMENSIONLESS BOUNDARY CORRESPONDENCE C DERIVATIVE AT SELECTED VALUES OF DIMENSIONLESS ARC LENGTH; C OUTPUT RESULTS FOR SUBSEQUENT GRAPH PLOTTING IF REQUIRED AND C TEST FOR SIGN CHANGES IN THIS DERIVATIVE. C SS=0E+0 SEND=0E+0 DO 60 IA=1,TNSUA NINTS=MAX(MININ,NINT(ARCLN(IA)*NXINT/TOTLN)) TINC=2E+0/NINTS DG=DGPOL(IA) IF (DG+1 .GT. MXCOF) THEN IER=50 RETURN ENDIF JT=JATYP(IA) AJT=ABS(JT) H0=H0VAL(AJT) BT=JACIN(AJT) AL=1E+0/(1E+0+BT) PT=PARNT(IA) MD=MIDPT(IA) HL=HALEN(IA) LOM=LOSUB(IA) LOD=(AJT-1)*NQPTS+1 IF (JT.GT.0) THEN CC=VTARG(IA)-VTARG(1) ELSE CC=VTARG(IA+1)-VTARG(1) ENDIF DO 30 I=1,DG+1 I1=I+LOM-1 JACOF(I)=SOLUN(I1) 30 CONTINUE SJT=SIGN(1E+0,REAL(JT)) DO 40 I=2,DG+1,2 JACOF(I)=SJT*JACOF(I) 40 CONTINUE TT=-1E+0 D=DSDT(TT) YMIN=YMAX IF (IA .EQ. 1) THEN XX=0E+0 IF (BT .LT. 0E+0) THEN YY=YMAX NASYM=NASYM+1 ASYMP(NASYM)=XX ELSE IF (BT .GT. 0E+0) THEN YY=0E+0 ELSE PHI=JACSUM(TT,DG,ACOEF(LOD),BCOEF(LOD),H0,JACOF) IF (D .EQ. 0E+0) THEN IER=51 RETURN ENDIF YY=TOTLN*PHI/D ENDIF IF (NINFD .EQ. 0E+0) THEN CPHCA=TUPI*ABS(YY)/TOTLN ENDIF IF (NZERD .EQ. 0E+0) THEN IF (YY .EQ. 0E+0) THEN CCAPH=YMAX ELSE CCAPH=TOTLN/TUPI/ABS(YY) ENDIF ENDIF WRITE(OUCH1,902) XX,YY YY=0E+0 WRITE(OUCH0,902) XX,YY CORXX(1)=0E+0 ENDIF C C ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS FOR INFINITE C DERIVATIVE CASES. C IF (BT .LT. 0E+0) THEN PHI=JACSUM(-1E+0,DG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) PHI=BCFSN(LOM)-2E+0*PHI COF=ABS(PHI)/D**(BT+1E+0) TERM=MCHEP**BT*COF IF (TERM .GT. CPHCA) THEN CPHCA=TERM COPHC=COF EXPHC=BT ENDIF ENDIF IF (BT .GT. 0E+0) THEN PHI=JACSUM(-1E+0,DG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) PHI=BCFSN(LOM)-2E+0*PHI IF (ABS(PHI) .EQ. 0E+0) THEN CCAPH=YMAX COCAP=YMAX EXCAP=AL-1E+0 ELSE COF=D/ABS(PHI)**AL TERM=MCHEP**(AL-1E+0)*COF IF (TERM .GT. CCAPH) THEN CCAPH=TERM COCAP=COF EXCAP=AL-1E+0 ENDIF ENDIF ENDIF C C "DO 50" LOOP FOR POINTS INTERIOR TO ARC NUMBER IA C DO 50 I=1,NINTS-1 TT=TT+TINC C C**** ARC LENGTH INCREASE BY GAUSS-LEGENDRE C SUM=0E+0 DO 45 K=1,NQPTS X=TT+5E-1*TINC*(QUPTS(K)-1E+0) SUM=SUM+QUWTS(K)*DSDT(X) 45 CONTINUE SINC=5E-1*TINC*SUM SS=SS+SINC XX=SS/TOTLN C C EVALUATE DIMENSIONLESS BCF DERIVATIVE *YY* C PHI=JACSUM(SJT*TT,DG,ACOEF(LOD),BCOEF(LOD),H0,JACOF) D=DSDT(TT) IF (D .EQ. 0E+0) THEN IER=51 RETURN ENDIF YY=TOTLN*(1E+0+SJT*TT)**BT*PHI/D WRITE(OUCH1,902) XX,YY YMIN=MIN(YY,YMIN) C C ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS FOR FINITE C DERIVATIVE CASES. C IF (NINFD .EQ. 0E+0) THEN CPHCA=MAX(CPHCA,TUPI*ABS(YY)/TOTLN) ENDIF IF (NZERD .EQ. 0E+0) THEN IF (YY .EQ. 0E+0) THEN CCAPH=YMAX ELSE CCAPH=MAX(CCAPH,TOTLN/TUPI/ABS(YY)) ENDIF ENDIF C C EVALUATE DIMENSIONLESS BCF *YY* C PHI=JACSUM(SJT*TT,DG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) PHI=BCFSN(LOM)-(1E+0-SJT*TT)*PHI YY=(CC+SJT*(1E+0+SJT*TT)**(1E+0+BT)*PHI)/TUPI WRITE(OUCH0,902) XX,YY 50 CONTINUE C C NEXT TAKE END POINT OF ARC NUMBER IA C TT=1E+0 D=DSDT(TT) SEND=SEND+ARCLN(IA) SS=SEND XX=SS/TOTLN C C EVALUATE DIMENSIONLESS BCF DERIVATIVE *YY* C IF (JT .LT. 0E+0) THEN IF (BT .LT. 0E+0) THEN YY=YMAX NASYM=NASYM+1 ASYMP(NASYM)=XX ELSE IF (BT .GT. 0E+0) THEN YY=0E+0 ELSE PHI=JACSUM(SJT*TT,DG,ACOEF(LOD),BCOEF(LOD),H0,JACOF) IF (D .EQ. 0E+0) THEN IER=51 RETURN ENDIF YY=TOTLN*PHI/D ENDIF ELSE PHI=JACSUM(TT,DG,ACOEF(LOD),BCOEF(LOD),H0,JACOF) IF (D .EQ. 0E+0) THEN IER=51 RETURN ENDIF YY=TOTLN*2E+0**BT*PHI/D ENDIF WRITE(OUCH1,902) XX,YY YMIN=MIN(YY,YMIN) IF (YMIN.LT.0E+0 .AND. (VTARG(IA+1) .GE. VTARG(IA))) THEN NPRVS=NPRVS+1 IPRVS(NPRVS)=IA BCDMN(NPRVS)=YMIN MAP11=.FALSE. ENDIF C C ESTIMATE FUNCTION EVALUATION CONDITION NUMBERS C IF (NINFD .EQ. 0E+0) THEN CPHCA=MAX(CPHCA,TUPI*ABS(YY)/TOTLN) ENDIF IF (NZERD .EQ. 0E+0) THEN IF (YY .EQ. 0E+0) THEN CCAPH=YMAX ELSE CCAPH=MAX(CCAPH,TOTLN/TUPI/ABS(YY)) ENDIF ENDIF IF (BT .LT. 0E+0) THEN PHI=JACSUM(-1E+0,DG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) PHI=BCFSN(LOM)-2E+0*PHI COF=ABS(PHI)/D**(BT+1E+0) TERM=MCHEP**BT*COF IF (TERM .GT. CPHCA) THEN CPHCA=TERM COPHC=COF EXPHC=BT ENDIF ENDIF IF (BT .GT. 0E+0) THEN PHI=JACSUM(-1E+0,DG-1,A1COF(LOD),B1COF(LOD),H1VAL(AJT), + BCFSN(LOM+1)) PHI=BCFSN(LOM)-2E+0*PHI IF (ABS(PHI) .EQ. 0E+0) THEN CCAPH=YMAX COCAP=YMAX EXCAP=AL-1E+0 ELSE COF=D/ABS(PHI)**AL TERM=MCHEP**(AL-1E+0)*COF IF (TERM .GT. CCAPH) THEN CCAPH=TERM COCAP=COF EXCAP=AL-1E+0 ENDIF ENDIF ENDIF C C EVALUATE DIMENSIONLESS BCF *YY* C YY=(VTARG(IA+1)-VTARG(1))/TUPI WRITE(OUCH0,902) XX,YY IF (JT.LT.0) THEN CORXX(PT+1)=XX ENDIF C 60 CONTINUE C 901 FORMAT(2E16.8,1X,A3) 902 FORMAT(2E16.8) C C NORMAL EXIT C IER=0 C END REAL FUNCTION DIAPHY(NARCS) INTEGER NARCS C C**** THE APPROXIMATE DIAMETER OF THE PHYSICAL DOMAIN. C C LOCAL VARIABLES C INTEGER I,IA,NH REAL A1,HH,T,DPD COMPLEX C1,CENTR,PARFUN PARAMETER (NH=5) EXTERNAL PARFUN C C**** GET ROUGH ESTIMATE OF CENTRE OF DOMAIN C CENTR=(0E+0,0E+0) DO 10 IA=1,NARCS C1=PARFUN(IA,(-1E+0,0E+0)) CENTR=CENTR+C1 C1=PARFUN(IA,(0E+0,0E+0)) CENTR=CENTR+C1 10 CONTINUE CENTR=CENTR/2E+0/NARCS C C**** GET ROUGH ESTIMATE OF MAXIMUM DISTANCE FROM CENTRE C DPD=0E+0 HH=2E+0/REAL(NH) DO 30 IA=1,NARCS T=-1E+0 DO 20 I=1,NH T=T+HH C1=PARFUN(IA,CMPLX(T))-CENTR A1=ABS(C1) DPD=MAX(DPD,A1) 20 CONTINUE 30 CONTINUE DIAPHY=2E+0*DPD C END SUBROUTINE DMCANP(NPTS,PHYPT,CANPT,INTER,CENTR,IGEOM,RGEOM, +ISNCA,RSNCA,ZSNCA,IQUCA,ZQUCA,WANTM,IER) C INTEGER NPTS,IER INTEGER IGEOM(*),ISNCA(*),IQUCA(*) REAL RGEOM(*),RSNCA(*) COMPLEX CENTR COMPLEX PHYPT(*),CANPT(*),ZSNCA(*),ZQUCA(*) LOGICAL INTER,WANTM C C ...................................................................... C C 1. DMCANP C DOMAIN MAPPING FOR THE CANONICAL --> PHYSICAL MAP. C C 2. PURPOSE C GIVEN A VECTOR OF ARBITRARY POINTS IN THE CANONICAL DOMAIN, C THIS ROUTINE COMPUTES THE VECTOR OF APPROXIMATE IMAGE POINTS C IN THE PHYSICAL DOMAIN. C C 3. CALLING SEQUENCE C CALL DMCANP(NPTS,PHYPT,CANPT,INTER,CENTR,IGEOM,RGEOM,ISNCA, C RSNCA,ZSNCA,IQUCA,ZQUCA,WANTM,IER) C C PARAMETERS C ON ENTRY C NPTS - INTEGER C THE NUMBER OF POINTS TO BE MAPPED. C C CANPT - COMPLEX ARRAY C A COMPLEX VECTOR OF SIZE AT LEAST NPTS. THIS IS C THE VECTOR OF GIVEN POINTS IN THE CANONICAL C DOMAIN. C C INTER - LOGICAL C TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE C OTHERWISE. (AS PREVIOUSLY USED IN JAPHYC, GQPHYC) C C CENTR - COMPLEX C THE POINT IN THE PHYSICAL PLANE THAT IS TO BE C MAPPED TO THE CENTRE OF THE UNIT DISC. FOR C EXTERIOR DOMAINS CENTR MUST BE SOME POINT IN THE C COMPLEMENTARY INTERIOR PHYSICAL DOMAIN. (AS PREV- C IOUSLY USED IN JAPHYC, GQPHYC) C C IGEOM - INTEGER ARRAY C THE INTEGER VECTOR IGEOM PREVIOUSLY SET UP BY C JAPHYC. C C RGEOM - REAL ARRAY C THE REAL VECTOR RGEOM PREVIOUSLY SET UP BY JAPHYC. C C C ISNCA - INTEGER ARRAY C THE INTEGER VECTOR PREVIOUSLY SET UP BY JACANP. C C RSNCA - REAL ARRAY C THE REAL VECTOR PREVIOUSLY SET UP BY JACANP. C C ZSNCA - COMPLEX ARRAY C THE COMPLEX VECTOR PREVIOUSLY SET UP BY JACANP. C C IQUCA - INTEGER ARRAY C THE INTEGER VECTOR PREVIOUSLY SET UP BY GQCANP. C C RQUCA - REAL ARRAY C THE REAL VECTOR PREVIOUSLY SET UP BY GQCANP. C C ZQUCA - COMPLEX ARRAY C THE COMPLEX VECTOR PREVIOUSLY SET UP BY GQCANP. C C WANTM - LOGICAL C IF WANTM IS TRUE THEN, ON AN ABNORMAL EXIT, AN C ERROR MESSAGE IS WRITTEN ON THE STANDARD OUTPUT C CHANNEL. IF WANTM IS FALSE THEN NO MESSAGE IS C WRITTEN. C ON EXIT C PHYPT - COMPLEX ARRAY C A COMPLEX VECTOR OF SIZE AT LEAST NPTS. PHYPT(K) C IS THE COMPUTED IMAGE IN THE PHYSICAL DOMAIN OF C THE GIVEN CANONICAL POINT CANPT(K), K=1,...,NPTS. C C IER - INTEGER C IF IER > 0 THEN AN ABNORMAL EXIT HAS OCCURRED; C A MESSAGE TO DESCRIBE THE ERROR IS WRITTEN ON C THE STANDARD OUTPUT CHANNEL IF WANTM IS TRUE. C IER=0 - NORMAL EXIT. C IER>0 - ABNORMAL EXIT; THE ERROR MESSAGE SHOULD C BE SELF EXPLANATORY. C C C 4. SUBROUTINES OR FUNCTIONS NEEDED C - THE CONFPACK LIBRARY. C - THE REAL FUNCTION R1MACH. C - THE USER SUPPLIED COMPLEX FUNCTIONS PARFUN AND DPARFN. C C C 5. FURTHER COMMENTS C - NOTE THAT THIS ROUTINE CAN ONLY BE USED A F T E R THE C ROUTINES JACANP AND GQCANP HAVE SUCCESSFULLY EXECUTED, C AND THAT MANY INPUT ARGUMENTS FOR DMCANP ARE OUTPUT VALUES C FROM JACANP AND GQCANP. C - THIS ROUTINE MAY BE USED FOR MAPPING POINTS ON THE UNIT C CIRCLE, BUT THE ROUTINE BMCANP WILL BE SOMEWHAT MORE C EFFICIENT FOR THIS CASE. C C ...................................................................... C AUTHOR: DAVID HOUGH, ETH, ZUERICH C LAST UPDATE: 3 JULY 1990 C ...................................................................... C C LOCAL VARAIBLES C INTEGER ACOFC,AICOC,BCOFC,BFSNC,BICOC,DGPOC,H0VLC,HALEN,HIVLC, +JAINC,JTYPC,LQSBG,LSUBC,MIDPT,MNCOF,MNSUA,MNSUC,MQUCA,NARCS,NJIND, +NPPQG,NQPTS,PARNT,PHPAS,PRNSA,QUWTC,QUPTC,SOLNC,TNGQP,TNSUC,VARGC, +VTARG,WPPQG,ZPPQG CHARACTER*6 IERTXT C EXTERNAL CATPH4,IERTXT C NARCS=ISNCA(1) NQPTS=ISNCA(2) TNSUC=ISNCA(3) MNSUC=ISNCA(5) MNCOF=ISNCA(6) MQUCA=IQUCA(4) MNSUA=IGEOM(4) C NJIND=NARCS+1 TNGQP=NJIND*NQPTS C C**** SET UP POINTERS TO IGEOM AND RGEOM, AS IN JAPHYC C PARNT=5 HALEN=3 MIDPT=MNSUA+3 VTARG=2*MNSUA+3 C C**** SET UP POINTERS TO ISNCA, RSNCA AND ZSNCA, AS IN JACANP C DGPOC=7 JTYPC=MNSUC+7 LSUBC=2*MNSUC+7 PRNSA=3*MNSUC+7 ACOFC=2 BCOFC=TNGQP+2 AICOC=2*TNGQP+2 BICOC=3*TNGQP+2 QUPTC=4*TNGQP+2 QUWTC=5*TNGQP+2 H0VLC=6*TNGQP+2 HIVLC=NJIND+6*TNGQP+2 JAINC=2*NJIND+6*TNGQP+2 PHPAS=4*NJIND+6*TNGQP+2 VARGC=MNSUC+4*NJIND+6*TNGQP+2 BFSNC=2 SOLNC=MNCOF+2 C C**** SET UP POINTERS TO IQUCA AND ZQUCA, AS IN GQCANP C LQSBG=5 NPPQG=MNSUC+5 WPPQG=2 ZPPQG=MQUCA+2 C C**** GET THE REQUIRED PHYSICAL POINTS C CALL CATPH4(NPTS,PHYPT,CANPT,NARCS,NQPTS,TNSUC,ISNCA(DGPOC), +ISNCA(JTYPC),ISNCA(LSUBC),IQUCA(LQSBG),IQUCA(NPPQG),IGEOM(PARNT), +ISNCA(PRNSA),RSNCA(AICOC),RSNCA(ACOFC),RSNCA(BICOC),RSNCA(BCOFC), +RSNCA(H0VLC),RSNCA(HIVLC),RGEOM(HALEN),RSNCA(JAINC),RGEOM(2), +RGEOM(MIDPT),RSNCA(PHPAS),RSNCA(QUPTC),RSNCA(QUWTC),RGEOM(VTARG), +RSNCA(VARGC),ZSNCA(BFSNC),CENTR,ZSNCA(1),ZSNCA(SOLNC), +ZQUCA(WPPQG),ZQUCA(ZPPQG),INTER,IER) C C**** SEND ERROR MESSAGE TO STANDARD OUTPUT OF NECESSARY C IF (IER.GT.0 .AND. WANTM) WRITE(*,*) IERTXT(IER) C END SUBROUTINE DMPHYC(NPTS,PHYPT,CANPT,INTER,CENTR,IGEOM,RGEOM,ISNPH, +RSNPH,IQUPH,RQUPH,ZQUPH,WANTM,IER) C INTEGER NPTS,IER INTEGER IGEOM(*),ISNPH(*),IQUPH(*) REAL RGEOM(*),RSNPH(*),RQUPH(*) COMPLEX CENTR COMPLEX PHYPT(*),CANPT(*),ZQUPH(*) LOGICAL INTER,WANTM C C ...................................................................... C C 1. DMPHYC C DOMAIN MAPPING FOR THE PHYSICAL --> CANONICAL MAP. C C 2. PURPOSE C GIVEN A VECTOR OF ARBITRARY POINTS IN THE PHYSICAL DOMAIN, C THIS ROUTINE COMPUTES THE VECTOR OF APPROXIMATE IMAGE POINTS C IN THE CANONICAL DOMAIN. C C C 3. CALLING SEQUENCE C CALL DMPHYC(NPTS,PHYPT,CANPT,INTER,CENTR,IGEOM,RGEOM,ISNPH, C RSNPH,IQUPH,RQUPH,ZQUPH,WANTM,IER) C C PARAMETERS C ON ENTRY C NPTS - INTEGER C THE NUMBER OF POINTS TO BE MAPPED. C C PHYPT - COMPLEX ARRAY C A COMPLEX VECTOR OF SIZE AT LEAST NPTS. THIS IS C THE VECTOR OF GIVEN POINTS IN THE PHYSICAL DOMAIN. C C INTER - LOGICAL C TRUE IF THE PHYSICAL DOMAIN IS INTERIOR, FALSE C OTHERWISE. (AS PREVIOUSLY USED IN JAPHYC, GQPHYC) C C CENTR - COMPLEX C THE POINT IN THE PHYSICAL PLANE THAT IS TO BE C MAPPED TO THE CENTRE OF THE UNI