SUBROUTINE TAR2EL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IXYZ) IMPLICIT NONE C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3) 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 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target 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 96.01.26 (BTD) Replaced IX,IY,IZ by IXYZ C 00.11.02 (BTD) modified to print composition to target.out C end history C C Copyright (C) 1993,1996,2000 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C 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 JXMAX=LMX1 JXMIN=LMX2 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 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=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 IXYZ(NAT+JA,1)=IXYZ(JA,1)+JXMAX IXYZ(NAT+JA,2)=IXYZ(JA,2) IXYZ(NAT+JA,3)=IXYZ(JA,3) 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 JX=1,NAT2 WRITE(IOSHP,FMT=9031)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO DO JX=NAT2+1,NAT WRITE(IOSHP,FMT=9032)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO 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 ICOMP(x,y,z)') 9031 FORMAT(I7,3I4,' 1 1 1') 9032 FORMAT(I7,3I4,' 2 2 2') END SUBROUTINE TAR2SP(A1,A2,A_1,B_1,A_2,B_2,PHI,PRINAX, & CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) IMPLICIT NONE C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) REAL A_1,A_2,B_1,B_2,PHI,PRINAX REAL A1(3),A2(3) C** Local variables: CHARACTER CMSGNM*70 INTEGER JA,JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2, & NAT1,NAT2 REAL AX2,AY2,AZ2,COSPHI, & R,RX2,RYZ2,RZ2,SINPHI,U,V, & X,XC1,XC2,Y,YC1,YC2,Z,ZC1,ZC2 C*********************************************************************** C Routine to construct two touching spheroids from "atoms" C First spheroid has symmetry axis in y direction. C Second spheroid is displaced in x direction, with symmetry C axis in yz plane, in direction ey*cos(phi)+ez*sin(phi). C Separation between centroids=(B_1+B_2)/2 C C Input: C A_1=length/d of 1st spheroid along symm.axis (d=lattice spacing) C A_2=length/d of 2nd spheroid along symm.axis C B_1=diameter/d of 1st spheroid C B_2=diameter/d of 2nd spheroid C PHI=angle (deg) specifying relative orientation of spheroids C PRINAX=0. to set A1=(1,0,0),A2=(0,1,0) in Target Framee C =1. to set A1,A2=principal axes of largest and second C largest moment of inertia 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target C ICOMP(1-NAT,1-3)=1 for sites within first spheroid C =2 for sites within second spheroid C CDESCR=description of target (up to 67 characters) C C B.T.Draine, Princeton Univ. Obs. C History: C 96.01.25 (BTD) First written. C 96.01.26 (BTD) Replaced IX,IY,IZ by IXYZ C 96.01.29 (BTD) Modified to add PRINAX to argument list, and C call to PRINAXIS to determine principal axes if C PRINAX=1. C 99.06.29 (BTD) Corrected to set ICOMP=2 in second spheroid C (previously it was incorrectly setting ICOMP=1 C in both spheroids) C 99.06.30 (BTD) Substantial modifications to change way two C spheroids are constructed. With new prescription, C should obtain a symmetric structure when called C with A_1,B_1 = A_2,B_2 and PHI=0 or 90 C 00.11.02 (BTD) Modified to write ICOMP to target.out C end history C C Copyright (C) 1996,1999,2000 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C*********************************************************************** C C Routine to construct target array representing two spheroids in C contact. C line connecting spheroid centers is parallel to x C symmetry axis of first spheroid is parallel to y C symmetry axis of second spheroid is in y-z plane, at angle PHI to y C C First spheroid: C Array sites contained within spheroidal surface defined by C ((X-XC1)/B_1)**2+((Y-YC1)/A_1)**2+((Z-ZC1)/B_1)**2=0.25 C Ideal volume V=(pi/6)*A_1*B_1*B_1 C C Second spheroid: C ((X-XC2)/B_2)**2+((U-UC1)/A_1)**2+((V-VC1)/B_1)**2=0.25 C where C U = Y*cos(PHI) + Z*sin(PHI) C V = -Y*sin(PHI) + Z*cos(PHI) C UC1 = YC1*cos(PHI) + ZC1*sin(PHI) C VC1 =-YC1*sin(PHI) + ZC1*cos(PHI) C C Dipoles are located at sites C (x,y,z)=(I,J,K)*d, I,J,K=integers C C Criterion for choosing XC1: C Point where spheroids contact each other should be midway between C lattice points. Thus take XC1+0.5*B_1 = 0.5 C or XC1 = 0.5(1.-B_1) C C Criterion for choosing YC1 and ZC1: C If XC1 is close to an integer, take YC1=ZC1=0 C If XC1 is close to a half-integer, take YC1=ZC1=0.5 C C Note that other criteria are obviously possible. For example, C could have chosen YC1 and ZC1 to maximize number of lattice C sites falling within spheroidal surfaces. C XC1=0.5*(1.-B_1) IF(XC1-INT(XC1).LT.0.25.OR.XC1-INT(XC1).GT.0.75)THEN YC1=0. ZC1=0. ELSE YC1=0.5 ZC1=0.5 ENDIF C NAT=0 C C Determine list of occupied sites in first spheroid C LMX1=INT(XC1-.5*B_1) LMX2=INT(XC1+.5*B_1) LMY1=INT(YC1-.5*A_1) LMY2=INT(YC1+.5*A_1) LMZ1=INT(ZC1-.5*B_1) LMZ2=INT(ZC1+.5*B_1) AX2=B_1*B_1 AY2=A_1*A_1 AZ2=B_1*B_1 DO 0300 JZ=LMZ1,LMZ2 Z=REAL(JZ)-ZC1 RZ2=Z*Z/AZ2 IF(RZ2.LT.0.25)THEN DO 0200 JY=LMY1,LMY2 Y=REAL(JY)-YC1 RYZ2=RZ2+Y*Y/AY2 IF(RYZ2.LT.0.25)THEN DO 0100 JX=LMX1,LMX2 X=REAL(JX)-XC1 R=RYZ2+X*X/AX2 IF(R.LT.0.25)THEN NAT=NAT+1 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDIF 0100 CONTINUE ENDIF 0200 CONTINUE ENDIF 0300 CONTINUE NAT1=NAT C C Determine occupied sites in second spheroid C XC2=XC1+.5*(B_1+B_2) YC2=YC1 ZC2=ZC1 LMX1=INT(XC2-.5*B_2) LMX2=INT(XC2+.5*B_2) C C This is not time-consuming, so no need to optimize choices of C LMY1,LMY2,LMZ1,LMZ2 C LMY1=INT(YC1-.5*MAX(A_2,B_2)) LMY2=INT(YC1+.5*MAX(A_2,B_2)) LMZ1=INT(ZC1-.5*MAX(A_2,B_2)) LMZ2=INT(ZC1+.5*MAX(A_2,B_2)) SINPHI=SIN(3.1415927*PHI/180.) COSPHI=COS(3.1415927*PHI/180.) C C transform (y,z) -> (u,v) with u=y*cos(phi)-z*sin(phi) C v=y*sin(phi)+z*cos(phi) C DO 0600 JX=LMX1,LMX2 X=REAL(JX)-XC2 RX2=(X/B_2)**2 IF(RX2.LT.0.25)THEN DO 0500 JY=LMY1,LMY2 Y=REAL(JY)-YC2 DO 0400 JZ=LMZ1,LMZ2 Z=REAL(JZ)-ZC2 U=(COSPHI*Y-SINPHI*Z)/A_2 V=(SINPHI*Y+COSPHI*Z)/B_2 R=RX2+U**2+V**2 IF(R.LT.0.25)THEN C Site is occupied: NAT=NAT+1 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDIF 0400 CONTINUE 0500 CONTINUE ENDIF 0600 CONTINUE NAT2=NAT-NAT1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TAR2SP',' NAT.GT.MXNAT ') ENDIF C DO 0900 JX=1,3 DO 0700 JA=1,NAT1 ICOMP(JA,JX)=1 0700 CONTINUE DO 0800 JA=NAT1+1,NAT ICOMP(JA,JX)=2 0800 CONTINUE 0900 CONTINUE C C Specify target axes A1 and A2 C If PRINAX=0. then C A1=(1,0,0) in target frame C A2=(0,1,0) in target frame C If PRINAX=1. then C A1,A2 are principal axes of largest,second largest moment of C inertia IF(PRINAX.LE.0.)THEN DO 0010 JX=1,3 A1(JX)=0. A2(JX)=0. 0010 CONTINUE A1(1)=1. A2(2)=1. ELSE CALL PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,A1,A2) ENDIF C C*********************************************************************** C Write target description into string CDESCR WRITE(CDESCR,FMT='(A,I7,A,I7,A,I7,A)') & ' Spheroids with',NAT1,' +',NAT2,'=',NAT,' dipoles' C*********************************************************************** CMSGNM=CDESCR CALL WRIMSG('TAR2SP',CMSGNM) IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)A_1,B_1,A_2,B_2,PHI,NAT,NAT1,NAT2,A1,A2 DO JX=1,NAT WRITE(IOSHP,FMT=9031)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3), & ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3) ENDDO CLOSE(UNIT=IOSHP) ENDIF RETURN 9020 FORMAT(' >TAR2SP: touching spheroids, A_1,B_2,A_2,B_2,PHI=',5F8.4, &/,3I7,' = NAT,NAT1,NAT2',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ ICOMP(x,y,z)') 9031 FORMAT(I7,3I4,3I2) END SUBROUTINE TAR3EL(A1,A2,AX,AY,AZ,CDESCR,IOSHP,MXNAT,NAT,IXYZ) IMPLICIT NONE C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3) 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 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target 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 96.01.26 (BTD) Replaced IX,IY,IZ by IXYZ C 00.11.02 (BTD) Modified to write ICOMP to target.out C end history C C Copyright (C) 1993,1996,2000 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. 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 JXMAX=LMX1 JXMIN=LMX2 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 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=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 IXYZ(NAT+JA,1)=IXYZ(JA,1)+JXMAX IXYZ(NAT+JA,2)=IXYZ(JA,2) IXYZ(NAT+JA,3)=IXYZ(JA,3) IXYZ(2*NAT+JA,1)=IXYZ(JA,1)+2*JXMAX IXYZ(2*NAT+JA,2)=IXYZ(JA,2) IXYZ(2*NAT+JA,3)=IXYZ(JA,3) 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 JX=1,NAT/3 WRITE(IOSHP,FMT=9031)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO DO JX=NAT/3+1,2*NAT/3 WRITE(IOSHP,FMT=9032)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO DO JX=2*NAT/3,NAT WRITE(IOSHP,FMT=9033)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO 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 ICOMP(x,y,z)') 9031 FORMAT(I7,3I4,' 1 1 1') 9032 FORMAT(I7,3I4,' 2 2 2') 9033 FORMAT(I7,3I4,' 3 3 3') END SUBROUTINE TARBLOCKS(A1,A2,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP, & CDESCR,MXNAT,NAT,IXYZ,ICOMP) IMPLICIT NONE C Arguments: INTEGER BLKSIZ,IPRINAX,NBLOCKS,IOSHP,MXNAT,NAT INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) INTEGER XYZB(3,100) REAL A1(3),A2(3) CHARACTER CDESCR*67 C Local variables: INTEGER JB,JD,JJX,JJY,JJZ,JOFF,JOFFX,JOFFY,JOFFZ,JX,JY,JZ C*********************************************************************** C Subroutine to construct a target out of cubic blocks C Given: C NBLOCKS =number of cubic blocks in target C BLKSIZ =(integer) ratio of block width/dipole spacing C (there will be BLKSIZ**3 dipoles per block) C XYZB(1-3,J) =x,y,z coordinates of block a, in units of block width C IPRINAX =0. to return A1=(1,0,0), A2=(0,1,0) C =1. to return A1,A2=principal axes with largest, 2nd largest C moment of inertia C IOSHP =device number for output file "target.out" C =-1 to suppress generation of file "target.out" C MXNAT =limit on largest allowed value of NAT=number of dipoles C C Returns: C NAT =number of dipoles in target C IXYZ(J,1-3) =x,y,z location of dipole J on lattice C ICOMP(J,1-3)=composition for dipole J; x,y,z directions C A1(1-3) =principal axis A1 in target frame (normalized) C A2(1-3) =principal axis A2 in target frame (normalized) C CDESCR =string describing target C C where A1 = principal axis with largest moment of inertia C A2 = principal axis with second-largest moment of inertia C C B.T. Draine, Princeton Univ. Observatory, 95.12.11 C History: C 96.01.02 (BTD) modified to compute eigenvalues and eigenvectors of C moment of inertia tensor, and to define vectors A1 and C A2 to be along eigenvectors with largest and second C largest eigenvalues, respectively. C 96.01.26 (BTD) Replaced IX,IY,IZ by IXYZ C Removed principal axis calculation to separate C subroutine PRINAXIS C 96.01.29 (BTD) Added variable IPRINAX to allow control over definition C of axes A1,A2 C 96.02.23 (BTD) Remove computation of principal axes and eigenvalues C as this is now carried out by routine PRINAXIS. C End history. C Copyright (C) 1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** c*** diagnostic write(0,*)'entered subroutine tarblock...' do jb=1,nblocks write(0,*)'jb=',jb,'xyz=',(xyzb(jx,jb),jx=1,3) enddo c*** C For the moment, let us assume that blocks do not overlap. JOFF=1-((BLKSIZ+1)/2) JD=0 DO JB=1,NBLOCKS JOFFX=JOFF+BLKSIZ*XYZB(1,JB) JOFFY=JOFF+BLKSIZ*XYZB(2,JB) JOFFZ=JOFF+BLKSIZ*XYZB(3,JB) DO JX=1,BLKSIZ JJX=JX+JOFFX DO JY=1,BLKSIZ JJY=JY+JOFFY DO JZ=1,BLKSIZ JJZ=JZ+JOFFZ JD=JD+1 IXYZ(JD,1)=JJX IXYZ(JD,2)=JJY IXYZ(JD,3)=JJZ C Homogeneous target: ICOMP(JD,1)=1 ICOMP(JD,2)=1 ICOMP(JD,3)=1 ENDDO ENDDO ENDDO ENDDO NAT=JD C IF(IPRINAX.NE.0)THEN C C Call PRINAXIS to (1) compute principal axes C (2) compute eigenvalues of moment of inertia tensor C c*** diagnostic write(0,*)'about to call prinaxis with' write(0,*)'mxnat=',mxnat write(0,*)'nat=',nat write(0,*)'icomp(3,3)=',icomp(3,3) write(0,*)'ixyz(3,3)=',ixyz(3,3) c*** CALL PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,A1,A2) c*** diagnostic write(0,*)'returned from prinaxis' c*** ELSE C C Do not compute principal axes: simply set C a1=(1,0,0) and a2=(0,1,0) in target frame. C A1(1)=1. A1(2)=0. A1(3)=0. A2(1)=0. A2(2)=1. A2(3)=0. ENDIF C C*********************************************************************** C Write target description into string CDESCR C Note: this will replace whatever CDESCR was read in from the C 'blocks.par' file. If you want to retain that description, C then comment out the following statement: C WRITE(CDESCR,FMT='(A,I7,A,I3,A)')'Target containing', C & NAT,' dipoles arranged in ',NBLOCKS,' cubic blocks' C*********************************************************************** C IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)NAT,A1,A2 DO JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3), & ICOMP(JX,1),ICOMP(JX,2),ICOMP(JX,3) ENDDO CLOSE(UNIT=IOSHP) ENDIF RETURN 9020 FORMAT(' >TARBLOCKS: cubic building blocks',/, &I7,' = NAT',/, &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 TARCEL(A1,A2,AX,AY,AZ,BX,BY,BZ,CDESCR,IOSHP,MXNAT,NAT, & IXYZ,ICOMP) IMPLICIT NONE C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3),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 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target 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 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ C 00.11.02 (BTD): remove XOFF,YOFF,ZOFF from target.out C cosmetic changes C end history C C Copyright (C) 1994,1996,2000 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. 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 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=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, & A1,A2 DO 40 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3), & 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,' = NAT,NIN',/, & 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 TARCYL(A1,A2,A,B,CDESCR,IOSHP,MXNAT,NAT,IXYZ) IMPLICIT NONE C Scalar arguments: CHARACTER CDESCR*67 INTEGER MXNAT,NAT REAL A,B C Array arguments: INTEGER*2 IXYZ(MXNAT,3) REAL A1(3),A2(3) C Local scalars: INTEGER IAXIS,IOSHP,JX,JY,JZ, & NB,NBD,NBMAX,NBMIN,NFAC,NLAY, & NX1,NX2,NY1,NY2,NZ1,NZ2 LOGICAL IN REAL ASPR,PI,REFF2,ROFF,X,XOFF,Y,YOFF,Z,ZOFF C*********************************************************************** C Purpose: to construct cylindrical target from "atoms", with cylinder C axis along x axis. C Input: 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms in target 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 ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,' 1 1 1') 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,IXYZ) IMPLICIT NONE C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3) 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 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target 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 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ C 00.11.02 (BTD): Write ICOMP to target.out C end history C C Copyright (C) 1993,1996,2000 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** 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 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. 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 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=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 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO 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 ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,' 1 1 1') END