SUBROUTINE TARGET(A1,A2,CSHAPE,CFLSHP,IDVSHP,IOSHP,CDESCR,MXNAT, & SHPAR,NAT,IXYZ,ICOMP,IDVOUT, & MXN3,NAT3) C C Arguments: CHARACTER CDESCR*67,CFLSHP*13,CSHAPE*6 INTEGER IDVOUT,IDVSHP,IOSHP,MXNAT,NAT,MXN3,NAT3 REAL A1(3),A2(3),SHPAR(6) INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3) EXTERNAL ERRMSG,REASHP,TARCYL,TARELL,TARHEX,TARREC,TARTET C C Local variables: INTEGER BLKSIZ,IPRINAX,J,J1,L,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 C IOSHP = device number for "target.out" file C = -1 to suppress printing of "target.out" C SHPAR(1-6) = up to 6 parameters defining target geometry C MXNAT = limit on largest allowed value of NAT C MXN3 = 3*MXNAT C IDVOUT = 1 to print out information on target shape C -1 to suppress printing information on target shape C CFLSHP = name of target shape file if CSHAPE='FRMFIL' C IDVSHP = device number to use to read target shape if C CSHAPE='FRMFIL' C Returns: C A1,A2 = two 3-vectors defining initial target orientation C CDESCR = string describing target C NAT = number of dipoles in target C NAT3 = 3*NAT C 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 end history C C Copyright (C) 1993,1994,1995,1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C******************************************************************** IF(CSHAPE.EQ.'RCTNGL')THEN CALL TARREC(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) C Homogeneous target: DO 100 J=1,NAT DO 95 l=1,3 ICOMP(J,L)=1 95 CONTINUE 100 CONTINUE ELSEIF(CSHAPE.EQ.'ELLIPS')THEN CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) C Homogeneous target: DO 110 J=1,NAT DO 105 L=1,3 ICOMP(J,L)=1 105 CONTINUE 110 CONTINUE ELSEIF(CSHAPE.EQ.'CYLNDR')THEN CALL TARCYL(A1,A2,SHPAR(1),SHPAR(2),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) C Homogeneous target: DO 120 J=1,NAT DO 115 L=1,3 ICOMP(J,L)=1 115 CONTINUE 120 CONTINUE ELSEIF(CSHAPE.EQ.'HEXGON')THEN CALL TARHEX(A1,A2,SHPAR(1),SHPAR(2),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) C Homogeneous target: DO 130 J=1,NAT DO 125 L=1,3 ICOMP(J,L)=1 125 CONTINUE 130 CONTINUE ELSEIF(CSHAPE.EQ.'TETRAH')THEN CALL TARTET(A1,A2,SHPAR(1),CDESCR,IOSHP,MXNAT,NAT,IXYZ) C Homogeneous target: DO 140 J=1,NAT DO 135 L=1,3 ICOMP(J,L)=1 135 CONTINUE 140 CONTINUE ELSEIF(CSHAPE.EQ.'UNICYL')THEN CALL TARCYL(A1,A2,SHPAR(1),SHPAR(2),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) 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 DO 150 J=1,NAT ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=2 150 CONTINUE C ELSEIF(CSHAPE.EQ.'ANIELL')THEN CALL TARELL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) C Dielectric constants 1,2,3 for E along directions x,y,z in Target Frame DO 160 J=1,NAT ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 160 CONTINUE C C TWOELL -> two touching identical ellipsoids, of materials 1 and 2 C second ellipsoid is displaced from first in x-direction ELSEIF(CSHAPE.EQ.'TWOELL')THEN CALL TAR2EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) J1=NAT/2 DO 165 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=1 ICOMP(J,3)=1 165 CONTINUE DO 166 J=J1+1,NAT ICOMP(J,1)=2 ICOMP(J,2)=2 ICOMP(J,3)=2 166 CONTINUE C C TWOAEL -> two touching identical ellipsoids with anisotropic C functions. C second ellipsoid is displaced from first in x-direction C 1st ellipsoid has dielectric fcns 1,2,3 in x,y,z dirs. C 2nd ellipsoid has dielectric fcns 4,5,6 in x,y,z dirs. ELSEIF(CSHAPE.EQ.'TWOAEL')THEN CALL TAR2EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) J1=NAT/2 DO 175 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 175 CONTINUE DO 176 J=J1+1,NAT ICOMP(J,1)=4 ICOMP(J,2)=5 ICOMP(J,3)=6 176 CONTINUE C C THRELL -> three touching identical ellipsoids, of materials 1,2,3 C second ellipsoid is displaced from first in +x-direction C third ellipsoid is displaced from second in +x-direction ELSEIF(CSHAPE.EQ.'THRELL')THEN CALL TAR3EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) J1=NAT/3 DO 185 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=1 ICOMP(J,3)=1 185 CONTINUE DO 186 J=J1+1,2*J1 ICOMP(J,1)=2 ICOMP(J,2)=2 ICOMP(J,3)=2 186 CONTINUE DO 187 J=2*J1+1,NAT ICOMP(J,1)=3 ICOMP(J,2)=3 ICOMP(J,3)=3 187 CONTINUE C C THRAEL -> three touching identical ellipsoids with anisotropic C functions. C second ellipsoid is displaced from first in x-direction C 1st ellipsoid has dielectric fcns 1,2,3 in x,y,z dirs. C 2nd ellipsoid has dielectric fcns 4,5,6 in x,y,z dirs. C 3rd ellipsoid has dielectric fcns 7,8.9 in x,y,z dirs. ELSEIF(CSHAPE.EQ.'THRAEL')THEN CALL TAR3EL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),CDESCR,IOSHP, & MXNAT,NAT,IXYZ) J1=NAT/3 DO 195 J=1,J1 ICOMP(J,1)=1 ICOMP(J,2)=2 ICOMP(J,3)=3 195 CONTINUE DO 196 J=J1+1,2*J1 ICOMP(J,1)=4 ICOMP(J,2)=5 ICOMP(J,3)=6 196 CONTINUE DO 197 J=2*J1+1,NAT ICOMP(J,1)=7 ICOMP(J,2)=8 ICOMP(J,3)=9 197 CONTINUE C C CONELL -> two concentric ellipsoids with isotropic dielectric functions C 1 (inner) and 2 (outer) ELSEIF(CSHAPE.EQ.'CONELL')THEN CALL TARCEL(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5), & SHPAR(6),CDESCR,IOSHP,MXNAT,NAT,IXYZ,ICOMP) C C TAR2SP -> two touching spheroids C ELSEIF(CSHAPE.EQ.'TWOSPH')THEN CALL TAR2SP(A1,A2,SHPAR(1),SHPAR(2),SHPAR(3),SHPAR(4),SHPAR(5), & SHPAR(6),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,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP, & CDESCR,MXNAT,NAT,IXYZ,ICOMP) C C Homogeneous target: C DO J=1,NAT DO L=1,3 ICOMP(J,L)=1 ENDDO ENDDO C C DW1996 -> the 13-cube target studied by Draine & Weingartner 1996 C SHPAR(1)=width per block in lattice units C ELSEIF(CSHAPE.EQ.'DW1996')THEN IPRINAX=1 NBLOCKS=13 BLKSIZ=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,NBLOCKS,BLKSIZ,XYZB,IPRINAX,IOSHP, & CDESCR,MXNAT,NAT,IXYZ,ICOMP) CDESCR='13-cube target used by Draine & Weingartner 1996' C Homogeneous target: DO J=1,NAT DO L=1,3 ICOMP(J,L)=1 ENDDO ENDDO C C add any other specific shape routines here C ELSEIF(CSHAPE.EQ.'FRMFIL')THEN CALL REASHP(CFLSHP,IDVOUT,IDVSHP,A1,A2,CDESCR,MXNAT,NAT, & IXYZ,ICOMP,MXN3,NAT3) 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 C Before returning, compute: C NAT3=3*NAT RETURN 9010 FORMAT(' >TARGET:',A,/, & 7X,I7,' = NAT0 = number of dipoles in target') END SUBROUTINE TARHEX(A1,A2,A,B,CDESCR,IOSHP,MXNAT,NAT,IXYZ) C Scalar arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3) REAL A,B C Array arguments: REAL A1(3),A2(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 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 end history C C Copyright (C) 1993,1995,1996,1997 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DO 0100 JA=1,3 A1(JA)=0. A2(JA)=0. 0100 CONTINUE A1(1)=1. A2(2)=1. C 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 1200 JZ=1,NMZ Z=REAL(JZ) DO 1100 JY=1,NMY Y=REAL(JY) DO 1000 JX=1,NMX X=REAL(JX) CALL TESTHEX(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 1000 CONTINUE 1100 CONTINUE 1200 CONTINUE C Now compute some geometric properties of pseudo-hexagonal prism C NLAY = number of layers = effective length of prism/d C NFAC = number of atoms in one hexagonal face NLAY=INT(A) NFAC=NAT/NLAY C BEFF = effective length of hexagonal size/d C computed for hexagon of area NFAC*d**2 C =(3/2)*SQRT(3)*BEFF**2 C Note: for hexagon, vertex-vertex diameter = 2*BEFF BEFF=.6204032*SQRT(FLOAT(NFAC)) C ASPR = effective aspect ratio of target = length/(2*BEFF) ASPR=.5*FLOAT(NLAY)/BEFF C Note: string CDESCR will be printed by subr. TARGET WRITE(CDESCR,FMT='(A,I7,A)')' Hexagonal prism of NAT=',NAT, & ' dipoles' C Here print any additional target info which is desired by using C subroutine WRIMSG WRITE(CMSGNM,FMT='(A,3F7.4,A)') & ' A1 = (',A1,') = hexagon symmetry axis' CALL WRIMSG('TARHEX',CMSGNM) WRITE(CMSGNM,FMT='(A,3F7.4,A)') & ' A2 = (',A2,') = center-vertex axis' CALL WRIMSG('TARHEX',CMSGNM) WRITE(CMSGNM,FMT='(A,I7,A,I4,A,I4,A,F7.4)') & ' NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY,' aspect ratio=',ASPR CALL WRIMSG('TARHEX',CMSGNM) WRITE(CDESCR,FMT='(A,I7,A,I4,A,I4,A,F7.4)') & ' Hex prism,NAT=',NAT,' NFAC=',NFAC,' NLAY=',NLAY, & ' asp.ratio=',ASPR C Now store atom coordinates in file IF(IOSHP.GT.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') C*** 3 lines of general description allowed: WRITE(IOSHP,9020)A,B,NAT,A1,A2 C*** DO 2000 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) 2000 CONTINUE 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',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) 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 TARREC(A1,A2,XV,YV,ZV,CDESCR,IOSHP,MXNAT,NAT,IXYZ) C*********************************************************************** C Routine to construct rectangular prism from "atoms" C Input: C XV=(x-length)/d (d=lattice spacing) C YV=(y-length)/d C ZV=(z-length)/d C IOSHP=device number for "target.out" file C =-1 to suppress printing of "target.out" C MXNAT=dimensioning information (max number of atoms) C Output: C A1(1-3)=unit vector (1,0,0) defining target axis 1 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 C end history C C Copyright (C) 1993,1995,1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments: REAL XV,YV,ZV INTEGER IOSHP,MXNAT,NAT CHARACTER CDESCR*67 INTEGER*2 IXYZ(MXNAT,3) REAL A1(3),A2(3) C C Local variables: INTEGER JA,JX,JY,JZ,NX,NY,NZ CHARACTER CMSGNM*70 C C External subroutines: EXTERNAL ERRMSG C C Intrinsic functions: INTEGER IFIX INTRINSIC IFIX C*********************************************************************** NX=IFIX(XV+0.5) NY=IFIX(YV+0.5) NZ=IFIX(ZV+0.5) WRITE(CMSGNM,FMT='(A,3I4)') & ' Rectangular prism; NX,NY,NZ=',NX,NY,NZ WRITE(CDESCR,FMT='(A,3I4)') & ' Rectangular prism; NX,NY,NZ=',NX,NY,NZ C C Specify target axes A1 and A2 C New convention: A1 will be (1,0,0) in target frame C A2 will be (0,1,0) in target frame C Convention: A1 will be along longest dimension C A2 will be along intermediate dimension DO 0100 JA=1,3 A1(JA)=0. A2(JA)=0. 0100 CONTINUE A1(1)=1. A2(2)=1. C C Now populate lattice: NAT=0 DO 30 JZ=1,NZ DO 20 JY=1,NY DO 10 JX=1,NX NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT') ENDIF IXYZ(NAT,1)=JX IXYZ(NAT,2)=JY IXYZ(NAT,3)=JZ 10 CONTINUE 20 CONTINUE 30 CONTINUE IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARREC',' NAT.GT.MXNAT ') ENDIF WRITE(CMSGNM,FMT='(A, I7, A)') & ' Rectangular prism of NAT=',NAT,' dipoles' IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)XV,YV,ZV,NAT,A1,A2 DO 40 JA=1,NAT WRITE(IOSHP,FMT=9030)JA,IXYZ(JA,1),IXYZ(JA,2),IXYZ(JA,3) 40 CONTINUE CLOSE (UNIT=IOSHP) ENDIF RETURN 9020 FORMAT (' >TARREC rectangular prism; AX,AY,AZ=',3F8.4,/, &I7,' = NAT',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) END SUBROUTINE TARTET(A1,A2,AX,CDESCR,IOSHP,MXNAT,NAT,IXYZ) C C** Arguments: CHARACTER CDESCR*67 INTEGER IOSHP,MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3) REAL AX REAL A1(3),A2(3) C** Local variables: INTEGER I,IMAX,IMIN,J,JMAX,JMIN,JX,K,KMAX,KMIN REAL FY,S,X,XOFF,Y,YMAX,YMIN,YOFF,Z,ZMAX,ZMAX0,ZOFF C*********************************************************************** C C Routine to construct regular tetrahedral target array C one of the faces is parallel to the y-z plane C one of the edges of this face is parallel to x-y plane C C Input: C AX = length of one edge of tetrahedron C IOSHP=device number for "target.out" file C =-1 to suppress printing of "target.out" C MXNAT = dimensioning information C Returns: C A1(1-3) = unit vector (1,0,0) (along one axis of tetrahedron) C A2(1-3) = unit vector (0,1,0) C CDESCR = string describing target (up to 67 chars.) C NAT = number of atoms in target C IXYZ(1-NAT,1-3)=x/d,y/d,z/d for atoms in target C C Occupied array sites are those within tetrahedral surface C Size: C S=length of each side of tetrahedron C Volume = S**3/(6*sqrt(2)) C Orientation: C Center of mass is at origin. C One face is parallel to yz plane. C Projection onto yz plane has vertex on y axis. C S=length of one side (in lattice units) C Vertices are at C A=( sqrt(3/8), 0, 0)*S C B=(-sqrt(1/24), sqrt(1/3), 0)*S C C=(-sqrt(1/24),-sqrt(1/12),-1/2)*S C D=(-sqrt(1/24),-sqrt(1/12), 1/2)*S C Length in x direction = S*sqrt(2/3) C y = S*sqrt(3)/2 C z = S C Angle(AOB)=arccos(-1/3)=109.4712 deg. C OA=OB=OC=OD=sqrt(3/8)*S C Occupied sites are assumed to be located at C (X,Y,Z)=(I+XOFF,J+YOFF,K+ZOFF) C where I,J,K are integers, and XOFF,YOFF,ZOFF are constants. C Program sets XOFF,YOFF,ZOFF depending on choice of parameter S. C C Criterion for choosing XOFF: C Base of tetrahedron is located at x = -sqrt(1/24)*S C Let IMIN be value of I for this plane C Choose XOFF so that IMIN+XOFF = -sqrt(1/24)*S + 0.5 C with -0.5 < XOFF < 0.5 C Criterion for choosing YOFF: C One edge of tetrahedron is located at y= -S/(2*sqrt(3)) C Let JMIN be value of J for this line C Choose YOFF so that JMIN+YOFF = -S/sqrt(12) + 0.5 C Criterion for choosing ZOFF: C One edge of tetrahedron is parallel to z axis C Choose ZOFF in order to have number of dipoles along C this edge as close as possible to S C e.g., if S=odd integer, then take ZOFF=0 to place dipoles C at integral locations C if S=even integer, then take ZOFF=1/2 to place dipoles C at half-integral locations C C B.T.Draine, Princeton Univ. Obs., 90.10.30 C History: C 90.11.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 end history C C Copyright (C) 1993,1995,1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** S=AX C Set XOFF (and IMIN,IMAX): IMIN=-INT(S*SQRT(1./24.)) XOFF=0.5-S*SQRT(1./24.)-IMIN IMAX=IMIN+INT(S*SQRT(2./3.)+0.5)-1 C Set YOFF (and JMIN,JMAX): JMIN=-INT(S/SQRT(12.)) YOFF=0.5-S/SQRT(12.)-JMIN JMAX=JMIN+INT(S*SQRT(.75)+0.5)-1 C Set ZOFF (and KMIN,KMAX): C Determine whether S is closest to even or odd integer. C (Temporarily let KMIN be integer which S is closest to) KMIN=INT(S+0.5) C If KMIN is even, then ZOFF=0.5 C If KMIN is odd, then ZOFF=0. ZOFF=0. IF(KMIN-2*(KMIN/2).EQ.0)ZOFF=0.5 KMIN=-INT(0.5*S+ZOFF) KMAX=KMIN+INT(S+0.5)-1 C C Determine list of occupied sites. NAT=0 DO 30 I=IMIN,IMAX X=REAL(I)+XOFF C YMAX=largest value of Y which can occur for this X value C YMIN=smallest value of Y which can occur for this X value C ZMAX0=largest value of Z which can occur for this X value YMAX=S*SQRT(3./16.)-X/SQRT(2.) YMIN=-S*SQRT(3./64.)+X/SQRT(8.) ZMAX0=3.*S/8.-X*SQRT(3./8.) DO 20 J=JMIN,JMAX Y=REAL(J)+YOFF IF(Y.GE.YMIN.AND.Y.LE.YMAX)THEN C ZMAX=largest value of Z which can occur for this (X,Y) FY=(Y-YMIN)/(YMAX-YMIN) ZMAX=(1.-FY)*ZMAX0 DO 10 K=KMIN,KMAX Z=REAL(K)+ZOFF IF(ABS(Z).LE.ZMAX)THEN C Site is occupied: NAT=NAT+1 IF(NAT.GT.MXNAT)THEN CALL ERRMSG('FATAL','TARTET', & ' NAT.GT.MXNAT ') ENDIF IXYZ(NAT,1)=I IXYZ(NAT,2)=J IXYZ(NAT,3)=K ENDIF 10 CONTINUE ENDIF 20 CONTINUE 30 CONTINUE C C*** Specify vectors A1 and A2 which will define target orientation C A1 is initially along x-axis C A2 is initially along y-axis A1(1)=1. A1(2)=0. A1(3)=0. A2(1)=0. A2(2)=1. A2(3)=0. C WRITE(CDESCR,FMT='(A,I7,A)')' Tetrahedron of NAT=',NAT,' dipoles' IF(IOSHP.GE.0)THEN OPEN(UNIT=IOSHP,FILE='target.out',STATUS='UNKNOWN') WRITE(IOSHP,FMT=9020)S,NAT,XOFF,YOFF,ZOFF,A1,A2 DO 40 JX=1,NAT WRITE(IOSHP,FMT=9030)JX,IXYZ(JX,1),IXYZ(JX,2),IXYZ(JX,3) 40 CONTINUE CLOSE(UNIT=IOSHP) ENDIF RETURN 9020 FORMAT(' >TARTET tetrahedral grain: S=',F8.3,/, &I7,3F8.4,' =NAT, XOFF,YOFF,ZOFF',/, &3F9.4,' = A_1 vector',/, &3F9.4,' = A_2 vector',/, &' JA IX IY IZ') 9030 FORMAT(I7,3I4) END SUBROUTINE TIMEIT(CMSGTM,DTIME) C C timeit_vms C 94.06.20 (PJF) add dtime to the formal parameters C 96.06.19 (PJF) modified to suppress printing if CMSGTM='noprint' C This version of TIMEIT is for VMS systems. C C Copyright (C) 1993,1995 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. c*********************************************************************** C .. Scalar Arguments .. CHARACTER CMSGTM* (*) C .. C .. Local Scalars .. REAL DTIME, OLDTIM, TIME INTEGER IAD CHARACTER CSTA*3, CMSGNM*70 C .. C .. External Subroutines .. EXTERNAL LIB$INIT_TIMER, LIB$SHOW_TIMER, WRIMSG C .. C .. Data statements .. DATA CSTA/'ONE'/ C .. IF(CSTA.EQ.'ONE')THEN CSTA='TWO' IAD=0 CALL LIB$INIT_TIMER(IAD) ELSEIF(CSTA.EQ.'TWO')THEN CSTA='ONE' IF(CMSGTM(1:7).NE.'NOPRINT'.AND.CMSGTM(1:7).NE.'noprint')THEN CALL LIB$SHOW_TIMER(IAD,2) WRITE(CMSGNM,FMT='(A,A)')' Timing results for: ',CMSGTM CALL WRIMSG('TIMEIT',CMSGNM) ENDIF ENDIF RETURN END SUBROUTINE VERSION(CSTAMP) C Arguments: CHARACTER CSTAMP*22 CSTAMP='DDSCAT.5a.8 [97.7.24]' 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 end history C C Copyright (C) 1993,1994,1995,1996,1997 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,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 C end history C Copyright (C) 1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C 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),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/ 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) c endif "if(cbinflag.eq.'USEBIN') then" END IF 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" END IF IF (IORI.EQ.0) THEN 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 weigths END IF RETURN END SUBROUTINE WRIMSG(CSUBRT,CMSGNM) C Standard procedure for writing messages C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C C 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,ICOMP,IXYZ0,MXNAT, & MXN3,NAT,NAT3,WAVEA,MXWAV,NWAV,AEFFA,MXRAD, & NRAD,ITHETA,IBETA,IPHI,THETA,MXTHET,NTHETA, & BETA,MXBETA,NBETA,PHI,MXPHI,NPHI,TIMERS, & MXTIMERS,NTIMERS,CBINFILE,IOBIN,CBINFLAG, & CNETFILE,IDNC,CNETFLAG,ICTHM,ILIN10,ILIN12, & IORI,IORTH,IPHIM,IRAD,IWAV,MXCOMP,MXSCA,NAT0, & 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,SMORI,THETAN, & CX1121,CXE01,CXE01R,CXE02,CXE02R,CXEPS,CXRFR, & CXF11,CXF21,CXF12,CXF22) C flatau note to Draine ---> there is a bug in PPOL (calculated from smori) C Bruce --- phin looks funny. Check. Cflatau ---> Bruce add writting qtable, qtable2, mtable in this subroutine? C History: C 96.11.14 (PJF) added binary write options (see also info in C "versions.f" about IDL postprocessing) C possible "chackpointing" 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 end history C 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,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,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),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) C .. Local Scalars .. COMPLEX CXVAR1 REAL CX01R,DEGRAD,PHIND,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 If CNETFLAG = 'NCCLOS' then call WRITESCA to C call WRITENET to C call NCCLOS to close netCDF files 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,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. IF(IORI.GT.0)THEN C*** Compute g= DO 12 J=1,IORTH G(J)=QSCAG(1,J)/QSCAT(J) 12 CONTINUE CALL NAMER(IWAV,IRAD,IORI,CFLSCA,CFLAVG) OPEN(UNIT=8,FILE=CFLSCA,STATUS='UNKNOWN') WRITE(8,FMT=9030)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0 WRITE(8,FMT=9032)AEFF,WAVE,XX WRITE(8,FMT=9031)(CXRFR(ND),CXEPS(ND),ND,ND=1,NCOMP) 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 13 I=1,17 QAV(I)=-999. 13 CONTINUE QEXT(2)=-999. QABS(2)=-999. QSCAT(2)=-999. QBKSCA(2)=-999. QPHA(2)=-999. G(2)=-999. DO 14 J=1,3 QSCAG(J,2)=-999. 14 CONTINUE C endif "IF(IORTH.EQ.1)THEN" END IF 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 16 J=1,3 QAV(8+J)=.5*(QSCAG(J,1)+QSCAG(J,2)) 16 CONTINUE IF(CMDTRQ.EQ.'DOTORQ')THEN DO 17 J=1,3 QAV(11+J)=.5*(QTRQAB(J,1)+QTRQAB(J,2)) QAV(14+J)=.5*(QTRQSC(J,1)+QTRQSC(J,2)) 17 CONTINUE 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) IF(IORTH.EQ.2)THEN WRITE(8,FMT=9155)(QSCAG(J,2),J=1,3),(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 CFLATAU --->DRAINE is j=1 set for IORTH=1 ???? in QTRQAB, QTRQSC ?? C add missing value code (set to -999.) DO 19 I=1,3 DO 21 J=1,2 QTRQAB(I,J)=-999. QTRQSC(I,J)=-999. 21 CONTINUE 19 CONTINUE C endif "IF(CMDTRQ.EQ.'DOTORQ')THEN" ENDIF IF(IORTH.EQ.1)THEN WRITE(8,FMT=9060) DO 22 ND=1,NSCAT C convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) RVAR1=CONJG(CXF11(ND))*CXF11(ND) RVAR2=CONJG(CXF21(ND))*CXF21(ND) CXVAR1=CONJG(CXF11(ND))*CXF21(ND) WRITE(8,FMT=9070)ND,THETND,PHIND,RVAR1,RVAR2,CXVAR1 22 CONTINUE ELSE IF(IORTH.EQ.2)THEN WRITE(8,FMT=9080) DO 27 ND=1,NSCAT C convert scattering angles to degrees PHIND=DEGRAD*PHIN(ND) THETND=DEGRAD*THETAN(ND) C C Compute PPOL = degree of linear polarization of scattered light C for incident unpolarized light C PPOL=SQRT(SM(ND,2,1)**2+SM(ND,3,1)**2)/ & SM(ND,1,1) WRITE(8,FMT=9090)THETND,PHIND,SM(ND,1,1),SM(ND,2,1), & SM(ND,3,1),SM(ND,4,1),SM(ND,1,2),SM(ND,1,3),PPOL 27 CONTINUE C endif "IF(IORTH.EQ.1)THEN" ENDIF CLOSE(8) C ELSE C This module writes out results from orientational averaging. C (arrive here when called with IORI=0) C OPEN(UNIT=8,FILE=CFLAVG,STATUS='UNKNOWN') WRITE(8,FMT=9030)CSTAMP,CDESCR,CMDFFT,CALPHA,CSHAPE,NAT0 WRITE(8,FMT=9032)AEFF,WAVE,XX WRITE(8,FMT=9031)(CXRFR(ND),CXEPS(ND),ND,ND=1,NCOMP) 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 C 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 28 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) 28 CONTINUE C C*** Write orientationally-averaged Q values (except Qpha) to 'qtable': OPEN(UNIT=10,FILE='qtable',STATUS='OLD') DO 31 ND=1,ILIN10 READ(10,FMT=*) 31 CONTINUE 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 32 ND=1,ILIN12 READ(12,FMT=*) 32 CONTINUE WRITE(12,FMT=9057)AEFF,WAVE,QPHSUM(1) ILIN12=ILIN12+1 CLOSE(12) ELSE IF(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 38 J=1,3 QAV(8+J)=.5*(QSCGSUM(J,1)+QSCGSUM(J,2)) 38 CONTINUE IF(CMDTRQ.EQ.'DOTORQ')THEN DO 41 J=1,3 QAV(11+J)=.5*(QTRQABSUM(J,1)+QTRQABSUM(J,2)) QAV(14+J)=.5*(QTRQSCSUM(J,1)+QTRQSCSUM(J,2)) 41 CONTINUE 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) WRITE(8,FMT=9155)(QSCGSUM(J,2),J=1,3),(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 C*** Write orientationally-averaged Q values (except Qpha) to 'qtable': OPEN(UNIT=10,FILE='qtable',STATUS='OLD') DO 46 ND=1,ILIN10 READ(10,FMT=*) 46 CONTINUE 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 47 ND=1,ILIN12 READ(12,FMT=*) 47 CONTINUE WRITE(12,FMT=9057)AEFF,WAVE,QAV(6),QAV(7),QAV(8) ILIN12=ILIN12+1 CLOSE(12) C WRITE(8,FMT=9080) DO 48 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,SMORI(ND,1,1), & SMORI(ND,2,1),SMORI(ND,3,1),SMORI(ND,4,1), & SMORI(ND,1,2),SMORI(ND,1,3),PPOL 48 CONTINUE 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,WAVEA,AEFFA,THETA,BETA,PHI,A1,A2, & ICOMP,IXYZ0,IWAV,IRAD,IORI,AEFF,WAVE,XX,AK1, & BETAD,THETAD,PHID,AKR,CX01R,CXE02R,CXRFR,CXEPS, & QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB,QTRQSC, & SM,PHIN,THETAN,TIMERS) C endif "if(cbinflag.ne.'NOTCDF') then" 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,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) C endif "if(cnetflag.ne.'NOTCDF') then " ENDIF RETURN 9030 FORMAT (' DDSCAT --- ',A,/,' TARGET ---',A,/,' ',A, & ' --- method of solution ',/,' ',A, & ' --- prescription for polarizabilies',/,' ',A, & ' --- shape ',/,' NAT0 = ',I7,' = number of dipoles') 9031 FORMAT (' n= (',F7.4,' , ',F7.4,'), epsilon= (',F8.4,' , ',F7.4, & ') for material',I2) 9032 FORMAT (' AEFF=',F10.5,' =',' effective radius (physical units)', & /,' WAVE=',F10.5,' = wavelength (physical units)',/, & ' K*A=',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',6X,'Qabs',7X,'Qsca',5X,'g=',4X,'Qbk',7X, & 'Qpha',/,1X,'JO=1: ',1P,3E10.3,E11.3,2E10.3) 9055 FORMAT (1X,'JO=2: ',1P,3E10.3,E11.3,2E10.3,/,1X,'mean: ',3E10.3, & E11.3,2E10.3,/,1X,'Qpol= ',E10.3,35X,'dQpha=',E10.3) 9056 FORMAT (1P,E10.4,E11.4,3E10.3,E11.3,E10.3) 9057 FORMAT (1P,E10.4,E11.4,3E11.3) 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 **',/, & 'theta phi S_11 S_21 S_31 S_41', & ' S_12 S_13 Pol.') 9090 FORMAT (F5.1,F6.1,1P,E10.3,1P,5E11.3,0P,F9.5) 9150 FORMAT (8X,'Qsca*g(1) Qsca*g(2) Qsca*g(3)',/,1X,'JO=1:',1P, & 3E11.3) 9155 FORMAT (1X,'JO=2:',1P,3E11.3,/,1X,'mean:',1P,3E11.3) 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