IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 AN1(0:200),BN1(0:200),AN2(0:200),BN2(0:200),M,CI, 1 T1,T2,T3 CHARACTER FQ*30,FS*30 DATA CI/(0.,1.)/ C C C CYL: CALCULATES EXTINCTION AND SCATTERING FROM AN INFINITE CYLINDER C Written by D. Mackowski dmckwski@eng.auburn.edu C C INPUT PARAMETERS: C X: SIZE PARAMETER C N,K: REAL AND IMAGINARY PARTS: M=N+IK C ALPHA: INCIDENT ANGLE OF RADIATION WITH RESPECT TO Z-AXIS (DEGREES). C NORMAL INCIDENCE: ALPHA=90 C NANG: NUMBER OF ANGLES AT WHICH SCATTERING IS CALCULATED C C OUTPUT: QE1,QS1: EXTINCTION AND SCATTERING EFFICIENCY FOR INCIDENT C RADIATION POLARIZED PARALLEL TO X-Z PLANE C QE2,QS2: EXTINCTION AND SCATTERING EFFICIENCY FOR INCIDENT C RADIATION POLARIZED PERPENDICULAR TO X-Z PLANE C T1, T2, T3: SCATTERING MATRIX ELEMENTS. REFER TO P. 202 OF C B&H FOR DETAILS. NORMALIZED WITH AVERAGE C SCATTERING EFFICIENCY. C PI=4.*DATAN(1.) 10 WRITE(*,'('' EXTINCTION OUTPUT FILE:'',\)') READ(*,'(A)') FQ IF(FQ.NE.' ') OPEN(1,FILE=FQ) WRITE(*,'('' NANG:'',\)') READ(*,*) NANG IF(NANG.NE.0) THEN WRITE(*,'('' SCATTERING OUTPUT FILE:'',\)') READ(*,'(A)') FS IF(FS.NE.' ') OPEN(2,FILE=FS) ENDIF 15 WRITE(*,'('' X,N,K,ALPHA:'',\)') READ (*,*) X,RN,RK,AD M=DCMPLX(RN,RK) CA=DCOS(AD*PI/180.) NT=NINT(X+4.*X**(1/3)+2) CALL CYLCOEF(X,M,CA,NT,AN1,BN1,AN2,BN2) QE1=.5*BN1(0) QS1=.5*BN1(0)*DCONJG(BN1(0)) QE2=.5*AN2(0) QS2=.5*AN2(0)*DCONJG(AN2(0)) DO 30 N=1,NT QE1=QE1+BN1(N) QE2=QE2+AN2(N) QS1=QS1+AN1(N)*DCONJG(AN1(N))+BN1(N)*DCONJG(BN1(N)) QS2=QS2+AN2(N)*DCONJG(AN2(N))+BN2(N)*DCONJG(BN2(N)) 30 CONTINUE QE1=4.*QE1/X QE2=4.*QE2/X QS1=4.*QS1/X QS2=4.*QS2/X QEA=.5*(QE1+QE2) QSA=.5*(QS1+QS2) WRITE(*,'('' QE,QSI,II:'',1P,4E11.3)') QE1,QS1,QE2,QS2 IF(FQ.NE.' ') WRITE(1,'(1P,5E11.3)') X,QE1,QS1,QE2,QS2 C IF(NANG.NE.0) THEN DO 60 I=0,NANG TD=I/DBLE(NANG)*180. T=PI*TD/180. T1=.5*BN1(0) T2=.5*AN2(0) T3=0. DO 40 N=1,NT CT=DCOS(N*T) ST=DSIN(N*T) T1=T1+CT*BN1(N) T2=T2+CT*AN2(N) T3=T3+ST*AN1(N) 40 CONTINUE S1=4*T1*DCONJG(T1)/QSA S2=4*T2*DCONJG(T2)/QSA S3=4*T3*DCONJG(T3)/QSA WRITE(*,'('' A,S1,S2,S3:'',F6.1,1P,3E11.3)') TD,S1,S2,S3 IF(FS.NE.' ') WRITE(2,'(1P,4E11.3)') TD,S1,S2,S3 60 CONTINUE ENDIF C WRITE(*,'('' MORE DATA(1),NEW FILE(2),OR QUIT(0):'',\)') READ(*,*) MORE IF(MORE.EQ.1) GOTO 15 IF(MORE.EQ.2) THEN IF(FQ.NE.' ')CLOSE(1) IF(NANG.NE.0) CLOSE(2) GOTO 10 ENDIF STOP END C SUBROUTINE CYLCOEF(X,M,CA,NT,AN1,BN1,AN2,BN2) IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 M,AN1(0:*),BN1(0:*),AN2(0:*),BN2(0:*),CI, 1 JN1(0:200),HN,HNP,DN1,ETA,FETA,AN,BN,CN,DN, 2 VN,WN,DEN REAL*8 JN(0:200),YN(0:200),JNP DATA CI/(0.,1.)/ C C CALCULATES THE MULTIPOLE COEFFICIENTS FOR EM SCATTERING BY C AN INFINITE CYLINDER, FOR S.P., R.I., AND COS(INC. ANG.) OF C X,M, AND CA. CALCULATES NT TERMS - COEF. RETURNED IN AN1, AN2, C BN1,BN2. REFER TO BOHREN & HUFFMAN FOR DETAILS. C SA=DSQRT((1.-CA)*(1.+CA)) ETA=X*CDSQRT(M*M-CA*CA) XI=X*SA FETA=CA*ETA*(XI*XI/ETA/ETA-1.) CALL JNREAL(XI,NT+1,JN,YN) CALL JNCOMP(ETA,NT+1,JN1) DO 40 N=0,NT DN1=(N*JN1(N)/ETA-JN1(N+1))/JN1(N) JNP=N*JN(N)/XI-JN(N+1) YNP=N*YN(N)/XI-YN(N+1) HN=JN(N)+CI*YN(N) HNP=JNP+CI*YNP DN=N*FETA*HN CN=N*FETA*JN(N) BN=XI*(M*M*XI*DN1*JN(N)-ETA*JNP) VN=XI*(M*M*XI*DN1*HN-ETA*HNP) WN=CI*XI*(ETA*HNP-XI*DN1*HN) AN=-CI*XI*(ETA*JNP-XI*DN1*JN(N)) DEN=WN*VN+CI*DN*DN AN1(N)=(CN*VN-BN*DN)/DEN BN1(N)=(WN*BN+CI*DN*CN)/DEN AN2(N)=(CI*CN*DN-AN*VN)/DEN BN2(N)=-CI*(CN*WN+AN*DN)/DEN 40 CONTINUE RETURN END C SUBROUTINE JNCOMP(Z,NS,JN) IMPLICIT REAL*8 (A-H,O-Z) COMPLEX*16 JN(0:*),A,Z C C CALCULATES THE INTEGER-ORDER BESSEL FUNCTION JN(Z) C N=0,1,...NS, FOR COMPLEX ARGUMENT Z. REFER TO BOHREN & HUFFMAN C FOR DETAILS. CODED 9/20/90 C ND=NINT((101+CDABS(Z))**.499+NS) ND=2*NINT(ND/2) JN(ND)=0. JN(ND-1)=1.D-32 A=0. DO 20 N=ND-1,3,-2 JN(N-1)=2.*N*JN(N)/Z-JN(N+1) JN(N-2)=2.*(N-1)*JN(N-1)/Z-JN(N) A=A+JN(N-1) 20 CONTINUE JN(0)=2.*JN(1)/Z-JN(2) A=JN(0)+2.*A DO 60 N=0,NS JN(N)=JN(N)/A 60 CONTINUE RETURN END C SUBROUTINE JNREAL(X,N,JN,YN) IMPLICIT REAL*8(A-H,O-Z) REAL*8 JN(0:*),YN(0:*) PARAMETER (IACC=40,BIGNO=1.D20,BIGNI=1.D-20) C JN(0)=BESSJ0(X) JN(1)=BESSJ1(X) YN(0)=BESSY0(X) YN(1)=BESSY1(X) C TOX=2./X DO 11 J=1,N-1 YN(J+1)=J*TOX*YN(J)-YN(J-1) 11 CONTINUE C IF(X.GT.FLOAT(N))THEN DO 21 J=1,N-1 JN(J+1)=J*TOX*JN(J)-JN(J-1) 21 CONTINUE ELSE M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2) JSUM=0 SUM=0. BJP=0. BJ=1. DO 22 J=M,1,-1 BJM=J*TOX*BJ-BJP BJP=BJ BJ=BJM IF(ABS(BJ).GT.BIGNO)THEN BJ=BJ*BIGNI BJP=BJP*BIGNI SUM=SUM*BIGNI DO 40 I=J+1,N JN(I)=JN(I)*BIGNI 40 CONTINUE ENDIF IF(JSUM.NE.0)SUM=SUM+BJ JSUM=1-JSUM IF((J.LE.N).AND.(J.GE.2)) JN(J)=BJP 22 CONTINUE SUM=2.*SUM-BJ DO 23 J=2,N JN(J)=JN(J)/SUM 23 CONTINUE ENDIF RETURN END C C THE FOLLOWING FUNCTIONS ARE FROM NUMERICAL RECIPES C REAL*8 FUNCTION BESSJ0(X) IMPLICIT REAL*8(A-H,O-Z) DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4, * -.2073370639D-5,.2093887211D-6/, Q1,Q2,Q3,Q4,Q5/-.1562499995D- *1, * .1430488765D-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/ DATA R1,R2,R3,R4,R5,R6/57568490574.D0,-13362590354.D0,651619640.7D *0, * -11214424.18D0,77392.33017D0,-184.9052456D0/, * S1,S2,S3,S4,S5,S6/57568490411.D0,1029532985.D0, * 9494680.718D0,59272.64853D0,267.8532712D0,1.D0/ IF(ABS(X).LT.8.)THEN Y=X**2 BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))) * /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) ELSE AX=ABS(X) Z=8./AX Y=Z**2 XX=AX-.785398164 BESSJ0=DSQRT(.636619772/AX)*(DCOS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y * *P5))))-Z*DSIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) ENDIF RETURN END C REAL*8 FUNCTION BESSY0(X) IMPLICIT REAL*8(A-H,O-Z) DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4, * -.2073370639D-5,.2093887211D-6/, Q1,Q2,Q3,Q4,Q5/-.1562499995D- *1, * .1430488765D-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/ DATA R1,R2,R3,R4,R5,R6/-2957821389.D0,7062834065.D0,-512359803.6D0 *, * 10879881.29D0,-86327.92757D0,228.4622733D0/, * S1,S2,S3,S4,S5,S6/40076544269.D0,745249964.8D0, * 7189466.438D0,47447.26470D0,226.1030244D0,1.D0/ IF(X.LT.8.)THEN Y=X**2 BESSY0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/(S1+Y*(S2+Y * *(S3+Y*(S4+Y*(S5+Y*S6)))))+.636619772*BESSJ0(X)*DLOG(X) ELSE Z=8./X Y=Z**2 XX=X-.785398164 BESSY0=DSQRT(.636619772/X)*(DSIN(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y* * P5))))+Z*DCOS(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) ENDIF RETURN END C REAL*8 FUNCTION BESSJ1(X) IMPLICIT REAL*8(A-H,O-Z) DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0,242396853.1D0 *, * -2972611.439D0,15704.48260D0,-30.16036606D0/, * S1,S2,S3,S4,S5,S6/144725228442.D0,2300535178.D0, * 18583304.74D0,99447.43394D0,376.9991397D0,1.D0/ DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5 *, * -.240337019D-6/, Q1,Q2,Q3,Q4,Q5/.04687499995D0,-.2002690873D-3 *, * .8449199096D-5,-.88228987D-6,.105787412D-6/ IF(ABS(X).LT.8.)THEN Y=X**2 BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))) * /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) ELSE AX=ABS(X) Z=8./AX Y=Z**2 XX=AX-2.356194491 BESSJ1=DSQRT(.636619772/AX)*(DCOS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y * *P5))))-Z*DSIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) ENDIF RETURN END C REAL*8 FUNCTION BESSY1(X) IMPLICIT REAL*8 (A-H,O-Z) DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5 *, * -.240337019D-6/, Q1,Q2,Q3,Q4,Q5/.04687499995D0,-.2002690873D-3 *, * .8449199096D-5,-.88228987D-6,.105787412D-6/ DATA R1,R2,R3,R4,R5,R6/-.4900604943D13,.1275274390D13,-.5153438139 *D11, * .7349264551D9,-.4237922726D7,.8511937935D4/, * S1,S2,S3,S4,S5,S6,S7/.2499580570D14,.4244419664D12, * .3733650367D10,.2245904002D8,.1020426050D6,.3549632885D3,1.D0/ IF(X.LT.8.)THEN Y=X**2 BESSY1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))/(S1+Y*(S2+Y* * (S3+Y*(S4+Y*(S5+Y*(S6+Y*S7))))))+.636619772 * *(BESSJ1(X)*DLOG(X)-1./X) ELSE Z=8./X Y=Z**2 XX=X-2.356194491 BESSY1=DSQRT(.636619772/X)*(DSIN(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y * *P5))))+Z*DCOS(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) ENDIF RETURN END