C +++ C C Source: src/trace/reflec.F C C ---------------------------------------------- C SHADOW C Center for X-ray Lithography C University of Wisconsin-Madison C 3731 Schneider Dr., Stoughton, WI, 53589 C ---------------------------------------------- C C Log: reflec.F C Revision 1.10 1991/07/06 20:07:38 khan C Grenoble and after. Minor changes C C Revision 1.9 91/04/05 15:05:59 cwelnak C changed quotes on #includes C C Revision 1.8 91/03/22 10:45:11 cwelnak C SUN version -- INCLUDE to #include C C Revision 1.7 91/02/20 13:02:16 cwelnak C dimensioning PIN might help... C C Revision 1.6 91/02/14 11:49:36 cwelnak C added variable PIN to make calling statement and subroutine consistent in C number of variables used. C C Revision 1.5 90/11/13 14:01:54 khan C Cleanup and SAVE statements C C Revision 1.4 90/10/24 12:03:12 khan C Added SAVE statement to make certain variables STATIC. C C Revision 1.3 90/07/19 21:37:57 khan C Put #ifdef's to make it work on BOTH VMS and Ultrix C C Revision 1.2 90/07/14 22:51:15 khan C All public include files (common.blk, etc) are now in ./../include/ dir. C C Revision 1.1 90/07/10 14:56:46 khan C Initial revision C C C --- #if defined(unix) || HAVE_F77_CPP # include #elif defined(vms) INCLUDE 'SHADOW$INC:HEADER.TXT/LIST' #endif C+++ C SUBROUTINE REFLEC C C PURPOSE To compute the local reflectivity of a mirror; C C FLAGS k_what: .lt. 0 --> initialization call. Reads C in data file. C .gt. 0 --> performs computations. C (see below) C C ARGUMENTS [ I ] PIN : (x,y,z) of the intercept C [ I ] wnum : wavenumber (cm-1) C [ I ] sin_ref : sine of angle from surface C [ I ] cos_pole : cosine of angle of normal from pole C [ O ] R_P : p-pol reflection coefficient C [ O ] R_S : s-pol " " C [ I ] ABSOR : film thickness C [ O ] ABSOR : absorption coefficient C C--- SUBROUTINE REFLEC (PIN,WNUM,SIN_REF,COS_POLE,R_P,R_S, $ PHASEP,PHASES,ABSOR,K_WHAT) #if defined(unix) || HAVE_F77_CPP c c This causes problems with F77 drivers, since can't use -I directive. c so I'll use the standard cpp directive instead. c c INCLUDE './../include/common.blk' c c # include #elif defined(vms) INCLUDE 'SHADOW$INC:COMMON.BLK/LIST' #endif DIMENSION ZF1(1000),ZF2(1000) DIMENSION PIN(3) dimension t_oe(200),gratio(200) C dimension do(200),de(200) dimension ener(300),delta_s(300),beta_s(300),delta_e(300) dimension beta_e(300),delta_o(300),beta_o(300) #if G77 logical ele(5),elo(5),els(5) #else logical*1 ele(5),elo(5),els(5) #endif character*80 file_grade dimension tspl (2,101,2,101),tx(101),ty(101),pds(6) dimension gspl (2,101,2,101),gx(101),gy(101) external dbcevl common /aaa/ t_oe,gratio common /bbb/ delo,beto,dele,bete,dels,bets C C SAVE the variables that need to be saved across subsequent invocations C of this subroutine. C SAVE QMIN, QMAX, QSTEP, DEPTH0, NREFL, TFILM, $ ZF1, ZF2, $ NIN, ENER, $ DELTA_S, BETA_S, $ NPAIR, $ DELTA_E,BETA_E, $ DELTA_O,BETA_O, $ TSPL,TX,TY,PDS, $ GSPL,GX,GY C C Initialization call. The ZF1,ZF2 values do NOT correspond to the F1,F2 C atomic scattering factors, as they contain a more complex form: C ZFi = Fi*RADIUS*4*PI*ATOMS C WNUM is the WAVENUMBER (cm-1) of the ray. C ALFA and GAMMA are the complex dielectric function C EPSILON = 1 - ALFA + i*GAMMA [ i=sqrt(-1) ] C and are dependent ONLY on the material, while the reflectivities C depend of course on the geometry too. C ALFA and GAMMA may be generated by PREREF.EXE in [CERRINA.ABS], C which is based on B.Henke data (see program header for complete reference). C C Two flags control the execution. C F_REFL = 0 ZF1,ZF2 are read in as arrays C = 1 ALFA and GAMMA are defined in the I/O session C and thus wavelength independent C = 2 d-spacings and optical constants for multilayer C are read in C K_WHAT = 0 Initialization C 1 Reflection case C 2 Absorption case C IF (K_WHAT.EQ.0) THEN IF (F_REFL.EQ.0) THEN #ifdef vms OPEN (23,FILE=FILE_REFL,STATUS='OLD', $ READONLY,FORM='UNFORMATTED') #else OPEN (23,FILE=FILE_REFL,STATUS='OLD', $ FORM='UNFORMATTED') #endif READ (23) QMIN,QMAX,QSTEP,DEPTH0 READ (23) NREFL READ (23) (ZF1(I),I=1,NREFL) READ (23) (ZF2(I),I=1,NREFL) CLOSE (23) TFILM = ABSOR D D DO ID=1,NREFL D DPHOT = (QMIN + QSTEP*(ID-1))/TWOPI*TOCM D WRITE (35,1010) DPHOT,ZF1(ID),ZF2(ID) D CONTINUE D1010 FORMAT (1X,3(E15.8,2X)) D RETURN ELSE IF (F_REFL.EQ.2) THEN c c this version allows specification of the individual c layer thicknesses. c c input parameters: c npair = no. of layer pairs (npair=0 means an ordinary mirror, elo is the mirror c xlam = wavelength (angstroms) c elo = odd layer material c ele = even layer material c els = substrate material c 1.0 - delo - i*beto = complex refractive index (odd) c 1.0 - dele - i*bete = complex refractive index (even) c t_oe = thickness t(odd)+t(even) in Angstroms of each layer pair c gratio = gamma ratio t(even)/(t(odd)+t(even)) of each layer pair c phr = grazing angle in radians c c iunit = 23 #ifdef vms open(unit=iunit,FILE=FILE_REFL,status='OLD',readonly) #else open(unit=iunit,FILE=FILE_REFL,status='OLD') #endif READ(iunit,*) NIN READ(iunit,*) (ENER(I), I = 1, NIN) DO 13 I=1,NIN READ(iunit,*) DELTA_S(I),BETA_S(I) 13 CONTINUE DO 23 I=1,NIN READ(iunit,*) DELTA_E(I),BETA_E(I) 23 CONTINUE DO 33 I=1,NIN READ(iunit,*) DELTA_O(I),BETA_O(I) 33 CONTINUE READ(iunit,*) NPAIR do 11 i = 1, npair read(iunit,*) t_oe(i),gratio(i) 11 CONTINUE c c Is the multilayer thickness graded ? c read (iunit,*) i_grade if (i_grade.eq.1) then read (iunit,101) file_grade 101 format (a80) OPEN (45, FILE=FILE_GRADE, STATUS='OLD', FORM='UNFORMATTED') READ (45) NTX, NTY READ (45) TX,TY DO 205 I = 1, NTX DO 205 J = 1, NTY READ (45) TSPL(1,I,1,J),TSPL(1,I,2,J), ! spline for t $ TSPL(2,I,1,J),TSPL(2,I,2,J) 205 CONTINUE C READ (45) (((TSPL(1,I,1,J),TSPL(1,I,2,J), ! spline for t C $ TSPL(2,I,1,J),TSPL(2,I,2,J)), J = 1,NTY), I = 1,NTX) READ (45) NGX, NGY READ (45) GX,GY DO 305 I = 1, NGX DO 305 J = 1, NGY READ (45) GSPL(1,I,1,J),GSPL(1,I,2,J), ! spline for gamma $ GSPL(2,I,1,J),GSPL(2,I,2,J) 305 CONTINUE C READ (45) (((GSPL(1,I,1,J),GSPL(1,I,2,J), ! spline for gamma C $ GSPL(2,I,1,J),GSPL(2,I,2,J)), J = 1,NGY), I = 1,NGX) CLOSE (45) end if close(unit=iunit) tfilm = absor RETURN END IF END IF C C This is the notmal calculation part; C C If F_REFL is 1, ALFA and GAMMA are defined during the input session C and are not modified anymore (single line or closely spaced lines case) C IF (F_REFL.EQ.0) THEN !Both absorp and normal !reflectivity INDEX = (WNUM - QMIN)/QSTEP + 1 IF (INDEX.LT.1) INDEX=1 C ('REFLEC','Photon energy below lower limit.',0) IF (INDEX.GT.NREFL) INDEX=NREFL-1 C ('REFLEC','Photon energy above upper limit.',0) IF (INDEX.EQ.NREFL) INDEX = INDEX - 1 WNUM0 = QSTEP*(INDEX-1) + QMIN DEL_X = WNUM - WNUM0 DEL_X = DEL_X/QSTEP ALFA = ZF1(INDEX) + (ZF1(INDEX+1)-ZF1(INDEX))*DEL_X GAMMA = ZF2(INDEX) + (ZF2(INDEX+1)-ZF2(INDEX))*DEL_X D WRITE (37,1020) WNUM,WNUM0,INDEX,DEL_X,ALFA,GAMMA D1020 FORMAT (1X,2(E15.8,1X),I4,3(E15.8,1X)) END IF IF (K_WHAT.EQ.1) THEN IF (F_REFL.NE.2) THEN C C Computes the optical coefficients. C COS_REF = SQRT(1.0D0 - SIN_REF**2) RHO = SIN_REF**2 - ALFA RHO = RHO + SQRT ((SIN_REF**2 - ALFA)**2 + GAMMA**2) RHO = SQRT(RHO/2) C C Computes now the reflectivities C RS1 = 4*(RHO**2)*(ABS(SIN_REF)-RHO)**2 + GAMMA**2 RS2 = 4*(RHO**2)*(ABS(SIN_REF)+RHO)**2 + GAMMA**2 R_S = RS1/RS2 C C Computes now the polarization ratio C RATIO1 = 4*RHO**2*(RHO*ABS(SIN_REF)-COS_REF**2)**2 + $ GAMMA**2*SIN_REF**2 RATIO2 = 4*RHO**2*(RHO*ABS(SIN_REF)+COS_REF**2)**2 + $ GAMMA**2*SIN_REF**2 RATIO = RATIO1/RATIO2 C C The reflectivity for p light will be C R_P = R_S*RATIO R_S = SQRT(R_S) R_P = SQRT(R_P) ELSE C C Multilayers reflectivity. C First interpolate for all the refractive indices. C XLAM = TWOPI/WNUM*1.0D8 ! Angstrom PHOT_ENER = WNUM/TWOPI*TOCM ! eV ELFACTOR = LOG10(1.0D04/30.0D0)/300.0D0 INDEX = LOG10(PHOT_ENER/ENER(1))/ELFACTOR + 1 C INDEX = 96.0*LOG10(PHOT_ENER/ENER(1)) + 1 IF (INDEX.LT.1) INDEX=1 C ('REFLEC','Photon energy too small.',2) IF (INDEX.GT.NIN) INDEX=NIN-1 C ('REFLEC','Photon energy too large.',2) DELS = DELTA_S(INDEX) + (DELTA_S(INDEX+1) - DELTA_S(INDEX)) $ *(PHOT_ENER - ENER(INDEX))/(ENER(INDEX+1) - ENER(INDEX)) BETS = BETA_S(INDEX) + (BETA_S(INDEX+1) - BETA_S(INDEX)) $ *(PHOT_ENER - ENER(INDEX))/(ENER(INDEX+1) - ENER(INDEX)) DELE = DELTA_E(INDEX) + (DELTA_E(INDEX+1) - DELTA_E(INDEX)) $ *(PHOT_ENER - ENER(INDEX))/(ENER(INDEX+1) - ENER(INDEX)) BETE = BETA_E(INDEX) + (BETA_E(INDEX+1) - BETA_E(INDEX)) $ *(PHOT_ENER - ENER(INDEX))/(ENER(INDEX+1) - ENER(INDEX)) DELO = DELTA_O(INDEX) + (DELTA_O(INDEX+1) - DELTA_O(INDEX)) $ *(PHOT_ENER - ENER(INDEX))/(ENER(INDEX+1) - ENER(INDEX)) BETO = BETA_O(INDEX) + (BETA_O(INDEX+1) - BETA_O(INDEX)) $ *(PHOT_ENER - ENER(INDEX))/(ENER(INDEX+1) - ENER(INDEX)) C C CALL FRESNEL (NPAIR,SIN_REF,COS_POLE,XLAM,R_S,R_P,PHASES,PHASEP) C c c If graded, compute the factor for t and gamma at the intercept PIN. c IF (I_GRADE.EQ.1) THEN XIN = PIN(1) YIN = PIN(2) CALL DBCEVL (TX,NTX,TY,NTY,TSPL,101,XIN,YIN,PDS,IER) IF (IER.NE.0) THEN CALL MSSG ('REFLEC','Spline error # ',IER) RETURN END IF TFACT = PDS(1) C CALL DBCEVL (GX,NGX,GY,NGY,GSPL,101,XIN,YIN,PDS,IER) IF (IER.NE.0) THEN CALL MSSG ('REFLEC','Spline error # ',IER) RETURN END IF GFACT = PDS(1) ELSE TFACT = 1.0D0 GFACT = 1.0D0 END IF CALL FRESNEL $ (TFACT,GFACT,NPAIR,SIN_REF,COS_POLE,XLAM,R_S,R_P,PHASES,PHASEP) END IF ELSE IF(K_WHAT.EQ.2) THEN C C This is the transmission case. SIN_REF is now the incidence angle C onto the filter. C C Computes now the penetration depth. SIN_REF is now the cosine of C the incidence angle of the ray on the screen. C c DEPTH = DEPTH0*GAMMA*SIN_REF/WNUM AB_COEFF = WNUM*GAMMA/ABS(SIN_REF) C C Computes the film absorption. The thickness is passed at the call with C K_WHAT = 0 C C ABSOR is the attenuation of the A vector. C ABSOR = EXP(-TFILM*AB_COEFF/2.0D0) END IF D D DPHOT = WNUM/TWOPI*TOCM D WRITE (34,1000) DPHOT,R_S,R_P,SIN_REF D1000 FORMAT (1X,4(E15.8,2X)) D RETURN END ****************************************************************************** subroutine FRESNEL (TFACT,GFACT,N,SIN_REF,COS_POLE,xlam,ans, $ anp,phaseS,PHASEP) c IMPLICIT REAL*8 (A-H,O-Z) DATA PI / 3.1415 92653 58979 32384 62643 D0 / DATA TODEG / 57.2957 79513 08232 08767 98155 D0 / complex*16 ci,fo,fe,fv,ffe,ffv,ffvp,ffo,ffep,ffop,re2, $ ro2,ao,ae,r,rp,fs,ffs,ffsp,rs2 dimension t_oe(200),gratio(200) common /aaa/ t_oe,gratio common /bbb/ delo,beto,dele,bete,dels,bets c ci=(0.0D0,1.0D0) ro2=(1.0D0-delo-ci*beto)**2 re2=(1.0D0-dele-ci*bete)**2 rs2=(1.0D0-dels-ci*bets)**2 SIN_REF2 = SIN_REF**2 COS_REF2 = 1.0D0 - SIN_REF2 C refo=(sin(phr))**2-2.0*delo C xmfo=-2.0*beto C fo=cmplx(refo,xmfo) fo = ro2 - COS_REF2 C refe=(sin(phr))**2-2.0*dele C xmfe=-2.0*bete C fe=cmplx(refe,xmfe) fe = re2 - COS_REF2 refv=SIN_REF2 xmfv=0.0D0 fv = Dcmplx(refv,xmfv) C refs=(sin(phr))**2-2.0*dels C xmfs=-2.0*bets C fs=cmplx(refs,xmfs) fs = rs2 - COS_REF2 C fo=cDsqrt(fo) fe=cDsqrt(fe) fv=cDsqrt(fv) fs=cDsqrt(fs) ffe=(fe-fo)/(fe+fo) ffo=-ffe ffv=(fv-fo)/(fv+fo) ffs=(fe-fs)/(fe+fs) ffep=(fe/re2-fo/ro2)/(fe/re2+fo/ro2) ffop=-ffep ffvp=(fv-fo/ro2)/(fv+fo/ro2) ffsp=(fe/re2-fs/rs2)/(fe/re2+fs/rs2) r=(0.0D0,0.0D0) rp=(0.0D0,0.0D0) do 1 j=1,n c c compute the thickness for the odd and even material : c THICK = T_OE(J) * TFACT GAMMA = GRATIO(J) * GFACT t_e = GAMMA * THICK t_o = (1.0D0-GAMMA) * THICK c ao=-ci*(pi*fo*t_o*cos_pole/xlam) ae=-ci*(pi*fe*t_e*cos_pole/xlam) ao=cDexp(ao) ae=cDexp(ae) if(j.eq.1)go to 6 r=(ae**4)*(r+ffe)/(r*ffe+1.0D0) rp=(ae**4)*(rp+ffep)/(rp*ffep+1.0D0) go to 7 6 r=(ae**4)*(r+ffs)/(r*ffs+1.0D0) rp=(ae**4)*(rp+ffsp)/(rp*ffsp+1.0D0) 7 r=(ao**4)*(r+ffo)/(r*ffo+1.0D0) rp=(ao**4)*(rp+ffop)/(rp*ffop+1.0D0) 1 continue r=(r+ffv)/(r*ffv+1.0D0) pp = Dimag(r) qq = Dreal(r) CALL ATAN_2 (PP,QQ,PHASES) ! S phase change in units of radians rp=(rp+ffvp)/(rp*ffvp+1.0D0) PP = DIMAG(RP) QQ = DREAL (RP) CALL ATAN_2 (PP,QQ,PHASEP) ! P phase change in units of radians anp=cDabs(rp) C anp=anp**2 ans=cDabs(r) C ans=ans**2 return end