/DOC/Derek Nielson Thesis/Research Documents/UAFWLIS-V1_01_FiniteSizeGround/OMEGA_complexZ/Che0.f
FORTRAN Legacy | 317 lines | 153 code | 31 blank | 133 comment | 0 complexity | 6fa1dcb9abd3db95817bb83142b20b6e MD5 | raw file
Possible License(s): GPL-3.0
- SUBROUTINE CHE0(A,Z,NUMA,MNUMA,PCKEXP1,PCKEXP2,SQASP11,SQASP12,
- & EXPAZ,H01,H02,H11,H12,HE01CX,HE02CX,HE01EX,
- & HE02EX,HE01SX,HE02SX,HE01D0,HE02D0,FLAG,RELERR)
-
- * *
- *******************************************************************************
- * REFERENCES: *
- * *
- * [2] S. L. Dvorak and E. F. Kuester, ``Numerical Computation of the *
- * Incomplete Lipschitz-Hankel Integral Je0(a,z),'' J. Comput. Phys., *
- * Vol. 87, No. 2, pp. 301--327, April 1990. *
- * [3] M. Abramowitz and I. E. Stegun, Handbook of Mathematical *
- * Functions with Formulas, Graphs, and Mathematical Tables, *
- * Washington, D. C.: U.S. Government Printing Office, 1972. *
- * [4] David L. Heckmann and S.L. Dvorak,"Numerical Computation of Hankel *
- * Functions of Integral Order for Complex-valued Arguments *
- * [6] S. L. Dvorak, ``Numerical Computation of the Incomplete *
- * Lipschitz-Hankel Integral Ye0(a,z),'' submitted to the SIAM *
- * Journal on Mathematical Analysis. *
- * [7] Z.ZHU, S.L.DVORAK AND DAVID HECKMAN,"Numerical Computation of *
- * Incomplete Lipschitz-Hankel Integrals of the Hankel Type for *
- * Complex-Valued Arguments *
- * *
- *******************************************************************************
- * PURPOSE: *
- * *
- * This subroutine determines which expansion will provide the most *
- * efficient computation of the incomplete Lipschitz-Hankel integrals *
- * AND OUTPUT THE NUMERICAL VALUE FOR: *
- * HE0*exp(az), CHEO*exp(az) and SCHE0*exp(az) *
- * HE0 means standard ILHI of Hankel type. *
- * CHE0 means complimentary ILHI of Hankel type *
- * SCHE0 means CHEO(a,*,z) *
- * *
- * where A and Z are complex numbers. and Real(Z)>=0
- * The integrals are evaluated to SD significant digits in most cases. *
- * The relative error is given by RELERR = .5*10.**(-SD). *
- *******************************************************************************
- * INPUTS: *
- * *
- * A = Factor in the exponential of the incomplete Lipschitz-Hankel *
- * integral. (Complex array, Size is set by MNUMA) *
- * Z = Upper limit of integration. (Complex, Real(Z)>0) *
- * NUMA = Number of actual terms in the A array. *
- * FLAG = 0; NEED TO COMPUTE BOTH HE01 AND HE02 *
- * = 1; ONLY NEED TO COMPUTE HE01 *
- * = 2; ONLY NEED TO COMPUTE HE02 *
- *************************************************************************************************
- * OUTPUTS: *
- * *
- * PCKEXP1/PCKEXP2 = Tells which expansion was used. *
- * HE01EX/HE02EX = Computed value for the incomplete Lipschitz-Hankel integral. *
- * HE01CX/HE02CX = Computed value for the complimentary incomplete Lipschitz-Hankel integral.*
- * HE01SX/HE02SX = Computed value for the incomplete Lipschitz-Hankel star integral. *
- * HE01D0/HE02D0 = RETURN HE01(A,DELTA,0)*EXPAZ/HE02(A,DELTA,0)*EXPAZ IF COMPUTED. *
- * SQASP11/SQASP12= SQRT(a^2+1) CONSIDERING THE BRANCH CUT WHEN USING ASYMPTOTIC OR *
- * QUASI-NEUMANN.OTHERWISE, ALWAYS POSITVE REAL PART. *
- *************************************************************************************************
- * SUBROUTINES AND FUNCTIONS CALLED: *
- * *
- * CHANKL = Computes Jn, Yn AND HN1, HN2 *
- * CFNDEXP1/CFNDEXP2 = Determine which expansion to use. *
- * CFDHE01/CFDHE02 = Compute HE01/HE02 using the proper expansion. *
- *******************************************************************************
- * VARIABLE DECLARATIONS:
- *
- $Declare
- PARAMETER (MAXNOJ=500)
-
- INTEGER,INTENT(IN):: NUMA,MNUMA,FLAG
- COMPLEX*16,INTENT(IN)::Z,A(MNUMA)
- COMPLEX*16,INTENT(OUT):: HE01D0(MNUMA),HE02D0(MNUMA)
- COMPLEX*16,INTENT(INOUT)::SQASP11(MNUMA),EXPAZ(MNUMA),
- & SQASP12(MNUMA)
- COMPLEX*16,INTENT(OUT):: HE01CX(MNUMA), HE02CX(MNUMA)
- COMPLEX*16,INTENT(OUT):: HE01EX(MNUMA), HE02EX(MNUMA)
- COMPLEX*16,INTENT(OUT):: HE01SX(MNUMA), HE02SX(MNUMA)
- COMPLEX*16,INTENT(OUT):: H01, H02, H11, H12
- INTEGER,INTENT(OUT):: PCKEXP1(MNUMA),PCKEXP2(MNUMA)
- DOUBLE PRECISION,INTENT(IN):: RELERR
-
- *********LOCAL VARIABLE************************************************
- DOUBLE PRECISION ZJASY,ZSP,COMPAR,ABASP1
- COMPLEX*16 JBESS(0:MAXNOJ),YBESS(0:MAXNOJ)
- COMPLEX*16 JBESSAZ(0:MAXNOJ),YBESSAZ(0:MAXNOJ)
- COMPLEX*16 HN1(0:MAXNOJ), HN2(0:MAXNOJ)
- COMPLEX*16 BFH0,BFH1,HN1AZ(0:MAXNOJ), HN2AZ(0:MAXNOJ)
- COMPLEX*16 ASP1,U(0:MAXNOJ),LOGA(MNUMA),LOGNA(MNUMA),FRONT(MNUMA)
- COMPLEX*16 EK(0:MAXNOJ),P(0:MAXNOJ), HE01, HE02,CE1
-
- INTEGER NOJ,KMAX_1(NUMA),INTZ,NOH1,NOH2,N,FWD_FLAG,NOY
- INTEGER NOJAZ,NOYAZ,NOH1AZ,NOH2AZ,KMAX_2(NUMA),K,FWD_FLAGAZ
-
- DOUBLE PRECISION APROX1, APROX2,AZ
- DOUBLE PRECISION IMZ
- LOGICAL NEWZ, NEEDJ,NEEDSTRU,NEEDAZ,NEEDUEP,NEEDJ_FLAG
- INTEGER FWD_FLAG1,FWD_FLAG2
- PARAMETER (E=2.71828)
- ****************************************************************************
- * NEEDJ_FLAG=0 ALL THE SERIES IS SECOND CONVERGENT, SO NO NEED TO COMPUTE Jn(Z)
- * =1 AT LEAST ONE OF THE SERIES NEED COMPUTE Jn(Z) BY BACKWARD
- * =2 ALL THE SERIES CAN USE FORWARD TO GET HN1 AND HN2, NO Jn(Z) NEEDED
- *******************************************************************************
- * COMMON BLOCKS
- *
- * FACTOR(K)=SQPI/(2.**(K+1)*GAMMA(K+1.5))
- DOUBLE PRECISION FACTOR(0:500)
- COMMON/FACT/FACTOR
- DOUBLE PRECISION PI, TWOPI
- COMPLEX*16 IM
- COMMON/CONST/ PI, TWOPI,IM
-
- INTEGER ZLO, ALO(2)
- COMMON /ZALO/ ZLO, ALO
-
- INTEGER SD
- COMMON /SIGDIG/ SD
- *******************************************************************************
- EXTERNAL CHANKL,CFDHEXP1,CFDHE01,CFDHEXP2,CFDHE02,CHAE01,CHAE02
- ****************************************************************************
- * Calculate an array of Hankel functions. This subroutine calculates
- * as many Hankel functions that are required for forward recurrence
- * to be stable for higher orders. The higher order Hankel functions
- * are calculated as needed.[4]
- NEEDJ=.TRUE.
-
- DO 10 N = 1, NUMA
- *SQASP1 is always choose the positive real part.
- SQASP12(N)=SQASP11(N)
- 10 CONTINUE
-
- ZSP=ABS(Z)
- INTZ=NINT(ZSP+2)
- ZJASY = REAL(SD+5)
-
- **************** Determine which expansion to use.
-
- IF(FLAG==2)THEN
- DO N=1,NUMA
- CALL CFDHEXP2(A(N),Z,SQASP12(N),RELERR,KMAX_2(N),PCKEXP2(N))
- END DO
- * Determine if the CHANKL.F need to be called, i.e. if need to compute Jn(Z) by backward
- * THOUGH PCKEXP=12 NEED NOT COMPUTE HN1/2(Z), HOWEVER, WE NEED RETURN HN1/2(Z) FIRST AND SECOND ORDER
- * SET DEFAULT TO USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
- NEEDJ_FLAG=.FALSE.
- DO N=1,NUMA
- IF (PCKEXP2(N)>=14.OR. ZSP<20.0) THEN
- * USE BACKWARD TO COMPUTE JBESS, HN1/2(0) AND HN1/2(1)
- NEEDJ_FLAG=.TRUE.
- END IF
- END DO
-
- IF(NEEDJ_FLAG==.FALSE.) THEN
- * USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
- CALL CHAE02(Z,RELERR*1.D-6,HN2(0),HN2(1))
- IF (REAL(Z)==0. .AND. IMAG(Z)<0.) THEN
- HN2(0)=IM*IMAG(HN2(0))
- HN2(1)=DBLE(HN2(1))
- END IF
- NOH2=1
- FWD_FLAG2=2
- ELSE
- * USE BEWARD TO COMPUTE Jn AND HN1/HN2
- CALL CHANKL(Z,RELERR,JBESS,NEEDJ,MAXNOJ,NOJ,HN1,HN2,
- & NOH1,NOH2,YBESS,NOY,FWD_FLAG)
- FWD_FLAG2=FWD_FLAG
- END IF
-
- ***********Compute the ILHIs using the proper expansion.
- NEEDSTRU=.TRUE.
- DO N = 1, NUMA
- NEEDAZ=.TRUE.
- NEEDUEP=.TRUE.
-
- IF(PCKEXP2(N)==14) THEN
- DO K = 0, MAXNOJ
- U(K) = 0.D0
- END DO
- END IF
-
- CALL CFDHE02(Z,A(N),RELERR,SQASP12(N),PCKEXP2(N),KMAX_2(N),
- & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
- & NEEDJ,HE02CX(N),HE02EX(N),HE02SX(N),HE02D0(N),U,
- & EK,P,NEEDUEP,ZJASY,FWD_FLAG2,BFH0,BFH1,NEEDSTRU,JBESSAZ,
- & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
- & FWD_FLAGAZ,CE1)
- END DO
-
- ELSE IF(FLAG==1)THEN
- DO N=1,NUMA
- CALL CFDHEXP1(A(N),Z,SQASP11(N),RELERR,KMAX_1(N),PCKEXP1(N))
- END DO
- * Determine if the CHANKL.F need to be called, i.e. if need to compute Jn(Z) by backward
- * THOUGH PCKEXP=12 NEED NOT COMPUTE HN1/2(Z), HOWEVER, WE NEED RETURN HN1/2(Z) FIRST AND SECOND ORDER
- * SET DEFAULT TO USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
- NEEDJ_FLAG=.FALSE.
- DO N=1,NUMA
- IF (PCKEXP1(N)>=14.OR. ZSP<20.0) THEN
- * USE BACKWARD TO COMPUTE JBESS, HN1/2(0) AND HN1/2(1)
- NEEDJ_FLAG=.TRUE.
- END IF
- END DO
-
- IF(NEEDJ_FLAG==.FALSE.) THEN
- * USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
- CALL CHAE01(Z,RELERR*1.D-6,HN1(0),HN1(1))
- IF(REAL(Z)==0. .AND. IMAG(Z)>0.) THEN
- HN1(0)=IM*IMAG(HN1(0))
- HN1(1)=DBLE(HN1(1))
- END IF
- NOH1=1
- FWD_FLAG1=1
- ELSE
- * USE BEWARD TO COMPUTE Jn AND HN1/HN2
- CALL CHANKL(Z,RELERR,JBESS,NEEDJ,MAXNOJ,NOJ,HN1,HN2,
- & NOH1,NOH2,YBESS,NOY,FWD_FLAG)
- FWD_FLAG1=FWD_FLAG
- END IF
-
- ***********Compute the ILHIs using the proper expansion.
- NEEDSTRU=.TRUE.
- DO N = 1, NUMA
- NEEDAZ=.TRUE.
- NEEDUEP=.TRUE.
-
- IF(PCKEXP1(N)==14) THEN
- DO K = 0, MAXNOJ
- U(K) = 0.D0
- END DO
- END IF
-
- CALL CFDHE01(Z,A(N),RELERR,SQASP11(N),PCKEXP1(N),KMAX_1(N),
- & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
- & NEEDJ,HE01CX(N),HE01EX(N),HE01SX(N),HE01D0(N),U,
- & EK,P,NEEDUEP,ZJASY,FWD_FLAG1,BFH0,BFH1,NEEDSTRU,JBESSAZ,
- & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
- & FWD_FLAGAZ,CE1)
-
- END DO
-
- ELSE
- DO N=1,NUMA
- CALL CFDHEXP1(A(N),Z,SQASP11(N),RELERR,KMAX_1(N),PCKEXP1(N))
- CALL CFDHEXP2(A(N),Z,SQASP12(N),RELERR,KMAX_2(N),PCKEXP2(N))
- END DO
- * Determine if the CHANKL.F need to be called, i.e. if need to compute Jn(Z) by backward
- * THOUGH PCKEXP=12 NEED NOT COMPUTE HN1/2(Z), HOWEVER, WE NEED RETURN HN1/2(Z) FIRST AND SECOND ORDER
- * SET DEFAULT TO USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
- NEEDJ_FLAG=.FALSE.
- DO N=1,NUMA
- IF (PCKEXP1(N)>=14 .OR. PCKEXP2(N)>=14
- & .OR. ZSP<20.0) THEN
- * USE BACKWARD TO COMPUTE JBESS, HN1/2(0) AND HN1/2(1)
- NEEDJ_FLAG=.TRUE.
- END IF
- END DO
-
- IF(NEEDJ_FLAG==.FALSE.) THEN
- * USE ASYMPTOTIC EXPRESION TO COMPUTE HN1/2(0) AND HN1/2(1)
- CALL CHAE01(Z,RELERR*1.D-6,HN1(0),HN1(1))
- CALL CHAE02(Z,RELERR*1.D-6,HN2(0),HN2(1))
- IF(REAL(Z)==0. .AND. IMAG(Z)>0.) THEN
- HN1(0)=IM*IMAG(HN1(0))
- HN1(1)=DBLE(HN1(1))
- ELSE IF (REAL(Z)==0. .AND. IMAG(Z)<0.) THEN
- HN2(0)=IM*IMAG(HN2(0))
- HN2(1)=DBLE(HN2(1))
- END IF
- NOH1=1
- NOH2=1
- FWD_FLAG1=1
- FWD_FLAG2=2
- ELSE
- * USE BEWARD TO COMPUTE Jn AND HN1/HN2
- CALL CHANKL(Z,RELERR,JBESS,NEEDJ,MAXNOJ,NOJ,HN1,HN2,
- & NOH1,NOH2,YBESS,NOY,FWD_FLAG)
- FWD_FLAG1=FWD_FLAG
- FWD_FLAG2=FWD_FLAG
- END IF
-
- ***********Compute the ILHIs using the proper expansion.
- NEEDSTRU=.TRUE.
- DO N = 1, NUMA
- NEEDAZ=.TRUE.
- NEEDUEP=.TRUE.
-
- IF(PCKEXP2(N)==14 .OR. PCKEXP1(N)==14) THEN
- DO K = 0, MAXNOJ
- U(K) = 0.D0
- END DO
- END IF
-
- CALL CFDHE01(Z,A(N),RELERR,SQASP11(N),PCKEXP1(N),KMAX_1(N),
- & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
- & NEEDJ,HE01CX(N),HE01EX(N),HE01SX(N),HE01D0(N),U,
- & EK,P,NEEDUEP,ZJASY,FWD_FLAG1,BFH0,BFH1,NEEDSTRU,JBESSAZ,
- & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
- & FWD_FLAGAZ,CE1)
-
- CALL CFDHE02(Z,A(N),RELERR,SQASP12(N),PCKEXP2(N),KMAX_2(N),
- & EXPAZ(N),JBESS,YBESS,MAXNOJ,NOJ,NOY,HN1,HN2,NOH1,NOH2,
- & NEEDJ,HE02CX(N),HE02EX(N),HE02SX(N),HE02D0(N),U,
- & EK,P,NEEDUEP,ZJASY,FWD_FLAG2,BFH0,BFH1,NEEDSTRU,JBESSAZ,
- & YBESSAZ,NOJAZ,NOYAZ,HN1AZ,HN2AZ,NOH1AZ,NOH2AZ,NEEDAZ,
- & FWD_FLAGAZ,CE1)
- END DO
-
- END IF
- H01=HN1(0)
- H11=HN1(1)
- H02=HN2(0)
- H12=HN2(1)
-
- RETURN
- END