C=====================================================================
C
C ***** STEVE SHENG'S FHT RESONATOR AND PROPAGATION PACKAGE *****
C ***** PULLED TOGETHER INTO ONE FILE BY AES, APRIL 1996 *****
C
C FURTHER DETAILS ON THIS CAN BE FOUND IN STEVE SHENG'S PH.D.
C DISSERTATION, S.-C. SHENG, STUDIES OF LASER RESONATORS AND
C BEAM PROPAGATION USING FAST TRANSFORM METHODS, DEPARTMENT OF
C APPLIED PHYSICS, STANFORD UNIVERSITY (MARCH 1980), AVAILABLE
C FROM UNIVERSITY MICROFILMS, ANN ARBOR, MICHIGAN.
C
C YOU *MAY* BE ABLE TO GET SOME ADDITIONAL INFORMATION FROM
C PROF. A. E. SIEGMAN, SIEGMAN@STANFORD.EDU, BUT THIS IS
C DEFINITELY *NOT* GUARANTEED!
C
C ON THE OTHER HAND, IF YOU USE THIS PACKAGE AND FIND IT USEFUL,
C I'D ALWAYS BE GLAD TO HEAR ABOUT THAT...
C
C=====================================================================
C
C
C THIS FILE CONTAINS:
C
C (1) A TEST DRIVER PROGRAM FOR TESTING SHENG'S FHT PROGRAM
C
C (2) AN INTERATIVE UNSTABLE-RESONATOR PROGRAM USING THE FHT
C
C (3) A ONE-STEP PROPAGATION-DIFFRACTION PROGRAM USING FHT
C
C (4) THE SUBROUTINES NEEDED BY THESE PROGRAMS
C
C ALL THE COMMON SUBROUTINES CALLED BY ALL OF THESE PROGRAMS
C ARE GROUPED IN A SINGLE BATCH AT THE END OF THESE THREE
C MAIN PROGRAMS
C
C SEARCH ON THE STRING " ***** " TO SEE WHAT'S IN THIS FILE
C NOTE: SOME OF THE ARRAYS IN SOME OF THE SUBROUTINES MAY HAVE
C TO BE REDIMENSIONED DIFFERENTLY IN ORDER TO WORK PROPERLY WITH
C THE DIFFERENT MAIN PROGRAMS -- CHECK THIS BEFORE USING
C
C=====================================================================
C
C ***** TEST DRIVER PROGRAM FOR STEVE SHENG'S FHT PROGRAM *****
C
C=====================================================================
C
C.....................................................................
C PROGRAM FHTNEW,TEST THE MODIFIED FAST HANKEL TRANSFORMATIONS
C WITH LOWER END CORRECTIONS(LEC).CALCULATIONS OF LB=0,1
C ARE INCLUDED & TWO CORRECTION TERMS ARE DONE.
C.....................................................................
C
C
DOUBLE PRECISION DSIN,DCOS,DSQRT,DEXP,DFLOAT
REAL*8 YR(130),YI(130),P1(256),P2(256),DALPH,PI,RC
REAL ITPOR,D(128),ALPH,RC1
REAL*8 RR,RJR,RJI,CF0,CF2
COMMON RR(256),RJR(512),RJI(512),CF0(256),CF2(256)
C
C YR,YI IS INPUT FUNCTION TO BE TRANSFORMED, THEY SHOULD BE
C DIMENSIONED AT LEAST TO (N+2). FOR LB=0, FUNCTION VALUE &
C ITS 2ND DERIVATIVE AT ORIGIN(A.O.) ARE STORED AT THE END
C OF THE ARRAYS.***** ASSUME INPUT 1ST,3RD DERIV A.O. ARE 0.
C FOR LB=1, 1ST & 3RD DERIV A.O. ARE STORED AT THE END OF
C THE ARRAYS.***** ASSUME INPUT FUNCTION VALUE AND 2ND DERIV
C A.O. ARE ZERO.
C IF 1ST,3RD DERIV. A.O. NOT ZERO FOR THE INPUT FUNCTION IN
C LB=0, THE OUTPUT 1ST,3RD DERIV. A.O. ARE STILL ALWAYS ZERO
C , BUT ITS FUNCTION VALUE & 2ND DERIV A.O. WILL INCLUDE TWO
C EXTRA TERMS COMING FROM THOSE NOZERO INPUT 1ST,3RD DERIV
C AT ORIGIN. THOSE TWO EXTRA TERMS ARE "NOT" INLUDED IN THIS
C PROGRAM. SEE NOTES IN SUBROUTINE FHTS.
C FOR LB=1, SIMILAR SITUATION HAPPENS, I.E. THIS PROGRAM DOES
C "NOT" INCLUDE THE EXTRA TWO TERMS FOR NOZERO FUNCTION VALUE
C AND 2ND DERIV A.O. ALSO SEE NOTES IN SUBROUTINE FHTS.
C
C THIS PROGRAM HAS ONLY TWO CORRECTION TERMS FOR LB=0,1.
C HIGHER BESSEL ORDER CORRECTION CAN BE INCLUDED WITH MINOR
C MODIFICATION OF CF0,CF2 AND THOSE TERMS IN SUBROUTINE FHTS.
C (NOTES WITH -*****-)
C
C CORRECTION ARRAYS CF0(N),CF2(N) CORRESPOND TO FOLLOWING:
C IN LB=0: CF0(FUNCTION VALUE A.O.), CF2(2ND DERIV A.O.)
C IN LB=1: CF0(1ST DERIV A.O.), CF2(2ND DERIV A.O.)
C
C THIS IS A COMPLETE REGENERATIVE PROGRAM,I.E. AFTER THE
C TRANSFORMATION, NEW FUNCTION VALUE & 2ND DERIV A.O.(FOR
C LB=0) OR 1ST & 2ND DERIV A.O.(FOR LB=1) ARE EVALUATED
C SO THAT FURTHER TRANSFORMATIONS ARE POSSIBLE, E.G.
C RESONATOR EIGENMODE PROBLEMS.
C
C ABOUT DETAIL MATHEMATICAL DESCRIPTION IN LEC, PLEASE SEE
C REFERENCE PAPER.
C
C RC CORRESPONDS TO X=-0.5 IN GARDNER SPACE, WHICH IS THE
C INTEGRATION STARTING POINT.
C
PI=3.1415926535898D0
C
N=128
NA2=N+2
NM=N-1
AK1=10.0
AK2=2.0
LB=0
NB=8
C
C NB IS BESSEL UPPER GENERATING ORDER.
C
CALL SETUPZ (N,LB,NB,AK1,AK2,DALPH,RRHO,BETAB,RC)
C
C WRITE (6,796)
C796 FORMAT ('1*** CORRECTION ARRAY CF0,CF2')
C DO 795 K=1,N
C CF0(K)=0.D0
C CF2(K)=0.D0
C WRITE(6,7962) RR(K),CF0(K),CF2(K)
C7962 FORMAT(1H0,3(D17.10,2X))
C795 CONTINUE
C
C THE FOLLOWING ARE INITIAL INPUT DATA. P1 IS USED TO SAVE
C THE INPUT DATA FOR LATER COMPARISON.
C
NZ=2
RC1=RR(1)
C
C-*****-
C FOR DIFFERENT LB,TEST FUNCTION INPUT CHANGE STARTS HERE.
C LB=0: YR(N+1),SFR0 ETC, VALUE AT ORIGIN,
C YR(N+2),SFR2 ETC, 2ND DERIV A.O.
C ASSUME INPUT 1ST,3RD DERIV A.O. ARE ZERO.
C LB=1: YR(N+1),SFR0 ETC, 1ST DERIV A.O.
C YR(N+2),SFR2 ETC, 3RD DERIV A.O.
C ASSUME INPUT FUNCTION VALUE & 2ND DERIV A.O. ARE ZERO.
C
DO 632 K=1,N
C IF LB=0, FOLLOWING TEST FUNCTION IS USED:
YR(K)=DEXP(-PI*RR(K)**2)
C IF LB=1, FOLLOWING TEST FUNCTION IS USED(LAGUERRE-GAUSSIAN
C L=1, NOT NORMALIZED)
C YR(K)=RR(K)*DEXP(-PI*RR(K)**2)
P1(K)=YR(K)
C P1 STORES THE INPUT FUNCTION FOR LATER COMPARISON PURPOSE.
632 YI(K)=0.D0
C FOLLOWING ARE INPUT FUNCTION VALUE & 2ND DERIV A.O. FOR LB=0:
YR(N+1)=1.D0
YR(N+2)=-2.D0*PI
YI(N+1)=0.D0
YI(N+2)=0.D0
C FOLLOWING ARE INPUT 1ST & 3RD DERIV A.O. FOR LB=1:
C YR(N+1)=1.D0
C YI(N+1)=0.D0
C YR(N+2)=-6.D0*PI
C YI(N+2)=0.D0
C
C FOR DIFFERENT LB, INPUT TEST FUNCTION CHANGE ENDS HERE.
C-*****-
C
C START DOING THE TRANSFORMATIONS.
C
DO 7964 KQ=1,NZ
C
C CALCULATE TOTAL POWER ITPOR.
C
SFR0=YR(N+1)
SFR2=YR(N+2)
SFI0=YI(N+1)
SFI2=YI(N+2)
ITPOR=(SFR0*RC1)**2+(SFR2*RC1**3)**2/12.D0+SFR0*SFR2*RC1**4/2.D0
ITPOR=ITPOR+(SFI0*RC1)**2+(SFI2*RC1**3)**2/12.D0
ITPOR=PI*(ITPOR+SFI0*SFI2*RC1**4/2.D0)
C
DO 633 K=1,NM
A1=RR(K+1)**2+RR(K+1)*RR(K)+RR(K)**2
AMP=A1*(YR(K+1)**2-YR(K)**2)*2.0*PI/3.0
AMP=AMP+A1*(YI(K+1)**2-YI(K)**2)*2.0*PI/3.0
A=-PI*(RR(K+1)+RR(K))*(YR(K+1)**2*RR(K)-YR(K)**2*RR(K+1))
A=A-PI*(RR(K+1)+RR(K))*(YI(K+1)**2*RR(K)-YI(K)**2*RR(K+1))
633 ITPOR=ITPOR+AMP+A
C
C PRINT SOME PARAMETERS.
C KT IS THE NUMBER OF TRANSFORMATION.
C
KT=KQ-1
WRITE (6,1991) KT,SFR0,SFR2,ITPOR,YR(1)
1991 FORMAT (1H0,'KT=',I4,1X,'SFR0=',D27.20,'SFR2=',
1 D27.20,'ITPOR=',E16.9,D27.20)
C
C P1 IS ORIGINAL GAUSSIAN INPUT,P2 IS THE DIFFERENCE BETWEEN
C THE INPUT AND OUTPUT OF FHT.
C
DO 965 KI=1,N
P2(KI)=YR(KI)-P1(KI)
WRITE(6,7965) RR(KI),P2(KI)
7965 FORMAT(1H0,2(D27.20,2X))
965 CONTINUE
C
C HERE IS THE TRANSFORMATION.
C
CALL FHTZ1(YR,YI,N,NA2,DALPH,RC,LB)
7964 CONTINUE
STOP
END
C......................................................................
C
C
C$DATA
//
<p>C
C=====================================================================
C
C ***** UNSTABLE RESONATOR PROGRAM USING THE FHT SUBROUTINES *****
C
C=====================================================================
C
***** STEVE SHENG'S ITERATIVE UNSTABLE-RESONATOR PROGRAM *****
<p>//FHTNEW JOB 'N71$D1,316,2.00,60','RESOIT',CLASS=E,REGION=512K
// EXEC FORTCG
//FORT.SYSIN DD *
C// EXEC WATFIV,LEVEL=5
C//GO.SYSIN DD *
C$WATFIV RUN=FREE
C WATFIV USES $DATA INSTEAD OF // AT THE END.
C// EXEC FORTCG
C//FORT.SYSIN DD *
C IF WE DON'T WANT TO PRINT OUT THE SOURCE PROGRAM, CHANGE
C THE JCL AS FOLLOWING:
C// EXEC FORTCG,PARM,FORT='NOSOURCE'
C.....................................................................
C PROGRAM "RESOIT",'RESONATOR EIGENMODE,USE TWO STEP TRANSFORM'.
C THE MOST RECENT USE OF THIS PROGRAM WAS ON AUG 21,1979.
C
C THE RESONATOR MUST BE CYLINDRICALLY SYMMETRIC, AND FHT (WITH
C LOWER END CORRECTION) IS EMPLOYED.
C NOTE: THE LEC IN FHT IS CORRECT FOR LB=0,1. BUT SLOPE CHANGE
C AT ORIGIN OF THE HUYGENS INTEGRAL IN THE SUBROUTINE PROP2A FOR
C E(R)*EXP(...R**2) IS CORRECT FOR LB=0 ONLY. ALSO ASSUME
C 1ST ORDER DERIVATIVE OF E(R) AT ORIGIN IS ZERO. FOR LB=1
C THE SLOPE CHANGES NEED TO BE MODIFIED.
C.....................................................................
C
C
DOUBLE PRECISION DSIN,DCOS,DSQRT,DEXP,DFLOAT
REAL*8 YR(2050),YI(2050),XM(2048),DALPH,PI,RC,AV
REAL*8 FN,A1,B1,D1,ZL,R1,R2,GAMMA1
REAL*8 T1,T2,T3,T4,GAMA,AFF,AFFQ
REAL*8 AXX,AXXQ,AXXP,XLAMDA,EPSILO,ED,AZ
REAL*8 RR,RJR,RJI,CF0,CF2
COMMON RR(2048),RJR(4096),RJI(4096),CF0(2048),CF2(2048)
C XM(2048) IS THE MIRROR REFLECTANCE FUNCTION.
C YR,YI IS THE FIELD FUNCTION WHICH PROPAGATES BACK AND FORTH
C INSIDE THE RESONATOR. THEY SHOULD BE DIMENSIONED AT LEAST TO
C (N+2), SINCE FUNCTION VALUE &amp; 2ND DERIVATIVES AT ORIGIN ARE
C STORED AT THE END OF THE ARRAYS.
C RC CORRESPONDS TO X=-0.5 IN GARDNER SPACE, WHICH IS THE
C INTEGRATION STARTING POINT.
C CF0,CF2 ARE CORRECTION ARRAYS DUE TO LOWER END CORRECTION
C FOR LB=0.
C
PI=3.1415926535898D0
C
N=1024
NA2=N+2
AK1=18.0
AK2=2.0
LB=0
NB=4
C
C NB IS BESSEL UPPER GENERATING ORDER.
C
CALL SETUPZ (N,LB,NB,AK1,AK2,DALPH,RRHO,BETAB,RC)
C
C
C GXX IS THE OUTPUT COUPLER RATIO= RR(N)/AXX, AXX IS THE OUTPUT
C COUPLER (MOST IN UNSTABLE RESONATOR) SIZE IN THE NUMERICAL SPACE.
C
GGX=6.0
C
C 6.0 IN THE GUARD BAND IS THE OVERALL WIDTH RATIO IN THE R DOMAIN,
C BUT THERE IS A MAGNIFICATION M=1.5, SO AXX, THE OUTPUT COUPLER
C WIDTH MUST BE SHRINKED BY M.***********
C
AXX=(RR(N)/(GGX*1.5))
AXX=AXX*2.0/1.5
C THE AXX IS CHANGED TO THE SLOWLY VARYING COS EDGE.
C
C FIND OUT KGG, S.T. RR(KGG)&lt;=AXX,RR(KGG+1)>AXX.
KGG=0
DO 80 K=1,N
IF(RR(K).GT.AXX) GOTO 80
KGG=KGG+1
80 CONTINUE
AXXP=AXX/2.D0
KGG1=0
DO 82 K=1,N
IF(RR(K).GT.AXXP) GOTO 82
KGG1=KGG1+1
82 CONTINUE
C
C FHT REALLY VIEWS THE OUTPUT COUPLER EDGE AS THE GEOMETRICAL MEAN
C OF RR(KGG) &amp; RR(KGG+1), SO AXX NEEDS TO BEADJUSTED AS FOLLOWING:
C AXX=RR(KGG)*DEXP(DALPH*0.5D0)
C AXXQ=AXX**2
C AXXQ IS (OUTPUT COUPLER EDGE)**2, IS USED LATER FOR PROPAGATIONS.
C
WRITE(6,8899)KGG,KGG1,RR(KGG),AXX,AXXP,GGX
8899 FORMAT(1H0,'*KGG=',I4,I4,'RR(KGG)=',D21.14,'AXX,COUPLER EDGE=',
1 D21.14,D13.5,'G=',E13.5)
C
KGGA=KGG+1
C WHEN REFLECTED FROM THE OUTPUT COUPLER,FIELD AT RR(K>=KGGA) ARE 0.
C
C FOLLOWING ARE INITIAL INPUT DATA.
C
DO 55 K=1,N
YR(K)=1.D0
YI(K)=1.D0
55 CONTINUE
C
C INPUT INITIAL WAVE IS ASSUMED TO BE REAL,PLANE WAVE WITH UNIFORM
C AMPLITUDE OF 1. SO THE FUNCTION VALUE AT ORIGIN IS 1 (REAL), 0
C (IMAG), 2ND DERIVATIVE AT ORIGIN IS ZERO FOR BOTH PARTS. (1ST
C DERIVATIVE AT ORIGIN IS ZERO ALL THE TIME OF THE RESONATOR
C CALCULATION DUE TO THE FACT LB=0.
C
YR(N+1)=1.D0
YR(N+2)=0.D0
YI(N+1)=0.D0
YI(N+2)=0.D0
C
C FOLLOWING ARE THE ROUND TRIP PARAMETERS. XLAMDA CAN BE GIVEN
C EITHER THE WAVELENGTH IN REAL DIMENSION, OR ANY NUMBER (SO THE
C RESONATOR IS SET IN THE NUMERICAL SPACE), AS LONG AS IT (WITH AXX
C AND A1,B1,D1) ENDS UP THE CORRECT COLLIMATED FRESNEL NUMBER AND
C EQUIVALENT FRESNEL NUMBER.
C
C ZL IS FREE SPACE PROPAGATION DISTANCE INSIDE THE RESONATOR,
C R1 IS THE OUTPUT COUPLER RADIUS, R2 IS THE RADIUS OF THE
C OTHER MIRROR (ASSUMED INFINITEIN SIZE-NO EDGE EFFECT).
C IN THIS EXAMPLE R1 IS CONVEX MIRROR, SO R1 IS NEGATIVE.
C R2 IS CONCAVE, SO R2 IS POSITIVE.
C
ZL=0.75D0
R1=-3.0D0
R2=4.5D0
C
C THE REFERENCE PLANE IS SET AT RIGHT BEFORE THE OUTPUT COUPLER.
C THE RAY VECTOR WILL START FROM THERE AND SEE MIRROR R1 FIRST,
C THEN THROUGH DISTANCE ZL, THEN MIRROR R2, THEN THROUGH DISTANCE
C ZL AGAIN TO RETURN TO THE STARTING REFERENCE PLANE.
C FOLLOWING ARE THE ROUND-TRIP A,B,C,D MATRIX:
C
A1=(1.D0-2.D0*ZL/R2)*(1.D0-2.D0*ZL/R1)-2.D0*ZL/R1
B1=ZL*(1.D0-2.D0*ZL/R2)+ZL
D1=1.D0-2.D0*ZL/R2
C
C THE WAVELENGTH IS SET TO 0.4163, COUPLED WITH AXX=B/(G*M)(B=B=8.46
C DUE TO N=1024, K1=18, K2=2, G=6.0) TO GIVE COLLIMATED FRESNEL
C NUMBER NC=M*AXX**2/(B1*XLAMDA)=2.547, EQUIVALENT FRESNEL NUMBER
C NEQ=(M**2-1)/(2*M**2)*NC=0.71. THE DESIGN GIVEN ABOVE IS A
C POSITIVE CONFOCAL UNSTABLE RESONATOR WITH M=1.5, SEE SIEGMAN'S PAPER
C AND THE PROGRAM DOCUMENTATIONS.
C
XLAMDA=0.416274D0
DO 84 K=1,KGG1
84 XM(K)=1.D0
KGG11=KGG1+1
DO 86 K=KGG11,KGG
AV=((RR(K)-RR(KGG11))/RR(KGG11))*PI/4.D0
XM(K)=DCOS(AV)
XM(K)=XM(K)**2
86 CONTINUE
DO 88 K=KGGA,N
88 XM(K)=0.D0
C
C THE ABOVE IS THE MIRROR FUNCTION.
C
C FOLLOWING IS ROUND-TRIP PROPAGATION.
C T2,T3 ARE THE ENERGY OF LAST ITERATION,PRESENT ITERATION OF THE
C FIELD AT RR(10), I.E.=YR(10)**2+YI(10)**2. NO SPECIAL REASON
C TO SELECT THE FIELD AT RR(10). AS LONG AS THE FIELD AT THIS POINT
C IS NOT ZERO (OR THE EIGENVALUE GAMA WILL BE OVERFLOW). THERE ARE
C SOME OTHER DELICATE WAY TO FIND GAMA, E.G. NORMALIZATION OR
C USING MORE THAN ONE POINT IN T2,T3.
C FOR INITIALIZATION, SET T2,T3 TO 1. NZ IS THE MAXIMUM NUMBER OF
C ITERATION. GAMMA1 IS LAST ITERATION EIGENVALUE. IF AZ=ABS(GAMA-
C GAMMA1)&lt;EPSILO, DESIRED CONVERGENCE IS REACHED AND ITERATION
C WILL STOP.
C
C
T2=1.D0
T3=1.D0
GAMMA1=1.D0
EPSILO=1.0D-4
NZ=60
NZ1=NZ
DO 180 NII=1,NZ
DO 89 KQ=1,N
YR(KQ)=YR(KQ)*XM(KQ)
89 YI(KQ)=YI(KQ)*XM(KQ)
C THE ABOVE IS THE MIRROR REFLECTION.
C
120 CALL PROP2A(N,NA2,DALPH,RC,LB,XLAMDA,A1,B1,D1,YR,YI)
IF(NII.LT.NZ1) GOTO 170
DO 65 KI=1,N
T1=YR(KI)**2+YI(KI)**2
WRITE(6,68) RR(KI),T1
65 CONTINUE
68 FORMAT(1H0,3(D16.9,2X))
C
170 T3=YR(10)**2+YI(10)**2
GAMA=DSQRT(T3/T2)
ED=GAMA-GAMMA1
AZ=DABS(ED)
WRITE (6,9992) NII,GAMA,T3,YR(10),YI(10),ED
9992 FORMAT (1H0,'KT=*',I4,1X,'GAMA',D17.10,'AMP**2',
1 D17.10,'YR(10)',D17.10,'YI(10)',D17.10,'ED=',D17.10)
IF(NZ1.LT.NZ) GOTO 60
IF(AZ.GT.EPSILO) GOTO 62
NZ1=1
C NZ1 IS SET TO NZ AT BEGINNING, ONCE THE PROGRAM DETECT THE
C AZ IS SMALLER THAN EPSILO, IT RESETS NZ1 TO 1, AND ENTER ONE
C MORE TIME ITERATION, SO THE INTERMEDIATE VALUE INSIDE THE
C RESONATOR (OTHER THAN THE REFERENCE PLANE) CAN BE RETRIEVED
C DURING THIS EXTRA TIME OF PROPAGATION (E.G. AFTER LINE 68,
C WE MAY PUT THE FOLLOWING COMMAND:
C
C IF(NZ1.GE.NZ) GOTO 170
C DO SOMETHING HERE
C AFTER THE CONVERGENCE
C IS DETECTED,
C E.G. IF I WANT TO KNOW
C THE EIGENMODE FIELD
C AT BEFORE MIRROR R2, I CAN
C DO A PROPAGATION FROM THE OUTPUT
C PLANE TO R1.(THIS EXTRA STEP IS
C NOT NEEDED WHEN WE DO THE ITERATION TO FIND EIGENMODE).
C
C AFTER THIS EXTRA ROUND TRIP PROPAGATION, THE PROGRAM STOPS.
C ** THE PROGRAM WILL STOP UNDER FOLLOWING TWO CONDITIONS:
C 1. NII>(NZ1=NZ), THE MAXIMUM ALLOWED ITERATION NUMBER IS REACHED.
C 2. NII>NZ1=1, WHEN ED&lt;EPSILO, CONVERGENCE REACHED.
62 GAMMA1=GAMA
T2=T3
C DO 190 KI=KGGA,N
C YR(KI)=0.D0
C190 YI(KI)=0.D0
C
C THIS IS OUTPUT COUPLER CUTOFF.
C
180 CONTINUE
60 STOP
END
C......................................................................
C
C
C
//
<p>C
C=====================================================================
C
C ***** DIFFRACTION PROGRAM USING STEVE SHENG'S FHT SUBROUTINES *****
C
C=====================================================================
C
C ***** STEVE SHENG'S ONE-STEP DIFFRACTION INTEGRAL PROGRAM *****
C
DOUBLE PRECISION DSIN,DCOS,DSQRT,DEXP,DFLOAT
REAL*8 DUM1,DUM2,CC1,CC2,DTHETA,PP,AX,BX,DX,SX,APX
REAL*8 WAL,XMAG1,XMAG2,THETA1,THETA2,EIG1,EIG2
REAL*8 XL1,XL2,XR1,XR2,XX,XXE,PHAR(1024),PHAI(1024)
REAL*8 PHDR(1024),PHDI(1024),TX
REAL*8 YR(1026),YI(1026),DALPH,PI,RC
REAL ITPOR,ALPH,RC1
REAL*8 RR,RJR,RJI,CF0,CF2
COMMON RR(1024),RJR(2048),RJI(2048),CF0(1024),CF2(1024)
C
C*****THIS PROGRAM IS DESIGNED FOR THE DIFFRACTION INTEGRAL
C AS SHOWN IN THE PAPER OR IN THE NOTE BOOK, WITH THE ONLY
C DIFFERENCE THAT THE EXTRA PHASE (T*J)**L IS NOT INCLUDED.
C ONE-STEP TRANSFORM IS APPLIED HERE(FHT). OTHER DIFFRACTION
C INTEGRAL FORM SUCH AS THE ONE STATED IN THE REVIEW PAPER
C BY PROF. A.E. SIEGMAN, WHERE VARIABLE TRANSFORMATIONS ARE
C CARRIED OUT WITH NC, M AS THE ONLY INDEPENDENT PARAMETERS,
C CAN ALSO BE CALCULATED BY THIS PROGRAM WITH MINOR MODIFI-
C CATION AND COORDINATES SCALING.
C*****
C
PI=3.1415926535898D0
C
N=1024
NA2=N+2
NM=N-1
AK1=10.0
AK2=2.0
LB=0
NB=3
C
C NB IS BESSEL UPPER GENERATING ORDER.
C
CALL SETUPZ (N,LB,NB,AK1,AK2,DALPH,RRHO,BETAB,RC)
C
WAL=1.15D-6
XL1=0.4057971008D0
XL2=XL1
XX=0.1739130433D0
XR1=2.D0*XX
XR2=2.D0*XL2-XR1
APX=1.D-3
C RR(NCUT+0.5)=XXE IS CUTOFF MIRROR EDGE.WAL IS OPERATING WAVE-
C LENGTH. XL1 IS LENGTH OF LEG1,XL2 IS LENGTH OF LEG2.
C XX IS THE CONFOCAL FOCUS LENGTH FROM MIRROR XR1. XR1 IS RADIUS
C OF CURVATURE OF MIRROR1, XR2 IS THAT OF MIRROR2.
C THE ABOVE GIVEN DATA IS SYMMETRIC CONFOCAL DESIGN:M=4/3,NEQ=-2.5
C TX IS SIGN OF BX,USED IN PHASE &amp; SLOPE MODIFICATIONS.
C
C FOLLOWING ARE GENERAL FORMULA FOR A,B,C,D MATRIX:
AX=(1.D0-2.D0*XL2/XR2)*(1.D0-2.D0*XL1/XR1)-2.D0*XL2/XR1
BX=(1.D0-2.D0*XL2/XR2)*XL1+XL2
DX=1.D0-2.D0*XL1/XR2
IF(BX.LT.0.D0)GOTO 5
TX=1.D0
C**HERE TX(SIGN OF BX) IS DETERMINED.
GOTO 7
5 TX=-1.D0
C
7 SX=APX/((DABS(BX)*WAL)**0.5)
C***** SX IS THE NORMALIZED RADIUS FOR FHT USE. LATER SX IS A
C SCALING FACTOR, WHICH HAS NO USE AT ALL FOR THIS (OR MOST)
C RESONATOR PROBLEM AND IS SET EQUAL TO ONE.
C*****
DO 6 K1=1,N
IF(RR(K1).GT.SX)GOTO 8
6 CONTINUE
8 NCUT=K1-1
NCC=K1
XXE=RR(NCUT)*DEXP(DALPH/2.D0)
WRITE(6,3990)NCUT
3990 FORMAT(' NCUT=',I4)
C
C***** RR(NCUT) IS THE CUTOFF EDGE IN NORMALIZED COORDINATES,
C XXE IS THE GEOMETRICAL MEAN OF RR(NCUT) &amp; RR(NCUT+1), WHICH
C THE CORRECT EDGE RADIUS FOR FRESNEL NUMBER CALCULATIONS.
C THIS RESULT CAME FROM THE PROPAGATION TESTS , IN WHICH
C TOPHAT INPUT FIELD IS TRANSFORMED BY FHT.
C*****
C
WRITE(6,2005)AX,BX,DX,XXE,SX
2005 FORMAT(' AX=',5(D22.15,1X))
C
C*****
C NO MATTER BX IS NEGATIVE OR POSITIVE,CALCULATION IN THIS
C PROGRAM IS CORRECT, BECAUSE IN THE RESONATOR EIGENMODE PART
C ONLY TX, SX(USE ABS(BX)) ARE USED, THE SIGN OF BX IS PROPERLY
C TAKEN CARE OF. IN THE FORWARD &amp; BACKWARD PROPAGATION PART
C THE RESCALING ONLY NEEDS ABS(BX).
C**WE MUST REMEMBER THE EXTRA PHASE (TJ)**L DEPENDS ON THE SIGN
C OF BX!!!
C*****
C
C*****FOLLOWING SX IS CHANGED TO 1.D0, WHICH IS A DUMMY(USELESS).
SX=1.D0
DO 10 K=1,N
DUM1=TX*SX**2*PI*AX*RR(K)**2
DUM2=TX*SX**2*PI*DX*RR(K)**2
PHAR(K)=DCOS(DUM1)
PHAI(K)=-DSIN(DUM1)
PHDR(K)=DCOS(DUM2)
10 PHDI(K)=-DSIN(DUM2)
C**FOLLOWING ARE INITIALIZATION OF EIGENVALUES &amp; INPUT FUNCTION:
XMAG1=1.D0
THETA1=0.D0
EIG1=1.D0
DO 12 K=1,N
YR(K)=100.D0
12 YI(K)=0.D0
YR(N+1)=YR(1)
YI(N+1)=0.D0
YR(N+2)=0.D0
YI(N+2)=0.D0
C CC1,CC2 USED IN FUNCTION VALUE &amp; SLOPES A.O. MODIFICATIONS.
PP=SX**2
CC1=TX*2.D0*PP*PI*AX
CC2=TX*2.D0*PP*PI*DX
C STARTS ITERATIONS!***
NZ=70
DO 100 KZ=1,NZ
DO 14 K=1,NCUT
DUM1=YR(K)
YR(K)=YR(K)*PHAR(K)-YI(K)*PHAI(K)
14 YI(K)=DUM1*PHAI(K)+YI(K)*PHAR(K)
DO 15 K=NCC,N
YR(K)=0.D0
15 YI(K)=0.D0
C
C FOLLOWING PRODUCTS ARE 2ND DERIV FOR LB=0:
C YR(N+1) ETC FUNCTION VALUE A.O. NOT CHANGE
YR(N+2)=YR(N+2)+CC1*YI(N+1)
YI(N+2)=YI(N+2)-CC1*YR(N+1)
C
C*****FOLLOWING PRODUCTS ARE 3RD DERIV FOR LB=1:
C YR(N+1) ETC 1ST DERIV A.O. NOT CHANGE.
C** YR(N+2)=YR(N+2)+3.D0*CC1*YI(N+1)
C** YI(N+2)=YI(N+2)-3.D0*CC1*YR(N+1)
CALL FHTZ1(YR,YI,N,NA2,DALPH,RC,LB)
C FOLLOWING ARE AFTER TRANSFORM PHASE &amp; VALUE,SLOPE MODIFICATION.
DO 16 K=1,N
DUM1=YR(K)
YR(K)=PP*(YR(K)*PHDR(K)-YI(K)*PHDI(K))
16 YI(K)=PP*(DUM1*PHDI(K)+YI(K)*PHDR(K))
C FOLLOWING PRODUCTS ARE 2ND DERIV FOR LB=0:
C YR(N+1) ETC FUNCTION VALUE A.O. ONLY*SX**2
YR(N+2)=PP*(YR(N+2)+CC2*YI(N+1))
YI(N+2)=PP*(YI(N+2)-CC2*YR(N+1))
YR(N+1)=PP*YR(N+1)
YI(N+1)=PP*YI(N+1)
C
C*****FOLLOWING PRODUCTS ARE 3RD DERIV FOR LB=1:
C YR(N+1) ETC 1ST DERIV A.O. ONLY*SX**2.
C** YR(N+2)=PP*(YR(N+2)+3.D0*CC2*YI(N+1))
C** YI(N+2)=PP*(YI(N+2)-3.D0*CC2*YR(N+1))
C** YR(N+1)=PP*YR(N+1)
C** YI(N+1)=PP*YI(N+1)
WRITE (6,1991) KZ,YR(670),YI(670),YR(671),YI(671)
1991 FORMAT (1H0,'KZ=',I4,1X,4(D25.18,1X))
C*****RR(670) IS PICKED FOR NO SPECIAL REASON, AS LONG AS
C THE FIELD IS NOT ZERO AT THIS POINT, IT IS ALL RIGHT
C TO CHECK THE EIGENVALUE &amp; CONVERGENCE THERE.
C*****
DUM1=DABS(YR(670))
IF(YI(670).GE.0.D0) GOTO 1881
DUM2=-1.D0
GOTO 1771
1881 DUM2=1.D0
1771 IF(DUM1.LE.1.0D-10) GOTO 18
THETA2=DATAN(YI(670)/YR(670))
IF(YR(670).GT.0.D0) GOTO 20
THETA2=THETA2+PI
GOTO 20
C IN 18 SIGN OF YI(670) DETERMINE THETA2:
18 THETA2=DUM2*PI/2.D0
20 IF(THETA2.GE.0.D0)GOTO 22
THETA2=THETA2+2.D0*PI
22 XMAG2=DSQRT(YR(670)**2+YI(670)**2)
EIG2=XMAG2/XMAG1
DTHETA=THETA2-THETA1
WRITE (6,1993) THETA2,XMAG2,EIG2,DTHETA
1993 FORMAT (1H0,'THETA2=',D25.18,'XMAG2=',
1 D18.11,'EIG2=',D18.11,'DTHETA=',D25.18)
DUM1=EIG1-EIG2
DUM1=DABS(DUM1)
WRITE(6,1995)DUM1
1995 FORMAT(' EIGENVALUE DIFFERENCE=',D25.18)
IF(DUM1.LE.1.0D-5) GOTO 30
EIG1=EIG2
THETA1=THETA2
XMAG1=XMAG2
100 CONTINUE
30 DO 32 K=1,N
PHDR(K)=YR(K)
32 PHDI(K)=YI(K)
DUM1=YR(N+1)
DUM2=YR(N+2)
XMAG1=YI(N+1)
XMAG2=YI(N+2)
SX=1.D0
WRITE(6,2003)SX
2003 FORMAT(' NEW SX=',D22.15)
C FOLLOWING TWO CASES,SX ARE THE SAME.
C FIRST RETREAT BY -XX: B=-XX IS NEGATIVE,SO TX=-1.D0
TX=-1.D0*TX
DO 34 K=1,N
PP=TX*SX**2*PI*1.D0*RR(K)**2*DABS(BX)/XX
PHAR(K)=DCOS(PP)
34 PHAI(K)=-DSIN(PP)
CC1=TX*2.D0*SX**2*PI*1.D0*DABS(BX)/XX
C** DABS(BX)/XX IS FOR PROPER RESCALING USE, SEE NOTES.
DO 36 K=1,N
EIG1=YR(K)
YR(K)=YR(K)*PHAR(K)-YI(K)*PHAI(K)
36 YI(K)=EIG1*PHAI(K)+YI(K)*PHAR(K)
C FOLLOWING FOR LB=0, 2ND DERIV A.O.:
YR(N+2)=YR(N+2)+CC1*YI(N+1)
YI(N+2)=YI(N+2)-CC1*YR(N+1)
C FOLLOWING FOR LB=1, 3RD DERIV A.O.:
C** YR(N+2)=YR(N+2)+3.D0*CC1*YI(N+1)
C** YI(N+2)=YI(N+2)-3.D0*CC1*YR(N+1)
CALL FHTZ1(YR,YI,N,NA2,DALPH,RC,LB)
PP=SX**4
DO 38 K=1,N
38 PHAR(K)=PP*(YR(K)**2+YI(K)**2)
DO 40 K=1,NCUT
YR(K)=PHDR(K)
40 YI(K)=PHDI(K)
YR(N+1)=DUM1
YR(N+2)=DUM2
YI(N+1)=XMAG1
YI(N+2)=XMAG2
DO 42 K=NCC,N
YR(K)=0.D0
42 YI(K)=0.D0
CALL FHTZ1(YR,YI,N,NA2,DALPH,RC,LB)
DO 44 K=1,N
44 PHAI(K)=PP*(YR(K)**2+YI(K)**2)
C PHAR(K) IS RETREAT -XX, PHAI(K) IS PROCEEDING TO COMING FOCUS.
DO 46 K=1,N
PP=PHDR(K)**2+PHDI(K)**2
WRITE(6,2001)RR(K),PHAR(K),PHAI(K),PP,PHDR(K),PHDI(K)
2001 FORMAT(4(D17.10,1X),2(D25.18,1X))
46 CONTINUE
WRITE(6,2007)DUM1,DUM2,XMAG1,XMAG2
2007 FORMAT(4(D25.18,1X))
STOP
END
C
C
//
<p>C
C=====================================================================
C
C ***** SETUPZ SUBROUTINE USED IN STEVE SHENG'S FHT PACKAGE *****
C
C=====================================================================
C
C
C
SUBROUTINE SETUPZ(N,LB,NB,AK1,AK2,DALPH,RRHO,BETAB,RC)
CCC---------------------------------------------------------------------
CCC INITIALIZES VARIOUS PARAMETERS INCLUDING COORDINATE ARRAY RR(N)
AND
CCC BESSEL ARRAY RJ(N4) NEEDED FOR FAST HANKEL TRANSFORM ROUTINES
CCC---------------------------------------------------------------------
C
REAL*8 RC,XC,XD,SCALE,WX
REAL*8 VD/0.0D0/,AD(1000),DBLE,DEXP,DALPH,DEALPH
REAL*8 RR,RJR,RJI,CF0,CF2,DUM1,DUM2,DUM3,DUM4,DUM5,DUM6
COMMON RR(1024),RJR(2048),RJI(2048),CF0(1024),CF2(1024)
C
C GENERATE FHT PARAMETER ALPH BY NEWTON'S METHOD
C
PI=3.1415926535898D0
AN = N
AK = AK1/AK2
Y = AN
DO 10 I=1,50
YP = (Y+AN)/(1.0+ALOG(AK*Y))
EPS = ABS((Y-YP)/Y)
IF (EPS.LT.1.0E-6) GO TO 20
10 Y = YP
WRITE (6,4010)
4010 FORMAT (' ----->WARNING - NEWTON METHOD FOR ALPH NOT CONVERGED')
20 ALPH = 1.0/Y
C
C GENERATE ADDITIONAL BASIC FHT PARAMETERS RRHO, BETAB, N2,N4
C
RRHO = AK2*ALPH/(AK1*AK1)
BETAB = 1.0/(AK2*ALPH)
R0 = SQRT(RRHO)
B = SQRT(BETAB)
N2 = 2*N
N4 = 4*N
C
C PRINT THE BASIC FHT PARAMETERS PLUS R0 AND B
C
WRITE (6,4020) N,LB,AK1,AK2,ALPH,RRHO,BETAB,R0,B
4020 FORMAT ('1----->SUBROUTINE SETUP EXECUTED:'/
# ' -----> N =',I5,2X,'LB =',I2,2X,'AK1 =',F5.2,2X,'AK2 =',F5.2/
# ' -----> ALPH =',1PE13.5,3X,'RRHO =',E13.5,3X,'BETAB =',
# E13.5/' -----> R0 =',1PE13.5,3X,'B =',1PE13.5)
C
C GENERATE THE RADIAL COORDINATE ARRAY RR(K) IN DOUBLE PRECISION.
C *** NOTE: THE PRECISION USED IN GENERATING RR(K) SHOULD BE THE
C SAME AS THAT IN RJR(K), OR THE SLIGHT SHIFT OF SAMPLING POINTS
C WILL CAUSE SIGNIFICANT ERROR.
C
RR(1) = SQRT(RRHO)
DALPH = DBLE(ALPH)
DEALPH = DEXP(DALPH)
RC=RR(1)/(DEALPH**0.5)
DO 30 K=2,N
30 RR(K) = DEALPH*RR(K-1)
RRNP1 = RR(N)*RR(N)/RR(N-1)
RR1=RR(1)
RRN=RR(N)
WRITE (6,4030) RR1,RRN,RRNP1
4030 FORMAT (' -----> RR(1) =',1PE13.5,3X,'RR(N) =',1PE13.5,
# 3X,'RRNP1 =',1PE13.5)
C
C GENERATE THE CORRECTION ARRAYS: CF0 IS EXACT, CF2 IS APPROXIMATE.
C
XD = 6.283185308D0*DBLE(RRHO)
XC=XD/(DEALPH**0.5)
C
C-*****-
C FOR DIFFERENT LB, MODIFICATIONS START HERE.
C NOW WE HAVE ONLY LB=0,1 UP TO TWO CORRECTION TERMS INCLUDED.
C
C FOLLOWING ARE CORRECTION ARRAYS FOR LB=0:
IF(LB.EQ.1) GOTO 42
LDIMA=1000
DO 45 K=1,N
CALL BESJ (XC,VD,NB,AD,LDIMA)
CF0(K)=AD(2)*RC
CF2(K)=PI*RC**4*RR(K)/2.D0-(PI*RC**2*RR(K))**3/3.D0
45 XC=XC*DEALPH
IF(LB.EQ.0) GOTO 47
C
C FOLLOWING ARE CORRECTION ARRAYS FOR LB=1:
42 DUM6=PI**2
DUM1=DUM6*RC**4/2.D0
DUM2=DUM6**2*RC**6/6.D0
DUM3=DUM6*RC**6/3.D0
DUM4=DUM6**2*RC**8/8.D0
C THE FORMULA OF CF0,CF2 FOR LB=1 ARE:
C CF0(K)=(PI*RR(K))**2*RC**4/2.D0-(PI*RR(K))**4*RC**6/6.D0
C CF2(K)=(PI*RR(K))**2*RC**6/3.D0-(PI*RR(K))**4*RC**8/8.D0
C
DO 44 K=1,N
DUM5=RR(K)**2
DUM6=RR(K)**4
CF0(K)=DUM1*DUM5-DUM2*DUM6
44 CF2(K)=DUM3*DUM5-DUM4*DUM6
C
C FOR DIFFERENT LB,MODIFICATIONS END HERE.
C-*****-
C
C GENERATE THE MODIFIED BESSEL FUNCTION ARRAY RJ(K)
C
47 DO 40 K=1,N2
X = SNGL(XD)
XO2PI=X/6.283185
LDIMA = MIN0(IFIX(1.25*(X+30.0)),1000)
CALL BESJ (XD,VD,NB,AD,LDIMA)
SCALE = DALPH*XD
RJR(K) = SCALE*AD(1+LB)/DFLOAT(N2)
RJI(K) = 0.0D0
40 XD = DEALPH*XD
C
C FOURIER TRANSFORM THE RJ(K) ARRAY
C
CALL DFFTS(RJR,RJI,N2,N2,N2,-1)
C
RETURN
C$PRINTON
END
C
C
&nbsp;
&nbsp;
&nbsp;
<p>C
C=====================================================================
C
C ***** FAST FOURIER TRANSFORM USED IN STEVE SHENG'S FHT PACKAGE *****
C
C=====================================================================
C
C
SUBROUTINE DFFTS(A,B,NTOT,N,NSPAN,ISN)
IMPLICIT REAL*8(A-H,O-Z)
C$PRINTOFF
C MULTIVARIATE COMPLEX FOURIER TRANSFORM, COMPUTED IN PLACE
C USING MIXED-RADIX FAST FOURIER TRANSFORM ALGORITHM.
C BY R. C. SINGLETON, STANFORD RESEARCH INSTITUTE, SEPT. 1968
C ARRAYS A AND B ORIGINALLY HOLD THE REAL AND IMAGINARY
C COMPONENTS OF THE DATA, AND RETURN THE REAL AND
C IMAGINARY COMPONENTS OF THE RESULTING FOURIER COEFFICIENTS.
C MULTIVARIATE DATA IS INDEXED ACCORDING TO THE FORTRAN
C ARRAY ELEMENT SUCCESSOR FUNCTION, WITHOUT LIMIT
C ON THE NUMBER OF IMPLIED MULTIPLE SUBSCRIPTS.
C THE SUBROUTINE IS CALLED ONCE FOR EACH VARIATE.
C THE CALLS FOR A MULTIVARIATE TRANSFORM MAY BE IN ANY ORDER.
C NTOT IS THE TOTAL NUMBER OF COMPLEX DATA VALUES.
C N IS THE DIMENSION OF THE CURRENT VARIABLE.
C NSPAN/N IS THE SPACING OF CONSECUTIVE DATA VALUES
C WHILE INDEXING THE CURRENT VARIABLE.
C THE SIGN OF ISN DETERMINES THE SIGN OF THE COMPLEX
C EXPONENTIAL, AND THE MAGNITUDE OF ISN IS NORMALLY ONE.
C A TRI-VARIATE TRANSFORM WITH A(N1,N2,N3), B(N1,N2,N3)
C IS COMPUTED BY
C CALL FFT(A,B,N1*N2*N3,N1,N1,1)
C CALL FFT(A,B,N1*N2*N3,N2,N1*N2,1)
C CALL FFT(A,B,N1*N2*N3,N3,N1*N2*N3,1)
C FOR A SINGLE-VARIATE TRANSFORM,
C NTOT = N = NSPAN = (NUMBER OF COMPLEX DATA VALUES), E.G.
C CALL FFT(A,B,N,N,N,1)
C THE DATA CAN ALTERNATIVELY BE STORED IN A SINGLE COMPLEX ARRAY C
C IN STANDARD FORTRAN FASHION, I.E., ALTERNATING REAL AND
C IMAGINARY PARTS. THEN WITH MOST FORTRAN COMPILERS, THE
C COMPLEX ARRAY C CAN BE EQUIVALENCED TO A REAL ARRAY A, THE
C MAGNITUDE OF ISN CHANGED TO TWO TO GIVE THE CORRECT INDEXING
C INCREMENT, AND A AND A(2) USED TO PASS THE INITIAL ADDRESSES
C FOR THE SEQUENCES OF REAL AND IMAGINARY VALUES, E.G.,
C COMPLEX C(NTOT)
C REAL A(2*NTOT)
C EQUIVALENCE (C(1),A(1))
C CALL FFT(A,A(2),NTOT,N,NSPAN,2)
C ARRAYS AT(MAXF), CK(MAXF), BT(MAXF), SK(MAXF), AND NP(MAXP)
C ARE USED FOR TEMPORARY STORAGE. IF THE AVAILABLE STORAGE
C IS INSUFFICIENT, THE PROGRAM IS TERMINATED BY A STOP.
C MAXF MUST BE .GE. THE MAXIMUM PRIME FACTOR OF N.
C MAXP MUST BE .GT. THE NUMBER OF PRIME FACTORS OF N.
C IN ADDITION, IF THE SQUARE-FREE PORTION K OF N HAS TWO OR
C MORE PRIME FACTORS, THEN MAXP MUST BE .GE. K-1.
DIMENSION A(1),B(1)
C ARRAY STORAGE IN NFAC FOR A MAXIMUM OF 15 PRIME FACTORS OF N.
C IF N HAS MORE THAN ONE SQUARE-FREE FACTOR, THE PRODUCT OF THE
C SQUARE-FREE FACTORS MUST BE .LE. 210
DIMENSION NFAC(11),NP(209)
C ARRAY STORAGE FOR MAXIMUM PRIME FACTOR OF 23
DIMENSION AT(23),CK(23),BT(23),SK(23)
EQUIVALENCE (I,II)
C THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH THE ARRAY DIMENSIONS.
MAXF=23
MAXP=209
IF(N .LT. 2) RETURN
INC=ISN
C72=0.30901699437494742D0
S72=0.95105651629515357D0
S120=0.86602540378443865D0
RAD=6.2831853071796D0
IF(ISN .GE. 0) GO TO 10
S72=-S72
S120=-S120
RAD=-RAD
INC=-INC
10 NT=INC*NTOT
KS=INC*NSPAN
KSPAN=KS
NN=NT-INC
JC=KS/N
RADF=RAD*DFLOAT(JC)*0.5D0
I=0
JF=0
C DETERMINE THE FACTORS OF N
M=0
K=N
GO TO 20
15 M=M+1
NFAC(M)=4
K=K/16
20 IF(K-(K/16)*16 .EQ. 0) GO TO 15
J=3
JJ=9
GO TO 30
25 M=M+1
NFAC(M)=J
K=K/JJ
30 IF(MOD(K,JJ) .EQ. 0) GO TO 25
J=J+2
JJ=J**2
IF(JJ .LE. K) GO TO 30
IF(K .GT. 4) GO TO 40
KT=M
NFAC(M+1)=K
IF(K .NE. 1) M=M+1
GO TO 80
40 IF(K-(K/4)*4 .NE. 0) GO TO 50
M=M+1
NFAC(M)=2
K=K/4
50 KT=M
J=2
60 IF(MOD(K,J) .NE. 0) GO TO 70
M=M+1
NFAC(M)=J
K=K/J
70 J=((J+1)/2)*2+1
IF(J .LE. K) GO TO 60
80 IF(KT .EQ. 0) GO TO 100
J=KT
90 M=M+1
NFAC(M)=NFAC(J)
J=J-1
IF(J .NE. 0) GO TO 90
C COMPUTE FOURIER TRANSFORM
100 SD=RADF/DFLOAT(KSPAN)
CD=2.0D0*DSIN(SD)**2
SD=DSIN(SD+SD)
KK=1
I=I+1
IF(NFAC(I) .NE. 2) GO TO 400
C TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)
KSPAN=KSPAN/2
K1=KSPAN+2
210 K2=KK+KSPAN
AK=A(K2)
BK=B(K2)
A(K2)=A(KK)-AK
B(K2)=B(KK)-BK
A(KK)=A(KK)+AK
B(KK)=B(KK)+BK
KK=K2+KSPAN
IF(KK .LE. NN) GO TO 210
KK=KK-NN
IF(KK .LE. JC) GO TO 210
IF(KK .GT. KSPAN) GO TO 800
220 C1=1.0D0-CD
S1=SD
230 K2=KK+KSPAN
AK=A(KK)-A(K2)
BK=B(KK)-B(K2)
A(KK)=A(KK)+A(K2)
B(KK)=B(KK)+B(K2)
A(K2)=C1*AK-S1*BK
B(K2)=S1*AK+C1*BK
KK=K2+KSPAN
IF(KK .LT. NT) GO TO 230
K2=KK-NT
C1=-C1
KK=K1-K2
IF(KK .GT. K2) GO TO 230
AK=C1-(CD*C1+SD*S1)
S1=(SD*C1-CD*S1)+S1
C1=2.0D0-(AK**2+S1**2)
S1=C1*S1
C1=C1*AK
KK=KK+JC
IF(KK .LT. K2) GO TO 230
K1=K1+INC+INC
KK=(K1-KSPAN)/2+JC
IF(KK .LE. JC+JC) GO TO 220
GO TO 100
C TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)
320 K1=KK+KSPAN
K2=K1+KSPAN
AK=A(KK)
BK=B(KK)
AJ=A(K1)+A(K2)
BJ=B(K1)+B(K2)
A(KK)=AK+AJ
B(KK)=BK+BJ
AK=-0.5D0*AJ+AK
BK=-0.5D0*BJ+BK
AJ=(A(K1)-A(K2))*S120
BJ=(B(K1)-B(K2))*S120
A(K1)=AK-BJ
B(K1)=BK+AJ
A(K2)=AK+BJ
B(K2)=BK-AJ
KK=K2+KSPAN
IF(KK .LT. NN) GO TO 320
KK=KK-NN
IF(KK .LE. KSPAN) GO TO 320
GO TO 700
C TRANSFORM FOR FACTOR OF 4
400 IF(NFAC(I) .NE. 4) GO TO 600
KSPNN=KSPAN
KSPAN=KSPAN/4
410 C1=1.0D0
S1=0
420 K1=KK+KSPAN
K2=K1+KSPAN
K3=K2+KSPAN
AKP=A(KK)+A(K2)
AKM=A(KK)-A(K2)
AJP=A(K1)+A(K3)
AJM=A(K1)-A(K3)
A(KK)=AKP+AJP
AJP=AKP-AJP
BKP=B(KK)+B(K2)
BKM=B(KK)-B(K2)
BJP=B(K1)+B(K3)
BJM=B(K1)-B(K3)
B(KK)=BKP+BJP
BJP=BKP-BJP
IF(ISN .LT. 0) GO TO 450
AKP=AKM-BJM
AKM=AKM+BJM
BKP=BKM+AJM
BKM=BKM-AJM
IF(S1 .EQ. 0) GO TO 460
430 A(K1)=AKP*C1-BKP*S1
B(K1)=AKP*S1+BKP*C1
A(K2)=AJP*C2-BJP*S2
B(K2)=AJP*S2+BJP*C2
A(K3)=AKM*C3-BKM*S3
B(K3)=AKM*S3+BKM*C3
KK=K3+KSPAN
IF(KK .LE. NT) GO TO 420
440 C2=C1-(CD*C1+SD*S1)
S1=(SD*C1-CD*S1)+S1
C1=2.0D0-(C2**2+S1**2)
S1=C1*S1
C1=C1*C2
C2=C1**2-S1**2
S2=2.0D0*C1*S1
C3=C2*C1-S2*S1
S3=C2*S1+S2*C1
KK=KK-NT+JC
IF(KK .LE. KSPAN) GO TO 420
KK=KK-KSPAN+INC
IF(KK .LE. JC) GO TO 410
IF(KSPAN .EQ. JC) GO TO 800
GO TO 100
450 AKP=AKM+BJM
AKM=AKM-BJM
BKP=BKM-AJM
BKM=BKM+AJM
IF(S1 .NE. 0) GO TO 430
460 A(K1)=AKP
B(K1)=BKP
A(K2)=AJP
B(K2)=BJP
A(K3)=AKM
B(K3)=BKM
KK=K3+KSPAN
IF(KK .LE. NT) GO TO 420
GO TO 440
C TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)
510 C2=C72**2-S72**2
S2=2.0D0*C72*S72
520 K1=KK+KSPAN
K2=K1+KSPAN
K3=K2+KSPAN
K4=K3+KSPAN
AKP=A(K1)+A(K4)
AKM=A(K1)-A(K4)
BKP=B(K1)+B(K4)
BKM=B(K1)-B(K4)
AJP=A(K2)+A(K3)
AJM=A(K2)-A(K3)
BJP=B(K2)+B(K3)
BJM=B(K2)-B(K3)
AA=A(KK)
BB=B(KK)
A(KK)=AA+AKP+AJP
B(KK)=BB+BKP+BJP
AK=AKP*C72+AJP*C2+AA
BK=BKP*C72+BJP*C2+BB
AJ=AKM*S72+AJM*S2
BJ=BKM*S72+BJM*S2
A(K1)=AK-BJ
A(K4)=AK+BJ
B(K1)=BK+AJ
B(K4)=BK-AJ
AK=AKP*C2+AJP*C72+AA
BK=BKP*C2+BJP*C72+BB
AJ=AKM*S2-AJM*S72
BJ=BKM*S2-BJM*S72
A(K2)=AK-BJ
A(K3)=AK+BJ
B(K2)=BK+AJ
B(K3)=BK-AJ
KK=K4+KSPAN
IF(KK .LT. NN) GO TO 520
KK=KK-NN
IF(KK .LE. KSPAN) GO TO 520
GO TO 700
C TRANSFORM FOR ODD FACTORS
600 K=NFAC(I)
KSPNN=KSPAN
KSPAN=KSPAN/K
IF(K .EQ. 3) GO TO 320
IF(K .EQ. 5) GO TO 510
IF(K .EQ. JF) GO TO 640
JF=K
S1=RAD/DFLOAT(K)
C1=DCOS(S1)
S1=DSIN(S1)
IF(JF .GT. MAXF) GO TO 998
CK(JF)=1.0
SK(JF)=0.0
J=1
630 CK(J)=CK(K)*C1+SK(K)*S1
SK(J)=CK(K)*S1-SK(K)*C1
K=K-1
CK(K)=CK(J)
SK(K)=-SK(J)
J=J+1
IF(J .LT. K) GO TO 630
640 K1=KK
K2=KK+KSPNN
AA=A(KK)
BB=B(KK)
AK=AA
BK=BB
J=1
K1=K1+KSPAN
650 K2=K2-KSPAN
J=J+1
AT(J)=A(K1)+A(K2)
AK=AT(J)+AK
BT(J)=B(K1)+B(K2)
BK=BT(J)+BK
J=J+1
AT(J)=A(K1)-A(K2)
BT(J)=B(K1)-B(K2)
K1=K1+KSPAN
IF(K1 .LT. K2) GO TO 650
A(KK)=AK
B(KK)=BK
K1=KK
K2=KK+KSPNN
J=1
660 K1=K1+KSPAN
K2=K2-KSPAN
JJ=J
AK=AA
BK=BB
AJ=0.0D0
BJ=0.0D0
K=1
670 K=K+1
AK=AT(K)*CK(JJ)+AK
BK=BT(K)*CK(JJ)+BK
K=K+1
AJ=AT(K)*SK(JJ)+AJ
BJ=BT(K)*SK(JJ)+BJ
JJ=JJ+J
IF(JJ .GT. JF) JJ=JJ-JF
IF(K .LT. JF) GO TO 670
K=JF-J
A(K1)=AK-BJ
B(K1)=BK+AJ
A(K2)=AK+BJ
B(K2)=BK-AJ
J=J+1
IF(J .LT. K) GO TO 660
KK=KK+KSPNN
IF(KK .LE. NN) GO TO 640
KK=KK-NN
IF(KK .LE. KSPAN) GO TO 640
C MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4)
700 IF(I .EQ. M) GO TO 800
KK=JC+1
710 C2=1.0D0-CD
S1=SD
720 C1=C2
S2=S1
KK=KK+KSPAN
730 AK=A(KK)
A(KK)=C2*AK-S2*B(KK)
B(KK)=S2*AK+C2*B(KK)
KK=KK+KSPNN
IF(KK .LE. NT) GO TO 730
AK=S1*S2
S2=S1*C2+C1*S2
C2=C1*C2-AK
KK=KK-NT+KSPAN
IF(KK .LE. KSPNN) GO TO 730
C2=C1-(CD*C1+SD*S1)
S1=S1+(SD*C1-CD*S1)
C1=2.0D0-(C2**2+S1**2)
S1=C1*S1
C2=C1*C2
KK=KK-KSPNN+JC
IF(KK .LE. KSPAN) GO TO 720
KK=KK-KSPAN+JC+INC
IF(KK .LE. JC+JC) GO TO 710
GO TO 100
C PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES
C PERMUTATION FOR SQUARE FACTORS OF N
800 NP(1)=KS
IF(KT .EQ. 0) GO TO 890
K=KT+KT+1
IF(M .LT. K) K=K-1
J=1
NP(K+1)=JC
810 NP(J+1)=NP(J)/NFAC(J)
NP(K)=NP(K+1)*NFAC(J)
J=J+1
K=K-1
IF(J .LT. K) GO TO 810
K3=NP(K+1)
KSPAN=NP(2)
KK=JC+1
K2=KSPAN+1
J=1
IF(N .NE. NTOT) GO TO 850
C PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE)
820 AK=A(KK)
A(KK)=A(K2)
A(K2)=AK
BK=B(KK)
B(KK)=B(K2)
B(K2)=BK
KK=KK+INC
K2=KSPAN+K2
IF(K2 .LT. KS) GO TO 820
830 K2=K2-NP(J)
J=J+1
K2=NP(J+1)+K2
IF(K2 .GT. NP(J)) GO TO 830
J=1
840 IF(KK .LT. K2) GO TO 820
KK=KK+INC
K2=KSPAN+K2
IF(K2 .LT. KS) GO TO 840
IF(KK .LT. KS) GO TO 830
JC=K3
GO TO 890
C PERMUTATION FOR MULTIVARIATE TRANSFORM
850 K=KK+JC
860 AK=A(KK)
A(KK)=A(K2)
A(K2)=AK
BK=B(KK)
B(KK)=B(K2)
B(K2)=BK
KK=KK+INC
K2=K2+INC
IF(KK .LT. K) GO TO 860
KK=KK+KS-JC
K2=K2+KS-JC
IF(KK .LT. NT) GO TO 850
K2=K2-NT+KSPAN
KK=KK-NT+JC
IF(K2 .LT. KS) GO TO 850
870 K2=K2-NP(J)
J=J+1
K2=NP(J+1)+K2
IF(K2 .GT. NP(J)) GO TO 870
J=1
880 IF(KK .LT. K2) GO TO 850
KK=KK+JC
K2=KSPAN+K2
IF(K2 .LT. KS) GO TO 880
IF(KK .LT. KS) GO TO 870
JC=K3
890 IF(2*KT+1 .GE. M) RETURN
KSPNN=NP(KT+1)
C PERMUTATION FOR SQUARE-FREE FACTORS OF N
J=M-KT
NFAC(J+1)=1
900 NFAC(J)=NFAC(J)*NFAC(J+1)
J=J-1
IF(J .NE. KT) GO TO 900
KT=KT+1
NN=NFAC(KT)-1
IF(NN .GT. MAXP) GO TO 998
JJ=0
J=0
GO TO 906
902 JJ=JJ-K2
K2=KK
K=K+1
KK=NFAC(K)
904 JJ=KK+JJ
IF(JJ .GE. K2) GO TO 902
NP(J)=JJ
906 K2=NFAC(KT)
K=KT+1
KK=NFAC(K)
J=J+1
IF(J .LE. NN) GO TO 904
C DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1
J=0
GO TO 914
910 K=KK
KK=NP(K)
NP(K)=-KK
IF(KK .NE. J) GO TO 910
K3=KK
914 J=J+1
KK=NP(J)
IF(KK .LT. 0) GO TO 914
IF(KK .NE. J) GO TO 910
NP(J)=-J
IF(J .NE. NN) GO TO 914
MAXF=INC*MAXF
C REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES
GO TO 950
924 J=J-1
IF(NP(J) .LT. 0) GO TO 924
JJ=JC
926 KSPAN=JJ
IF(JJ .GT. MAXF) KSPAN=MAXF
JJ=JJ-KSPAN
K=NP(J)
KK=JC*K+II+JJ
K1=KK+KSPAN
K2=0
928 K2=K2+1
AT(K2)=A(K1)
BT(K2)=B(K1)
K1=K1-INC
IF(K1 .NE. KK) GO TO 928
932 K1=KK+KSPAN
K2=K1-JC*(K+NP(K))
K=-NP(K)
936 A(K1)=A(K2)
B(K1)=B(K2)
K1=K1-INC
K2=K2-INC
IF(K1 .NE. KK) GO TO 936
KK=K2
IF(K .NE. J) GO TO 932
K1=KK+KSPAN
K2=0
940 K2=K2+1
A(K1)=AT(K2)
B(K1)=BT(K2)
K1=K1-INC
IF(K1 .NE. KK) GO TO 940
IF(JJ .NE. 0) GO TO 926
IF(J .NE. 1) GO TO 924
950 J=K3+1
NT=NT-KSPNN
II=NT-INC+1
IF(NT .GE. 0) GO TO 924
RETURN
C ERROR FINISH, INSUFFICIENT ARRAY STORAGE
998 ISN=0
PRINT 999
STOP
999 FORMAT(44H0ARRAY BOUNDS EXCEEDED WITHIN SUBROUTINE FFT)
C$PRINTON
END
<p>C
C=====================================================================
C
C ***** SUBROUTINE BESJ USED IN STEVE SHENG'S FHT PACKAGE *****
C
C=====================================================================
C
C
SUBROUTINE BESJ(X,V,N,A,LDIMA)
CCC------------------------------------------------------------------
CCC GENERATES BESSEL FUNCTIONS
CCC------------------------------------------------------------------
C$PRINTOFF
CC
CC BESSEL FUNCTIONS OF THE FIRST (J) OR SECOND (Y) KIND (ANL C370S)
CC (DOUBLE PRECISION) (CSD)
CC SHARE CODE: 20.50.90 TYPE: SUBPROGRAM
CC LANGUAGE - SOURCE: FORTRAN
CC CALLABLE FROM: PL/I
CC SOURCE - LIBRARY: SYS4.CSD.ANL.SRC
CC CARDS: 314
CC LOAD MODULE LIBRARY: NONE
CC CORE - MINIMUM REQUIRED: 0
CC AVAILABLE ON: CAMPUS/168
CC DOCUMENTATION: Argonne National Laboratory (ANL) Program
CC LIBRARY, Volumes 1, 2, and 3
CC
CC ENTRY BESJ OF SUBROUTINE CALCULATES THE J-BESSEL
CC FUNCTION OF X FOR ORDERS V,1+V,2+V,...,N+V
CC R(K) IS STORED IN A(K+2) (TEMPORARILY)
CC J(V+K) IS STORED INTO A(K+1)
CC X>0
CC 0&lt;=V&lt;1.D0
CC N>=0
CC DIMENSION OF ARRAY A MUST BE AT LEAST MAX(X,N)+16
CC LDIMA IS THE DIMENSION OF A SUPPLIED BY THE
CC USER.
CC
CC E R R O R R E T U R N S
CC FOR X, V, OR N OUT OF RANGE OR IF DIMENSION
CC OF ARRAY FURNISHED IS TOO SMALL.
CC
CC ERROR RETURN INFORMATION INCLUDES
CC X, V, N, MU, AND LDIMA WHERE MU IS THE
CC SIZE-2 OF THE ARRAY A NEEDED (MU IS
CC MEANINGLESS IF ANOTHER PARAMETER IS
CC OUT OF RANGE, E.G.,X=-5.)
CC
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 X,A,V
DIMENSION A(1)
DIMENSION GA(8),GB(7)
DATA PI/3.141592653589793D0/
DATA GA/5.563696434080738D-1,-1.549702804727049D1,
X-3.626839829961286D0,-1.127104620330990D1,
X8.701641788922069D0,-2.691911585444793D0,
X1.341319873241111D-1,-6.260629961830687D0/
DATA GB/1.137393573796284D1,4.798545220728675D1,
X6.087645511180565D0,1.055597915497223D2,
X-3.623508015808819D0,6.463221025538290D0,
X4.816923395585939D-1 /
DATA PI2/1.570796326794897D0/
EQUIVALENCE (GO,GOLRGV,GOSMLV)
DATA THRSHV/.789D-3/
DATA EC /Z4093C467E37DB0C8/
DATA S2/1.644934066848226/,
1 S2D2SQ/Z40AD2BF5ABCE72C3/,
2 S5/1.03692775514337D0/,
3 S33/.400685634386531/,
4 S4/1.082323233711138D0/
DATA PI4D72/Z4115A57EB579CE55/
DATA THSJ/.12D0/
DIMENSION GOVSMC(4)
ASSIGN 510 TO JY1
ASSIGN 110 TO JY4
80 CONTINUE
IF(X.LE.0.OR.V.LT.0.OR.V.GE.1.0.OR.N.LT.0) GO TO 995
C
C SET OVERFLOW INDICATOR OFF
C
C 87 CALL OVERFL(MU)
VP1=V+1.0D0
IX=X
LMB=IX+1
MU=LMB
GO TO JY4,(110,120)
110 CONTINUE
IF(LMB.LT.N) MU=N
120 CONTINUE
C
C LMB IS THE VALUE OF I FOR WHICH WE ARE ASSURED
C THAT J(V+I,X) IS ON THE TAIL OF THE FUNCTION,
C I.E.,LMB=IFIX(X+1)
C LET LL=MAX(LMB,N)...
C DEFINE MU TO BE THE POINT FROM WHICH WE
C MUST RECUR TO ASSURE THAT J(V+LL,X)
C IS ACCURATELY DETERMINED
C
C FOR I=LL,LL+1,...,MU
C WE STORE R(I)=J(I+1,X)/J(I,X) INTO A(I+2).
C THIS AVOIDS THE PROBLEM OF OVERFLOW
C INHERENT IN THE GROWTH OF THE J-FUNCTION
C WHEN RECURRING BACKWARD ON ITS TAIL
C NOTE THAT R(I-1)=X/(2*(V+I)-X*R(I))
TWODX=2.0D0/X
DR=TWODX*(V+DFLOAT(MU))
FKP1=1.D0
FK=0.D0
180 CONTINUE
C
C ITERATE UNTIL FKP1 IS GREATER THAN REGISTER ACCURACY
C
MU=1+MU
DR=DR+TWODX
FKM1=FK
FK=FKP1
FKP1=DR*FK-FKM1
IF ((FKP1+1.D0).NE.FKP1) GO TO 180
C
C THE VALUE OF MU IS NOW WELL DETERMINED
C
MU=MU+1
M=MU
IF((M+2).GT.LDIMA) GO TO 995
A(M+2)=0.D0
200 CONTINUE
IF (M.EQ.LMB) GO TO 250
DR=2*(M+V)
M=M-1
A(M+2)=X/(DR-X*A(M+3))
GO TO 200
C
C STORE JBAR(I) INTO A(I+1) FOR I=0,1,2,...,LMB+1
C
250 CONTINUE
A( M +1)=1.0D0/A( M +2)
A( M +2)=1.0D0
C
C RECUR BACKWARD TO GET JBAR'S
C
280 CONTINUE
DR=2*(M+V)/X
A(M)=DR*A(M+1)-A(M+2)
M=M-1
IF (M.GT.0) GO TO 280
C
C CLEAR R'S UNDERFLOW MAY OCCUR HERE
C
LMB=LMB+1
DO 290 M=LMB,MU
A(M+2)=A(M+2)*A(M+1)
C
C 287 CALL OVERFL(I)
C 288 IF (I.EQ.3) A(M+2)=0.D0
290 CONTINUE
C
C NORMALIZE SEQUENCE OF JBARS BY SUMMATION
C
IF (V.EQ.0.D0) GO TO 305
VX=V+2.0D0
BGAM=GA(1)+GB(1)/(VX+
X GA(2)+GB(2)/(VX+
X GA(3)+GB(3)/(VX+
X GA(4)+GB(4)/(VX+
X GA(5)+GB(5)/(VX+
X GA(6)+GB(6)/(VX+
X GA(7)+GB(7)/(VX+GA(8))))))))
BGAM=BGAM/VP1
BGAMSQ=BGAM*BGAM
Q2DXPV=(2.0D0/X)**V
305 CONTINUE
C
C SUMMATION
ALPHA=A(1)
PHI=2.0D0
IF (V.EQ.0.0D0) GO TO 320
D1=1.0D0
D2=V
EN2=V
EN1=V+2.0D0
PHI=Q2DXPV*BGAM
ALPHA=PHI*ALPHA
320 CONTINUE
DO 350 M=1, MU,2
IF (V.EQ.0.D0) GO TO 330
PHI=PHI*(EN2/D2)*(EN1/D1)
D2=D2+2.0D0
EN1=EN1+2.0D0
D1=1.0D0+D1
EN2=EN2+1.0D0
330 CONTINUE
ALPHA=ALPHA+PHI*A(M+2)
350 CONTINUE
C A(1),A(2),...,A(LMB+2) CONTAIN JBAR(0),JBAR(1),...,JBAR(LMB+1)
C A(LMB+3),...,A(MU+2) CONTAIN R(LMB+1),...,R(MU)
C
C NORMALIZE JBARS
C
M=MU+2
DO 460 I=1,M
457 A(I)=A(I)/ALPHA
460 CONTINUE
500 CONTINUE
GO TO JY1,(510,1100)
510 CONTINUE
RETURN
995 PRINT 1,X,V,N,MU,LDIMA
RETURN
1 FORMAT('0ERROR IN BESL X=',E14.5,' V=',E14.5,' N=',I5,' MU=',I5,
X' LDIMA=',I5)
ENTRY BESY(X,V,N,A,LDIMA)
ASSIGN 1100 TO JY1
ASSIGN 120 TO JY4
GO TO 80
1100 CONTINUE
VX=V
TJ=DABS(A(1))-THSJ
C
C COMPUTE GO AND G1
IF (V.LT.THRSHV) GO TO 1200
VX=V-1.0D0
IF ((VX+THRSHV).GT.0.D0) GO TO 1190
VX=V
GOLRGV=DCOTAN(V*PI)-(Q2DXPV**2/PI)*BGAMSQ/V
GO TO 1220
1190 CONTINUE
C
C FUNCTIONS OF V USED IN COMPUTING GO AND G1
C MUST BE TRANSFORMED TO FUNCTIONS OF VX=1-V
C
Q2DXPV=Q2DXPV/TWODX
BGAMSQ=BGAMSQ/V**2
1200 CONTINUE
C
C COMPUTE GO USING EXPANSION
Z=EC+DLOG(X/2.0D0)
GO= Z/PI2
IF (V.NE.0.0D0) GO TO 1210
G1=2.0D0/PI2
BGAMSQ=1.0D0
Q2DXPV=BGAMSQ
GO TO 1230
1210 CONTINUE
ASQ=Z*Z
ACUB=Z*ASQ
AFORTH=Z*ACUB
AFIFTH=AFORTH*Z/15.D0
G1=-VX
GOSMLV=
X (((((2.0D0*(AFIFTH+ACUB*S2/3.D0+ASQ*S33+Z*S2D2SQ)
X +Z*S4/2.D0+S2*S33+S5)*G1+
X 2.0D0*Z*S33+S2D2SQ+ASQ*S2
X +AFORTH/3.0D0+PI4D72)*G1+
X S33+Z*S2+ACUB/1.5D0 )*G1+
X 1.5D0*S2+ASQ )*G1+Z)/PI2
1220 CONTINUE
C COMPUTE G1
G1=(Q2DXPV**2/PI2)*BGAMSQ*(2.0D0+VX)/(1.0D0-VX)
1230 CONTINUE
C COMPUTE YO FROM SUM(J'S) FORM
C
EN3=VX+1.0D0
EN2=VX+EN3
EN1=VX+4.0D0
D1=2.0D0
D2=D1-VX
D3=D1+VX
TJE =0
IF(TJ.GE.0.D0.AND.VX.GE.0.D0) GO TO 1232
TJE=1
C
C Y(VX+1,X) MUST ALSO BE COMPUTED BY A SUM
C
THVDX=3.0D0*VX/X
PSIZ=-BGAMSQ*Q2DXPV**2/(PI2*X)
PSI1=GO-.5D0*G1
1232 CONTINUE
IF (VX.LT.0.D0) GO TO 1233
M=3
YV=GO*A(1)
IF (TJ.GE.0) GO TO 1238
YVP1=PSIZ*A(1)+PSI1*A(2)
GO TO 1238
1233 CONTINUE
Z=TWODX*V*A(1)-A(2)
YV=GO*Z
M=2
YVP1=PSIZ*Z+PSI1*A(1)
1238 CONTINUE
DO 1250 I=M,MU,2
YV=G1*A(I)+YV
G=G1
G1=-G1*(EN1/D1)*(EN2/D2)*(EN3/D3)
EN1=EN1+2.0D0
EN2=EN2+1.0D0
EN3=EN3+1.0D0
D1=D1+1.0D0
D2=1.0D0+D2
D3=D3+2.0D0
IF (TJE) 1240,1250,1240
1240 CONTINUE
YVP1=YVP1+THVDX*G*A(I)+.5*(G-G1)*A(I+1)
1250 CONTINUE
IF (VX.GE.0.D0) GO TO 1260
Z=YVP1
YVP1=V*Z*TWODX-YV
YV=Z
GO TO 1400
1260 CONTINUE
IF (TJ.LT.0.D0) GO TO 1400
1270 CONTINUE
C
C NOW COMPUTE Y(V+1)
C WRONSKIAN PROVIDED NOT NEAR A ZERO OF J
C
YVP1=(YV*A(2)-1.0D0/(X*PI2))/A(1)
1400 CONTINUE
C
C RECUR FORWARD TO GET Y'S (WISE?)
C
A(1)=YV
A(2)=YVP1
G=V*TWODX
DO 1500 I=2,N
C
G=G+TWODX
C OVERFLOW MAY OCCUR HERE
A(I+1)=G*A(I)-A(I-1)
C1487 CALL OVERFL(LMB)
C1488 IF (LMB.EQ.1) A(I+1)=1.D78
1500 CONTINUE
C$PRINTON
RETURN
END
<p>C
C=====================================================================
C
C ***** SUBROUTINE PROP2A USED IN STEVE SHENG'S FHT PACKAGE *****
C
C=====================================================================
C
C
SUBROUTINE PROP2A(N,NA2,DALPH,RC,LB,XLAMDA,A,B,D,YR,YI)
REAL*8 DALPH,RC,XLAMDA,A,B,D,YR(NA2),YI(NA2),DUM1,DUM2
REAL*8 PI,V1,V2,T1,T2,T3
REAL*8 RR,RJR,RJI,CF0,CF2
COMMON RR(2048),RJR(4096),RJI(4096),CF0(2048),CF2(2048)
PI=3.1415926535898D0
C THE FOLLOWING DERIVATIVE CHANGE AT ORIGIN IS TRUE ONLY FOR LB=0.
C
C MULTIPLY THE INPUT FIELD 2ND DERIVATIVE AT ORIGIN WITH A FACTOR.
C
DUM2=PI*(A-1.D0)/(B*XLAMDA)
YR(N+2)=YR(N+2)+2.D0*DUM2*YI(N+1)
YI(N+2)=YI(N+2)-2.D0*DUM2*YR(N+1)
C MODIFY INPUT FIELD.
DO 10 K=1,N
DUM1=RR(K)**2*DUM2
V1=DCOS(DUM1)
V2=DSIN(DUM1)
T1=YR(K)
YR(K)=YR(K)*V1+YI(K)*V2
YI(K)=-T1*V2+YI(K)*V1
10 CONTINUE
C
CALL FHTZ1(YR,YI,N,NA2,DALPH,RC,LB)
C MULTIPLY A CHIRP IN THE FREQUENCY DOMAIN.
DUM2=PI*B*XLAMDA
DO 20 K=1,N
DUM1=DUM2*RR(K)**2
V1=DCOS(DUM1)
V2=DSIN(DUM1)
T1=YR(K)
YR(K)=YR(K)*V1-YI(K)*V2
YI(K)=T1*V2+YI(K)*V1
20 CONTINUE
C MODIFY THE 2ND DERIVATIVE AT ORIGIN IN THE FREQUENCY DOMAIN, DUE
C TO THE CHIRP.
YR(N+2)=YR(N+2)-2.D0*DUM2*YI(N+1)
YI(N+2)=YI(N+2)+2.D0*DUM2*YR(N+1)
CALL FHTZ1(YR,YI,N,NA2,DALPH,RC,LB)
C
C MULTIPLY THE OUTPUT FIELD 2ND DERIVATIVE AT ORIGIN WITH A FACTOR.
DUM2=PI*(D-1.D0)/(B*XLAMDA)
YR(N+2)=YR(N+2)+2.D0*DUM2*YI(N+1)
YI(N+2)=YI(N+2)-2.D0*DUM2*YR(N+1)
C MODIFY OUTPUT FIELD.
DO 30 K=1,N
DUM1=RR(K)**2*DUM2
V1=DCOS(DUM1)
V2=DSIN(DUM1)
T1=YR(K)
YR(K)=YR(K)*V1+YI(K)*V2
YI(K)=-T1*V2+YI(K)*V1
30 CONTINUE
C
RETURN
END
C
<p>C
C=====================================================================
C
C ***** SUBROUTINE FHTZ1 USED IN STEVE SHENG'S FHT PACKAGE *****
C
C=====================================================================
C
C
SUBROUTINE FHTZ1(HRR,HII,N,NA2,DALPH,RC,LB)
C--------------------------------------------------------------------
C FAST-HANKEL-TRANSFORMS ARRAY HRR(N),HII(N) IN PLACE
C--------------------------------------------------------------------
C
C HANKEL-TRANSFORMS HRR,HII OVER K=1 TO N -- TRUNCATES TO ZEROS
C OVER K = N+1 TO N2 ON INPUT, BUT NOT ON OUTPUT -- INCLUDES
C ZERO-ORDER "LOWER-END CORRECTION"
C
C USE BESSEL-FUNCTION ARRAY RJR,RJI RADIAL COORDINATE ARRAY RR(N)
C
C USES SUBROUTINE DFFTS-- REQUIRES SUBROUTINE SETUPX TO BE CALLED TO
C INITIALIZE PARAMETERS BEFORE CALLING THIS SUBROUTINE
C
C ****-----**** HR,HI SHOULD BE DIMENSIONED AT LEAST 2*N.
C
REAL*8 PI/3.1415926535898D0/
REAL*8 HRR(NA2),HII(NA2),HR(2048),HI(2048),RC,DUM,DALPH
REAL*8 SFR01,SFR21,SFI01,SFI21,SFR0,SFR2,SFI0,SFI2
REAL*8 RR,RJR,RJI,CF0,CF2,Q
COMMON RR(1024),RJR(2048),RJI(2048),CF0(1024),CF2(1024)
C
C HR,HI SHOULD HAVE DIMENSION OF 2N.
C
N2=N*2
C
C TAKE THOSE DERIVATIVES OUT OF HRR,HII, AND CALCULATE
C THE TRANSFORMED DERIVATIVES AT THE ORIGIN.
C
C-*****-
C DIFFERENT LB,MODIFICATIONS OF "TRANSFORMED VALUE &amp; DERIV A.O."
C START HERE.
C
C IN THIS SUBROUTINE ONLY LB=0,1 ARE INCLUDED.
C FOR LB=0, SFR"0" ETC CORRESPOND TO FUNCTION VALUE A.O.,
C SFR"2" ETC CORRESPOND TO 2ND DERIV A.O.
C*****(FOR LB=0, IF INPUT FUNCTION DOES HAVE 1ST,3RD DERIV A.O.
C THE TRANSFORMED 1ST,3RD DERIV A.O. ALWAYS VANISH, BUT
C FUNCTION VALUE &amp; 2ND DERIV A.O. WILL CONTAIN TWO EXTRA
C TERMS FROM THOSE INPUT 1ST,3RD DERIV A.O., WHICH ARE
C "NOT" INCLUDED IN THE FOLLOWING CALCULATIONS)
C FOR LB=1: SFR"0" ETC CORRESPOND TO 1ST DERIV A.O.
C SFR"2" ETC CORRESPOND TO 3RD DERIV A.O.
C THE TRANSFORMED FUNCTION VALUE &amp; 2ND DERIV A.O. VANISH.
C*****(FOR LB=1, IF INPUT FUNCTION DOES HAVE NONZERO FUNCTION
C VALUE &amp; 2ND DERIV A.O., THE TRANSFORMED FUNCTION VALUE
C &amp; 2ND DERIV A.O. ALWAYS VANISH, BUT 1ST &amp; 3RD DERIV A.O.
C WILL INCLUDE TWO EXTRA TERMS FROM THOSE NONZERO INPUT
C FUNCTION VALUE &amp; 2ND DERIV A.O.)
C
SFR0=HRR(N+1)
SFI0=HII(N+1)
SFR2=HRR(N+2)
SFI2=HII(N+2)
SFR01=0.D0
SFI01=0.D0
SFR21=0.D0
SFI21=0.D0
C
IF(LB.EQ.1) GOTO 35
C
DO 622 K=1,N
SFR01=SFR01+RR(K)**2*HRR(K)
SFI01=SFI01+RR(K)**2*HII(K)
SFR21=SFR21+RR(K)**4*HRR(K)
SFI21=SFI21+RR(K)**4*HII(K)
622 CONTINUE
SFR01=2.D0*PI*(DALPH*SFR01+RC**2*SFR0/2.D0+RC**4*SFR2/8.D0)
SFI01=2.D0*PI*(DALPH*SFI01+RC**2*SFI0/2.D0+RC**4*SFI2/8.D0)
SFR21=-(2.D0*PI)**3*(DALPH*SFR21+RC**4*SFR0/4.D0+RC**6*SFR2/
1 12.D0)
SFR21=SFR21/2.D0
SFI21=-(2.D0*PI)**3*(DALPH*SFI21+RC**4*SFI0/4.D0+RC**6*SFI2/
1 12.D0)
SFI21=SFI21/2.D0
C
IF(LB.EQ.0) GOTO 37
35 DO 623 K=1,N
SFR01=SFR01+RR(K)**3*HRR(K)
SFI01=SFI01+RR(K)**3*HII(K)
SFR21=SFR21+RR(K)**5*HRR(K)
SFI21=SFI21+RR(K)**5*HII(K)
623 CONTINUE
Q=PI**2
SFR01=2.D0*Q*(DALPH*SFR01+RC**4*SFR0/4.D0+RC**6*SFR2/36.D0)
SFI01=2.D0*Q*(DALPH*SFI01+RC**4*SFI0/4.D0+RC**6*SFI2/36.D0)
SFR21=2.D0*PI**4*(DALPH*SFR21+RC**6*SFR0/6.D0+RC**8*SFR2/
1 48.D0)
SFR21=-3.D0*SFR21
SFI21=2.D0*PI**4*(DALPH*SFI21+RC**6*SFI0/6.D0+RC**8*SFI2/
1 48.D0)
SFI21=-3.D0*SFI21
C
C DIFFERENT LB, MODIFICATIONS END HERE.
C-*****-
C
C PATCHING HR,HI WITH ZEROS FROM N+1 TO 2N.
C
37 NP1=N+1
DO 10 K=NP1,N2
HR(K)=0.D0
10 HI(K)=0.D0
C
C
DO 5 K=1,N
HR(K)=RR(K)*HRR(K)
5 HI(K)=RR(K)*HII(K)
C
C FAST-HANKEL-TRANSFORM INPUT FUNCTION RF(K) USING RJ(K)
C
CALL DFFTS(HR,HI,N2,N2,N2,+1)
DO 30 K=1,N2
DUM = HR(K)
HR(K) = (DUM*RJR(K) - HI(K)*RJI(K))
30 HI(K) = (DUM*RJI(K)+HI(K)*RJR(K))
CALL DFFTS(HR,HI,N2,N2,N2,+1)
C
C-*****-
C DIFFERENT LB, MODIFICATION STARTS HERE.
C LOWER-END CORRECTIONS ARE TAKEN IN:
C CONTRIBUTION ONLY CORRECT FOR LB=0,1.
C
DUM=4.D0*DFLOAT(LB)
DO 60 K=1,N
HRR(K)=(HR(K)+CF0(K)*SFR0+CF2(K)*SFR2/(2.D0+DUM))/RR(K)
60 HII(K)=(HI(K)+CF0(K)*SFI0+CF2(K)*SFI2/(2.D0+DUM))/RR(K)
C
C FOR DIFFERENT LB, MODIFICATION ENDS HERE.
C-*****-
C
C
HRR(N+1)=SFR01
HII(N+1)=SFI01
HRR(N+2)=SFR21
HII(N+2)=SFI21
C
C END OF SUBROUTINE FHT
C
C 99 WRITE (6,901) N
C901 FORMAT (' ----->SUBROUTINE FHT EXECUTED: N =',I6)
RETURN
END
C