SUBROUTINE SCAT(AK,AKS,EM1,EM2,E02, & CBKSCA,CSCA,CSCAG,CXE00,CXF1L,CXF2L,CXP,CXSCR1,CXSCR2,CXSCR3, & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR, & PHIN,SCRRS1,THETAN,IXYZ) C C Arguments: REAL CBKSCA,CSCA,CSCAG,E02 INTEGER ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR INTEGER*2 IXYZ(MXNAT,3) COMPLEX CXE00(3),CXF1L(MXSCA),CXF2L(MXSCA),CXP(NAT,3), & CXSCR1(MXN3),CXSCR2(MXN3),CXSCR3(MXN3) REAL AK(3),AKS(3,MXSCA),EM1(3,MXSCA),EM2(3,MXSCA), $ PHIN(MXSCA),SCRRS1(MXN3),THETAN(MXSCA) C C Local variables: COMPLEX CXI REAL AK2,AK3,AKK,CONST,COSPHI,COSTH,DCOSTH,DPHI,EINC,ES2, & PI,PHI,RRR,SINPHI,SINTH,TERM,W,WTH INTEGER ICOSTH,IPHI,IPHM,J,K,ND COMPLEX CXES(3) REAL AKS0(3),AKS1(3),AKS2(3) C C Intrinsic functions: INTRINSIC COS,EXP,MIN,SIN,SQRT,CONJG,REAL C C Data statements: DATA CXI/(0.,1.)/ C*********************************************************************** C C SCAT computes energy radiated by array of oscillating dipoles C and corresponding scattering properties if oscillations are in C response to incident E field. C C Given: C AK(1-3) = KX,KY,KZ = components of incident k-vector C AKS(3,MXSCA)=scattered k vectors C EM1(3,MXSCA)=pol.vector 1 for each scattering direction C EM2(3,MXSCA)=pol.vector 2 for each scattering direction C E02 = |E_0|^2 , where E_0 = incident complex E field C NAT = number of dipoles C NAT3 = 3*NAT C XYZ(NAT,1-3) = X,Y,Z for dipoles 1-NAT C CXE00(1-3) = (complex) reference polarization vector C CXP(1-NAT,1-3) = components of polarization of each dipole C ICTHM = number of THETA values for calculation of CSCA and G C (ICTHM must be odd; recommend ICTHM > 5*(2*pi*a/lambda) C for accuracy) C IPHIM = number of PHI values for calculation of CSCA and G C NDIR = number of directions at which scattering matrix elements C are to be computed C THETAN(1-NDIR) : scattering angle (radians); C cos(THETAN) = n0 dot n C PHIN(1-NDIR) : scattering angle (radians); C PHIN is azimuthal angle of n from plane defined C by n0 and Re(CXE00) C C Returns: C CBKSCA=differential scattering cross section for theta=pi C (lattice units) C CSCA = scattering cross section (lattice units) C G = , where theta=scattering angle C F1L(NDIR) = f_{1L} for direction NDIR C F2L(NDIR) = f_{2L} for direction NDIR C where f_{1L} connects incident polarization L to C outgoing theta polarization (E in scattering plane) C f_{2L} connects incident polarization L to outgoing phi C polarization (E perpendicular to scattering plane) C This is same f_{ml} notation as used by Draine (1988; C Ap.J., 333, 848). C C Notes: C Scratch vectors introduced to permit vectorization: C SCRRS1 = real vector of length.GE.NAT C SCRCS1,SCRCS2,SCRCS3 = complex vectors of length.GE.NAT C C B.T.Draine, Princeton Univ. Obs., 87/1/4 C History: C 88.06.29 (BTD): modified C 89.11.28 (BTD): modified C 90.11.02 (BTD): modified to enable use with vacuum sites. C 90.11.20 (BTD): corrected error in calculation of f_ml (missing C factor of k). C 90.12.03 (BTD): changed ordering of XYZ0 and CXP C eliminated cxe0 from argument list C 90.12.10 (BTD): replaced XYZ0 with IXYZ0 to reduce memory use C 91.08.14 (BTD): add CBKSCA to argument list for SCAT C add code to eliminate sum over phi when theta=0 or pi C add code to calculate CBKSCA C 94.04.04 (BTD): removed superfluous CMPLX( ) operation from 2 lines C in DO 200 loop. Also added new variable TERM to move C some computation out of last DO loop. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** PI=4.*ATAN(1.) C C ICTHM,IPHIM determine angular grid used for calculation of total C scattering cross section and G= C ICTHM=number of theta values for which scattering is evaluated C (0 to pi) C IPHIM=number of phi values for which scattering is evaluated C (0 to 2*pi) C Note: ICTHM must be odd as Simpson's rule integration is employed. C Check whether ICTHM is odd, and increase by 1 if not. C IF(ICTHM.EQ.2*(ICTHM/2))ICTHM=ICTHM+1 C C Incident k vector AK defines one axis (for spherical integration) C Use (real part of) reference polarization vector CXE00 to define C second coordinate axis. C Construct third coordinate axis by taking cross product of first C two axes. C AKS1, AKS2 = vectors of length equal to AK lying along second and C third coordinate axes. C AK2=AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3) AKK=SQRT(AK2) AK3=AKK*AK2 AKS1(1)=REAL(CXE00(1)) AKS1(2)=REAL(CXE00(2)) AKS1(3)=REAL(CXE00(3)) C EINC = magnitude of (real part of) reference electric field EINC=SQRT(AKS1(1)**2+AKS1(2)**2+AKS1(3)**2) IF(EINC.LE.0.)CALL ERRMSG('FATAL','SCAT', & 'CXE00 IS PURE IMAGINARY!') RRR=AKK/EINC AKS1(1)=RRR*AKS1(1) AKS1(2)=RRR*AKS1(2) AKS1(3)=RRR*AKS1(3) AKS2(1)=(AK(2)*AKS1(3)-AK(3)*AKS1(2))/AKK AKS2(2)=(AK(3)*AKS1(1)-AK(1)*AKS1(3))/AKK AKS2(3)=(AK(1)*AKS1(2)-AK(2)*AKS1(1))/AKK C DCOSTH=2./REAL(ICTHM-1) DPHI=2.*PI/REAL(IPHIM) CSCA=0. CSCAG=0. C C CONST = multiplicative factor entering weights C factor 2/(ICTHM-1) = d(cos(theta)) C factor 2*pi/IPHIM = d(phi) C factor 3 in denom from Simpson's rule (later mult. by 1,4,2..) C Use Simpson's rule integration over d(cos(theta)) CONST=4.*PI/ (3.*REAL((ICTHM-1)*IPHIM)) C Note that for cos(theta)=-1 or 1 we do not need to sum over phi, C so treat these two directions separately DO 80 ICOSTH=2,ICTHM-1 WTH=4. IF(ICOSTH.GT.2*(ICOSTH/2)) WTH=2. C IF(ICOSTH.EQ.1) WTH=1. C IF(ICOSTH.EQ.ICTHM) WTH=1. COSTH=1.-REAL(ICOSTH-1)*DCOSTH SINTH=SQRT(1.-MIN(COSTH*COSTH,1.)) IPHM=IPHIM C Now evaluate total weight factor W W=CONST*WTH DO 70 IPHI=1,IPHM PHI=DPHI*(REAL(IPHI)-0.5) SINPHI=SIN(PHI) COSPHI=COS(PHI) C Now compute scattered k vector C Initialize CXES(1-3) = scattered electric field DO 10 J=1,3 AKS0(J)=COSTH*AK(J)+SINTH* & (COSPHI*AKS1(J)+SINPHI*AKS2(J)) CXES(J)=(0.,0.) 10 CONTINUE C C Now compute CXES(1-3) = (r/k*k)*exp(-ikr+iwt)*E(ks,r,t) C where E(ks,r,t) is complex electric field vector propagating C with wave vector ks, at distance r from origin C DO 50 J=1,NAT CXSCR1(J)=(0.,0.) C Compute SCRCS1(J)=ks dot p(j) / ks**2 DO 20 K=1,3 CXSCR1(J)=CXSCR1(J)+AKS0(K)*CXP(J,K) 20 CONTINUE CXSCR1(J)=CXSCR1(J)/AK2 SCRRS1(J)=0. C SCRRS1(J) is ks dot r(j) DO 30 K=1,3 SCRRS1(J)=SCRRS1(J)+IXYZ(J,K)*AKS0(K) 30 CONTINUE C CXSCR2(J) is exp(-i*ks*r(j)) factor for atom j CXSCR2(J)=EXP(-CXI*SCRRS1(J)) DO 40 K=1,3 CXES(K)=CXES(K)+(CXP(J,K)-AKS0(K)*CXSCR1(J))* & CXSCR2(J) 40 CONTINUE 50 CONTINUE C C Have completed computation of vector ES(1-3) for scattering direction C (THETA,PHI). C ES2=0.0 DO 60 K=1, 3 ES2=ES2+CXES(K)*CONJG(CXES(K)) 60 CONTINUE CSCA=CSCA+W*ES2 CSCAG=CSCAG+W*COSTH*ES2 70 CONTINUE 80 CONTINUE C*** C Now do special case theta=0 and theta=pi W=CONST*IPHM DO 180 ICOSTH=-1,1,2 C Compute scattered k vector C Initialize CXES(1-3) = scattered electric field DO 110 J=1,3 AKS0(J)=ICOSTH*AK(J) CXES(J)=(0.,0.) 110 CONTINUE C C Now compute CXES(1-3) = (r/k*k)*exp(-ikr+iwt)*E(ks,r,t) C where E(ks,r,t) is complex electric field vector propagating C with wave vector ks, at distance r from origin C DO 150 J=1,NAT CXSCR1(J)=(0.,0.) C Compute SCRCS1(J)=ks dot p(j) / ks**2 DO 120 K=1,3 CXSCR1(J)=CXSCR1(J)+AKS0(K)*CXP(J,K) 120 CONTINUE CXSCR1(J)=CXSCR1(J)/AK2 SCRRS1(J)=0. C SCRRS1(J) is ks dot r(j) DO 130 K=1,3 SCRRS1(J)=SCRRS1(J)+IXYZ(J,K)*AKS0(K) 130 CONTINUE C CXSCR2(J) is exp(-i*ks*r(j)) factor for atom j CXSCR2(J)=EXP(-CXI*SCRRS1(J)) DO 140 K=1,3 CXES(K)=CXES(K)+(CXP(J,K)-AKS0(K)*CXSCR1(J))*CXSCR2(J) 140 CONTINUE 150 CONTINUE C C Have completed computation of vector ES(1-3) for scattering direction C theta=-1 or 1 C ES2=0.0 DO 160 K=1, 3 ES2=ES2+CXES(K)*CONJG(CXES(K)) 160 CONTINUE CSCA=CSCA+W*ES2 CSCAG=CSCAG+W*ICOSTH*ES2 C Save backscattering cross section IF(ICOSTH.EQ.-1)CBKSCA=AK2*AK2*ES2/E02 180 CONTINUE C*** CSCA=AK2*AK2*CSCA/E02 CSCAG=AK2*AK2*CSCAG/E02 IF(NDIR.LE.0)RETURN C C*** Calculate the 2 matrix elements FM1 for selected directions C F1L(NDIR) is to the polarization state e1 (E in scattering plane) C F2L(NDIR) is to the polarization state e2 (E perp. to scatt. plane) C IMPORTANT NOTE: The following calculation assumes that the incident C wave is linearly polarized, with E00 having no imaginary component. C TERM=AK3/SQRT(E02) DO 230 ND=1,NDIR CXF1L(ND)=(0.0,0.0) CXF2L(ND)=(0.0,0.0) DO 210 J=1,NAT CXSCR2(J)=(0.0,0.0) CXSCR3(J)=(0.0,0.0) DO 200 K=1, 3 CXSCR2(J)=CXSCR2(J)+CXP(J,K)*EM1(K,ND) CXSCR3(J)=CXSCR3(J)+CXP(J,K)*EM2(K,ND) 200 CONTINUE 210 CONTINUE DO 220 J=1, NAT CXSCR1(J)=EXP(-CXI*(AKS(1,ND)*IXYZ(J,1)+ & AKS(2,ND)*IXYZ(J,2)+ & AKS(3,ND)*IXYZ(J,3))) CXF1L(ND)=CXF1L(ND)+CXSCR2(J)*CXSCR1(J) CXF2L(ND)=CXF2L(ND)+CXSCR3(J)*CXSCR1(J) 220 CONTINUE C C Units: Dipole polarizations are in units of [E]*d**3, where d=lattice C spacing, and [E]=dimensions of E field C Therefore, up to this point CXF1L, CXF2L are in units C of [E]*d**3 C Now multiply by AK3/SQRT(E02) (units of 1/([E]d**3)) to get C dimensionless result. C CXF1L(ND)=CXF1L(ND)*TERM CXF2L(ND)=CXF2L(ND)*TERM 230 CONTINUE RETURN C END SUBROUTINE SCAVEC(MXSCA,NSCAT,THETAN,PHIN,CXE01,ENSC,EM1,EM2) C*********************************************************************** C Given: C MXSCA=dimension information for ENSC,PHIN,THETAN C NSCAT=number of scattering directions C THETAN(1-NSCAT)=scattering angles theta C PHIN(1-NSCAT)=scattering angles phi C CXE01(1-3)=complex polarization vector 1 (phi=0 direction C is defined by Real component of CXE01) C Returns: C ENSC(1-3,1-NSCAT)=scattering vectors in Lab Frame C EM1(1-3,1-NSCAT)=scattered pol vectors parallel to scat. plane C in Lab Frame C EM2(1-3,1-NSCAT)=scattered pol vectors perp. to scat. plane C in Lab Frame C It is assumed that incident propagation vector is in x-direction C in Lab Frame C C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C Arguments: INTEGER MXSCA,NSCAT COMPLEX CXE01(3) REAL EM1(3,MXSCA), EM2(3,MXSCA), ENSC(3,MXSCA), PHIN(MXSCA), & THETAN(MXSCA) C C Local variables: INTEGER I REAL COSPHI,COSTHE,DENOM,PHI0,SINPHI,SINTHE C C*********************************************************************** C*** First determine orientation of phi=0 direction DENOM=SQRT((REAL(CXE01(1)))**2 + (REAL(CXE01(2)))**2 + & (REAL(CXE01(3)))**2 ) PHI0=ACOS(REAL(CXE01(2))/DENOM) DO 1000 I=1,NSCAT COSPHI=COS(PHIN(I)+PHI0) SINPHI=SIN(PHIN(I)+PHI0) COSTHE=COS(THETAN(I)) SINTHE=SIN(THETAN(I)) ENSC(1,I)=COSTHE ENSC(2,I)=SINTHE*COSPHI ENSC(3,I)=SINTHE*SINPHI EM1(1,I)=-SINTHE EM1(2,I)=COSTHE*COSPHI EM1(3,I)=COSTHE*SINPHI EM2(1,I)=0. EM2(2,I)=-SINPHI EM2(3,I)=COSPHI 1000 CONTINUE RETURN END SUBROUTINE TARCYL(A1,A2,A,B,AXIS,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) C*********************************************************************** C Purpose: to construct cylindrical target from "atoms" C Input: C AXIS = 1.,2.,3. for symmetry axis along x,y,z axis C A = cylinder length (in units of lattice spacing d) C B = cylinder diameter (in units of lattice spacing d) C IOSHP=device number for "target.out" file C =-1 to suppress printing of "target.out" C MXNAT = dimensioning information C Returns: C A1(1-3) = unit vector along cylinder axis C A2(1-3) = unit vector perpendicular to cylinder axis C CDESCR = string describing target (up to 67 chars.) C NAT = number of atoms in target C IX(1-NAT)=x/d for atoms in target C IY(1-NAT)=y/d C IZ(1-NAT)=z/d C C B.T.Draine, Princeton Univ. Obs., 87.04.03 C History: C 89.12.21 (BTD): Modified for use by ddscat 89.12.21 C 90.12.03 (BTD): Modified for conformity with other target routines C 91.01.05 (BTD): Changed I4 -> I7 when printing NAT C 91.04.22 (BTD): Added XV,YV,ZV and A1,A2 to WRITE(IOSHP,FMT=9020) C 93.03.12 (BTD): Changed CDESCR*(*) -> CDESCR*67 C 95.04.07 (BTD): Corrected error in code which affected target C construction for 4.5TARELL cylindrical target: Length=',F8.4, &' Diameter=',F8.4,' IAXIS=',I1,/, &I7,' = NAT',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) END SUBROUTINE TESTCYL(X,Y,Z,A,B,IAXIS,IN) C Scalar arguments: INTEGER IAXIS LOGICAL IN REAL A,B,X,Y,Z C Local scalars: REAL R2,S IN=.FALSE. IF(IAXIS.EQ.1)THEN S=X R2=Y*Y+Z*Z ELSEIF(IAXIS.EQ.2)THEN S=Y R2=X*X+Z*Z ELSEIF(IAXIS.EQ.3)THEN S=Z R2=X*X+Y*Y ENDIF IF(R2.GT.0.25*B*B)RETURN IF(S.LT.1.)RETURN IF(S.GT.INT(A+.5))RETURN IN=.TRUE. RETURN END SUBROUTINE TARELL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) C*********************************************************************** C Routine to construct ellipsoid from "atoms" C Input: C AX=(x-length)/d (d=lattice spacing) C AY=(y-length)/d C AZ=(z-length)/d C IOSHP=device number for "target.out" file C =-1 to suppress printing of "target.out" C MXNAT=dimensioning information (max number of atoms) C Output: C A1(1-3)=(1,0,0)=unit vector defining target axis 1 in Target Frame C A2(1-3)=(0,1,0)=unit vector defining target axis 2 in Target Frame C NAT=number of atoms in target C IX(1-NAT)=x/d for atoms of target C IY(1-NAT)=y/d C IZ(1-NAT)=z/d C CDESCR=description of target (up to 67 characters) C B.T.Draine, Princeton Univ. Obs. C History: C 91/01/05 (BTD): Changed I4 -> I7 when printing NAT C 91/04/22 (BTD): Added AX,AY,AZ and A1,A2 to WRITE(IOSHP,FMT=9020) C 91/05/01 (BTD): TEMPORARY hack to be able to reproduce all arrays C used by Draine 88 (in which XOFF=YOFF=ZOFF=0.5 always). C 92/10/26 (BTD): Added code to set string CDESCR C Added one call to WRIMSG C 93/01/01 (BTD): Changed method for determining target axes A1, A2 C [previous criteria -- A1 = longest dimension, C A2 = intermediate dimension -- created confusion] C new rule: A1=(1,0,0), A2=(0,1,0) C 93/03/12 (BTD): Changed CDESCR*(*) -> CDESCR*67 C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT) REAL AX,AY,AZ REAL A1(3),A2(3) C** Local variables: INTEGER JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2 REAL AX2,AY2,AZ2,R,RYZ2,RZ2,X,XOFF,Y,YOFF,Z,ZOFF C*********************************************************************** C C Routine to construct pseudo-ellipsoidal target aray. C With occupied array sites contained within ellipsoidal surface C defined by (X/AX)**2+(Y/AY)**2+(Z/AZ)**2=0.25 C Ideal volume V=(pi/6)*AX*AY*AZ C C Dipoles are located at sites C (x,y,z)=(I+XOFF,J+YOFF,Z+KOFF), I,J,K=integers C XOFF,YOFF,ZOFF=constants C C For sphere: call with AX=AY=AZ C For spheroid: call with AX=AY (or AY=AZ or AX=AZ) C B.T.Draine, Princeton Univ. Obs., 88/8/12 C C C Criterion for choosing XOFF: C If AX is close to an even integer, take XOFF=1/2 C If AX is close to an odd integer, take XOFF=0 JX=INT(AX+.5) IF(2*(JX/2).EQ.JX)THEN XOFF=0.5 ELSE XOFF=0. ENDIF C Similar criterion for YOFF: JY=INT(AY+.5) IF(2*(JY/2).EQ.JY)THEN YOFF=0.5 ELSE YOFF=0. ENDIF C Similar criterion for ZOFF: JZ=INT(AZ+.5) IF(2*(JZ/2).EQ.JZ)THEN ZOFF=0.5 ELSE ZOFF=0. ENDIF C C*** temporary hack for experimental purposes: C Note: this was commented out 93/03/11 by BTD -- cannot remember why C or when this experimental change was made C XOFF=0.5 C YOFF=0.5 C ZOFF=0.5 C LMX1=-INT(.5*AX+.5) LMX2=INT(.5*AX-.25) LMY1=-INT(.5*AY+.5) LMY2=INT(.5*AY-.25) LMZ1=-INT(.5*AZ+.5) LMZ2=INT(.5*AZ-.25) AX2=AX*AX AY2=AY*AY AZ2=AZ*AZ NAT=0 C C Specify target axes A1 and A2 C Previous convention: A1 along longest dimension C A2 along intermediate dimension C New convention: A1=(1,0,0) in target frame C A2=(0,1,0) in target frame DO 0010 JX=1,3 A1(JX)=0. A2(JX)=0. 0010 CONTINUE C Old code: C IF(AX.GE.AY.AND.AX.GE.AZ)THEN C A1(1)=1. C IF(AY.GE.AZ)THEN C A2(2)=1. C ELSE C A2(3)=1. C ENDIF C ELSEIF(AY.GT.AX.AND.AY.GE.AZ)THEN C A1(2)=1. C IF(AX.GE.AZ)THEN C A2(1)=1. C ELSE C A2(3)=1. C ENDIF C ELSEIF(AZ.GT.AX.AND.AZ.GT.AY)THEN C A1(3)=1. C IF(AX.GE.AY)THEN C A2(1)=1. C ELSE C A2(2)=1. C ENDIF C ENDIF C New code: A1(1)=1. A2(2)=1. C C Determine list of occupied sites. DO 0300 JZ=LMZ1,LMZ2 Z=REAL(JZ)+ZOFF RZ2=Z*Z/AZ2 IF(RZ2.LT.0.25)THEN DO 0200 JY=LMY1,LMY2 Y=REAL(JY)+YOFF RYZ2=RZ2+Y*Y/AY2 IF(RYZ2.LT.0.25)THEN DO 0100 JX=LMX1,LMX2 X=REAL(JX)+XOFF R=RYZ2+X*X/AX2 IF(R.LT.0.25)THEN C Site is occupied: NAT=NAT+1 IX(NAT)=JX IY(NAT)=JY IZ(NAT)=JZ ENDIF 0100 CONTINUE ENDIF 0200 CONTINUE ENDIF 0300 CONTINUE C IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ') ENDIF C*********************************************************************** C Write target description into string CDESCR IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN WRITE(CDESCR,FMT='(A,I7,A)')' Spherical target containing', & NAT,' dipoles' ELSE WRITE(CDESCR,FMT='(A,I7,A)')' Ellipsoidal target containing', & NAT,' dipoles' ENDIF C*********************************************************************** IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)AX,AY,AZ,NAT,XOFF,YOFF,ZOFF,A1,A2 DO 40 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IX(JX),IY(JX),IZ(JX) 40 CONTINUE CLOSE(UNIT=IOSHP) ENDIF RETURN 9020 FORMAT(' >TARELL ellipsoidal grain; AX,AY,AZ=',3F8.4,/, &I7,3F8.4,' = NAT,XOFF,YOFF,ZOFF',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) END SUBROUTINE TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,NAT,IX,IY,IZ,ICOMP,IDVOUT, & MXN3,NAT3) C C Arguments: CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6 INTEGER IDVOUT,IDVSHP,IOSHP,MXNAT,NAT,MXN3,NAT3 REAL A1(3),A2(3),SHPAR(6) INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT),ICOMP(MXNAT,3) EXTERNAL ERRMSG,REASHP,TARCYL,TARELL,TARHEX,TARREC,TARTET C C Local variables: INTEGER J,J1,L,L1,L2,L3,LC C C*********************************************************************** C Given: C CSHAPE = string describing target geometry; will determine which C "shape" routine is invoked. Current options are: C RCTNGL, ELLIPS, CYLNDR, HEXGON, TETRAH, UNICYL, ANIELL, C TWOELL, TWOAEL, FRMFIL, THRELL, THRAEL, CONELL C IOSHP = device number for "target.out" file C = -1 to suppress printing of "target.out" C SHPAR(1-6) = up to 6 parameters defining target geometry C MXNAT = limit on largest allowed value of NAT C MXN3 = 3*MXNAT C IDVOUT = 1 to print out information on target shape C -1 to suppress printing information on target shape C CFLSHP = name of target shape file if CSHAPE='FRMFIL' C IDVSHP = device number to use to read target shape if C CSHAPE='FRMFIL' C Returns: C A1,A2 = two 3-vectors defining initial target orientation C CDESCR = string describing target C NAT = number of dipoles in target C NAT3 = 3*NAT C IX(1-NAT),IY(1-NAT),IZ(1-NAT)=integers fixing x,y,z coordinates of C occupied sites in target C ICOMP(1-NAT,1-3)=composition for E in x,y,z direction at each site C C P.J.Flatau, Colorado State University C B.T.Draine, Princeton Univ. Observatory C History: C 90.11.08 (BTD): Changed DO...ENDDO to DO #...# CONTINUE C 90.12.02 (BTD): Changed ordering of data in ICOMP C 90.12.14 (BTD): Changed option MKRECT -> RCTNGL C 91.05.02 (BTD): Added IDVOUT to argument list C Added IDVOUT to argument list for subr. REASHP C 91.05.23 (BTD): Added A1,A2 to arg. list for subr. REASHP C 91.11.12 (BTD): Added IDVSHP to arg. list for TARGET and REASHP C 92.09.21 (BTD): Added option UNIELL (uniaxial ellipsoid) C 93.01.07 (BTD): Added option TWOELL (two touching ellipsoids, C consisting of materials 1 and 2) C 92.01.07 (BTD): Added option TWOAEL (two touching anisotropic C ellipsoids, 1st with diel.fnc.1,2,3 along x,y,z C 2nd with diel.fnc.4,5,6 along x,y,z C 92.01.19 (BTD): Added option THRELL (three touching anisotropic C ellipsoids, 1st with diel.fnc.1,2,3 along x,y,z C 2nd with diel.fnc.4,5,6 along x,y,z C 3rd with diel.fnc.7,8,9 along x,y,z C 93.03.12 (BTD): Changed CDESCR*60 -> CDESCR*67 C Moved WRITE(IDVOUT,9010) from DDSCAT to TARGET C 93.12.14 (BTD): Corrected handling of THRELL (had assigned C ICOMP=1 to 50% of sites, ICOMP=2 to other 50%) C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) to C allow up to 6 shape parameters C Added call to TARCEL to generate two concentric C ellipsoids C 94.03.21 (BTD): Changed "CONCEL" to "CONELL" for consistency with C REAPAR C C Copyright (C) 1993,1994 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C******************************************************************** IF(CSHAPE.EQ.'RCTNGL') THEN CALL TARREC(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C Homogeneous target: DO 100 J=1,NAT DO 95 l=1,3 ICOMP(J,L)=1 95 CONTINUE 100 CONTINUE ELSEIF(CSHAPE.EQ.'ELLIPS')THEN CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C Homogeneous target: DO 110 J=1,NAT DO 105 L=1,3 ICOMP(J,L)=1 105 CONTINUE 110 CONTINUE ELSEIF(CSHAPE.EQ.'CYLNDR')THEN CALL TARCYL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C Homogeneous target: DO 120 J=1,NAT DO 115 L=1,3 ICOMP(J,L)=1 115 CONTINUE 120 CONTINUE ELSEIF(CSHAPE.EQ.'HEXGON')THEN CALL TARHEX(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C Homogeneous target: DO 130 J=1,NAT DO 125 L=1,3 ICOMP(J,L)=1 125 CONTINUE 130 CONTINUE ELSEIF(CSHAPE.EQ.'TETRAH')THEN CALL TARTET(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C Homogeneous target: DO 140 J=1,NAT DO 135 L=1,3 ICOMP(J,L)=1 135 CONTINUE 140 CONTINUE ELSEIF(CSHAPE.EQ.'UNICYL')THEN CALL TARCYL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C SHPAR(3)=1,2,3 for cylinder axis in x,y,z directions C Assume c axis to be parallel to cylinder axis C Assume component 1 to be dielectric const. for E par. c axis C Assume component 2 to be dielectric const. for E perp. c axis LC=INT(SHPAR(3)) DO 150 J=1,NAT DO 145 L=1,3 ICOMP(J,L)=2 145 CONTINUE ICOMP(J,LC)=1 150 CONTINUE C ELSEIF(CSHAPE.EQ.'ANIELL')THEN CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT C Dielectric constants 1,2,3 for E along directions a1,a2,a3=a1xa2 C L1=1,2,3 for axis a1 in x,y,z directions (in Target Frame) C L2=1,2,3 for axis a2 in x,y,z directions (in Target Frame) C L3=1,2,3 for axis a3=a1xa2 in x,y,z directions (in Target Frame) L1=A1(1)+2*A1(2)+3*A1(3)+1.E-6 L2=A2(1)+2*A2(2)+3*A2(3)+1.E-6 L3=(A1(2)*A2(3)-A1(3)*A2(2))+2*(A1(3)*A2(1)-A1(1)*A2(3)) & +3*(A1(1)*A2(2)-A1(2)*A2(1))+1.E-6 DO 160 J=1,NAT ICOMP(J,L1)=1 ICOMP(J,L2)=2 ICOMP(J,L3)=3 160 CONTINUE C C TWOELL -> two touching identical ellipsoids, of materials 1 and 2 C second ellipsoid is displaced from first in x-direction ELSEIF(CSHAPE.EQ.'TWOELL')THEN CALL TAR2EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT J1=NAT/2 DO 165 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=1 ICOMP(J,3)=1 165 CONTINUE DO 166 J=J1+1,NAT ICOMP(J,1)=2 ICOMP(J,2)=2 ICOMP(J,3)=2 166 CONTINUE C C TWOAEL -> two touching identical ellipsoids with anisotropic C functions. C second ellipsoid is displaced from first in x-direction C 1st ellipsoid has dielectric fcns 1,2,3 in x,y,z dirs. C 2nd ellipsoid has dielectric fcns 4,5,6 in x,y,z dirs. ELSEIF(CSHAPE.EQ.'TWOAEL')THEN CALL TAR2EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT J1=NAT/2 DO 175 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 175 CONTINUE DO 176 J=J1+1,NAT ICOMP(J,1)=4 ICOMP(J,2)=5 ICOMP(J,3)=6 176 CONTINUE C C THRELL -> three touching identical ellipsoids, of materials 1,2,3 C second ellipsoid is displaced from first in +x-direction C third ellipsoid is displaced from second in +x-direction ELSEIF(CSHAPE.EQ.'THRELL')THEN CALL TAR3EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT J1=NAT/3 DO 185 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=1 ICOMP(J,3)=1 185 CONTINUE DO 186 J=J1+1,2*J1 ICOMP(J,1)=2 ICOMP(J,2)=2 ICOMP(J,3)=2 186 CONTINUE DO 187 J=2*J1+1,NAT ICOMP(J,1)=3 ICOMP(J,2)=3 ICOMP(J,3)=3 187 CONTINUE C C THRAEL -> three touching identical ellipsoids with anisotropic C functions. C second ellipsoid is displaced from first in x-direction C 1st ellipsoid has dielectric fcns 1,2,3 in x,y,z dirs. C 2nd ellipsoid has dielectric fcns 4,5,6 in x,y,z dirs. C 3rd ellipsoid has dielectric fcns 7,8.9 in x,y,z dirs. ELSEIF(CSHAPE.EQ.'THRAEL')THEN CALL TAR3EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IX,IY,IZ) NAT3=3*NAT J1=NAT/3 DO 195 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 195 CONTINUE DO 196 J=J1+1,2*J1 ICOMP(J,1)=4 ICOMP(J,2)=5 ICOMP(J,3)=6 196 CONTINUE DO 197 J=2*J1+1,NAT ICOMP(J,1)=7 ICOMP(J,2)=8 ICOMP(J,3)=9 197 CONTINUE C C CONELL -> two concentric ellipsoids with isotropic dielectric functions C 1 (inner) and 2 (outer) ELSEIF(CSHAPE.EQ.'CONELL')THEN CALL TARCEL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5), & SHPAR(6),CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ,ICOMP) NAT3=3*NAT C C add any other specific shape routines here C ELSEIF(CSHAPE.EQ.'FRMFIL') THEN CALL REASHP(CFLSHP,IDVOUT,IDVSHP,A1,A2,CDESCR,MXNAT,NAT, $ IX,IY,IZ,ICOMP,MXN3,NAT3) ENDIF C Check to see that arrays are large enough IF(NAT.GT.MXNAT)CALL ERRMSG('FATAL','TARGET', & ' MXNAT < NAT; must increase MXNAT in DDSCAT') C Write out target description WRITE(IDVOUT,FMT=9010)CDESCR,NAT C RETURN 9010 FORMAT(' >TARGET:',A,/, & 7X,I7,' = NAT0 = number of dipoles in target') END SUBROUTINE TARHEX(A1,A2,A,B,AORI,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) C Scalar arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT) REAL A,AORI,B C Array arguments: REAL A1(3),A2(3) C Local scalars: CHARACTER CMSGNM*70 INTEGER IAXIS,IFACE,IORI,JA,JX,JY,JZ, $ NA,NB,NC,NFAC,NLAY,NMX,NMY,NMZ LOGICAL IN REAL ACM,ASPR,BCM,BEFF,CCM,PI,X,Y,Z C*********************************************************************** C Routine to construct hexagonal prism from "atoms" C Input: C A=prism length/d C B=width of one hexagonal edge/d C AORI=1.,2.,3.,4.,5.,6. to specify orientation of prism C IOSHP=device number for "target.out" file C =-1 to suppress printing of "target.out" C MXNAT=dimensioning information (max number of atoms) C Output: C A1(1-3)=unit vector defining target axis 1 (parallel to prism) C A2(1-3)=unit vector defining target axis 2 (perp. to A1) C NAT=number of atoms in target C IX(1-NAT)=x/d for atoms of target C IY(1-NAT)=y/d C IZ(1-NAT)=z/d C CDESCR=string describing target (up to 67 chars.) C C B.T.Draine, Princeton Univ. Obs., 87/4/3 C C History: C 90/12/02 (BTD): IOSHP is now output device for target.out C (and is now closed before returning) C 90/12/14 (BTD): Rewritten so that atoms are always on integral sites C but CM coords may be either integral or half-integral C 91/01/05 (BTD): Change I4 -> I7 when printing NAT. C 91/04/22 (BTD): Added XV,YV,ZV and A1,A2 to WRITE(IOSHP,FMT=9020) C 92/09/09 (BTD+PJF): in output, NLAY was too large by factor 2 and NFAC C too small by factor of 2: now fixed. C 93/03/12 (BTD): Changed CDESCR*(*) -> CDESCR*67 C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** PI=4.*ATAN(1.) IORI=AORI+1.e-6 C IORI =1 for symmetry axis along X axis, face normal to Y axis. C 2 X Z C 3 Y Z C 4 Y X C 5 Z X C 6 Z Y C Determine target axis vectors A1 and A2 C Convention: A1 will be along prism symmetry axis C A2 will be perp. to A1 along vertex direction DO 0100 JA=1,3 A1(JA)=0. A2(JA)=0. 0100 CONTINUE IF(IORI.EQ.1)THEN IAXIS=1 IFACE=2 ELSEIF(IORI.EQ.2)THEN IAXIS=1 IFACE=3 ELSEIF(IORI.EQ.3)THEN IAXIS=2 IFACE=3 ELSEIF(IORI.EQ.4)THEN IAXIS=2 IFACE=1 ELSEIF(IORI.EQ.5)THEN IAXIS=3 IFACE=1 ELSEIF(IORI.EQ.6)THEN IAXIS=3 IFACE=2 ENDIF a1(1)=1. a2(2)=1. C C IAXIS=1,2,3 for symmetry axis along x,y,z axis C IFACE=1,2,3 for "face" normal to x,y,z axis C A = prism length C B = prism "diameter" C B/2 = length of one side of hexagon C B= distance between opposite vertices=max.diameter of hexagon C for hexagon, area=(3/8)*SQRT(3)*B**2 C for hex prism, volume=A*area=(3/8)*SQRT(3)*A*B**2 C Now determine limits for testing x,y,z values C Along axis, run from 1 to NA=INT(A+0.5) C Perp. to axis, run from 1 to NB=INT(B+0.5) NA=INT(A+0.5) NB=INT(B+0.5) NC=INT(.5*SQRT(3.)*B+0.5) C ACM="center" of figure in axial direction C BCM="center" of figure in B direction C CCM="center" of figure in C direction ACM=REAL(NA)/2.+0.5 BCM=REAL(NB)/2.+0.5 CCM=REAL(NC)/2.+0.5 IF(IAXIS.EQ.1)THEN NMX=NA NMY=NB NMZ=NB IF(IFACE.EQ.2)NMY=NC IF(IFACE.EQ.3)NMZ=NC ELSEIF(IAXIS.EQ.2)THEN NMY=NA NMX=NB NMZ=NB IF(IFACE.EQ.1)NMX=NC IF(IFACE.EQ.3)NMZ=NC ELSEIF(IAXIS.EQ.3)THEN NMZ=NA NMX=NB NMY=NB IF(IFACE.EQ.1)NMX=NC IF(IFACE.EQ.2)NMY=NC ENDIF NAT=0 write(0,*)'nmx,nmy,nmz=',nmx,nmy,nmz write(0,*)'acm,bcm,ccm=',acm,bcm,ccm write(0,*)'iaxis,iface=',iaxis,iface DO 1200 JZ=1,NMZ Z=REAL(JZ) DO 1100 JY=1,NMY Y=REAL(JY) DO 1000 JX=1,NMX X=REAL(JX) CALL TESTHEX(IAXIS,IFACE,ACM,BCM,CCM,A,B,X,Y,Z,IN) IF(IN)THEN NAT=NAT+1 IX(NAT)=JX IY(NAT)=JY IZ(NAT)=JZ ENDIF 1000 CONTINUE 1100 CONTINUE 1200 CONTINUE C Now compute some geometric properties of pseudo-hexagonal prism C NLAY = number of layers = effective length of prism/d C NFAC = number of atoms in one hexagonal face NLAY=INT(A) NFAC=NAT/NLAY C BEFF = effective length of hexagonal size/d C computed for hexagon of area NFAC*d**2 C =(3/2)*SQRT(3)*BEFF**2 C Note: for hexagon, vertex-vertex diameter = 2*BEFF BEFF=.6204032*SQRT(FLOAT(NFAC)) C ASPR = effective aspect ratio of target = length/(2*BEFF) ASPR=.5*FLOAT(NLAY)/BEFF C Note: string CDESCR will be printed by subr. TARGET WRITE(CDESCR,FMT='(A,I7,A)')' Hexagonal prism of NAT=',NAT, & ' dipoles' C Here print any additional target info which is desired by using C subroutine WRIMSG WRITE(CMSGNM,FMT='(A,3F7.4,A)') $ ' A1 = (',A1,') = hexagon symmetry axis' CALL WRIMSG('TARHEX',CMSGNM) WRITE(CMSGNM,FMT='(A,3F7.4,A)') $ ' A2 = (',A2,') = center-vertex axis' CALL WRIMSG('TARHEX',CMSGNM) WRITE(CMSGNM,FMT='(A,I7,A,I4,A,I4,A,F7.4)') $ ' NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY,' aspect ratio=',ASPR CALL WRIMSG('TARHEX',CMSGNM) WRITE(CDESCR,FMT='(A,I7,A,I4,A,I4,A,F7.4)') $ ' Hex prism,NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY, $ ' asp.ratio=',ASPR C NOW STORE ATOM COORDINATES IN FILE IF(IOSHP.GT.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') C*** 3 lines of general description allowed: WRITE(IOSHP,9020)A,B,NAT,IAXIS,IFACE, A1, A2 C*** DO 2000 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IX(JX),IY(JX),IZ(JX) 2000 CONTINUE CLOSE(UNIT=IOSHP) ENDIF RETURN 9020 FORMAT(' Hexagonal Prism: A=',F7.4,' B=',F7.4,/, &I7,'= NAT. Orientation: IAXIS=',I1,' IFACE=',I2,/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) END SUBROUTINE TESTHEX(IAXIS,IFACE,ACM,BCM,CCM,A,B,X,Y,Z,IN) C Arguments: INTEGER IAXIS,IFACE REAL ACM,BCM,CCM,A,B,X,Y,Z LOGICAL IN C Local variables: REAL AX,Q,U,V C C*** ACM,BCM,CCM=location of center in A,B,C directions C A direction is along axis C B direction is vertex to vertex of hexagon C C direction is face to face of hexagon C IF(IAXIS.EQ.1)THEN AX=X IF(IFACE.EQ.2)THEN U=Y V=Z ELSE U=Z V=Y ENDIF ELSEIF(IAXIS.EQ.2)THEN AX=Y IF(IFACE.EQ.1)THEN U=X V=Z ELSE U=Z V=X ENDIF ELSEIF(IAXIS.EQ.3)THEN AX=Z IF(IFACE.EQ.1)THEN U=X V=Y ELSE U=Y V=X ENDIF ENDIF IN=.FALSE. C*** Test along axis: IF(2.*ABS(AX-ACM).GT.A)RETURN C*** Test along C direction C .4330127=SQRT(3)/4 IF(ABS(U-CCM).GT.0.4330127*B)RETURN C*** Now test whether closer to CM than line bounding edge C .5773503=1/SQRT(3) Q=0.5*ABS(U-CCM)+0.8660254*ABS(V-BCM) IF(Q.GT.0.4330127*B)RETURN IN=.TRUE. RETURN END