SUBROUTINE TAR2EL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) C*********************************************************************** C Routine to construct two touching ellipsoids from "atoms" C Input: C AX=(x-length of one ellipsoid)/d (d=lattice spacing) C AY=(y-length of one ellipsoid)/d C AZ=(z-length of one ellipsoid)/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 C Note: atoms 1 - NAT/2 are in first ellipsoid C NAT/2+1 - NAT are in second ellipsoid C C B.T.Draine, Princeton Univ. Obs. C History: C 92/01/07 (BTD) adapted from tarell.f 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: CHARACTER CMSGNM*70 INTEGER JA,JX,JXMAX,JXMIN,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NAT2 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 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 C 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 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 A1(1)=1. A2(2)=1. C C Determine list of occupied sites in first ellipsoid C 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 IF(NAT.GT.1)THEN IF(JX.GT.JXMAX)JXMAX=JX IF(JX.LT.JXMIN)JXMIN=JX ELSE JXMAX=JX JXMIN=JX ENDIF ENDIF 0100 CONTINUE ENDIF 0200 CONTINUE ENDIF 0300 CONTINUE C C Now create duplicate second ellipsoid JXMAX=JXMAX-JXMIN+1 DO 0400 JA=1,NAT IX(NAT+JA)=IX(JA)+JXMAX IY(NAT+JA)=IY(JA) IZ(NAT+JA)=IZ(JA) 0400 CONTINUE NAT=2*NAT IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ') ENDIF C*********************************************************************** C Write target description into string CDESCR NAT2=NAT/2 IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN WRITE(CDESCR,FMT='(A,I7,A)')' Two spheres, each containing', & NAT2,' dipoles' ELSE WRITE(CDESCR,FMT='(A,I7,A)')' Two ellipsoids, each containing', & NAT2,' dipoles' ENDIF C*********************************************************************** CMSGNM=CDESCR CALL WRIMSG('TARELL',CMSGNM) 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 TAR3EL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) C*********************************************************************** C Routine to construct three collinear touching ellipsoids from "atoms" C Input: C AX=(x-length of one ellipsoid)/d (d=lattice spacing) C AY=(y-length of one ellipsoid)/d C AZ=(z-length of one ellipsoid)/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 Note: atoms 1 - NAT/3 are in first ellipsoid C NAT/3+1 - 2*NAT/3 are in second ellipsoid C 2*NAT/3+1 - NAT are in third ellipsoid C Ellipsoids are displaced from one another in x-direction C First ellipsoid is at smallest x values, third at largest C B.T.Draine, Princeton Univ. Obs. C C History: C 93/01/20 (BTD) adapted from tar2el.f 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: CHARACTER CMSGNM*70 INTEGER JA,JX,JXMAX,JXMIN,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NAT2 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 C 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 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 A1(1)=1. A2(2)=1. C C Determine list of occupied sites in first ellipsoid C 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 IF(NAT.GT.1)THEN IF(JX.GT.JXMAX)JXMAX=JX IF(JX.LT.JXMIN)JXMIN=JX ELSE JXMAX=JX JXMIN=JX ENDIF ENDIF 0100 CONTINUE ENDIF 0200 CONTINUE ENDIF 0300 CONTINUE C C Now create duplicate second and third ellipsoids JXMAX=JXMAX-JXMIN+1 DO 0400 JA=1,NAT IX(NAT+JA)=IX(JA)+JXMAX IY(NAT+JA)=IY(JA) IZ(NAT+JA)=IZ(JA) IX(2*NAT+JA)=IX(JA)+2*JXMAX IY(2*NAT+JA)=IY(JA) IZ(2*NAT+JA)=IZ(JA) 0400 CONTINUE NAT=3*NAT IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARELL',' NAT.GT.MXNAT ') ENDIF C*********************************************************************** C Write target description into string CDESCR NAT2=NAT/3 IF(AX.EQ.AY.AND.AX.EQ.AZ)THEN WRITE(CDESCR,FMT='(A,I7,A)')' Three spheres, each containing', & NAT2,' dipoles' ELSE WRITE(CDESCR,FMT='(A,I7,A)')' Three ellipsoids, each ', & 'containing',NAT2,' dipoles' ENDIF C*********************************************************************** CMSGNM=CDESCR CALL WRIMSG('TARELL',CMSGNM) 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 TARCEL(A1,A2,AX,AY,AZ,BX,BY,BZ,CDESCR,IOSHP,MXNAT,NAT, & IX,IY,IZ,ICOMP) C*********************************************************************** C Routine to construct target consisting of two materials, with outer C surface an allipsoid of dimensions AX,AY,AZ, C and core/mantle interface an ellipsoid of dimensions BX,BY,BZ C Input: C AX=(x-length)/d of outer ellipsoid (d=lattice spacing) C AY=(y-length)/d " " " C AZ=(z-length)/d " " " C BX=(x-length)/d of inner ellipsoid C BY=(y-length)/d " " " C BZ=(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 ICOMP(1-NAT,1-3)=1 for sites within inner ellipsoid, C =2 for sites between inner and outer ellipsoids C CDESCR=description of target (up to 67 characters) C C B.T.Draine, Princeton Univ. Obs. C History: C 94.01.26 (AP) : Adapted from TARELL by Antonio Peimbert, Princeton U. C 94.01.27 (BTD): Modified for compatibility with TARGET, using C AX,AY,AZ,BX,BY,BZ to transfer shape parameters read by C subroutine REAPAR C 94.03.21 (BTD): Corrected format statement C C Copyright (C) 1994 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),ICOMP(MXNAT,3) REAL AX,AY,AZ,BX,BY,BZ REAL A1(3),A2(3) C** Local variables: INTEGER JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NIN REAL AX2,AY2,AZ2,BX2,BY2,BZ2,R,RR,RYZ2,RZ2,X,XOFF,Y,YOFF,Z,ZOFF C*********************************************************************** C C Inner ellipsoid is defined by C (X/AX)**2+(Y/AY)**2+(Z/AZ)**2 = 0.25 C Outer ellipsoid is defined by C (X/BX)**2+(Y/BY)**2+(Z/BZ)**2 = 0.25 C Material within inner ellipsoid is of composition 1 C Material between inner and outer ellipsoids is of composition 2 C C Ideal volume V=(pi/6)*BX*BY*BZ 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 concentric spheres: call with AX=AY=AZ, BX=BY=BZ C For spheroids: call with AX=AY (or AY=AZ or AX=AZ) and BX=BY (or BY=BZ C or BX=BZ) C C*********************************************************************** C C Check that input parameters have "inner" ellipsoid smaller than C "outer" ellipsoid: IF(AX.LT.BX)CALL ERRMSG('FATAL','TARCEL',' AX < BX ') IF(AY.LT.BY)CALL ERRMSG('FATAL','TARCEL',' AY < BY ') IF(AZ.LT.BZ)CALL ERRMSG('FATAL','TARCEL',' AZ < BZ ') C Criterion for choosing XOFF: try to optimize outer ellipsoid surface 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 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 BX2=BX*BX BY2=BY*BY BZ2=BZ*BZ NAT=0 NIN=0 C C Specify target axes A1 and A2 C 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 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 RR=Z*Z/BZ2+Y*Y/BY2+X*X/BX2 IF(RR.LT.0.25)THEN C within inner ellipsoid: ICOMP(NAT,1)=1 ICOMP(NAT,2)=1 ICOMP(NAT,3)=1 NIN=NIN+1 ELSE C between inner and outer ellipsoids: ICOMP(NAT,1)=2 ICOMP(NAT,2)=2 ICOMP(NAT,3)=2 ENDIF 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,I6,A,I7,A)')' Concentric ellipsoids with', & NIN,' inner ',NAT,' total dipoles' ENDIF C*********************************************************************** IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)AX,AY,AZ,BX,BY,BZ,NAT,NIN,XOFF,YOFF,ZOFF, & A1,A2 DO 40 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IX(JX),IY(JX),IZ(JX),ICOMP(JX,1), & ICOMP(JX,2),ICOMP(JX,3) 40 CONTINUE CLOSE(UNIT=IOSHP) ENDIF RETURN 9020 FORMAT(' >TARCEL: concentric ellipsoids; AX,AY,AZ=',3F8.4, & ' BX,BY,BZ=',3F8.4,/, & I7,I7,3F8.4,' = NAT,NIN,XOFF,YOFF,ZOFF',/, & 3F9.4,' = A_1 vector',/, & 3F9.4,' = A_2 vector',/, & ' JA IX IY IZ ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,3I2) END SUBROUTINE TARREC(A1,A2,XV,YV,ZV,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) C*********************************************************************** C Routine to construct rectangular prism from "atoms" C Input: C XV=(x-length)/d (d=lattice spacing) C YV=(y-length)/d C ZV=(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)=unit vector (1,0,0) defining target axis 1 C A2(1-3)=unit vector (0,1,0) defining target axis 2 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 C B.T.Draine, Princeton Univ. Obs. C History: 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 93.03.12 (BTD): Specified CDESCR*67, and now use it. C 95.07.12 (BTD): Changed definition of axes A1,A2. They used to be C longest,next longest dimention. They are now C fixed to always be (1,0,0) and (0,1,0). 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 Arguments: REAL XV,YV,ZV INTEGER IOSHP,MXNAT,NAT CHARACTER CDESCR*67 INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT) REAL A1(3),A2(3) C C Local variables: INTEGER JA,JX,JY,JZ,NX,NY,NZ CHARACTER CMSGNM*70 C C External subroutines: EXTERNAL ERRMSG C C Intrinsic functions: INTRINSIC IFIX C*********************************************************************** NX=IFIX(XV+0.5) NY=IFIX(YV+0.5) NZ=IFIX(ZV+0.5) WRITE(CMSGNM,FMT='(A,3I4)') & ' Rectangular prism; NX,NY,NZ=',NX,NY,NZ WRITE(CDESCR,FMT='(A,3I4)') & ' Rectangular prism; NX,NY,NZ=',NX,NY,NZ C C Specify target axes A1 and A2 C New convention: A1 will be (1,0,0) in target frame C A2 will be (0,1,0) in target frame C Convention: A1 will be along longest dimension C A2 will be along intermediate dimension DO 0100 JA=1,3 A1(JA)=0. A2(JA)=0. 0100 CONTINUE A1(1)=1. A2(2)=1. C C Now populate lattice: NAT=0 DO 30 JZ=1,NZ DO 20 JY=1,NY DO 10 JX=1,NX NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT') ENDIF IX(NAT)=JX IY(NAT)=JY IZ(NAT)=JZ 10 CONTINUE 20 CONTINUE 30 CONTINUE IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT ') ENDIF WRITE(CMSGNM,FMT='(A, I7, A)') & ' Rectangular prism of NAT=',NAT,' dipoles' IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)XV,YV,ZV,NAT,A1,A2 DO 40 JA=1,NAT WRITE(IOSHP,FMT=9030)JA,IX(JA),IY(JA),IZ(JA) 40 CONTINUE CLOSE (UNIT=IOSHP) ENDIF RETURN 9020 FORMAT (' >TARREC rectangular prism; AX,AY,AZ=',3F8.4,/, &I7,' = NAT',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) END SUBROUTINE TARTET(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IX,IY,IZ) 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 I,IMAX,IMIN,J,JMAX,JMIN,JX,K,KMAX,KMIN REAL FY,S,X,XOFF,Y,YMAX,YMIN,YOFF,Z,ZMAX,ZMAX0,ZOFF C*********************************************************************** C C Routine to construct regular tetrahedral target array C one of the faces is parallel to the y-z plane C one of the edges of this face is parallel to x-y plane C C Input: C AX = length of one edge of tetrahedron 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 (1,0,0) (along one axis of tetrahedron) C A2(1-3) = unit vector (0,1,0) 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 Occupied array sites are those within tetrahedral surface C Size: C S=length of each side of tetrahedron C Volume = S**3/(6*sqrt(2)) C Orientation: C Center of mass is at origin. C One face is parallel to yz plane. C Projection onto yz plane has vertex on y axis. C S=length of one side (in lattice units) C Vertices are at C A=( sqrt(3/8), 0, 0)*S C B=(-sqrt(1/24), sqrt(1/3), 0)*S C C=(-sqrt(1/24),-sqrt(1/12),-1/2)*S C D=(-sqrt(1/24),-sqrt(1/12), 1/2)*S C Length in x direction = S*sqrt(2/3) C y = S*sqrt(3)/2 C z = S C Angle(AOB)=arccos(-1/3)=109.4712 deg. C OA=OB=OC=OD=sqrt(3/8)*S C Occupied sites are assumed to be located at C (X,Y,Z)=(I+XOFF,J+YOFF,K+ZOFF) C where I,J,K are integers, and XOFF,YOFF,ZOFF are constants. C Program sets XOFF,YOFF,ZOFF depending on choice of parameter S. C C Criterion for choosing XOFF: C Base of tetrahedron is located at x = -sqrt(1/24)*S C Let IMIN be value of I for this plane C Choose XOFF so that IMIN+XOFF = -sqrt(1/24)*S + 0.5 C with -0.5 < XOFF < 0.5 C Criterion for choosing YOFF: C One edge of tetrahedron is located at y= -S/(2*sqrt(3)) C Let JMIN be value of J for this line C Choose YOFF so that JMIN+YOFF = -S/sqrt(12) + 0.5 C Criterion for choosing ZOFF: C One edge of tetrahedron is parallel to z axis C Choose ZOFF in order to have number of dipoles along C this edge as close as possible to S C e.g., if S=odd integer, then take ZOFF=0 to place dipoles C at integral locations C if S=even integer, then take ZOFF=1/2 to place dipoles C at half-integral locations C C B.T.Draine, Princeton Univ. Obs., 90/10/30 C History: C 90/11/ 8 (BTD): Changed DO...ENDDO to DO #...# CONTINUE C 90/11/ 9 (BTD): Corrected error in location of center of mass. 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 C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** S=AX C Set XOFF (and IMIN,IMAX): IMIN=-INT(S*SQRT(1./24.)) XOFF=0.5-S*SQRT(1./24.)-IMIN IMAX=IMIN+INT(S*SQRT(2./3.)+0.5)-1 C Set YOFF (and JMIN,JMAX): JMIN=-INT(S/SQRT(12.)) YOFF=0.5-S/SQRT(12.)-JMIN JMAX=JMIN+INT(S*SQRT(.75)+0.5)-1 C Set ZOFF (and KMIN,KMAX): C Determine whether S is closest to even or odd integer. C (Temporarily let KMIN be integer which S is closest to) KMIN=INT(S+0.5) C If KMIN is even, then ZOFF=0.5 C If KMIN is odd, then ZOFF=0. ZOFF=0. IF(KMIN-2*(KMIN/2).EQ.0)ZOFF=0.5 KMIN=-INT(0.5*S+ZOFF) KMAX=KMIN+INT(S+0.5)-1 C C Determine list of occupied sites. NAT=0 DO 30 I=IMIN,IMAX X=REAL(I)+XOFF C YMAX=largest value of Y which can occur for this X value C YMIN=smallest value of Y which can occur for this X value C ZMAX0=largest value of Z which can occur for this X value YMAX=S*SQRT(3./16.)-X/SQRT(2.) YMIN=-S*SQRT(3./64.)+X/SQRT(8.) ZMAX0=3.*S/8.-X*SQRT(3./8.) DO 20 J=JMIN,JMAX Y=REAL(J)+YOFF IF(Y.GE.YMIN.AND.Y.LE.YMAX)THEN C ZMAX=largest value of Z which can occur for this (X,Y) FY=(Y-YMIN)/(YMAX-YMIN) ZMAX=(1.-FY)*ZMAX0 DO 10 K=KMIN,KMAX Z=REAL(K)+ZOFF IF(ABS(Z).LE.ZMAX)THEN C Site is occupied: NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARTET', & ' NAT.GT.MXNAT ') ENDIF IX(NAT)=I IY(NAT)=J IZ(NAT)=K ENDIF 10 CONTINUE ENDIF 20 CONTINUE 30 CONTINUE C C*** Specify vectors A1 and A2 which will define target orientation C A1 is initially along x-axis C A2 is initially along y-axis A1(1)=1. A1(2)=0. A1(3)=0. A2(1)=0. A2(2)=1. A2(3)=0. C WRITE(CDESCR,FMT='(A,I7,A)')' Tetrahedron of NAT=',NAT,' dipoles' IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)S,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(' >TARTET tetrahedral grain: S=',F8.3,/, &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 TIMEIT(CMSGTM,DTIME) C C timeit_vms C 94/06/20 (PJF) add dtime to the formal parameters C C This version of TIMEIT is for VMS systems. 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 .. Scalar Arguments .. CHARACTER CMSGTM* (*) C .. C .. Local Scalars .. REAL DTIME, OLDTIM, TIME INTEGER IAD CHARACTER CSTA*3, CMSGNM*70 C .. C .. External Subroutines .. EXTERNAL LIB$INIT_TIMER, LIB$SHOW_TIMER, WRIMSG C .. C .. Data statements .. DATA CSTA/'ONE'/ C .. IF (CSTA .EQ. 'ONE') THEN CSTA = 'TWO' IAD = 0 CALL LIB$INIT_TIMER(IAD) ELSE IF (CSTA .EQ. 'TWO') THEN WRITE (CMSGNM,FMT='(A, A)') 'TIMING RESULTS OF: ', CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) CSTA = 'ONE' CALL LIB$SHOW_TIMER(IAD,2) END IF RETURN END SUBROUTINE VERSION(CSTAMP) C Arguments: CHARACTER CSTAMP*11 C******************************************************************************* C History of major changes to entire DDSCAT package: C 90.12.21 (BTD): Added code to DDSCAT to create output files C 'qtable' (summary of Q values) C 'mtable' (summary of diel. func. for material 1) C 91.08.14 (BTD): Add QBKSCA to arg. list for GETFML C Modify code to print out Q_bk C 91.08.15 (BTD): Provide different headings for qtable depending on C whether IORTH=1 or 2 (add statement 9043 to DDSCAT). C 91.09.17 (BTD): Remove calls to ALPHA and EVALA in DDSCAT (needed to C move these to subroutine GETFML since alpha from Lattice C Dispersion Relation depends on propagation direction C and polarization state) C Add CALPHA,CXEPS,ICOMP,MXCOMP to argument list for C GETFML (information required by ALPHA) C 91.11.12 (BTD): Added variable IDVSHP to argument list for TARGET C (device number for REASHP to use in reading shape file) C 93.01.15 (BTD): Modified DDSCAT to divided output file qtable into 2 files: C qtable (containing Q_ext,Q_abs,Q_sca,g,Q_bk) and C qtable2 (containing Q_pha, Q_pol, Q_cpol) C 93.01.16 (BTD): Modify DDSCAT so that qtable, qtable2, and mtable are closed C after writing, and reopened for each new write C 93.01.20 (BTD): Add MXNX, MXNY, MXNZ to argument list for subroutine C EXTEND to permit checking of target size against C maximum allowed dimensions C 93.03.11 (BTD): Deleted all code associated with unused variable C CMACHN (originally included to identify machine/OS) C 93.03.12 (BTD): Changed CDESCR*60 -> CDESCR*67 C Moved WRITE(IDVOUT,9010) to subr. TARGET C 93.06.02 (BTD): Corrected error in FORMAT statement 9045 in DDSCAT C 93.06.03 (BTD): Corrected error in computation of angle-averaged C backscattering cross section QBKSUM(JO) in DDSCAT (had C neglected to reset sum to zero for each new target). C 93.09.28 (BTD): Add "fix" in DDSCAT.f to work around Sun compiler/OS C bug (see description below. C 93.12.15 (BTD): Corrected calculation of PPOL in DDSCAT.f Added C comments on correspondence between our notation and C elements of 2x2 amplitude scattering matrix and 4x4 C Mueller scattering matrix. C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) to C allow up to 6 shape parameters in DDSCAT, REAPAR, C TARGET. C 94.06.20 (PJF): Version 4c.1: C Modifications to various routines to support use of C GPFA FFT code of Temperton by specifying option C NEWTMP in ddscat.par. Changes made in C cxfft3 C eself C extend C reapar C Add FFT timing code TSTFFT to distribution. C Change made to timeit_xxx routines for compatibility C with TSTFFT C 94.12.20 (BTD): Add subroutine VERSION to consolidate log of C significant changes to DDSCAT package. C 95.04.07 (BTD): REASHP: Changed CDESCR*60->CDESCR*67 for compatibility C with calling routine TARGET C TARCYL: corrected error in code C 95.05.15 (BTD): REASHP: Corrected READ statement C added module for inhomogeneous/anisotropic C targets. C 95.05.26 (BTD): REAPAR: corrected construction of orthogonal C polarization when E01 is complex (e.g., C circular polarization). C 95.06.14 (BTD): REAPAR: added code to check that user is not requesting C computation of f_ml when using other than C linearly polarized E01 C C Copyright (C) 1993,1994,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** CSTAMP='DDSCAT.4c.1' RETURN END SUBROUTINE WRIMSG(CSUBRT,CMSGNM) C Standard procedure for writing messages 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 .. Scalar Arguments .. CHARACTER CMSGNM* (*), CSUBRT* (*) C .. C .. Local Scalars .. INTEGER IOOUT C .. C .. External Subroutines .. EXTERNAL GETSET C .. CALL GETSET('GET','IOOUT',IOOUT) WRITE (IOOUT,FMT=9000) CSUBRT, CMSGNM 9000 FORMAT (' >',A,' ',A) RETURN END SUBROUTINE XNORM3(A,XN) C auxilary for "rotation" routines 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 .. Scalar Arguments .. REAL XN C .. C .. Array Arguments .. REAL A(3) C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. XN = SQRT(A(1)**2+A(2)**2+A(3)**2) RETURN END