PROGRAM ASMAA C C THIS PROGRAM CALCULATES THE SCATTERING COEFFICIENTS , SCATTERED C INTENSITY,SCATTERING AND EXTINCTION CROSS SECTIONS , AND THE C MUELLER SCATTERING MATRIX ELEMENTS FOR THE SCATTERING OF ELECTRO- C MAGNETIC PLANE WAVES FROM A CYLINDER AT OBLIQUE INCIDENCE. C THE CYLINDER IS INFINITELY LONG AND OF CIRCULAR CROSS SECTION. C THE REFRACTIVE INDEX AND THE RADIUS OF THE CYLINDER ARE ARBITRARY. C THIS PROGRAM IS AN IMPROVED VERSION OF THE ONE THAT WAS PUBLISHED IN C COMPUTER PHYSICS COMMUNICATIONS, VOLUME 69, PP. 406-414, 1992 BY C HASHIM YOUSIF AND EDWARD BOUTROS. THE IMSL LIBRARY IS NO LONGER NEEDED C TO RUN THE PROGRAM, AND ALSO THE MEDIUM OUTSIDE THE CYLINDER IS NOLONGER C A FREE SPACE (ANY NONCONDUCTING MATERIAL). ANY PROBLEMS OR SUGGESTIONS C SHOULD BE DIRECTED TO C HASHIM YOUSIF C UNIVERSITY OF PITTSBURGH AT BRADFORD C BRADFORD, PA 16701 C YOUSIF@PITT.EDU C INTEGER ISIZE,NSIZE,IFF,IGG C C****************************************************************** PARAMETER (ISIZE=40000) PARAMETER (NSIZE=500) PARAMETER (IFF=1000,IGG=1000) C C ISIZE IS THE LEADING DIMENSION FOR BESSEL FUNCTIONS; NSIZE C IS THAT FOR THE SCATTERING ANGLE ******************************************************************* INTEGER NUMBER,ICOE,ISIJ,NBG,I REAL PI,WAVEL,RAD,REF1,REFM1,RPERB1 REAL WAVEL0,REF0,REM0 REAL DTHETA,ANGLEI,ANGLEF,THETA,WANO REAL AWANO,V,AFY,PERB,PERB1,EPSO,EO,HO REAL BJV(ISIZE),BYV(ISIZE),CPARA,CPERP REAL EXTM,EXTE,ANORM REAL I1(0:NSIZE),I2(0:NSIZE),I3(0:NSIZE) REAL I4(0:NSIZE),S11(0:NSIZE),ANGLE(0:NSIZE) COMPLEX CCC,WAN1,AWAN1,U,BJU(ISIZE) COMPLEX HV(ISIZE),BJUD(ISIZE),HVD(ISIZE) complex G(0:isize) COMPLEX AM(-ISIZE:ISIZE) COMPLEX BM(-ISIZE:ISIZE) COMPLEX CM(-ISIZE:ISIZE) COMPLEX DM(-ISIZE:ISIZE) COMPLEX T1(0:NSIZE),T4(0:NSIZE) COMPLEX T2(0:NSIZE),T3(0:NSIZE) COMPLEX TTM0(0:NSIZE),TTE0(0:NSIZE) C OPEN(UNIT=6,FILE='OBLI2.OUT',STATUS='NEW') C PI=4.0*ATAN(1.0) CCC=CMPLX(0.0,1.0) C C------------------------------------------------------ C-------THE INPUT PARAMETERS ARE READ FROM A CONSOLE--- C------------------------------------------------------ C WRITE(5,*)'HIT CR AFTER ENTERING EACH INPUT DATA' C C REFMED IS THE REFRACTIVE INDEX OF THE MEDIUM OUSTIDE THE CYLINDER, THIS C PARAMETER IS REAL WRITE(5,*)'WHAT IS THE REFRACTIVE INDEX OF THE OUTSIDE MEDIUM?' READ(5,*)REFMED C C WAVEL0 IS THE FREE SPACE WAVELENGTH OF THE INCIDENT RADIATION C IN METERS OR THE SAME UNITS AS RADIUS (RAD) C WRITE(5,*)'WHAT IS THE FREE-SPACE WAVELENGTH?' READ(5,*)WAVEL0 C C SCALING THE WAVELENGTH IS DONE IN THE NEXT STATEMENT WAVEL=WAVEL0/REFMED C C RAD IS THE RADIUS OF THE CYLINDER . RAD HAS C THE SAME UNIT AS WAVEL C WRITE(5,*)'WHAT IS THE RADIUS?' READ(5,*)RAD C C REF0 IS THE REAL PART OF THE REFRACTIVE C INDEX OF THE CYLINDER C WRITE(5,*)'WHAT IS THE REAL REFRACTIVE INDEX?' READ(5,*)REF0 C SCALING THE REAL REFRACTIVE INDEX REF1=REF0/REFMED C C REFM0 IS THE IMAGINARY PART OF THE REFRACTIVE C INDEX OF THE CYLINDER. IT MUST BE A POSITIVE C REAL NUMBER C WRITE(5,*)'WHAT IS THE IMAGINARY PART OF THE INDEX?' READ(5,*)REFM0 C C SCALING THE IMIGINARY REFRACTIVE INDEX OF THE CYLINDER C REFM1=REFM0/REFMED C C RPERB1 IS THE RELATIVE PERMEABILITY OF THE C CYLINDER WITH RESPECT TO THE OUTSIDE MEDIUM. C PERB1=1 MEANS THAT THE PERMEABILITY OF THE CYLINDER IS THE C SAME AS THE OUTSIDE MEDIUM. C THIS PARAMETER MUST BE A POSITIVE REAL NUMBER C WRITE(5,*)'WHAT IS THE RELATIVE PERMEABILITY?' READ(5,*)RPERB1 C C DTHETA IS THE ANGLE BETWEEN THE INCIDENT FIELD C AND THE AXIS OF THE CYLINDER IN DEGREES (SEE FIG. 1) C WRITE(5,*)'WHAT IS THE TILT ANGLE?' READ(5,*)DTHETA C C ANGLEI IS THE INITIAL SCATTERING ANGLE IN DEGREES. C THE FORWARD SCATTERING ANGLE CORRESPONDS TO ANGLEI=0. C WRITE(5,*)'WHAT IS THE INITIAL SCATTERING ANGLE?' READ(5,*)ANGLEI C C ANGLEF IS THE FINAL SCATTERING ANGLE IN DEGREES. C WRITE(5,*)'WHAT IS THE FINAL SCATTERING ANGLE?' READ(5,*)ANGLEF C C NUMBER IS THE NUMBER OF INTERVALS BETWEEN ANGLEI C AND ANGLEF C WRITE(5,*)'WHAT IS THE NUMBER OF INTERVALS?' READ(5,*)NUMBER C C ICOE IS AN INTEGER, IF IT IS SET TO ZERO THE SCATTERING C COEFFICIENTS WILL NOT BE PRINTED C WRITE(5,*)'IF YOU DO NOT WANT THE SCATTERING' WRITE(5,*)'COEFFICIENTS TO BE PRINTED ENTER 0' WRITE(5,*)'OTHERWISE ENTER ANY NUMBER' READ(5,*)ICOE C C ISIJ IS AN INTEGER, IF IT IS SET TO ZERO SIJ AND THE MAGNITUDE C SQUARED OF THE SCATTERED AMPLITUDES WILL NOT BE PRINTED C WRITE(5,*)'IF YOU DO NOT WANT SIJ AND THE MAGNITUDE' WRITE(5,*)'SQUARED OF THE SCATTERED AMPLITUDES TO BE' WRITE(5,*)'PRINTED ENTER 0 OTHERWISE ENTER ANY NUMBER' READ(5,*)ISIJ C_________________________________________________________ C_______THE END OF THE INPUT PARAMETERS___________________ C C THETA IS THE SAME AS DTHETA BUT IN RADIANS C THETA=DTHETA*PI/180.0 C C WANO IS THE WAVE NUMBER IN THE MEDIUM OUTSIDE THE CYLINDER. C WANO=2.0*PI/WAVEL C C AWANO IS THE LEFT HAND SIDE OF EQ.(4) C AWANO=WANO*SIN(THETA) C C V IS GIVEN BY EQ.(14) C V=AWANO*RAD C C WAN1 IS THE WAVE NUMBER INSIDE THE CYLINDER C GIVEN BY EQ.(11) C WAN1=WANO*CMPLX(REF1,-REFM1) C C AWAN1 IS THE AXIAL WAVE NUMBER INSIDE THE CYLINDER C AWAN1=SQRT(WAN1**2-(WANO*COS(THETA))**2) C C U IS GIVEN BY EQ.(13) C U=AWAN1*RAD C C AFY IS THE ANGULAR FREQUENCY OF THE INCIDENT C RADIATION C AFY=2.0*PI*3.0 E+08/WAVEL C C PERB IS THE PERMEABILITY OF FREE SPACE C PERB=PI*4.0 E-07 C C PERB1 IS THE PERMEABILITY OF THE CYLINDER C PERB1=RPERB1*PERB C C EPSO IS THE PERMITTIVITY OF FREE SPACE C EPSO=1.0 E-09/(36.0*PI) C C NBG IS THE NUMBER OF TERMS WHICH IS NECESSARY C FOR THE SERIES SOLUTION TO CONVERGE GIVEN BY EQ.(31) C NBG=V+4.0*(V)**0.334+2 C C EO IS THE AMPLITUDE OF THE INCIDENT ELECTRIC FIELD C PARALLEL TO THE INCIDENT PLANE. C EO=1.0 C C HO IS THE AMPLITUDE OF THE INCIDENT MAGNETIC FIELD C PARALLEL TO THE INCIDENT PLANE. C HO=1.0 C C IT SHOULD BE EMPHASISED HERE THAT THE SCATTERING COEFFICIENTS ARE C CALCULATED UNDER THE ASSUMPTION THAT BOTH EO AND HO ARE 1. HOWEVER, C THE SCATTERED AMPLITUDES ARE CALCULATED UNDER THE ASSUMPTION THAT BOTH C THE INCIDENT ELECTRIC FIELDS PARALLEL AND PERPENDICULAR TO THE INCIDENT C PLANE IS 1. MOREOVER, THE CROSS SECTIONS ARE CALCULATED UNDER THE USUAL C ASSUMPTION THAT THE INCIDENT INTENSITY IS ONE FOR EACH POLARIZATION C C THE FOLLWING CALL TO SUBROUTINE BESSEL_START DETRMINE THE INITIAL C INDEX (MST) TO START THE RECURRENCE RELATION C CALL BESSEL_START(V,NBG+2,MST) C C THE FOLLOWING CALL TO SUBROUTINE LOGDER_START IS TO DETRMINE THE C STARTING INDEX (NMX) TO START RECURRENCE RELATION C CALL LOGDER_START(U,NBG,NMX) C THE FOLLOWING CALL ON SUBROUTINE CHECK ENSURES THAT ENOUGH SPACES C ARE ALLOCATED TO THE ARRAYS C CALL CHECK(ISIZE,NSIZE,NBG,NUMBER,MST,IFF,IGG,NMX) C CALL BESSEL(v,bjv,byv,isize,nbg+2,IFF,MST) C CALL DRLOG(U,NBG,G,ISIZE,IGG,NMX) C------------------------------------------------------------- WRITE(6,10) 10 FORMAT(10X,'PROGRAM OUTPUT',/) WRITE(6,20)RAD 20 FORMAT(5X,'RADIUS=',E13.7,/) WRITE(6,30)WAVEL0 30 FORMAT(5X,'WAVELENGTH=',E13.7,/) WRITE(6,31)REFMED 31 FORMAT(5X,'REFRACTIVE INDEX OF THE OUSIDE MEDIUM=',E13.7,/) WRITE(6,40)DTHETA 40 FORMAT(5X,'TILT ANGLE=',E13.7,/) WRITE(6,50)WANO*RAD 50 FORMAT(5X,'SIZE PARAMETER=',E13.7,/) WRITE(6,60)REF0,REFM0 60 FORMAT(5X,'REFRACTIVE INDEX=',E13.6,1X,',-',E13.6,/) WRITE(6,70)RPERB1,NBG 70 FORMAT(5X,'RELATIVE PERMEABILITY=',E13.7,2X,'NBG=',I10,//) C C CONSTRUCTING HANKEL FUNCTION OF THE SECOND KIND (HV) USING C EQ.(27) C DO 80 I=1,NBG+2 HV(I)=CMPLX(BJV(I),-BYV(I)) 80 CONTINUE C C THE FOLLOWING CALL TO SUBROUTINE CDERV CALCULATES THE C FIRST DERIVATIVES OF HANKEL FUNCTIONS OF THE SECOND C KIND (HVD) AND OF THE REAL ARGUMENT (V) C CALL CDERV(HVD,HV,NBG,ISIZE) C C THE SCATTERING COEFFICIENTS FOR THE TM POLARIZATION C (AM AND BM) ARE CALCULATED USING THE FOLLOWING CALL C TO SUBROUTINE COEFF CALL COEFF(AM,BM,V,U,WANO,WAN1,THETA,NBG,AFY 6 ,BJV,HV,HVD,G,EPSO,PERB,PERB1,EO,1,ISIZE) C C THE SCATTERING COEFFICIENTS FOR THE TE POLARIZATION C (CM AND DM) ARE CALCULATED USING THE FOLLOWING CALL C TO SUBROUTINE COEFF C CALL COEFF(CM,DM,V,U,WANO,WAN1,THETA,NBG,AFY 6 ,BJV,HV,HVD,G,EPSO,PERB,PERB1,HO,2,ISIZE) C IF(ICOE.NE.0)THEN WRITE(6,90) 90 FORMAT(7X,'ORDER',15X,'AM',30X,'BM') DO 100 I=-NBG,NBG WRITE(6,120)I,AM(I),BM(I) 100 CONTINUE WRITE(6,110) 110 FORMAT(/,7X,'ORDER',15X,'CM',30X,'DM') DO 130 I=-NBG,NBG WRITE(6,120)I,CM(I),DM(I) 120 FORMAT(3X,I7,5X,2(E13.7,2X),5X,2(E13.7,2X)) 130 CONTINUE ENDIF C C THE FOUR FOLLOWING CALLS TO SUBROUTINE AMPLITUDE C CALCULATES T1,T4,T2, AND T3 ,RESPECTIVELY (SEE C EQS.(16,17,23, AND 24)). C CALL AMPLITUDE(T1,NBG,1.0,AM 6 ,ANGLEI,ANGLEF,NUMBER,THETA,ISIZE,NSIZE) CALL AMPLITUDE(T4,NBG,PERB*AFY/WANO,BM 6 ,ANGLEI,ANGLEF,NUMBER,THETA,ISIZE,NSIZE) CALL AMPLITUDE(T2,NBG,1.0,CM 6 ,ANGLEI,ANGLEF,NUMBER,THETA,ISIZE,NSIZE) CALL AMPLITUDE(T3,NBG,WANO/(PERB*AFY),DM 6 ,ANGLEI,ANGLEF,NUMBER,THETA,ISIZE,NSIZE) C C THE FOLLOWING CALL TO SUBROUTINE CROSS CALCULATES C THE SCATTERING CROSS SECTION PER UNIT LENGTH FOR C THE INCIDENT TM POLARIZATION (CPARA) USING EQ.(18) C CALL CROSS(CPARA,AM,BM,1,NBG,PERB,EPSO,AWANO 6 ,AFY,ISIZE) C C THE FOLLOWING CALL TO SUBROUTINE CROSS CALCULATES C THE SCATTERING CROSS SECTION PER UNIT LENGTH FOR THE C INCIDENT TE POLARIZATION (CPERP) USING EQ.(25) C CALL CROSS(CPERP,DM,CM,2,NBG,PERB,EPSO,AWANO 6 ,AFY,ISIZE) C C THE FOLLOWING CALL TO SUBROUTINE AMPLITUDE CALCULATES THE C COMPLEX SCATTERED AMPLITUDE IN THE FORWARD DIRECTION FOR C THE TM INCIDENT POLARIZATION (TTM0). C CALL AMPLITUDE(TTM0,NBG,1.0,AM,0.0,180.0,0 6 ,THETA,ISIZE,NSIZE) C C CALCULATING THE EXTINCTION CROSS SECTION PER UNIT LENGTH C FOR THE INCIDENT TM POLARIZATION (EXTM) USING THE OPTICAL C THEOREM USING EQ.(19) C EXTM=-4.0/WANO*REAL(TTM0(0)) C C THE FOLLOWING CALL TO SUBROUTINE AMPLITUDE CALCULATES THE C COMPLEX SCATTERED AMPLITUDE IN THE FORWARD DIRECTION FOR C THE INCIDENT TE POLARIZATION (TTE0) C CALL AMPLITUDE(TTE0,NBG,1.0,CM,0.0,180.0,0 6 ,THETA,ISIZE,NSIZE) C C THE EXTINCTION CROSS SECTION PER UNIT LENGTH FOR THE C TE INCIDENT POLARIZATION (EXTE) IS CALCULATED USING C THE OPTICAL THEOREM C EXTE=-4.0/WANO*REAL(TTE0(0)) C C EFFICIENCY IS CROSS SECTION PER UNIT AREA C WRITE(6,140)CPARA/(2.0*RAD),EXTM/(2.0*RAD) 140 FORMAT(3X,'TM SCATTERING AND EXTINCTION EFFICIENCIES=', 6 2(E13.7,2X)) WRITE(6,150)CPERP/(2.*RAD),EXTE/(2.0*RAD) 150 FORMAT(3X,'TE SCATTERING AND EXTINCTION EFFICIENCIES=', 6 2(E13.7,2X),//) C C ANORM IS S11 IN THE FORWARD DIRECTION; THIS QUANTITY CALCULATED C SEPARATELY BECAUSE THE INITIAL SCATTERING ANGLE MAY OR MAY NOT C BE ZERO. C ANORM=0.5*(TTM0(0)*CONJG(TTM0(0))+TTE0(0)*CONJG(TTE0(0))) C C THE FOLLOWING CALL TO SUBROUTINE SIJ CALCULATES THE MAGNITUDE C SQUARED OF THE SCATTERED AMPLITUDES AND THE MUELLER SCATTERING C MATRIX ELEMENTS (SIJ) IF THE PARAMETER ISIJ IS SET TO ANY NUMBER C OTHER THAN ZERO. C IF(ISIJ.NE.0)THEN CALL SIJ(T1,T2,T3,T4,I1,I2,I3,I4,ANGLEI,ANGLEF, 6 NUMBER,ANORM,NSIZE,ANGLE,S11) ENDIF STOP END c SUBROUTINE BESSEL_START(X,NORDER,NST) C C THIS SUBROUTINE DETRMINES THE STATING INDEX OF THE DOWNWARD C RECURRENCE RELATION, SEE THE FOLLOWING REFERENCE FOR MORE DETAILS: C COMPUTER PHYSICS COMMUNICATIONS, VOL. 106 PP. 199-206, 1997 REAL X INTEGER NORDER,N1,N2,NST C IF(NORDER.NE.1)THEN N1=10+X+8.41*(X)**0.393 N2=20+NORDER+4.7*NORDER**0.24 NST=MAX(N1,N2) ELSE NST=10+X+8.41*X**0.393 ENDIF RETURN END C SUBROUTINE LOGDER_START(ARG,NN,NMX) C THIS SUBROUTINE DETRMINES THE STATING INDEX OF THE DOWNWARD C RECURRENCE RELATION, SEE THE FOLLOWING REFERENCE FOR MORE DETAILS: C COMPUTER PHYSICS COMMUNICATIONS, VOL. 106 PP. 199-206, 1997 C COMPLEX ARG INTEGER NBG,NN N1=CABS(ARG) NMX=1.1*MAX(NN,N1)+15 RETURN END C SUBROUTINE CDERV(XOUT,XIN,NSTOP,ISIZE) C C THIS SUBROUTINE CALCULATES THE FIRST DERIVATIVES C OF ANY BESSEL FUNCTIONS .IT USES EQ.(28) FOR THAT. C NSTOP IS THE HIGHEST ORDER OF BESSEL FUNCTION FOR C WHICH THE FIRST DERIVATIVE IS COMPUTED. C INTEGER ISIZE,NSTOP,I COMPLEX XOUT(ISIZE),XIN(ISIZE) C DO 10 I=1,NSTOP+1 IF(I.EQ.1)XOUT(I)=-XIN(I+1) IF(I.NE.1)THEN XOUT(I)=(XIN(I-1)-XIN(I+1))/2.0 END IF 10 CONTINUE RETURN END C SUBROUTINE COEFF(AS,BS,V,U,WANO,WAN1,THETA,NBG,AFY 6 ,BJV,HV,HVD,G,EPSO,PERB,PERB1,AMPL,MODE,ISIZE) C C THIS SUBROUTINE CALCULATES THE SCATTERING COEFFICIENTS FOR THE C TM POLARIZATION IN THE FIRST CALL AND TE SCATTERING COEFFICIENTS C IN THE SECOND CALL. C INTEGER ISIZE,M,NBG,JJ,MODE REAL PI,V,WANO,BJV(ISIZE),PERB1,PERB REAL AMPL,THETA,EPSO,AFY COMPLEX CCC,QM,U,G(0:ISIZE) COMPLEX SM,HVD(ISIZE),HV(ISIZE) COMPLEX R,N,WAN1,X,AZ,D COMPLEX AS(-ISIZE:ISIZE) COMPLEX BS(-ISIZE:ISIZE) C CCC=CMPLX(0.0,1.0) PI=4.0*ATAN(1.0) C DO 10 M=0,NBG JJ=M+1 C C QM AND SM ARE GIVEN BY EQ.(9) C QM=G(M)/U SM=HVD(JJ)/(V*HV(JJ)) C C R IS GIVEN BY EQ.(15) C R=(1.0/U)**2-(1.0/V)**2 C C N IS THE REFRACTIVE INDEX OF THE CYLINDER GIVEN C BY EQ.(12) C N=WAN1/WANO C C IF MODE EQUALS ONE, THE TM SCATTERING COEFFICIENTS ARE C CALCULATED. IF IT IS TWO, THE TE COEFFICIENTS ARE C CALCULATED C IF(MODE.EQ.1)X=PERB1/PERB IF(MODE.EQ.2)X=PERB/PERB1*(WAN1/WANO)**2 C C D IS THE EXPRESSION GIVEN BY EQ.(10) C D=(SM-X*QM)*(SM-N**2*QM/X) IF(M.EQ.0)GO TO 5 D=D-(R*M*COS(THETA))**2 5 D=D*PI*(V*HV(JJ))**2/2.0 AZ=CCC**M*AMPL*SIN(THETA) C C AS AND BS ARE THE SCATTERING COEFFICIENTS GIVEN BY C EQS.(7 AND 8) FOR THE TM INCIDENT POLARIZATION C AS(M)=-AZ*(BJV(JJ)/HV(JJ)+CCC*(SM-X*QM)/D) IF(MODE.EQ.1.AND.THETA.NE.PI/2)THEN BS(M)=WANO/(PERB*AFY)*AZ*R*M*COS(THETA)/D ENDIF IF(MODE.EQ.2.AND.THETA.NE.PI/2)THEN BS(M)=-WANO/(EPSO*AFY)*AZ*R*M*COS(THETA)/D ENDIF IF(THETA.EQ.PI/2)BS(M)=0.0 10 CONTINUE C C USING THE SYMMETRY OF THE SCATTERING COEFFICIENTS (EQ.29)) DO 20 M=1,NBG AS(-M)=(-1)**M*AS(M) BS(-M)=-(-1)**M*BS(M) 20 CONTINUE RETURN END C SUBROUTINE AMPLITUDE(TT,NBG,XFACTOR,BB1,ANGLEI 6 ,ANGLEF,NUMBER,THETA,ISIZE,NSIZE) C C THIS SUBROUTINE CALCULATES THE AMPLITUDE SACTTERING MATRIX C ELEMENTS T1,T2,T3, AND T4 GIVEN BY EQS.(16,17,23,AND 24) C INTEGER I,NUMBER,NBG,M,ISIZE,NSIZE REAL PI,ANGLE,ANGLEI,ANGLEF,PHI,XFACTOR,THETA COMPLEX CCC,SUMT,EMPH,BB1(-ISIZE:ISIZE) COMPLEX TT(0:NSIZE) C CCC=CMPLX(0.0,1.0) PI=4.0*ATAN(1.0) C C THE PROGRAM CONSIDERS THE FORWARD SCATTERING ANGLE AS ZERO. C I.E , THE SCATTERING ANGLE IN EQS.(16,17,23,AND 24) IS C SHIFTED BY 180 DEGREES. C DO 10 I=0,NUMBER IF(NUMBER.NE.0)THEN ANGLE=ANGLEI+FLOAT(I)*(ANGLEF-ANGLEI)/NUMBER END IF IF(NUMBER.EQ.0)ANGLE=ANGLEI PHI=(ANGLE+180)*PI/180. SUMT=CMPLX(0.0,0.0) DO 20 M=-NBG,NBG EMPH=CMPLX(COS(M*PHI),-SIN(M*PHI)) SUMT=SUMT+BB1(M)*EMPH*CCC**M 20 CONTINUE TT(I)=XFACTOR*SUMT/SIN(THETA) 10 CONTINUE RETURN END C SUBROUTINE SIJ(T1,T2,T3,T4,I1,I2,I3,I4,ANGLEI 6 ,ANGLEF,NUMBER,ANORM,NSIZE,ANGLE,S11) C C THIS SUBROUTINE CALCULATES THE MUELLER SCATTERING C MATRIX ELEMENTS(SIJ) C INTEGER NSIZE,I,NUMBER REAL S11(0:NSIZE),S12 REAL S13,S14,S22,S23 REAL S24,S33,S34,S44 REAL ANGLE(0:NSIZE),ANGLEI,ANGLEF,ANORM REAL I1(0:NSIZE),I2(0:NSIZE) REAL I3(0:NSIZE),I4(0:NSIZE) COMPLEX T1(0:NSIZE),T2(0:NSIZE) COMPLEX T3(0:NSIZE),T4(0:NSIZE) C WRITE(6,10) 10 FORMAT(/,20X,'NORMALIZED SIJ') WRITE(6,20) 20 FORMAT(/,5X,'ANGLE',17X,'S11',12X,'S12') C C I1,I2,I3,AND I4 ARE THE MAGNITUDE SQUARED OF T1,T2,T3,AND T4, C RESPECTIVELY. C C THE MUELLER SCATTERING MATRIX ELEMENTS (SIJ) ARE CALCULATED C USING EQ.(26) . SIJ EXCEPT S11 ARE NORMALIZED TO S11 AT ANY C ANGLE. S11 IS NORMALIZED TO S11(0) C DO 30 I=0,NUMBER C IF(NUMBER.NE.0)THEN ANGLE(I)=ANGLEI+FLOAT(I)*(ANGLEF-ANGLEI)/NUMBER END IF IF(NUMBER.EQ.0)ANGLE(I)=ANGLEI I1(I)=T1(I)*CONJG(T1(I)) I2(I)=T2(I)*CONJG(T2(I)) I3(I)=T3(I)*CONJG(T3(I)) I4(I)=T4(I)*CONJG(T4(I)) S11(I)=(I1(I)+I2(I)+I3(I)+I4(I))/2.0 S12=(I1(I)-I2(I)-I3(I)+I4(I))/2.0 WRITE(6,120)ANGLE(I),S11(I)/ANORM,S12/S11(I) C 30 CONTINUE C WRITE(6,40) 40 FORMAT(/,5X,'ANGLE',17X,'S13',12X,'S14') C DO 50 I=0,NUMBER C S13=REAL(T1(I)*CONJG(T3(I))+T4(I)*CONJG(T2(I))) S14=AIMAG(T3(I)*CONJG(T1(I))+T2(I)*CONJG(T4(I))) WRITE(6,120)ANGLE(I),S13/S11(I),S14/S11(I) C 50 CONTINUE C WRITE(6,60) 60 FORMAT(/,5X,'ANGLE',17X,'S22',12X,'S23') C DO 70 I=0,NUMBER S22=(I1(I)+I2(I)-I4(I)-I3(I))/2.0 S23=REAL(T1(I)*CONJG(T3(I))-T2(I)*CONJG(T4(I))) WRITE(6,120)ANGLE(I),S22/S11(I),S23/S11(I) C 70 CONTINUE C WRITE(6,80) 80 FORMAT(/,5X,'ANGLE',17X,'S24',12X,'S33') C DO 90 I=0,NUMBER C S24=AIMAG(T3(I)*CONJG(T1(I))+T4(I)*CONJG(T2(I))) S33=REAL(T2(I)*CONJG(T1(I))+T3(I)*CONJG(T4(I))) WRITE(6,120)ANGLE(I),S24/S11(I),S33/S11(I) 90 CONTINUE C WRITE(6,100) 100 FORMAT(/,5X,'ANGLE',17X,'S34',12X,'S44') C DO 110 I=0,NUMBER S34=AIMAG(T2(I)*CONJG(T1(I))+T3(I)*CONJG(T4(I))) S44=REAL(T1(I)*CONJG(T2(I))-T4(I)*CONJG(T3(I))) WRITE(6,120)ANGLE(I),S34/S11(I),S44/S11(I) 110 CONTINUE C WRITE(6,130) 130 FORMAT(/,5X,'MAGNITUDE SQUARED OF SCATTERED AMPLITUDES') C WRITE(6,140) 140 FORMAT(/,5X,'ANGLE',17X,'I1',12X,'I2') C DO 150 I=0,NUMBER WRITE(6,120)ANGLE(I),I1(I),I2(I) 150 CONTINUE C WRITE(6,160) 160 FORMAT(/,5X,'ANGLE',17X,'I3',12X,'I4') C DO 170 I=0,NUMBER WRITE(6,120)ANGLE(I),I3(I),I4(I) 170 CONTINUE C 120 FORMAT(3X,E13.7,6X,2(E13.7,2X)) RETURN END C SUBROUTINE CROSS(CR,AS,BS,MODE,NBG,PERB,EPSO,AWANO 6 ,AFY,ISIZE) C C THIS SUBROUTINE IS DESIGNED TO CALCULATED THE C SCATTERING CROSS SECTION PER UNIT LENGTH FOR C THE CURRENT POLARIZATION C INTEGER ISIZE,I,NBG,MODE REAL PI,ETA,SUMA,SUMB,XFACTOR REAL AFY,AWANO,CR COMPLEX AS(-ISIZE:ISIZE),BS(-ISIZE:ISIZE) C PI=4.0*ATAN(1.0) C C ETA IS THE IMPEDANCE OF FREE SPACE C ETA=120.*PI SUMA=0.0 SUMB=0.0 DO 10 I=-NBG,NBG SUMA=SUMA+AS(I)*CONJG(AS(I)) SUMB=SUMB+BS(I)*CONJG(BS(I)) 10 CONTINUE SUMA=SUMA*EPSO SUMB=SUMB*PERB C C IF MODE=1,THE CURRENT POLARIZATION IS TM. HOWEVER, C MODE=2 INDICATES THE CURRENT ONE IS TE. C IF(MODE.EQ.1)XFACTOR=4.0*ETA*AFY/AWANO**2 IF(MODE.EQ.2)XFACTOR=4.0*AFY/(ETA*AWANO**2) C C CR IS THE SCATTERING CROSS SECTION PER UNIT LENGTH FOR C THE CURRENT POLARIZATION GIVEN BY EQ.(18) FOR THE INCIDENT C TM AND BY EQ.(25) FOR THE TE C CR=XFACTOR*(SUMA+SUMB) RETURN END C SUBROUTINE CHECK(ISIZE,NSIZE,NBG,NUMBER,MST,IFF,IGG,NMX) C C THIS SUBROUTINE IS TO ENSURE THAT NO ARRAY C IS OVERWRITTEN USING EQ.(30) C INTEGER ISIZE,NSIZE,NBG,NUMBER IERROR=0 JERROR=0 KERROR=0 LERROR=0 C IF((IGG+ISIZE).LT.NMX)THEN CALL MESSAGE IERROR=1 WRITE(6,17)NMX-ISIZE 17 FORMAT('ALLOCATED SPACE FOR PARAMETER IGG IS NOT 6 ENOUGH, IGG MUST BE AT LEAST',I7) WRITE(6,*) ENDIF C IF((IFF+ISIZE).LT.MST)THEN CALL MESSAGE JERROR=1 WRITE(6,11)MST-IFF 11 FORMAT(3X,'ALLOCATED SPACE FOR PARAMETER IFF IS NOT 6 ENOUGH, IFF MUST BE AT LEAST', I7) WRITE(6,*) END IF C IF(ISIZE.LE.NBG+2)THEN CALL MESSAGE KERROR=1 WRITE(6,10)NBG+2 10 FORMAT(3X,'ALLOCATED SPACE FOR ISIZE IS NOT 6 ENOUGH, ISIZE MUST BE AT LEAST' ,I7) WRITE(6,*) END IF C IF(NSIZE.LT.NUMBER+2)THEN CALL MESSAGE LERROR=1 WRITE(6,20)NUMBER+2 20 FORMAT(3X,'ALLOCATED SPACE FOR NSIZE IS NOT ENOUGH, 6 NSIZE MUST BE AT LEAST' ,I7) WRITE(6,*) END IF C IF(IERROR.EQ.1.OR.JERROR.EQ.1)STOP IF(KERROR.EQ.1.OR.LERROR.EQ.1)STOP RETURN END SUBROUTINE DRLOG(ARG,NN,G,ISIZE,IGG,NMX) C C THIS SUBROUTINE CALCULATES THE LOGARITHMIC DERAVITATIVES C OF BESSEL FUNCTIONS OF THE COMPLEX ARGUMENT. IT USES THE C THE DOWNWARD RECURSION FORMULAE. SETTING THE C FUNCTION TO ZERO AT NMX (G(NMX)=0) C INTEGER NMX,NN,I,ISIZE COMPLEX G(0:IGG+ISIZE) COMPLEX ARG C C NN IS THE DESIRED HIGHEST ORDER C N1=CABS(ARG) c c C NMX=1.1*MAX(NN,N1)+15 G(NMX)=CMPLX(0.0,0.0) DO 10 I=NMX,1,-1 G(I-1)=FLOAT((I-1))/ARG-1.0/((FLOAT(I)/ARG)+G(I)) 10 CONTINUE RETURN END c subroutine BESSEL(x,bj,by,isize,nc,IFF,MST) C THIS SUBROUTINE CALCULATES THE BESSEL FUNCTION OF THE C FIRST KIND (BJ) AND THE SECOND KIND (BY) OF THE REAL ARGUMENT (X) C THE DETAILS ARE GIVEN BY THE BOOK OF BOHREN AND HUFFMAN. C real x real bj(isize),by(isize),f(isize+IFF),xstop,alpha,rn real*8 gamma,twodivpi integer nstop,nmx,j,n,nn,k,mst,m1,l,ml,m2 integer ll,m3,m4,kk,ns gamma=0.57721566490153d0 twodivpi=0.63661977236758d0 c c Bessel functions j(n) are computed by downward recurrence c beginning at n = MST c Bessel functions y(n) are computed by upward recurrence, c bj(n+1) = j(n), by(n+1) = y(n). c nstop=nc mst=(mst/2)*2 f(mst+1)=0.0 f(mst)=1.0e-32 m1=mst-1 do 201 l=1,m1 ml=mst-l f(ml)=2.0*float(ml)*f(ml+1)/x-f(ml+2) 201 continue alpha=f(1) m2=mst-2 do 202 ll=2,m2,2 alpha=alpha+2.0d0*f(ll+1) 202 continue do 203 n=1,m1 bj(n)=f(n)/alpha 203 continue by(1)=bj(1)*(alog(x/2.0)+gamma) m4=mst/2-1 do 204 l=1,m4 by(1)=by(1)-2.0*((-1.0)**l)*bj(2*l+1)/float(l) 204 continue by(1)=twodivpi*by(1) by(2)=(bj(2)*by(1)-twodivpi/x)/bj(1) ns=nstop-1 do 205 kk=1,ns by(kk+2)=2.0*float(kk)*by(kk+1)/x-by(kk) 205 continue c return end C SUBROUTINE MESSAGE WRITE(6,10) 10 FORMAT(3X,'THE ALLOCATED SPACE OF THE PARAMETER') WRITE(6,20) 20 FORMAT('INDICATED BELOW IS NOT SUFFICIENT, THIS IS CORRECTED BY') WRITE(6,30) 30 FORMAT('LOCATING THE RIGHT SPACE FOR THE PARAMETER') WRITE(6,40) 40 FORMAT('IN THE STATEMENT OF THE MAIN PROGRAM') RETURN END