SUBROUTINE TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,DX,NAT,IXYZ,ICOMP,IDVOUT, & MXN3,NAT3) IMPLICIT NONE C C Arguments: C CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6 INTEGER IDVOUT,IDVSHP,IOSHP,MXNAT,NAT,MXN3,NAT3 REAL A1(3),A2(3),DX(3),SHPAR(6) INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3) EXTERNAL ERRMSG,REASHP,TARCYL,TARELL,TARHEX,TARREC,TARTET C C Local variables: C INTEGER BLKSIZ,IPRINAX,J,J1,NBLOCKS INTEGER XYZB(3,100) 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, BLOCKS, C DW1996, TWOSPH, LYRSLB, NSPHER, PRISM3, SPHARM 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 DX(1-3)= d_x/d, d_y/d, d_z/d for lattice, where d=(d_x*d_y*d_z)**(1/3) C is the effective lattice spacing, and d_x,d_y,d_z are lattice C spacings in x,y,z directions 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 IXYZ(1-NAT,1-3)=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 95.12.11 (BTD): Modified to handle "BLOCKS" as target option, using C subroutine TARBLOCKS C 96.01.04 (BTD): Modified to handle "DW1996" as target option C 96.01.25 (BTD): Modified to handle "TWOSPH" as target option C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ C 96.01.29 (BTD): Added variable IPRINAX to control choice of axes A1,A2 C in option BLOCKS C 97.06.04 (BTD): Corrected bug in handling of option "UNICYL" -- was C not setting composition correctly after changing REAPAR C so that no longer had choice of cylinder axis. C Cylinder axis is always in x direction in target frame. C Thanks for Mike Wolff for calling attention to the C bug. C 97.11.01 (BTD): add DX to argument list to allow use of noncubic lattice C add DX to argument list of TARELL C 97.12.26 (BTD): add DX to argument list of TARREC, TARCYL, TARHEX, C TARTET, TAR2EL, TAR3EL, TARCEL, TAR2SP C 98.04.27 (BTD): restored missing comma in argument list of TAR2EL C 98.12.29 (BTD): added new option 'LYRSLB' for layered slab C 00.06.12 (BTD): modified to support option TARNSP for multisphere C target, including variable lattice spacing DX C all data type conversions now explicit C 00.07.18 (BTD): corrected bug in calls to tarblocks C 00.11.02 (BTD): added ICOMP to arg. list for TAR2EL,TAR3EL,TARCYL, C TARELL,TARHEX,TARREC,RCTNGL,TARTET C 02.02.12 (BTD): added option PRISM3 and call to routine TARPRSM C 03.05.21 (BTD): added option SPHARM = irregular target C generated using spherical harmonics C target generated by TARSPHARM C 03.05.22 (BTD): bug fix: define CFLSHP='spharm.par' when calling C TARSPHARM C end history C C Copyright (C) 1993,1994,1995,1996,1997,1998,2000,2002,2003 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C******************************************************************* C----------------------------------------------------------------------- C TARREC -> homogeneous, isotropic, rectangular target IF(CSHAPE.EQ.'RCTNGL')THEN CALL TARREC(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C ELLIPS -> homogeneous, isotropic, ellipsoidal target ELSEIF(CSHAPE.EQ.'ELLIPS')THEN CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C CYLNDR -> homogeneous, isotropic cylinder ELSEIF(CSHAPE.EQ.'CYLNDR')THEN CALL TARCYL(A1,A2,SHPAR(1),SHPAR(2),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C HEXGON -> homogeneous, isotropic hexagonal prism ELSEIF(CSHAPE.EQ.'HEXGON')THEN CALL TARHEX(A1,A2,SHPAR(1),SHPAR(2),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C TETRAH -> homogeneous, isotropic tetrahedron ELSEIF(CSHAPE.EQ.'TETRAH')THEN CALL TARTET(A1,A2,SHPAR(1),DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ, & ICOMP) C----------------------------------------------------------------------- C UNICYL -> uniaxial cylinder C Cylinder axis is in x direction in Target Frame. 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 ELSEIF(CSHAPE.EQ.'UNICYL')THEN CALL TARCYL(A1,A2,SHPAR(1),SHPAR(2),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C Set composition DO J=1,NAT ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=2 ENDDO C----------------------------------------------------------------------- C ANIELL -> anisotropic ellipsoid ELSEIF(CSHAPE.EQ.'ANIELL')THEN CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C Dielectric constants 1,2,3 for E along directions x,y,z in Target Frame DO J=1,NAT ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 ENDDO 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),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) 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),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) J1=NAT/2 DO J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 ENDDO DO J=J1+1,NAT ICOMP(J,1)=4 ICOMP(J,2)=5 ICOMP(J,3)=6 ENDDO 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 C 1st,2nd,3rd ellipsoid has ICOMP=1,2,3 (set in TAR3EL) ELSEIF(CSHAPE.EQ.'THRELL')THEN CALL TAR3EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C THRAEL -> three touching identical ellipsoids with anisotropic C functions. Need to override TAR3EL. 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),DX,CDESCR,IOSHP, & MXNAT,NAT,IXYZ,ICOMP) J1=NAT/3 DO J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 ENDDO DO J=J1+1,2*J1 ICOMP(J,1)=4 ICOMP(J,2)=5 ICOMP(J,3)=6 ENDDO DO J=2*J1+1,NAT ICOMP(J,1)=7 ICOMP(J,2)=8 ICOMP(J,3)=9 ENDDO 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),DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C TAR2SP -> two touching spheroids ELSEIF(CSHAPE.EQ.'TWOSPH')THEN CALL TAR2SP(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5), & SHPAR(6),DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C BLOCKS -> target is constructed from identical cubic blocks C locations of blocks, and size of blocks (lattice units) C is specified in file 'blocks.par' C file structure: C target description C IPRINAX = 0 (a1=100,a2=010) or 1 (a1,a2=principal axes) C N = number of blocks C SIZE = width of one block, in lattice units C X Y Z = location of block 1 (in units of block width) C ... C X Y Z = " " " N " " " " " C ELSEIF(CSHAPE.EQ.'BLOCKS')THEN OPEN(UNIT=IDVSHP,FILE='blocks.par',STATUS='OLD', & FORM='FORMATTED') READ(IDVSHP,FMT='(A67)')CDESCR READ(IDVSHP,*)IPRINAX READ(IDVSHP,*)NBLOCKS READ(IDVSHP,*)BLKSIZ DO J=1,NBLOCKS READ(IDVSHP,*)XYZB(1,J),XYZB(2,J),XYZB(3,J) ENDDO CLOSE(IDVSHP) CALL TARBLOCKS(A1,A2,DX,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP, & CDESCR,MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C DW1996 -> the 13-cube target studied by Draine & Weingartner 1996 C SHPAR(1)=width per block in lattice units ELSEIF(CSHAPE.EQ.'DW1996')THEN IPRINAX=1 NBLOCKS=13 BLKSIZ=NINT(SHPAR(1)) XYZB(1,1)=0 XYZB(2,1)=1 XYZB(3,1)=0 XYZB(1,2)=1 XYZB(2,2)=1 XYZB(3,2)=0 XYZB(1,3)=0 XYZB(2,3)=2 XYZB(3,3)=0 XYZB(1,4)=1 XYZB(2,4)=2 XYZB(3,4)=0 XYZB(1,5)=0 XYZB(2,5)=1 XYZB(3,5)=1 XYZB(1,6)=1 XYZB(2,6)=1 XYZB(3,6)=1 XYZB(1,7)=0 XYZB(2,7)=2 XYZB(3,7)=1 XYZB(1,8)=1 XYZB(2,8)=2 XYZB(3,8)=1 XYZB(1,9)=2 XYZB(2,9)=1 XYZB(3,9)=0 XYZB(1,10)=2 XYZB(2,10)=2 XYZB(3,10)=0 XYZB(1,11)=0 XYZB(2,11)=0 XYZB(3,11)=1 XYZB(1,12)=0 XYZB(2,12)=0 XYZB(3,12)=2 XYZB(1,13)=0 XYZB(2,13)=1 XYZB(3,13)=2 CALL TARBLOCKS(A1,A2,DX,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP, & CDESCR,MXNAT,NAT,IXYZ,ICOMP) CDESCR='13-cube target used by Draine & Weingartner 1996' C----------------------------------------------------------------------- C TARSLB -> Layered slab: ELSEIF(CSHAPE.EQ.'LYRSLB')THEN CALL TARSLB(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5), & SHPAR(6),DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C C slab has dimensions a,b,c in x,y,z dimensions C current version allows for up to 4 layers C slab is layered in x direction C 0 < x < x1 is composition 1 C x1 < x < x2 is composition 2 C x2 < x < x3 is composition 3 C x3 < x < a is composition 4 C for fewer than 4 layers, let x2 or x3=0 C SHPAR(1-3)= a/d , b/d , c/d where d=(dx*dy*dz)**{1/3} C SHPAR(4) = x1/a C SHPAR(5) = x2/a C SHPAR(6) = x3/a C----------------------------------------------------------------------- C NSPHER -> Multisphere target: ELSEIF(CSHAPE.EQ.'NSPHER')THEN CALL TARNSP(A1,A2,SHPAR(1),SHPAR(2),DX,CFLSHP,CDESCR, & IDVSHP,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C PRISM3 -> triangular prism: ELSEIF(CSHAPE.EQ.'PRISM3')THEN CALL TARPRSM(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),DX, & CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C SPHARM -> spherical harmonic target C SHPAR(1)=a_eff/d C SHPAR(2)=0. to use A1=(1,0,0),A2=(0,1,0) in TF C 1. to use A1,A2 = principal axes with largest,2nd largest C moment of inertia C a_lm values are input from file 'spharm.par' ELSEIF(CSHAPE.EQ.'SPHARM')THEN CFLSHP='spharm.par' CALL TARSPHARM(SHPAR(1),SHPAR(2),A1,A2,CFLSHP, & CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C----------------------------------------------------------------------- C add any other specific shape routines here C----------------------------------------------------------------------- C FRMFIL -> read list of occupied sites and compositions from file ELSEIF(CSHAPE.EQ.'FRMFIL')THEN CALL REASHP(CFLSHP,IDVOUT,IDVSHP,A1,A2,DX,CDESCR,MXNAT,NAT, & IXYZ,ICOMP,MXN3,NAT3) C----------------------------------------------------------------------- C Reach here only in event of error: ELSE CALL ERRMSG('FATAL','TARGET', & ' Uknown shape requested: check ddscat.par') 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 Before returning, compute: NAT3=3*NAT RETURN 9010 FORMAT(' >TARGET:',A,/, & 7X,I7,' = NAT0 = number of dipoles in target') END SUBROUTINE TARHEX(A1,A2,A,B,DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) IMPLICIT NONE C Scalar arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) REAL A,B C Array arguments: REAL A1(3),A2(3),DX(3) C Local scalars: CHARACTER CMSGNM*70 INTEGER JA,JX,JY,JZ,NA,NB,NC,NFAC,NLAY,NMX,NMY,NMZ LOGICAL IN REAL ACM,ASPR,BCM,BEFF,CCM,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 DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z C directions, and d=(dx*dy*dz)**(1/3)=effective lattice C spacing 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 (prism axis) C A2(1-3)=unit vector (0,1,0) defining target axis 2 (normal to C one of the rectangular faces. 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=string describing target (up to 67 chars.) C C B.T.Draine, Princeton Univ. Obs., 87.04.03 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 95.07.13 (BTD): Simplify: restrict to case where prism axis is in C x direction and y axis is normal to a rectangular C face. C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ C 96.02.23 (BTD): Delete PI (unused) C 97.04.28 (BTD): Initialize NAT=0 (failure to do so reported by C Michael Wolff). C 97.12.26 (BTD): Added DX(3) to argument list to support nonuniform C lattice spacing. C *** Note: additional code modifications required C to support noncubic lattice! C 98.03.06 (BTD): If called with noncubic lattice (not yet supported) C call ERRMSG and halt with fatal error. C 98.03.07 (BTD): Modify to write DX to file "target.out" for C compatibility with REASHP C 00.11.02 (BTD): Added ICOMP to argument list C set ICOMP=1 for occupied sites C write ICOMP to target.out C end history C C Copyright (C) 1993,1995,1996,1997,1998,2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Current version of TARHEX is restricted to cubic lattices IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN CALL ERRMSG('FATAL','TARHEX', & ' tarhex does not support noncubic lattice') ENDIF DO JA=1,3 A1(JA)=0. A2(JA)=0. ENDDO A1(1)=1. A2(2)=1. 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 NMX=NA NMY=NC NMZ=NB NAT=0 DO JZ=1,NMZ Z=REAL(JZ) DO JY=1,NMY Y=REAL(JY) DO JX=1,NMX X=REAL(JX) CALL TESTHEX(ACM,BCM,CCM,A,B,X,Y,Z,IN) IF(IN)THEN NAT=NAT+1 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDIF ENDDO ENDDO ENDDO C Set composition DO JZ=1,3 DO JX=1,NAT ICOMP(JX,JZ)=1 ENDDO ENDDO 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') WRITE(IOSHP,9020)A,B,NAT,A1,A2,DX 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(' Hexagonal Prism: A=',F7.4,' B=',F7.4,/, &I7,'= NAT. Prism axis along x, one rect. face. normal to y',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/, &' JA IX IY IZ ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,3I2) END SUBROUTINE TESTHEX(ACM,BCM,CCM,A,B,X,Y,Z,IN) C Arguments: 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 AX=X U=Y-CCM V=Z-BCM 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).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)+0.8660254*ABS(V) IF(Q.GT.0.4330127*B)RETURN IN=.TRUE. RETURN END SUBROUTINE TARNSP(A1,A2,DIAMX,PRINAX,DX,CFLSHP,CDESCR, & IDVSHP,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C** Arguments: IMPLICIT NONE CHARACTER CDESCR*67,CFLSHP*13 INTEGER IDVSHP,IOSHP,MXNAT,NAT INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) REAL DIAMX,PRINAX REAL A1(3),A2(3),DX(3) C** Local variables: INTEGER NSPHMX PARAMETER(NSPHMX=100) CHARACTER CMSGNM*70 LOGICAL OCC INTEGER JA,JX,JY,JZ,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2,NSPH REAL R2,SCALE,X,XMAX,XMIN,Y,YMAX,YMIN,Z,ZMAX,ZMIN REAL AS(NSPHMX),XS(NSPHMX),YS(NSPHMX),ZS(NSPHMX) C*********************************************************************** C Routine to construct multisphere target C C Input: C DIAMX =max extent in X direction/d C PRINAX=0 to set A1=(1,0,0), A2=(0,1,0) C =1 to use principal axes for vectors A1 and A2 C DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z C directions, and d=(dx*dy*dz)**(1/3)=effective lattice C spacing C CFLSHP= name of file containing locations and radii of spheres 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 C and, from input file CFLSHP: C C NSPH = number of spheres C one descriptive line which will be ignored (may be blank) C XS(1) YS(1) ZS(1) AS(1) C XS(2) YS(2) ZS(2) AS(2) C ... C XS(NSPH) YS(NSPH) ZS(NSPH) AS(NSPH) C C where XS(J),YS(J),ZS(J) = x,y,z coordinates in target frame, in C arbitrary units, of center of sphere J C AS(J) = radius of sphere J (same arbitrary units) C C units used for XS,YS,ZS, and AS are really arbitrary: actual size of C overall target, in units of lattice spacing d, is controlled by C parameter SHPAR(1) in ddscat.par C 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)=composition identifier (currently 1 1 1) C CDESCR=description of target (up to 67 characters) C C B.T.Draine, Princeton Univ. Obs., 2000.06.12 C C History: C 00.06.12 (BTD) adapted to version 6.0 (introduce DX) C 00.11.02 (BTD) write ICOMP to target.out C 01.01.13 (BTD) corrected comments C 01.04.21 (BTD) corrected comments C 01.04.21 (BTD) corrected typo: changed C YMAX=YS(1)+YS(1) C to YMAX=YS(1)+AS(1) C Thanks to Ivan O. Sosa Perez (UNAM) for discovering C this error. C 03.06.06 (BTD) corrected error in setting LMX1,...,LMZ2 C end history C C Copyright (C) 2000,2001,2003 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Read parameters of spheres: OPEN(UNIT=IDVSHP,FILE=CFLSHP) READ(IDVSHP,*)NSPH READ(IDVSHP,*) DO JX=1,NSPH READ(IDVSHP,*)XS(JX),YS(JX),ZS(JX),AS(JX) ENDDO CLOSE(IDVSHP) C Dipoles are located at sites C (x,y,z)=(I,J,K), I,J,K=integers C Determine max extent in X direction: XMIN=XS(1)-AS(1) XMAX=XS(1)+AS(1) YMIN=YS(1)-AS(1) YMAX=YS(1)+AS(1) ZMIN=ZS(1)-AS(1) ZMAX=ZS(1)+AS(1) IF(NSPH.GT.1)THEN DO JX=2,NSPH IF(XS(JX)-AS(JX).LT.XMIN)XMIN=XS(JX)-AS(JX) IF(XS(JX)+AS(JX).GT.XMAX)XMAX=XS(JX)+AS(JX) IF(YS(JX)-AS(JX).LT.YMIN)YMIN=YS(JX)-AS(JX) IF(YS(JX)+AS(JX).GT.YMAX)YMAX=YS(JX)+AS(JX) IF(ZS(JX)-AS(JX).LT.ZMIN)ZMIN=ZS(JX)-AS(JX) IF(ZS(JX)+AS(JX).GT.ZMAX)ZMAX=ZS(JX)+AS(JX) ENDDO ENDIF SCALE=DIAMX/(XMAX-XMIN) C Now determine min,max values of I,J,K: LMX1=NINT(SCALE*XMIN/DX(1)-0.01) LMX2=NINT(SCALE*XMAX/DX(1)+0.01) LMY1=NINT(SCALE*YMIN/DX(2)-0.01) LMY2=NINT(SCALE*YMAX/DX(2)+0.01) LMZ1=NINT(SCALE*ZMIN/DX(3)-0.01) LMZ2=NINT(SCALE*ZMAX/DX(3)+0.01) C Convert AS to squared radius: DO JX=1,NSPH AS(JX)=(SCALE*AS(JX))**2 XS(JX)=SCALE*XS(JX) YS(JX)=SCALE*YS(JX) ZS(JX)=SCALE*ZS(JX) ENDDO C Determine list of occupied sites NAT=0 DO JZ=LMZ1,LMZ2 Z=REAL(JZ)*DX(3) DO JY=LMY1,LMY2 Y=REAL(JY)*DX(2) DO JX=LMX1,LMX2 X=REAL(JX)*DX(1) OCC=.FALSE. DO JA=1,NSPH R2=(X-XS(JA))**2+(Y-YS(JA))**2+(Z-ZS(JA))**2 IF(R2.LT.AS(JA))OCC=.TRUE. ENDDO IF(OCC)THEN C Site is occupied: NAT=NAT+1 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDIF ENDDO ENDDO ENDDO IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARNSP',' NAT.GT.MXNAT ') ENDIF C Homogeneous target: DO JA=1,NAT DO JX=1,3 ICOMP(JA,JX)=1 ENDDO ENDDO 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 C IF(PRINAX.LE.0.)THEN DO JX=1,3 A1(JX)=0. A2(JX)=0. ENDDO A1(1)=1. A2(2)=1. ELSE CALL PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,DX,A1,A2) ENDIF C*********************************************************************** C Write target description into string CDESCR WRITE(CDESCR,FMT='(A,I7,A)')' Multisphere cluster containing', & NAT,' dipoles' C*********************************************************************** CMSGNM=CDESCR CALL WRIMSG('TARNSP',CMSGNM) IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)NSPH,DIAMX,NAT,A1,A2,DX 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(' >TARNSP multisphere target composed of ',I3,' spheres, ', & 'DIAMX=',F8.4,/, & I7,' = NAT dipoles',/, & 3F9.4,' = A_1 vector',/, & 3F9.4,' = A_2 vector',/, & 3F9.4,' = lattice spacings (d_x,d_y,d_z)/d',/, &' JA IX IY IZ ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,3I2) END SUBROUTINE TARPRSM(A1,A2,A,BA,CA,LA,DX,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,BA,CA,LA C Array arguments: REAL A1(3),A2(3),DX(3) C Local scalars: CHARACTER CMSGNM*70 INTEGER JA,JX,JY,JZ,NFAC,NLAY,NX,NY,NZ LOGICAL IN REAL & B,BETA,C,COTBETA,COTGAMMA,GAMMA,L,X, & XMAX,XMIN,Y,YMAX,YMIN,Z,ZMAX,ZMIN C*********************************************************************** C Routine to construct triangular prism from "atoms". C triangle side lengths = a,b,c ; prism length = L C prism axis is assumed to be in x direction C normal to prism face of width a : (0,1,0) C normal to prism face of width b : (0,-cos(gamma),sin(gamma)) C normal to prism face of width c : (0,-cos(beta),-sin(beta)) C angles alpha,beta,gamma are opposite sides a,b,c C first layer of target array is at x=0 (target boundary at x=-0.5) C layer along side a is at int(b*sin(gamma)) C boundary is at ymax=int(b*sin(gamma))+0.5 C Input: C A = a/d C BA = b/a C CA = c/a C LA = L/a C DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z C directions, and d=(dx*dy*dz)**(1/3)=effective lattice C spacing 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 (prism axis) C A2(1-3)=unit vector (0,1,0) defining target axis 2 (normal to C rectangular faces with sides a,L 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)=dielectric function identifier for dipole C locations and 3 directions in space C CDESCR=string describing target (up to 67 chars.) C C B.T.Draine, Princeton Univ. Obs., 2002.02.12 C C History: C 02.02.12 (BTD): adapted from tarprsm.f for ver5a10 C 04.02.25 (BTD): corrected typographical error found by In-Ho Lee C end history C C Copyright (C) 2002,2004 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C----------------------------------------------------------------------- C Current version of TARPRSM is restricted to cubic lattices IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN CALL ERRMSG('FATAL','TARPRSM', & ' tarprsm does not support noncubic lattice') ENDIF DO JA=1,3 A1(JA)=0. A2(JA)=0. ENDDO A1(1)=1. A2(2)=1. C A = a/d C B = b/d C C = c/d C L = L/d B=BA*A C=CA*A L=LA*A c sanity check: IF((A.GE.B+C).OR.(B.GE.A+C).OR.(C.GE.A+B)) & CALL ERRMSG('FATAL','TARPRSM', &'prism sides a,b,c do not satisfy triangle inequality:check input' &) BETA=ACOS((A*A+C*C-B*B)/(2.*A*C)) GAMMA=ACOS((A*A+B*B-C*C)/(2.*A*B)) COTBETA=COS(BETA)/SIN(BETA) COTGAMMA=COS(GAMMA)/SIN(GAMMA) c ideal prism extent: c x = xmin to x=xmax c triangle vertices at (0,ymax,zmin),(0,ymax,zmax),(0,ymin,za) c where xmin=-0.5 c xmax=xmin+L c ymax=int[b*sin(gamma)]+0.5 c ymin=ymax-b*sin(gamma) c zmax=a/2 c zmin=-a/2 c za=-a/2+c*cos(beta) XMIN=-0.5 XMAX=XMIN+L YMAX=INT(B*SIN(GAMMA))+0.5 YMIN=YMAX-B*SIN(GAMMA) ZMAX=0.5*A ZMIN=-ZMAX C Now determine limits for testing x,y,z values C Along axis (x), run from 0 to NX=INT(XMAX) C Along y direction, run from 0 to NY=INT(YMAX) C Along z direction, run from -NZ to NZ=INT(ZMAX) NX=INT(XMAX) NY=INT(YMAX) NZ=INT(ZMAX) NAT=0 DO JZ=-NZ,NZ Z=REAL(JZ) DO JY=0,NY Y=REAL(JY) IF((Z.LE.ZMAX+(Y-YMAX)*COTGAMMA) & .AND. & (Z.GE.ZMIN-(Y-YMAX)*COTBETA))THEN DO JX=0,NX NAT=NAT+1 IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDDO ENDIF ENDDO ENDDO NLAY=NX+1 NFAC=NAT/NLAY DO JX=1,3 DO JA=1,NAT ICOMP(JA,JX)=1 ENDDO ENDDO 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,') = prism axis direction' CALL WRIMSG('TARPRSM',CMSGNM) WRITE(CMSGNM,FMT='(A,3F7.4,A)') & ' A2 = (',A2,') = normal to face of width a' CALL WRIMSG('TARPRSM',CMSGNM) WRITE(CMSGNM,FMT='(A,I7,A,I4,A,I4)') & ' NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY CALL WRIMSG('TARPRSM',CMSGNM) WRITE(CDESCR,FMT='(A,F7.4,A,F8.4,A,F7.4,A,F8.4)') & ' tri.prism: a/d=',A, & ' b/a=',BA,' c/a=',CA,' L/a=',LA 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,A1,A2 C*** 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( & ' Triangular Prism: a/d=',F7.4,' b/a=',F7.4,' c/a=',F7.4, & ' L/a=',F7.4,/, & I7,'= NAT. Prism axis along x, rect. face a normal to y',/, & 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,DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ, & ICOMP) IMPLICIT NONE C Arguments: REAL XV,YV,ZV INTEGER IOSHP,MXNAT,NAT CHARACTER CDESCR*67 INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) REAL A1(3),A2(3),DX(3) C Local variables: INTEGER JA,JX,JY,JZ,NX,NY,NZ CHARACTER CMSGNM*70 C External subroutines: EXTERNAL ERRMSG C Intrinsic functions: INTEGER IFIX INTRINSIC IFIX 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 DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z C directions, and d=(dx*dy*dz)**(1/3)=effective lattice C spacing 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 in Target Frame C A2(1-3)=unit vector (0,1,0) 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 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 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ C 97.12.26 (BTD): Added DX(1-3) to argument list to support nonuniform C lattice spacing. C 98.02.11 (BTD): Modify to generate target on noncubic lattice. C 00.11.02 (BTD): Add ICOMP to argument list, set ICOMP=1, C write ICOMP to target.out C end history C C Copyright (C) 1993,1995,1996,1997,1998,2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C*********************************************************************** NX=IFIX(XV/DX(1)+0.5) NY=IFIX(YV/DX(2)+0.5) NZ=IFIX(ZV/DX(3)+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 Specify target axes A1 and A2 C Convention: A1 will be (1,0,0) in target frame C A2 will be (0,1,0) in target frame DO JA=1,3 A1(JA)=0. A2(JA)=0. ENDDO A1(1)=1. A2(2)=1. C C Now populate lattice: NAT=0 DO JZ=1,NZ DO JY=1,NY DO JX=1,NX NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDDO ENDDO ENDDO IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT ') ENDIF C homogeneous, isotropic composition: DO JX=1,3 DO JA=1,NAT ICOMP(JA,JX)=1 ENDDO ENDDO 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,DX DO JA=1,NAT WRITE(IOSHP,FMT=9030)JA,IXYZ(JA,1),IXYZ(JA,2),IXYZ(JA,3), & ICOMP(JA,1),ICOMP(JA,2),ICOMP(JA,3) ENDDO 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',/, &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/, &' JA IX IY IZ ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,3I2) END SUBROUTINE TARSLB(A1,A2,XV,YV,ZV,F1,F2,F3, & DX,CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) IMPLICIT NONE C*********************************************************************** C Routine to construct rectangular prism from "atoms" C with layered structure C slab has dimensions a,b,c in x,y,z dimensions C current version allows for up to 4 layers C slab is layered in x direction C 0 < x < f1*a is composition 1 C f1*a < x < f2*a is composition 2 C f2*a < x < f3*a is composition 3 C f3*a < x < a is composition 4 C C for fewer than 4 layers, let x2=x3=0 (for two layers) C or x3=0 (for three layers) C C Input: C XV=(x-length)/d (d=lattice spacing) C YV=(y-length)/d C ZV=(z-length)/d C F1,F2,F3 = interface locators (see above) C DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z C directions, and d=(dx*dy*dz)**(1/3)=effective lattice C spacing 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 in Target Frame C A2(1-3)=unit vector (0,1,0) 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)=composition C C B.T.Draine, Princeton Univ. Obs. C History: C 98.12.29 (BTD): Created using TARREC as starting material. C 00.11.02 (BTD): Write ICOMP to target.out C end history C C Copyright (C) 1998,2000 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments: REAL F1,F2,F3,XV,YV,ZV INTEGER IOSHP,MXNAT,NAT CHARACTER CDESCR*67 INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) REAL A1(3),A2(3),DX(3) C C Local variables: INTEGER JA,JX,JY,JZ,NCOMP,NX,NX1,NX2,NX3,NY,NZ CHARACTER CMSGNM*70 C C External subroutines: EXTERNAL ERRMSG C C Intrinsic functions: INTEGER IFIX INTRINSIC IFIX C*********************************************************************** c*** diagnostic write(0,*)'entered tarslb' c*** IF(F2.LE.0.)THEN NCOMP=2 ELSE IF(F3.LE.0.)THEN NCOMP=3 ELSE NCOMP=4 ENDIF ENDIF NX=IFIX(XV/DX(1)+0.5) NY=IFIX(YV/DX(2)+0.5) NZ=IFIX(ZV/DX(3)+0.5) WRITE(CMSGNM,FMT='(A,3I4)') & ' Layered slab; NX,NY,NZ=',NX,NY,NZ WRITE(CDESCR,FMT='(A,3I4)') & ' Layered slab; NX,NY,NZ=',NX,NY,NZ C NX1=IFIX(F1*NX+0.5) NX2=NX NX3=NX IF(NCOMP.GE.3)NX2=IFIX(F2*NX+0.5) IF(NCOMP.EQ.4)NX3=IFIX(F3*NX+0.5) C C Specify target axes A1 and A2 C Convention: A1 will be (1,0,0) in target frame C A2 will be (0,1,0) in target frame C 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 JZ=1,NZ DO JY=1,NY DO JX=1,NX1 NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARSLB',' NAT.GT.MXNAT') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ICOMP(NAT,1)=1 ICOMP(NAT,2)=1 ICOMP(NAT,3)=1 ENDDO DO JX=NX1+1,NX2 NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARSLB',' NAT.GT.MXNAT') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ICOMP(NAT,1)=2 ICOMP(NAT,2)=2 ICOMP(NAT,3)=2 ENDDO ENDDO ENDDO IF(NCOMP.GE.3)THEN DO JZ=1,NZ DO JY=1,NY DO JX=NX2+1,NX3 NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARSLB',' NAT.GT.MXNAT') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ICOMP(NAT,1)=3 ICOMP(NAT,2)=3 ICOMP(NAT,3)=3 ENDDO ENDDO ENDDO ENDIF IF(NCOMP.EQ.4)THEN DO JZ=1,NZ DO JY=1,NY DO JX=NX3+1,NX NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARSLB',' NAT.GT.MXNAT') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ICOMP(NAT,1)=4 ICOMP(NAT,2)=4 ICOMP(NAT,3)=4 ENDDO ENDDO ENDDO ENDIF 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,DX DO JA=1,NAT WRITE(IOSHP,FMT=9030)JA,IXYZ(JA,1),IXYZ(JA,2),IXYZ(JA,3), & ICOMP(JA,1),ICOMP(JA,2),ICOMP(JA,3) ENDDO 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',/, &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/, &' JA IX IY IZ ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,3I2) END SUBROUTINE TARSPHARM(AEFF,PRINAX,A1,A2,CFLSHP, & CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) IMPLICIT NONE C Parameters: C note: this must agree with corresponding parameter in FUNCTION RSPHARM C below INTEGER NMAX PARAMETER(NMAX=10) C Arguments: CHARACTER CDESCR*67,CFLSHP*(*) INTEGER IOSHP,MXNAT,NAT INTEGER*2 ICOMP(MXNAT,3),IXYZ(MXNAT,3) REAL AEFF,PRINAX REAL A1(3),A2(3) REAL AL0(0:NMAX) COMPLEX ALM(1:NMAX,1:NMAX) C Local variables INTEGER NRES PARAMETER(NRES=20) CHARACTER CMSGNM*70 INTEGER & JA,JPHI,JTH,JX,JY,JZ, & L,LMAX,LMX1,LMX2,LMY1,LMY2,LMZ1,LMZ2, & M,NFAC,NLAY,NMX,NMY,NMZ,NX,NY,NZ REAL & COSPH,COSTH,PHI,PHIARG,PI,R,R2,RAD,RYZ,RYZ2,SINPH,SINTH,SUM, & THETA,TWOPI,X,XMAX,XMIN,XX,Y,YMAX,YMIN,YY,Z,ZMAX,ZMIN,ZZ REAL & PHIS(0:NRES*NMAX,3), & RADS(0:NRES*NMAX,3), & THETS(0:NRES*NMAX,3) COMPLEX CXI C External functions REAL RSPHARM EXTERNAL RSPHARM C*********************************************************************** C Routine to construct irregular grain using spherical harmonics C to define surface. C C Input: C AEFF=effective radius/d (d=lattice spacing) C CDESCR=name of file containing spherical harmonic amplitudes 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms of target C CDESCR=string describing target (up to 67 chars.) C C B.T.Draine, Princeton Univ. Obs., 02.11.12 C C History: C 02.11.12 (BTD): Start to write this routine C 02.11.13 (BTD): further modification C 03.05.21 (BTD): bug fix (set NAT=0) C bug fix to suppress cross section info when IOSHP < 1 C end history C C Copyright (C) 2002 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C----------------------------------------------------------------------- CXI=(0.,1.) c*** diagnostic c write(0,*)'entered tarspharm with aeff=',aeff c*** c*** read parameters from file PI=4.*ATAN(1.D0) TWOPI=2.*PI OPEN(UNIT=17,FILE=CFLSHP) READ(17,*)LMAX IF(LMAX.LE.1)THEN WRITE(0,*)'FATAL ERROR: LMAX=0 does not make sense for ', & 'TARSPHARM -- why sphere?' STOP ENDIF c on input: c c to facilitate playing with amplitudes and phases of the a_lm: c we input a_lm via |a_lm|^2 and phase information as follows: c for l = 0 : AL0 = |a00|^2 c for l > 0, c m = 0 : AL0 = +/- |a_l0|^2 where al0 = +/- |a_l0| c m > 0 : c Re(ALM) = |a_lm|^2 c Im(ALM) = zeta/2*pi where alm=|a_lm|exp(i*zeta) READ(17,*)AL0(0) DO L=1,LMAX READ(17,*)AL0(L),(ALM(L,M),M=1,L) ENDDO CLOSE(17) c now convert to standard representation of a_lm as complex numbers: AL0(0)=SQRT(AL0(0)) DO L=1,LMAX AL0(L)=AL0(L)/SQRT(ABS(AL0(L))) DO M=1,L XX=REAL(ALM(L,M)) YY=AIMAG(ALM(L,M)) ALM(L,M)=SQRT(XX)*EXP(CXI*TWOPI*YY) ENDDO ENDDO c if necessary, renormalize the a_lm : SUM=AL0(0)**2 DO L=1,LMAX SUM=SUM+AL0(L)**2 DO M=1,L SUM=SUM+2.*ABS(ALM(L,M))**2 ENDDO ENDDO SUM=SQRT(SUM) IF(ABS(SUM-1.).GT.1.D-3)THEN WRITE(0,*)'SUM=',SUM,' : need to renormalize the a_lm' AL0(0)=AL0(0)/SUM DO L=1,LMAX AL0(L)=AL0(L)/SUM DO M=1,L ALM(L,M)=ALM(L,M)/SUM ENDDO ENDDO ENDIF DO JA=1,3 A1(JA)=0. A2(JA)=0. ENDDO A1(1)=1. A2(2)=1. c if IOSHP > 0 , calculate three cross-sectional slices c polar slice with phi=0 IF(IOSHP.GT.0)THEN PHI=0. DO JTH=0,NRES*LMAX THETA=TWOPI*REAL(JTH)/(REAL(NRES*LMAX)) COSTH=COS(THETA) THETS(JTH,1)=THETA PHIS(JTH,1)=PHI PHIARG=PHI IF(THETA.GT.PI)PHIARG=PHI+PI RADS(JTH,1)=RSPHARM(COSTH,PHIARG,LMAX,AL0,ALM) ENDDO c polar slice with phi=pi/2 PHI=PI/2.D0 DO JTH=0,NRES*LMAX THETA=TWOPI*REAL(JTH)/(REAL(NRES*LMAX)) COSTH=COS(THETA) THETS(JTH,2)=THETA PHIS(JTH,2)=PHI PHIARG=PHI IF(THETA.GT.PI)PHIARG=PHI+PI RADS(JTH,2)=RSPHARM(COSTH,PHIARG,LMAX,AL0,ALM) ENDDO c equatorial slice THETA=PI/2.D0 COSTH=0. DO JPHI=0,NRES*LMAX PHI=TWOPI*REAL(JPHI)/(REAL(NRES*LMAX)) THETS(JPHI,3)=THETA PHIS(JPHI,3)=PHI RADS(JPHI,3)=RSPHARM(COSTH,PHI,LMAX,AL0,ALM) ENDDO OPEN(UNIT=IOSHP,FILE='spharm.out') WRITE(IOSHP,7600)LMAX DO JTH=0,NRES*LMAX WRITE(IOSHP,7700) & (THETS(JTH,L),PHIS(JTH,L),RADS(JTH,L),L=1,3) ENDDO CLOSE(UNIT=IOSHP) ENDIF c sample surface to determine total extent of target in x,y,z directions XMIN=0. XMAX=0. YMIN=0. YMAX=0. ZMIN=0. ZMAX=0. DO JTH=1,NRES*LMAX THETA=PI*REAL(JTH-0.5)/REAL(NRES*LMAX) COSTH=COS(THETA) SINTH=SIN(THETA) DO JPHI=1,NRES*LMAX PHI=TWOPI*REAL(JPHI)/REAL(NRES*LMAX) RAD=AEFF*RSPHARM(COSTH,PHI,LMAX,AL0,ALM) XX=COSTH*RAD YY=SINTH*COS(PHI)*RAD ZZ=SINTH*SIN(PHI)*RAD IF(XX.LT.XMIN)XMIN=XX IF(XX.GT.XMAX)XMAX=XX IF(YY.LT.YMIN)YMIN=YY IF(YY.GT.YMAX)YMAX=YY IF(ZZ.LT.ZMIN)ZMIN=ZZ IF(ZZ.GT.ZMAX)ZMAX=ZZ ENDDO ENDDO LMX1=-INT(-XMIN) LMY1=-INT(-YMIN) LMZ1=-INT(-ZMIN) LMX2=INT(XMAX) LMY2=INT(YMAX) LMZ2=INT(ZMAX) c construct list of occupied sites NAT=0 DO JZ=LMZ1,LMZ2 Z=REAL(JZ) DO JY=LMY1,LMY2 Y=REAL(JY) RYZ2=Y**2+Z**2 RYZ=SQRT(RYZ2) COSPH=Y/RYZ SINPH=Z/RYZ DO JX=LMX1,LMX2 X=REAL(JX) R=SQRT(RYZ2+X**2) COSTH=X/R IF(COSPH.GE.0.)THEN PHI=ASIN(SINPH) ELSEIF(COSPH.LT.0)THEN PHI=PI-ASIN(SINPH) ENDIF RAD=AEFF*RSPHARM(COSTH,PHI,LMAX,AL0,ALM) IF(RAD.GE.R)THEN C Site is occupied: NAT=NAT+1 IF(NAT.GT.MXNAT)THEN c*** c write(0,*)'nat=',nat c write(0,*)'mxnat=',mxnat c*** CALL ERRMSG('FATAL','TARSPHARM',' NAT.GT.MXNAT ') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ ENDIF ENDDO ENDDO ENDDO C Homogeneous target: DO JA=1,NAT DO JX=1,3 ICOMP(JA,JX)=1 ENDDO ENDDO 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 C of inertia IF(PRINAX.LE.0.)THEN DO JX=1,3 A1(JX)=0. A2(JX)=0. ENDDO A1(1)=1. A2(2)=1. ELSE c*** diagnostic write(0,*)'about to call prinaxis...' CALL PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,A1,A2) write(0,*)'...returned from prinaxis' ENDIF C----------------------------------------------------------------------- C Write target description into string CDESCR C Note: string CDESCR will be printed by subr. TARGET WRITE(CDESCR,FMT='(A,I7,A)') & ' Spheroidal harmonic target of NAT=',NAT,' dipoles' C----------------------------------------------------------------------- C Here print any additional information which is desired by using C subroutine WRIMSG WRITE(CMSGNM,FMT='(A,I7,A)') & ' Target generated using spherical harmonics: NAT=',NAT, & ' dipoles' CALL WRIMSG('TARSPHARM',CMSGNM) IF(IOSHP.GT.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') C*** 3 lines of general description allowed: WRITE(IOSHP,9020)AEFF,NAT,A1,A2 C*** DO JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) ENDDO CLOSE(UNIT=IOSHP) ENDIF RETURN 7600 FORMAT( & 'three cross sectional slices for spherical harmonic grain',/, & I3,' = LMAX (highest order for Y_lm)',/, & '------ slice 1 ------ ------ slice 2 ------ ', & '------ slice 3 ------',/, & ' theta phi r/a_eff theta phi r/a_eff ', & 'theta phi r/a_eff') 7700 FORMAT(F6.4,F7.4,F8.4,2(2F7.4,F8.4)) 9020 FORMAT(' Spherical harmonic grain: AEFF=',F7.4,/, &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 FUNCTION RSPHARM(COSTH,PHI,LMAX,AL0,ALM) IMPLICIT NONE C Parameters C note: this must agree with corresponding parameter in C SUBROUTINE TARSPHARM above INTEGER NMAX PARAMETER(NMAX=10) C Arguments: INTEGER LMAX REAL COSTH,PHI,RSPHARM REAL AL0(0:NMAX) COMPLEX ALM(1:NMAX,1:NMAX) C Local variables: INTEGER L,M REAL AX,FAC,FAC2,FOURPI,SIGN,SUM,Y_L0 COMPLEX CXI,CXTERM c External functions REAL P_LM EXTERNAL P_LM C----------------------------------------------------------------------- C FUNCTION RSPHARM C Given: C COSTH = cos(theta) C PHI = phi (radians) C LMAX = maximum order of spherical harmonics C AL0 = real array of spherical harmonic amplitudes a_l0 C ALM = complex array of spherical harmonic amplitudes a_lm, m>0 C Returns C RSPHARM = distance from "center" to surface for target with C a_eff=1. C Note: we assume that a_{l,-m} = (-1)^m conjg(a_{l,m}) C so that sum over spherical harmonics gives real function C C History C 02.11.12 (BTD) first written C end history C C Copyright (C) 2002 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C----------------------------------------------------------------------- FOURPI=16.*ATAN(1.) CXI=(0.,1.) SUM=AL0(0)/SQRT(FOURPI) DO L=1,LMAX FAC=SQRT((2*L+1)/FOURPI) Y_L0=FAC*P_LM(L,0,COSTH) c*** diagnostic c write(0,*)'Y_L0=',Y_L0,' L=',L,' COSTH=',COSTH c*** SUM=SUM+AL0(L)*Y_L0 SIGN=1. DO M=1,L FAC=FAC/SQRT(REAL((L+M)*(L+1-M))) SIGN=-SIGN CXTERM=ALM(L,M)*EXP(CXI*M*PHI) c CXSUM=CXSUM+FAC*P_LM(L,M,COSTH)*(CXTERM+SIGN*CONJG(CXTERM)) SUM=SUM+FAC*P_LM(L,M,COSTH)*2.*REAL(CXTERM) ENDDO ENDDO RSPHARM=(FOURPI*SUM**2)**(1./3.) RETURN END FUNCTION P_LM(L,M,X) c arguments INTEGER L,M REAL P_LM,X c local variables INTEGER I,LL REAL FACT,PLL,PMM,PMMP1,SINTH c----------------------------------------------------------------------- c FUNCTION P_LM = Associated Legendre polynomial P_{lm}(x=cos(theta)) c given: c L,M c X c returns: c P_LM = P_lm(x) c associated Legendre polynomial c this routine follows Numerical Recipes c history c 02.11.12 (BTD) first written c end history c----------------------------------------------------------------------- PMM=1. IF(M.GT.0)THEN SINTH=SQRT(1.-X*X) FACT=1. DO I=1,M PMM=PMM*FACT*SINTH FACT=FACT+2. ENDDO ENDIF IF(L.EQ.M)THEN P_LM=PMM ELSE PMMP1=X*(2*M+1)*PMM IF(L.EQ.M+1)THEN P_LM=PMMP1 ELSE DO LL=M+2,L PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M) PMM=PMMP1 PMMP1=PLL ENDDO P_LM=PLL ENDIF ENDIF RETURN END SUBROUTINE TARTET(A1,A2,AX,DX,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 AX REAL A1(3),A2(3),DX(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 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 Input: C AX = length of one edge of tetrahedron C DX(1-3)=(dx,dy,dz)/d where dx,dy,dz=lattice spacings in x,y,z C directions, and d=(dx*dy*dz)**(1/3)=effective lattice C spacing 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 IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms in target C ICOMP(1-NAT,1-3)=1 (composition identifier) 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.08 (BTD): Changed DO...ENDDO to DO #...# CONTINUE C 90.11.09 (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 95.07.13 (BTD): Eliminated unused arguments AY,AZ C 96.01.26 (BTD): Replaced IX,IY,IZ by IXYZ C 97.12.26 (BTD): Added DX(3) to argument list to support noncubic C lattice. C *** Note: additional code modifications required C to support noncubic lattice! C 98.03.06 (BTD): If called with noncubic lattice (not yet supported) C call ERRMSG and halt with fatal error. C 98.03.07 (BTD): Modify to write DX to file "target.out" for C compatibility with REASHP C 00.11.02 (BTD): Add ICOMP to argument list C set ICOMP=1 at occupied sites C write ICOMP to target.out C end history C C Copyright (C) 1993,1995,1996,1997,1998,2000 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Current version of TARTET is restricted to cubic lattices IF(DX(1).NE.1..OR.DX(2).NE.1.)THEN CALL ERRMSG('FATAL','TARTET', & ' tartet does not support noncubic lattice') ENDIF 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 Determine list of occupied sites. NAT=0 DO 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 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 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 IXYZ(NAT,1)=I IXYZ(NAT,2)=J IXYZ(NAT,3)=K ENDIF ENDDO ENDIF ENDDO ENDDO C Initialize composition: DO K=1,3 DO I=1,NAT ICOMP(I,K)=1 ENDDO ENDDO 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. 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,DX 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(' >TARTET tetrahedral grain: S=',F8.3,/, &I7,3F8.4,' =NAT, XOFF,YOFF,ZOFF',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &3F9.6,' = lattice spacings (d_x,d_y,d_z)/d',/, &' JA IX IY IZ ICOMP(x,y,z)') 9030 FORMAT(I7,3I4,3I2) END SUBROUTINE VERSION(CSTAMP) IMPLICIT NONE C Arguments: CHARACTER CSTAMP*22 CSTAMP='DDSCAT.6.0 [04.02.25]' C******************************************************************************* C History of major changes to 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, eself, extend, 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 95.06.14 (BTD): Version 5a.1: first development C Differs from version 4d in computing: C full radiation pressure force (transverse components, C in additional to longitudinal component computed C previously); C Torque due to scattering and absorption of radiation C as measured by "vector torque cross section". C Changes made in: C DDSCAT C GETFML (including argument list) C SCAT (including argument list) C VERSION (CSTAMP is now CHARACTER*22, so that a date C can be attached to version number) C 95.06.15 (BTD): Added CXE to argument list of SCAT C Added new control parameters C CMDSOL (to choose solution method) C CMDTRQ (to choose whether or not to compute torques) C 95.06.20 (PJF+: Major modification of DDSCAT and GETFML to support C BTD): use of modular iterative solvers (e.g., CCGPACK and C PIM routines). C : add CMDSOL to argument list to specify method C of solution, with "supported" options C PETRKP = Petravic&Kuo-Petravic CCG routine (from C CCGPACK, coded by P.J. Flatau based on C earlier code by B.T. Draine) C PBCGST = Preconditioned Biconjugate Gradient with C Stabilization (from PIM package, coded by C R. Dias da Cunha and T. Hopkins) C : add CMDTRQ to argument list of REAPAR and GETFML C 95.07.21 (BTD): Major revision to subroutine SCAT C Added additional scratch arrays to DDSCAT,GETFML,SCAT C 95.07.24 (BTD): Corrected errors in SCAT C 95.08.14 (BTD): Modified GETFML and PROGRESS to include C COMMON/NORMRES/ERRSCAL to allow renormalization of C error |Ax-b| relative to |b|. C 96.01.05 (BTD): Version 5a.2: C Differs from version 5a.1 in including C shape options BLOCKS and DW1996, including C automatic calculation of principal axes A1 and A2 C with largest and second-largest moment of inertia C eigenvalue and eigenvector problem is solved using C routines BALANC, BALBAK, EIGV, ELMHS0, ELTRN0, HQR2, C IPMPAR, and SPMPAR taken from the Naval Surface C Warfare Center Library of Mathematics Subroutines. C 96.01.25 (BTD): Added shape option TWOSPH (routine TAR2SP, and C modifications to TARGET,REAPAR). C 96.02.23 (BTD): Subroutine PRINAXIS now carries out calculation of C principal axes and eigenvalues of moment of inertia C for arbitary grain shape (if invoked by appropriate C "target generation" subroutine), and prints out C 3 principal axes and eigenvalues C 96.10.17 (BTD): Version 5a.3: C Added NAT0 to argument list for routine ALPHA. C This required modification of ALPHA and GETFML. C 96.10.18 (BTD): Changed NEWTMP to GPFAFT (for Generalized Prime C Factor Algorithm for Fourier Transform). C 96.11.06 (BTD): Version 5a.4: C Removed code from DDSCAT and moved to subroutines C GETMUELLER and WRITESCA C Changed definition of scattering angle phi. C Previously, phi was measured from plane containing C incident k vector and Re(CXE01) C Henceforth, phi is measured from Lab x,y plane C (incident radiation is along x axis). C This involved changes to SCAVEC C Modified GETMUELLER so that for IORTH=2 it C automatically computes the complete Mueller scattering C matrix for each target orientation, and for the C orientational average. C Modified WRITESCA so that it writes out selected C elements of the Mueller matrix. C 96.11.14 (PJF): Add timers(mxtimers), ntimers to DDSCAT and getfml C timers contains information about CPU times of C various parts of the code and informations needed C for timming (e.g. number of iterations) C Further changes of "WRITESCA" to include binary C file option. This should be of use, for example, C for users of IDL. C Added "cbinflag" to DDSCA, reapar. Added C "cbinfile", "iobin" and "open" unformatted file in DDSCAT C "close(iobin) in DDSCAT. C Removed "getset" routine. Hardwired "ioerr" in C "errmsg", and "ioout" in "wrimsg". C Combined all "rotation" code in rot.2 C Added dd.pro package for DDSCAT output postprocessing. C The package includes IDL reading code of binary C DDSCAT file and Bohren and Huffman Mie routine coded in IDL. C Added 10,000 lines of LAPACK code to solve 3x3 C eigenvalue problem. Removed NSWC l(eispack) C driver because (a) it contains ald r1mach C tedious to support code (b) because LAPACK is C more modern. C Add several *.ps files to doc and BibTeX C bibliography of recent DDA references. C 96.11.15 (PJF): Added NetCDF binary option to "writesca.f" C One needs "include" in writesca.f which is non-standard C in FORTRAN77 but will probably work on all architectures. C Also -lnetcdf flag is needed during the loading step C if NetCDF is being used. Otherwise "empty" NetCDF routines C should be substittuted to satisfy loader (c.f. dummynetcdf.f) C "ncclose" is executed from DDSCAT.f instead of writesca.f C NetCDF is available from http://www.unidata.ucar.edu/ C for most popular computers. C Added "writebin.f" and "writecdf.f" for 5a7 C 96.11.21 (BTD): moved NCCLOS to WRITENET so that all dependence on C netCDF is through subroutine WRITENET C This required changes to DDSCAT,WRITESCA, and C WRITENET, and adding additional call to WRITESCA C from DDSCAT, and additional call to WRITENET from C WRITESCA C added variable IDVOUT2 to COMMON/NORMERR to permit C device number for standard output to be passed to C subroutine PROGRESS C 96.12.02 (BTD): Corrected error in DDSCAT.f : SMORI was not C initialized, and hence was incorrect except (possibly) C for first wavelength and radius. C Also eliminated some superfluous variables from C DDSCAT.f C 97.04.28 (BTD): Corrected bug in TARHEX: NAT was not initialized to 0 C 97.06.04 (BTD): Corrected bug in TARGET, affecting only target option C "UNICYL" (target composition was not being set C correctly). C 97.07.24 (BTD): Corrected bug in GETMUELLER: was using PHI rather C than PHIN in computation of amplitude scattering C matrix and Mueller matrix from f_ml computed by C GETFML. Required change in GETMUELLER argument, C list, with corresponding change required in DDSCAT.f C Also corrected bugs in GETMUELLER pertaining to C computation of Mueller matrix elements S_14 and C S_31. Note that because S_31 was incorrect, PPOL C (evaluated in WRITESCA) was numerically incorrect. C --------------- begin creation of version 6.0 ----------------------- C continuing to support "public" version 5a7 while C carrying on private development of version 6.0 C --------------------------------------------------------------------- C 98.01.01 (BTD): Version 6.0, with capability of using noncubic C rectangular lattices. Extensive changes required, C because polarizability tensor will now be nondiagonal C even when dielectric tensor is isotropic. C Changes include following modules: C DDSCAT C ALPHA C CMATVEC (in cprod.f) C CPROD C ESELF C EVALA C EVALE C EVALQ C GETFML C MATVEC (in cprod.f) C PRINAXIS C REASHP C SCAT C TARELL C TARGET C TARBLOCKS C TARCEL C TARCYL C TARHEX C TARREC C TARTET C TAR2EL C TAR3EL C TAR2SP C and additional module C NONCUBIC (called by subroutine ALPHA) C 98.01.20 (BTD): Modified subroutine ALPHA to use (and print out) C constant C4 in evaluation of first-order corrections C to alpha. C 98.03.10 (BTD): Changes to: C ERRMSG C REAPAR C TAR2EL C TAR2SP C TAR3EL C TARBLOCKS C TARCEL C TARELL C TARHEX C TARTET C 98.04.27 (BTD): Minor changes made in C ALPHA C CGCOMMON C EVALQ C GETFML C NONCUBIC C PRINAXIS C RESTORE C TARCYL C TARGET C 98.10.07 (BTD) Rewrote cgcommon.f, function SCNRM2, C because g77 optimization generates bad code C complicated code replaced with simple code C -------------- [separate creation of version 5a9 at this time] ---- C 98.10.08 (BTD) modified cgcommon.f, subroutine PROGRESS, C to keep track of minimum error, and to C note when new error minimum is attained. This is C to observe convergence behavior of PBCGST. C 98.10.26 (BTD) modified subroutine PROGRESS in cgcommon.f C to also print information on C rate of convergence (variable RATE). C 98.10.27 (BTD) modified subroutine SMACHCONS in cgcommon.f C because Solaris compiler complained about value of C OVERFLOW parameter. Reduced by factor 2 C 98.11.16 (BTD) Corrected error in value of CXE01R passed to WRITEBIN C for binary write option. C 98.12.07 (BTD) Found and corrected bug reported by Timo Nousiainen C (tpnousia@lumi.meteo.helsinki.fi) which caused Mueller C matrix elements S_12 and S_13 to be computed C due to inconsistency in definition of f_ml. Removed C code in GETFML which had been inserted to change C f_ml to refer to "usual" incident pol states l=1,2 C whereas GETFML assumed that incident pol states C l=1,2 were those specified by the user in ddscat.par C Stick with this latter convention. C Modified GETMUELLER so that Mueller matrix elements C are computed correctly even when user-specified C incident polarization state CXE01 is not linearly C polarized. C 98.12.07 (BTD) Edited UserGuide to give more detailed (and now correct) C discussion of computation of Mueller matrix elements C from the f_ml. C 98.12.07 (BTD) Inserted (commented out) code into getfml to allow C verification of final error in approximate solution C to confirm reliability of error reported by iteration C machinery. After confirming that these errors appear C to be reasonable, code was commented out, but can be C reactivated at future date if repetition of such tests C is desired. C 98.12.16 (BTD) In getfml.f, increase upper limit on iterations from C 300 to 10000 C 99.01.26 (BTD) Modified DDSCAT.f and reapar.f to allow user to specify C first IWAV,IRAD,IORI (normally 0 0 0). This allows C one to restart calculational runs which abort after C doing only some of the desired wavelengths, radii, C and orientations. Note: structure of ddscat.par is C changed as a result! C 99.03.01 (BTD) rewrote subroutine orient.f to make it handle odd C NTHETA as described in the UserGuide C 99.03.05 (BTD) corrected errors in new orient.f C 99.04.30 (BTD) corrected new errors in orient.f (introduced in C 99.03.05 revision) C 99.06.30 (BTD) modifications to reapar (was not reading correct C number of SHPAR values for case TWOSPH) and tar2sp C in order to deal correctly with two spheroid case. C (Although have not yet implemented changes necessary C to do TWOSPH with noncubic lattice.) C------------------------------------------------------------------------ C 00.06.12 (BTD) [Separate release of version 5a10 at this time] C------------------------------------------------------------------------ C 00.06.12 (BTD) New target option NSPHER supported by routine TARNSP C Replaced old LAPACK source code with latest version C of routines from www.netlib.org C New LAPACK routines required by PRINAXIS are now C collected in two files, lapackblas.f and lapacksubs.f C with new LAPACK routines, PRINAXIS now executes properly C when compiled with g77 C To support TARNSP, it was necessary to modify argument C lists for REAPAR and TARGET to include character C variable CFLSHP to pass name of file containing C description of multisphere target C Minor "tidying" of a number of files to make C all data type conversion explicit. C All variables in CXFFT3 now declared explicitly. C 01.04.21 (BTD) copied over changes made to REASHP in version 5a10 C modified all target routines so that in all cases C the target.out file created when using CALLTARGET C will contain composition information, so that there C is now a standard format for shape.dat input files C Modified as-delivered REASHP to read ICOMP. C 01.04.21 (BTD) copied over changes made to TARNSP to correct bug C in version 5a10 C 02.02.12 (BTD) added support for target option PRISM3, with C new routine TARPRSM and concomitant changes to C TARGET and REAPAR C 02.11.15 (BTD) working with Matthew Collinge to add MPI capability to C code. C DDSCAT.f now calls C MPI_INIT C MPI_COMM_RANK C MPI_COMM_SIZE C MPI_FINALIZE C created new module mpi_subs.f with following routines: C SUBROUTINE COLSUM, which calls MPI_REDUCE C SUBROUTINE SHARE1, which calls MPI_BCAST C SUBROUTINE SHARE2, which calls MPI_BCAST C for use on systems where MPI is not installed, C created new module mpi_fake.f with dummy routines C SUBROUTINE MPI_INIT C SUBROUTINE MPI_COMM_RANK C SUBROUTINE MPI_COMM_SIZE C SUBROUTINE MPI_BCAST C SUBROUTINE MPI_REDUCE C SUBROUTINE MPI_FINALIZE C so that calls to these subroutines can remain in C the code. C 03.02.13 (BTD) added 1 more digit of precision to qtable and qtable2 C outputs C 03.04.13 (BTD) working with Matthew Collinge to add MPI capability to C added character variable CFFLOG to contain name of C running output file C added routine NAMID to generate a unique name for output C file printed by each MPI-spawned process. C output to logfile is no longer unbuffered (if required C for debugging purposes, this can be reactivated -- see C comments in DDSCAT.f C added variables ITNUM(2) and MXITER to argument lists C for GETFLM and WRITESCA C 03.07.13 (BTD) added new variable ITNUMMX(2) to DDSCAT so that C waarbbori.avg file will contain *maximum* number of C iterations required by any of the orientations being C averaged over C 03.10.30 (BTD) corrected typo in WRITENET C 04.02.25 (BTD) corrected typo in TARPRSM C end history C C Copyright (C) 1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003, C 2004 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** RETURN END SUBROUTINE WRITEBIN(CBINFLAG,CBINFILE,IOBIN,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,DX,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & IORI,AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR, & CXE01R,CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G, & QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,SM,PHIN,THETAN, & TIMERS) C*********************************************************************** C Subroutine "writebin" writes FORTRAN unformatted binary file C to file cbinfile ('dd.bin') C C Original version created by P.J.Flatau, University of California, SD C History: C 96.11.15 (PJF): First version. C 98.01.11 (BTD): Added DX(3) to argument list, and write DX to C binary output upon entry (DX doesn't change). C 99.06.30 (BTD): Added SAVE IENTRY statement. C end history C Copyright (C) 1996,1998,1999 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C IMPLICIT NONE C .. Scalar Arguments .. REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD, & THTMID,THTMXD,TOL,WAVE,XX INTEGER ICTHM,IOBIN,IORI,IORTH,IPHIM,IRAD,IWAV,MXBETA,MXCOMP,MXN3, & MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA,MXTHET,MXTIMERS, & MXWAV,NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI,NRAD,NSCAT, & NTHETA,NTIMERS,NWAV,NX,NY,NZ CHARACTER CALPHA*6,CBINFLAG*6,CMDFFT*6,CMDTRQ*6,CSHAPE*6, & CBINFILE*14,CSTAMP*22,CDESCR*67 C .. Array Arguments .. COMPLEX CXE01R(3),CXE02R(3),CXEPS(MXCOMP),CXRFR(MXCOMP) REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),DX(3),G(2), & PHI(MXPHI),PHIN(MXSCA), & QABS(2),QAV(17),QBKSCA(2),QEXT(2),QPHA(2), & QSCAT(2),QTRQAB(3,2),QTRQSC(3,2),SM(MXSCA,4,4),THETA(MXTHET), & THETAN(MXSCA),TIMERS(MXTIMERS),WAVEA(MXWAV) INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3) C .. Local Scalars .. INTEGER I,IENTRY,IS,J C .. External Subroutines .. EXTERNAL WRIMSG C .. Data statements .. DATA IENTRY/0/ SAVE IENTRY IF(IENTRY.EQ.0) THEN CALL WRIMSG('WRITEBIN','First call to binary write') IENTRY=1 C General run information: WRITE(IOBIN)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,CMDTRQ C mostly dimensions: WRITE(IOBIN)IORTH,NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI, & NWAV,NRAD,NTHETA,NPHI,NBETA,NTIMERS C maximum dimension (not really needed) WRITE(IOBIN)MXNX,MXNY,MXNZ,MXCOMP,MXSCA WRITE(IOBIN)BETMID,BETMXD,THTMID,THTMXD,PHIMID,PHIMXD C dump scattering angles WRITE(IOBIN)(PHIN(I),I=1,NSCAT),(THETAN(I),I=1,NSCAT) C misc run information WRITE(IOBIN)ICTHM,IPHIM,TOL WRITE(IOBIN)(WAVEA(I),I=1,NWAV),(AEFFA(I),I=1,NRAD), & (THETA(I),I=1,NTHETA),(BETA(I),I=1,NBETA), & (PHI(I),I=1,NPHI),(A1(I),I=1,3),(A2(I),I=1,3), & (ICOMP(I),I=1,NAT3),(IXYZ0(I,1),I=1,NAT), & (IXYZ0(I,2),I=1,NAT),(IXYZ0(I,3),I=1,NAT),(DX(I),I=1,3) c endif "if(cbinflag.eq.'USEBIN') then" ENDIF C C binary write of for particular "iwav, irad, itheta, ibeta iphi" C Let us try to keep only one version of the binary file. Thus, in C case of IORTH=1 or 'DOTORQ' some values are "missing" C flatau -->draine do we need smori for "ppol" do we need cxf_ij ? C IF(IORI.NE.0.AND.CBINFLAG.EQ.'ALLBIN')THEN WRITE(IOBIN)IWAV,IRAD,IORI,(TIMERS(I),I=1,NTIMERS) WRITE(IOBIN)AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR,CXE01R, & CXE02R,(CXRFR(I),I=1,NCOMP), (CXEPS(I),I=1,NCOMP) WRITE(IOBIN) QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB,QTRQSC C write sm as 16 independent matrices (IDL handles 7 subscripts only) DO 10 I=1,4 DO 10 J=1,4 WRITE(IOBIN)(SM(IS,I,J),IS=1,NSCAT) 10 CONTINUE C endif "if(iori. ne. 0 .and. cbinflag.eq.'ALLCDF') then" ENDIF IF(IORI.EQ.0)THEN C C BTD 98.05.01 add CONTINUE statement in lieu of actual code C CONTINUE C C iori=0 case Cflatau ---> note. I need to add code for orientational averages C but may be I can recalculate it from sm for iwav, irad, itheta, C ---> dump weights C ENDIF RETURN END SUBROUTINE WRIMSG(CSUBRT,CMSGNM) IMPLICIT NONE 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 History: C 96.11.14 (PJF) Remove "getset" and hardwire" ioout" C 96.11.20 (BTD) change IOOUT to IDVOUT C*********************************************************************** C .. Scalar Arguments .. CHARACTER CMSGNM* (*), CSUBRT* (*) C .. C .. Local Scalars .. INTEGER IDVOUT SAVE IDVOUT C .. C Note: on some systems (e.g., Solaris) IDVOUT=0 generates unbuffered C output to "standard output". On other systems IDVOUT=0 may not be C valid; then set IDVOUT=6 to get "standard output", probably buffered. C DATA IDVOUT/0/ WRITE (IDVOUT,FMT=9000) CSUBRT, CMSGNM 9000 FORMAT (' >',A,' ',A) RETURN END SUBROUTINE WRITESCA(MXNX,MXNY,MXNZ,NX,NY,NZ, & DX,ICOMP,IXYZ0,MXNAT,MXN3,NAT,NAT3, & WAVEA,MXWAV,NWAV,AEFFA,MXRAD,NRAD, & ITHETA,IBETA,IPHI,THETA,MXTHET,NTHETA, & BETA,MXBETA,NBETA,PHI,MXPHI,NPHI,NSMELTS, & TIMERS,MXTIMERS,NTIMERS,CBINFILE,IOBIN, & CBINFLAG,CNETFILE,IDNC,CNETFLAG, & ICTHM,ILIN10,ILIN12,IORI,IORTH, & IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & ITNUM,MXITER, & NCOMP,NORI,NSCAT, & CALPHA,CDESCR,CFLAVG,CFLSCA, & CMDFFT,CMDTRQ,CSHAPE,CSTAMP, & AEFF,AK1,AKR,BETAD,BETMID,BETMXD, & PHID,PHIMID,PHIMXD,THETAD,THTMID, & THTMXD,TOL,WAVE,XX, & A1,A2,PHIN,QABS,QABSUM,QBKSCA, & QBKSUM,QEXSUM,QEXT,QPHA,QPHSUM, & QSCAG,QSCAT,QSCGSUM,QSCSUM,QTRQAB, & QTRQABSUM,QTRQSC,QTRQSCSUM, & S1111,S2121,SM,SMIND1,SMIND2,SMORI, & THETAN,CX1121,CXE01,CXE01R,CXE02, & CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) IMPLICIT NONE C Scalar Arguments .. REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD, & THTMID,THTMXD,TOL,WAVE,XX INTEGER IBETA,ICTHM,IDNC,ILIN10,ILIN12,IOBIN,IORI,IORTH,IPHI, & IPHIM,IRAD,ITHETA,IWAV,MXBETA,MXCOMP,MXITER,MXN3,MXNAT, & MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA,MXTHET,MXTIMERS,MXWAV, & NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI,NRAD,NSCAT,NSMELTS, & NTHETA,NTIMERS,NWAV,NX,NY,NZ INTEGER ITNUM(2) CHARACTER CALPHA*6,CBINFLAG*6,CMDFFT*6,CMDTRQ*6,CNETFLAG*6, & CSHAPE*6,CFLAVG*13,CBINFILE*14,CFLSCA*14,CNETFILE*14, & CSTAMP*22,CDESCR*67 C Array Arguments .. COMPLEX CX1121(MXSCA),CXE01(3),CXE01R(3),CXE02(3),CXE02R(3), & CXEPS(MXCOMP),CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA), & CXF22(MXSCA),CXRFR(MXCOMP) REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),DX(3), & PHI(MXPHI),PHIN(MXSCA),QABS(2), & QABSUM(2),QBKSCA(2),QBKSUM(2),QEXSUM(2),QEXT(2),QPHA(2), & QPHSUM(2),QSCAG(3,2),QSCAT(2),QSCGSUM(3,2),QSCSUM(2), & QTRQAB(3,2),QTRQABSUM(3,2),QTRQSC(3,2),QTRQSCSUM(3,2), & S1111(MXSCA),S2121(MXSCA),SM(MXSCA,4,4),SMORI(MXSCA,4,4), & THETA(MXTHET),THETAN(MXSCA),TIMERS(MXTIMERS),WAVEA(MXWAV) INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3) INTEGER SMIND1(9),SMIND2(9) C Local Scalars .. COMPLEX CXVAR1 REAL DEGRAD,MKD,PHIND,PI,PPOL,RVAR1,RVAR2,RVAR3,THETND INTEGER I,J,ND C Local Arrays .. REAL G(2),QAV(17) C External Subroutines .. EXTERNAL NAMER,WRITEBIN,WRITENET C Intrinsic Functions .. INTRINSIC CONJG,SQRT C Save statement .. SAVE DEGRAD C Data statements .. DATA DEGRAD/57.29578/ C*********************************************************************** C Purpose of this module is to collect a lot of the output code into C a single module. C If CNETFLAG .eq 'NCCLOS' then close NETCDF compliant output file and C return with no other action taken C If CNETFLAC .ne.'NCCLOS' then: C If IORI = 0, then output only orientational averages C If IORI > 0, then output for specific orientation C If CNETFLAG .ne.'NOTCDF' use WRITENET to write NETCDF compliant C output file c .eq.'NOTCDF' then do not write NETCDF compliant output file C If CBINFLAG .ne.'NOTBIN' use WRITEBIN to write binary output file C .eq.'NOTBIN' then do not write binary output file C----------------------------------------------------------------------- C History: C 96.11.14 (PJF) added binary write options (see also info in C "versions.f" about IDL postprocessing) C possible "checkpointing" could be done here "close(iobin)" C open(iobin,file=cbinfile,form='unformatted',access='append') C (which is non-standard Fortran 77 (but standard in F90) C open(iobin,file=cbinfile,form='unformatted',status='old') C 96.11.21 (BTD) added call to WRITENET if called with CNETFLAG='NCCLOS' C in order to get call to NCCLOS out of DDSCAT.f C 98.01.11 (BTD) added DX(3) to argument list to pass information on C nonuniform lattice spacing [development of DDSCAT.6.0] C 98.11.16 (BTD) corrected CX01R -> CXE01R in argument list of WRITEBIN C 98.12.12 (BTD) Introduce NSMELTS,SMIND1,SMIND2 to argument list to C allow selection of Muller matrix elements to be C printed out. Add code to allow from 1 to 9 matrix C elements to be printed. C 98.12.13 (BTD) Modify to print out validity criterion |m|kd C for each composition C 03.01.30 (BTD) disabled output line to print out dx,dy,dz C 03.02.13 (BTD) added one more digit of precision to output Q values C with appropriate modifications to format statements C 03.04.14 (BTD) added ITNUM(2) and MXITER to argument list C added ITNUM(2) = number of iterations C and MXITER = max number of iterations allowed C to output statements 9150 and 9155 C 03.09.01 (BTD) cosmetic change to output statement: g -> g(1) C 04.01.04 (BTD) add and modify comments C end history C C Copyright (C) 1996,1998,2003 C B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C If CNETFLAG = 'NCCLOS' then call WRITENET to close netCDF files C IF(CNETFLAG.EQ.'NCCLOS')THEN CALL WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,DX,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & ITHETA,IBETA,IPHI, & IORI,AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR, & CXE01R,CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G, & QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,SM,SMORI, & PHIN,THETAN,TIMERS) RETURN ENDIF C If IORI > 0, then write out scattering properties for this specific C target orientation. PI=4.*ATAN(1.) IF(IORI.GT.0)THEN C*** Compute g= DO J=1,IORTH G(J)=QSCAG(1,J)/QSCAT(J) ENDDO CALL NAMER(IWAV,IRAD,IORI,CFLSCA,CFLAVG) OPEN(UNIT=8,FILE=CFLSCA,STATUS='UNKNOWN') WRITE(8,FMT=9030)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0 C following code to be enabled for noncubic treatment: C & ,DX WRITE(8,FMT=9032)AEFF,WAVE,XX DO J=1,NCOMP MKD=(4.*PI/(3.*NAT0))**(1./3.)* & SQRT(REAL(CXRFR(J)*CONJG(CXRFR(J))))*XX WRITE(8,FMT=9031)CXRFR(J),CXEPS(J),MKD,J ENDDO WRITE(8,FMT=9033)TOL,ICTHM,IPHIM,A1,A2 WRITE(8,FMT=9035)AKR,CXE01R,CXE02R WRITE(8,FMT=9040)BETAD,THETAD,PHID WRITE(8,FMT=9050)QEXT(1),QABS(1),QSCAT(1),G(1),QBKSCA(1), & QPHA(1) IF(IORTH.EQ.1)THEN Cflatau assign "missing value" value if there is only one polarization DO I=1,17 QAV(I)=-999. ENDDO QEXT(2)=-999. QABS(2)=-999. QSCAT(2)=-999. QBKSCA(2)=-999. QPHA(2)=-999. G(2)=-999. DO J=1,3 QSCAG(J,2)=-999. ENDDO ENDIF IF(IORTH.EQ.2)THEN QAV(1)=.5*(QEXT(1)+QEXT(2)) QAV(2)=.5*(QABS(1)+QABS(2)) QAV(3)=.5*(QSCAT(1)+QSCAT(2)) QAV(4)=.5*(QSCAG(1,1)+QSCAG(1,2))/QAV(3) QAV(5)=.5*(QBKSCA(1)+QBKSCA(2)) QAV(6)=.5*(QPHA(1)+QPHA(2)) QAV(7)=QEXT(1)-QEXT(2) QAV(8)=QPHA(1)-QPHA(2) DO J=1,3 QAV(8+J)=.5*(QSCAG(J,1)+QSCAG(J,2)) ENDDO IF(CMDTRQ.EQ.'DOTORQ')THEN DO J=1,3 QAV(11+J)=.5*(QTRQAB(J,1)+QTRQAB(J,2)) QAV(14+J)=.5*(QTRQSC(J,1)+QTRQSC(J,2)) ENDDO ENDIF WRITE(8,FMT=9055)QEXT(2),QABS(2),QSCAT(2),G(2), & QBKSCA(2),QPHA(2),(QAV(J),J=1,8) ENDIF WRITE(8,FMT=9150)(QSCAG(J,1),J=1,3),ITNUM(1),MXITER IF(IORTH.EQ.2)THEN WRITE(8,FMT=9155)(QSCAG(J,2),J=1,3),ITNUM(2),MXITER, & (QAV(J),J=9,11) ENDIF IF(CMDTRQ.EQ.'DOTORQ')THEN WRITE(8,FMT=9250)(QTRQAB(J,1),J=1,3),(QTRQSC(J,1),J=1,3) IF(IORTH.EQ.2)WRITE(8,FMT=9255)(QTRQAB(J,2),J=1,3), & (QTRQSC(J,2),J=1,3),(QAV(J),J=12,17) ELSE C add missing value code (set to -999.) DO I=1,3 DO J=1,2 QTRQAB(I,J)=-999. QTRQSC(I,J)=-999. ENDDO ENDDO ENDIF IF(IORTH.EQ.1)THEN WRITE(8,FMT=9060) DO ND=1,NSCAT C convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) RVAR1=REAL(CONJG(CXF11(ND))*CXF11(ND)) RVAR2=REAL(CONJG(CXF21(ND))*CXF21(ND)) CXVAR1=CONJG(CXF11(ND))*CXF21(ND) WRITE(8,FMT=9070)ND,THETND,PHIND,RVAR1,RVAR2,CXVAR1 ENDDO ELSEIF(IORTH.EQ.2)THEN WRITE(8,FMT=9080) IF(NSMELTS.EQ.1)THEN WRITE(8,FMT=9081)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.2)THEN WRITE(8,FMT=9082)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.3)THEN WRITE(8,FMT=9083)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.4)THEN WRITE(8,FMT=9084)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.5)THEN WRITE(8,FMT=9085)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.6)THEN WRITE(8,FMT=9086)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.7)THEN WRITE(8,FMT=9087)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.8)THEN WRITE(8,FMT=9088)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.9)THEN WRITE(8,FMT=9089)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ENDIF DO ND=1,NSCAT C convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) C Compute PPOL = degree of linear polarization of scattered light C for incident unpolarized light PPOL=SQRT(SM(ND,2,1)**2+SM(ND,3,1)**2)/SM(ND,1,1) WRITE(8,FMT=9090)THETND,PHIND,PPOL, & (SMORI(ND,SMIND1(J),SMIND2(J)),J=1,NSMELTS) ENDDO ENDIF CLOSE(8) ELSE C This module writes out results from orientational averaging. C (arrive here when called with IORI=0) OPEN(UNIT=8,FILE=CFLAVG,STATUS='UNKNOWN') WRITE(8,FMT=9030)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0 C following code to be enabled for noncubic treatment: C & ,DX WRITE(8,FMT=9032)AEFF,WAVE,XX DO J=1,NCOMP MKD=(4.*PI/(3.*NAT0))**(1./3.)* & SQRT(REAL(CXRFR(J)*CONJG(CXRFR(J))))*XX WRITE(8,FMT=9031)CXRFR(J),CXEPS(J),MKD,J ENDDO WRITE(8,FMT=9033)TOL,ICTHM,IPHIM,A1,A2 RVAR1=AK1 RVAR2=0. RVAR3=0. WRITE(8,FMT=9037)RVAR1,RVAR2,RVAR3,CXE01,CXE02 WRITE(8,FMT=9042)BETMID,BETMXD,NBETA,THTMID,THTMXD,NTHETA, & PHIMID,PHIMXD,NPHI,NORI,IORTH RVAR1=QSCGSUM(1,1)/QSCSUM(1) WRITE(8,FMT=9050)QEXSUM(1),QABSUM(1),QSCSUM(1),RVAR1, & QBKSUM(1),QPHSUM(1) IF(IORTH.EQ.1)THEN WRITE(8,FMT=9150)(QSCGSUM(J,1),J=1,3) IF(CMDTRQ.EQ.'DOTORQ')WRITE(8,FMT=9250) & (QTRQABSUM(J,1),J=1,3),(QTRQSCSUM(J,1),J=1,3) WRITE(8,FMT=9060) DO ND=1,NSCAT C Convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) WRITE(8,FMT=9070)ND,THETND,PHIND,S1111(ND),S2121(ND), & CX1121(ND) ENDDO C*** Write orientationally-averaged Q values (except Qpha) to 'qtable': OPEN(UNIT=10,FILE='qtable',STATUS='OLD') DO ND=1,ILIN10 READ(10,FMT=*) ENDDO RVAR1=QSCGSUM(1,1)/QSCSUM(1) WRITE(10,FMT=9056)AEFF,WAVE,QEXSUM(1),QABSUM(1),QSCSUM(1), & RVAR1,QBKSUM(1) ILIN10=ILIN10+1 CLOSE(10) OPEN(UNIT=12,FILE='qtable2',STATUS='OLD') DO ND=1,ILIN12 READ(12,FMT=*) ENDDO WRITE(12,FMT=9057)AEFF,WAVE,QPHSUM(1) ILIN12=ILIN12+1 CLOSE(12) ELSEIF(IORTH.EQ.2)THEN RVAR1=QSCGSUM(1,2)/QSCSUM(2) QAV(1)=.5*(QEXSUM(1)+QEXSUM(2)) QAV(2)=.5*(QABSUM(1)+QABSUM(2)) QAV(3)=.5*(QSCSUM(1)+QSCSUM(2)) QAV(4)=.5*(QSCGSUM(1,1)+QSCGSUM(1,2))/QAV(3) QAV(5)=.5*(QBKSUM(1)+QBKSUM(2)) QAV(6)=.5*(QPHSUM(1)+QPHSUM(2)) QAV(7)=QEXSUM(1)- QEXSUM(2) QAV(8)=QPHSUM(1)- QPHSUM(2) DO J=1,3 QAV(8+J)=.5*(QSCGSUM(J,1)+QSCGSUM(J,2)) ENDDO IF(CMDTRQ.EQ.'DOTORQ')THEN DO J=1,3 QAV(11+J)=.5*(QTRQABSUM(J,1)+QTRQABSUM(J,2)) QAV(14+J)=.5*(QTRQSCSUM(J,1)+QTRQSCSUM(J,2)) ENDDO ENDIF WRITE(8,FMT=9055)QEXSUM(2),QABSUM(2),QSCSUM(2),RVAR1, & QBKSUM(2),QPHSUM(2),(QAV(J),J=1,8) WRITE(8,FMT=9150)(QSCGSUM(J,1),J=1,3),ITNUM(1),MXITER WRITE(8,FMT=9155)(QSCGSUM(J,2),J=1,3),ITNUM(2),MXITER, & (QAV(J),J=9,11) IF(CMDTRQ.EQ.'DOTORQ')THEN WRITE(8,FMT=9250)(QTRQABSUM(J,1),J=1,3), & (QTRQSCSUM(J,1),J=1,3) WRITE(8,FMT=9255)(QTRQABSUM(J,2),J=1,3), & (QTRQSCSUM(J,2),J=1,3),(QAV(J),J=12,17) ENDIF C*** Write orientationally-averaged Q values (except Qpha) to 'qtable': OPEN(UNIT=10,FILE='qtable',STATUS='OLD') DO ND=1,ILIN10 READ(10,FMT=*) ENDDO WRITE(10,FMT=9056)AEFF,WAVE,QAV(1),QAV(2),QAV(3),QAV(4), & QAV(5) ILIN10=ILIN10+1 CLOSE(10) OPEN(UNIT=12,FILE='qtable2',STATUS='OLD') DO ND=1,ILIN12 READ(12,FMT=*) ENDDO WRITE(12,FMT=9057)AEFF,WAVE,QAV(6),QAV(7),QAV(8) ILIN12=ILIN12+1 CLOSE(12) WRITE(8,FMT=9080) IF(NSMELTS.EQ.1)THEN WRITE(8,FMT=9081)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.2)THEN WRITE(8,FMT=9082)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.3)THEN WRITE(8,FMT=9083)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.4)THEN WRITE(8,FMT=9084)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.5)THEN WRITE(8,FMT=9085)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.6)THEN WRITE(8,FMT=9086)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.7)THEN WRITE(8,FMT=9087)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.8)THEN WRITE(8,FMT=9088)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ELSEIF(NSMELTS.EQ.9)THEN WRITE(8,FMT=9089)(SMIND1(J),SMIND2(J),J=1,NSMELTS) ENDIF DO ND=1,NSCAT C Convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) C Compute PPOL = degree of linear polarization of scattered light C for incident unpolarized light (Bohren & Huffman p. 382) PPOL=SQRT(SMORI(ND,2,1)**2+SMORI(ND,3,1)**2)/ & SMORI(ND,1,1) WRITE(8,FMT=9090)THETND,PHIND,PPOL, & (SMORI(ND,SMIND1(J),SMIND2(J)),J=1,NSMELTS) ENDDO ENDIF CLOSE(8) ENDIF C write FORTRAN unformatted file IF(CBINFLAG.NE.'NOTBIN')THEN CALL WRITEBIN(CBINFLAG,CBINFILE,IOBIN,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,DX,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & IORI,AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR, & CXE01R,CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G, & QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,SM, & PHIN,THETAN,TIMERS) ENDIF IF(CNETFLAG.NE.'NOTCDF')THEN CALL WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,DX,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & ITHETA,IBETA,IPHI, & IORI,AEFF,WAVE,XX,AK1,BETAD,THETAD,PHID,AKR, & CXE01R,CXE02R,CXRFR,CXEPS,QEXT,QABS,QSCAT,G, & QBKSCA,QPHA,QAV,QTRQAB,QTRQSC,SM,SMORI, & PHIN,THETAN,TIMERS) ENDIF RETURN 9030 FORMAT(' DDSCAT --- ',A,/,' TARGET ---',A,/,' ',A, & ' --- method of solution ',/,' ',A, & ' --- prescription for polarizabilies',/,' ',A, & ' --- shape ',/,I7,' = NAT0 = number of dipoles',/, & /) c following code to be enabled for noncubic treatment C & 3F8.5,' = normalized lattice spacings dx,dy,dz') 9031 FORMAT('n= (',F7.4,' , ',F7.4,'), eps.= (',F8.4,' , ',F7.4, & ') |m|kd=',0PF8.4,' for subs.',I2) 9032 FORMAT(' AEFF=',F10.5,' =',' effective radius (physical units)', & /,' WAVE=',F10.5,' = wavelength (physical units)',/, & 'K*AEFF=',F10.5,' = 2*pi*aeff/lambda') 9033 FORMAT(' TOL=',1P,E10.3,' = error tolerance for CCG method',/, & ' ICTHM=',I3,' = theta values used in comp. of Qsca,g',/, & ' IPHIM=',I3,' = phi values used in comp. of Qsca,g',/,1X, & 0P,'(',F8.5,2F9.5,') = target axis A1 in Target Frame',/, & 1X,'(',F8.5,2F9.5,') = target axis A2 in Target Frame') 9035 FORMAT(1X,'(',F8.5,2F9.5,') = k vector (latt. units) in TF',/,1X, & 3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 1 in TF',/,1X, & 3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 2 in TF') 9037 FORMAT(1X,'(',F8.5,2F9.5, & ') = k vector (latt. units) in Lab Frame',/,1X, & 3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 1 in LF',/,1X, & 3 ('(',F8.5,',',F8.5,')'),'=inc.pol.vec. 2 in LF') 9040 FORMAT(' BETA =',F7.3,' = rotation of target around A1',/, & ' THETA=',F7.3,' = angle between A1 and k',/,' PHI =', & F7.3,' = rotation of A1 around k') 9042 FORMAT(2F8.3,' = beta_min, beta_max ; NBETA =',I2,/,2F8.3, & ' = theta_min, theta_max; NTHETA=',I2,/,2F8.3, & ' = phi_min, phi_max ; NPHI =',I2,/, & ' Results averaged over ',I3,' target orientations',/, & ' and ',I3,' incident polarizations') 9050 FORMAT(10X,'Qext',7X,'Qabs',8X,'Qsca',5X,'g(1)=',3X,'Qbk',8X, & 'Qpha',/,1X,'JO=1: ',1P3E11.4,1PE12.4,1P2E11.4) 9055 FORMAT(1X,'JO=2: ',1P,3E11.4,E12.4,2E11.4,/,1X,'mean: ',3E11.4, & E12.4,2E11.4,/,1X,'Qpol= ',E11.4,39X,'dQpha=',E11.4) 9056 FORMAT(1PE10.4,1PE11.4,1P4E12.4,1PE11.4) 9057 FORMAT(1PE10.4,1PE11.4,1P3E12.4) 9060 FORMAT('**** Selected scattering directions ', & ' [note: incident pol state 1 is FIXED!]',/, & ' ND THETA PHI <|f11|^2> <|f21|^2> Re', & ' Im') 9070 FORMAT(I3,2F6.1,1P,2E10.3,2E11.3) C 9080 FORMAT(' ***************** Selected scattering directions *****', C & '************',/, C & ' ND THETA PHI <|f11|^2> <|f21|^2> <|f12|^2> ', C & '<|f22|^2> DEPOLR PPOL') 9080 FORMAT(12X,'** Mueller matrix elements for selected scattering ', & 'directions **') 9081 FORMAT('theta phi Pol. S_',I1,I1) 9082 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1) 9083 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1) 9084 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1) 9085 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1) 9086 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1) 9087 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1) 9088 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1) 9089 FORMAT('theta phi Pol. S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1,7X, & 'S_',I1,I1,7X,'S_',I1,I1,7X,'S_',I1,I1) 9090 FORMAT(F5.1,F6.1,F9.5,1P9E11.3) 9150 FORMAT(9X,'Qsca*g(1) Qsca*g(2) Qsca*g(3) iter mxiter',/, & 1X,'JO=1:',1P3E12.4,2I7) 9155 FORMAT(1X,'JO=2:',1P3E12.4,2I7,/,1X,'mean:',1P3E12.4) 9250 FORMAT(8X,'Qtrqab(1) Qtrqab(2) Qtrqab(3) ', & 'Qtrqsc(1) Qtrqsc(2) Qtrqsc(3)',/,1X,'JO=1:',1P,6E11.3) 9255 FORMAT(1X,'JO=2:',1P,6E11.3,/,1X,'mean:',1P,6E11.3) END