SUBROUTINE INTERP(X,Y,Z,XA,YA,ZA,IDVOUT,MXTAB,NTAB) C Arguments: INTEGER IDVOUT,MXTAB,NTAB REAL X,Y,Z REAL XA(MXTAB),YA(MXTAB),ZA(MXTAB) C Local variables: INTEGER I1,I2,I3,I4,INC REAL SGN,X1,X2,X3,X4,Y1,Y2,Y3,Y4 C*********************************************************************** C Given array of tabulated values XA(1-NTAB), YA(1-NTAB), ZA(1-NTAB), C and given independent variable X, this routine interpolates for two C dependent variables: Y and Z. C Uses only parabolic interpolation for safety. C When called with X outside range of tabulated values XA, returns C extrapolated values Y and Z but issues warning statement. C Note: XA values must be either monotonically increasing OR C monotonically decreasing. C B.T.Draine, Princeton Univ. Observatory C 89/12/19 (BTD): Revised. C 91/05/02 (BTD): Added IDVOUT to argument list. C Changed WRITE(0 -> WRITE(IDVOUT C 92/03/24 (BTD): Added test IF(I2.GE.NTAB) in case INTERP is used for C data sets with different lengths. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DATA I2/1/ IF(I2.GE.NTAB)I2=NTAB-1 INC=0 SGN=1. C*** Check whether X is increasing or decreasing: IF(XA(1).GE.XA(NTAB))SGN=-1. C*** Check whether outside table limits IF((SGN*(X-XA(1))).LT.0.)THEN I2=2 WRITE(IDVOUT,6990)X GOTO 4600 ENDIF IF(SGN*(X-XA(NTAB)).GT.0.)THEN I2=NTAB-1 WRITE(IDVOUT,6990)X GOTO 4600 ENDIF C*** X is within table limits. Find I2 1000 IF(SGN*(XA(I2)-X))2000,3000,4000 2000 IF(INC)2100,2200,2200 2100 IF(I2+2-NTAB)4700,4700,2500 2200 INC=1 I2=I2+1 IF(I2+1.LE.NTAB)GO TO 1000 2500 I2=NTAB-1 GO TO 4600 3000 Y=YA(I2) Z=ZA(I2) RETURN 4000 IF(INC)4200,4200,4100 4100 I2=I2-1 IF(I2-2)4500,4700,4700 4200 INC=-1 I2=I2-1 IF(I2.GE.2)GOTO 1000 4500 I2=2 4600 I1=I2-1 I3=I2+1 X1=XA(I1) X2=XA(I2) X3=XA(I3) Y1=YA(I1) Y2=YA(I2) Y3=YA(I3) CALL PARAB3(X,Y,X1,X2,X3,Y1,Y2,Y3) Y1=ZA(I1) Y2=ZA(I2) Y3=ZA(I3) CALL PARAB3(X,Z,X1,X2,X3,Y1,Y2,Y3) RETURN 4700 I1=I2-1 I3=I2+1 I4=I2+2 X1=XA(I1) X2=XA(I2) X3=XA(I3) X4=XA(I4) Y1=YA(I1) Y2=YA(I2) Y3=YA(I3) Y4=YA(I4) CALL PARAB4(X,Y,X1,X2,X3,X4,Y1,Y2,Y3,Y4) Y1=ZA(I1) Y2=ZA(I2) Y3=ZA(I3) Y4=ZA(I4) CALL PARAB4(X,Z,X1,X2,X3,X4,Y1,Y2,Y3,Y4) RETURN 6990 FORMAT('Warning from INTERP: outside table limits for X=',1PE12.5) END SUBROUTINE PARAB3(X,Y,X1,X2,X3,Y1,Y2,Y3) REAL A,B,X,X1,X2,X3,Y,Y1,Y2,Y3 C*********************************************************************** C Subroutine PARAB3 does parabolic interpolation, with parabola C constrained to fit (x1,y1),(x2,y2),(x3,y3) exactly. C B.T.Draine, Institute for Advanced Study, March 1980. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** A=(Y3-Y2-(X3-X2)*(Y1-Y2)/(X1-X2))/((X3-X2)*(X3-X1)) B=(Y1-Y2)/(X1-X2)-A*(X1-X2) Y=(A*(X-X2)+B)*(X-X2)+Y2 RETURN END SUBROUTINE PARAB4(X,Y,X1,X2,X3,X4,Y1,Y2,Y3,Y4) REAL A,B,X,X1,X2,X3,X4,Y,Y1,Y2,Y3,Y4 C*********************************************************************** C Subroutine PARAB4 is designed to do parabolic interpolation, with C parabola constrained to match (x2,y2) and (x3,y3) exactly, and to C minimize sum of squared deviations from (x1,y1) and (x4,y4). C It is assumed that x1.lt.x2.le.x.lt.x3.lt.x4 C Fit parabola Y=A*(X-X2)**2+B*(X-X2)+Y2 C B.T.Draine, Institute for Advanced Study, March 1980. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. *********************************************************************** A=((X1-X2)*(Y3-Y2)/(X3-X2)+Y2-Y1)*(X1-X2)*(X1-X3) & +((X4-X2)*(Y3-Y2)/(X3-X2)+Y2-Y4)*(X4-X2)*(X4-X3) A=-A/(((X1-X2)*(X1-X3))**2+((X4-X2)*(X4-X3))**2) B=(Y3-Y2)/(X3-X2)-A*(X3-X2) Y=(A*(X-X2)+B)*(X-X2)+Y2 RETURN END SUBROUTINE NAMER(IWAV,IRAD,IDIR,CFLSCA,CFLAVG) C*********************************************************************** C Purpose: to generate file names for specific (wavelength, C size, C incident direction) C Present version allows up to 100 wavelengths, (00-99) C 100 sizes (00-99) C 1000 directions (000-999) C Given: C IWAV=integer (0-99) identifying wavelength C IRAD=integer (0-99) identifying size C IDIR=integer (0-999) identifying direction C Returns: C CFLSCA=name for output file containing Qext,Qabs,Qpha,Qsca,g C and fml values for selected directions C CFLAVG=name for output file for given wavelength and size, C containing scattering etc. averaged over inc.dir. C C B.T.Draine, Princeton Univ. Observatory, 1988 C History: C 90/11/21 (BTD): Modified to allow 3 digits to identify direction C and to change indices to begin from w00r00k000.sca C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments: INTEGER IDIR,IDIR0,IRAD,IRAD0,IWAV,IWAV0 CHARACTER CFLSCA*14,CFLAVG*13 C C Local variables: CHARACTER A1*1,A2*1,A3*1,EXT*4,ZSCA*14,ZAVG*13 CHARACTER N(10)*1,V1(2)*1,V2(2)*1,V3(3)*1,W(14)*1 C C Equivalences: EQUIVALENCE (A1,W(1)),(V1,W(2)),(A2,W(4)),(V2,W(5)),(A3,W(7)), $ (V3,W(8)),(EXT,W(11)) EQUIVALENCE (W,ZSCA) EQUIVALENCE (W,ZAVG) C C DATA statements: DATA N/'0','1','2','3','4','5','6','7','8','9'/ C*********************************************************************** C*** Filename CFLSCA will be of the form w01r01k001 C (for IWAV=2,IRAD=2,IDIR=2) C A1='w' A2='r' A3='k' C*** Set IWAV0=IWAV-1 to run over range 0-99 C Set IRAD0=IRAD-1 to run over range 0-99 C Set IDIR0=IDIR-1 to run over range 0-999 IWAV0=IWAV-1 IRAD0=IRAD-1 IDIR0=IDIR-1 V1(1)=N(IWAV0/10+1) V1(2)=N(IWAV0-10*(IWAV0/10)+1) V2(1)=N(IRAD0/10+1) V2(2)=N(IRAD0-10*(IRAD0/10)+1) V3(1)=N(IDIR0/100+1) V3(2)=N((IDIR0-100*(IDIR0/100))/10+1) V3(3)=N(IDIR0-10*(IDIR0/10)+1) EXT='.sca' CFLSCA=ZSCA C*** CFLAVG is of slightly different form, C e.g., w00r00ori.avg A3='o' V3(1)='r' V3(2)='i' V3(3)='.' EXT='avg ' CFLAVG=ZAVG RETURN C END SUBROUTINE NULLER(CXVEC,IOCC,MXNAT,MXN3,NAT) C Arguments: INTEGER MXNAT,MXN3,NAT INTEGER*2 IOCC(MXNAT) COMPLEX*8 CXVEC(NAT,3) C Local variables: INTEGER J,L C*********************************************************************** C This routine takes as input a complex vector with NAT3 elements C and sets to zero those elements corresponding to "vacuum" sites. C This is accomplished by multiplying by vector IOCC(J), whose C elements are either 1 or 0 depending on whether site is occupied C or "vacuum". C Note: CXVEC should be given dimension CXVEC(NAT,3) here C so that first 3*NAT elements of CXVEC are employed. C C B.T.Draine, Princeton Univ. Obs., 90/11/1 C History: C 90/12/15 (BTD): Corrected error in dimensioning of CXVEC. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DO 110 L=1,3 DO 100 J=1,NAT CXVEC(J,L)=IOCC(J)*CXVEC(J,L) 100 CONTINUE 110 CONTINUE RETURN END SUBROUTINE ORIENT(BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX, & MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI, & BETA,THETA,PHI,WGTA,WGTB) C Arguments: INTEGER MXBETA,MXTHET,MXPHI,NBETA,NTHETA,NPHI REAL BETAMI,BETAMX,THETMI,THETMX,PHIMIN,PHIMAX REAL BETA(MXBETA),THETA(MXTHET),PHI(MXPHI), & WGTA(MXTHET,MXPHI),WGTB(MXBETA) C Local variables: INTEGER I,J REAL DELTA,WG C*********************************************************************** C Given: C BETAMI=minimum value of beta (radians) C BETAMX=maximum value of beta (radians) C THETMI=minimum value of theta (radians) C THETMX=maximum value of theta (radians) C PHIMIN=minimum value of phi (radians) C PHIMAX=maximum value of phi (radians) C MXBETA,MXTHET,MXPHI=parameters for dimensioning of arrays C BETA,THETA,PHI C NBETA=desired number of values of beta C NTHETA=desired number of values of theta C NPHI=desired number of values of PHI C Returns: C BETA(1-NBETA)=beta values (radians) C THETA(1-NTHETA)=theta values (radians) C PHI(1-NPHI)=phi values (radians) C WGTA(1-NTHETA,1-NPHI)=weighting of each orientation of target C axis a1 in Lab Frame C (sum of weights = 1) C WGTB(1-NBETA)=weighting of each rotation of target around a1 C (sum of weights = 1) C Note: it is assumed that target orientation weight function C can be factored into WGTA*WGTB -- i.e., that rotations C around a1 are decoupled from orientation of a1. C C Purpose: to generate list of desired target orientations C C Definitions: beta=angle of rotation of target around target axis a1 C (beta=0 is defined to be such that a2 lies in a1-x C plane, with a2 in direction of increasing theta) C theta=angle between target axis a1 and x axis C (0.le.theta.le.pi) C phi=angle of rotation of axis a1 around axis C (phi=0 is defined to be such that a1 lies in xy plane) C Present version assumes: C beta to be uniformly distributed between BETAMI and BETAMX C cos(theta) to be uniformly distributed between cos(THETMI) and C cos(THETMX) C phi to be uniformly distributed between PHIMIN and PHIMAX C C Values assigned to BETA and PHI are midpoints of uniform intervals in C beta and phi C Note that if NBETA=1, BETA is set to midpoint of range in phi C Likewise, if NPHI=1, PHI is set to midpoint of range in phi C C If NTHETA=1: specify a single value of THETA C cos(theta)=0.5*(cos(thetami)+cos(thetamx)) C C If NTHETA>1: C If NTHETA is even: C divide range of cos(theta) into NTHETA equal intervals C set THETA to midpoints of these intervals C give intervals equal weights C If NTHETA is odd: C divide range of cos(theta) into NTHETA-1 equal intervals C set THETA to endpoints of these intervals; first value of C THETA is THETAMI and last value of THETA is THETAMX C C Note: values chosen for beta, phi are at midpoints of C uniform intervals C C B.T.Draine, Princeton Univ. Obs., 89.11.20 C History: C 92.04.01 (BTD) Modify to allow Simpson's rule or trapezoidal rule C for integration over cos(theta). C 99.03.01 (BTD) Modified to properly assign THETA values when NPHI>1 C Code up to this time was not consistent with C description in UserGuide. This problem was brought C to my attention by Miroslav Kocifaj C (Miroslav.Kocifaj@swh.sk). C 99.03.05 (BTD) Corrected errors in assignment of THETA values. C Thanks to Miroslav Kocifaj for detecting these. C 99.04.30 (BTD) Corrected errors in assignment of BETA and C PHI values when BETAMI.ne.0 or PHIMIN.ne.0 C These errors (inadvertent omission of BETAMI and C PHIMIN from expressions assigning BETA(J) and PHI(J)) C were introduced in 99.03.05 revision. C Thanks to Henriette Lemke for calling attention to C these errors. C Copyright (C) 1993,1999 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C Assign values to BETA: C DELTA=(BETAMX-BETAMI)/REAL(NBETA) DO J=1,NBETA BETA(J)=BETAMI+DELTA*(REAL(J)-0.5) ENDDO C C Assign values to THETA: C IF(NTHETA.EQ.2*(NTHETA/2))THEN C C THETA is even: C DELTA=(COS(THETMX)-COS(THETMI))/REAL(NTHETA) DO J=1,NTHETA THETA(J)=ACOS(COS(THETMI)+DELTA*(REAL(J)-0.5)) ENDDO ELSE C C THETA is odd: C IF(NTHETA.EQ.1)THEN THETA(1)=ACOS(0.5*(COS(THETMI)+COS(THETMX))) ELSE DELTA=(COS(THETMX)-COS(THETMI))/REAL(NTHETA-1) THETA(1)=THETMI THETA(NTHETA)=THETMX DO J=2,NTHETA-1 THETA(J)=ACOS(COS(THETMI)+DELTA*REAL(J-1)) ENDDO ENDIF ENDIF C C Assign values to PHI: C DELTA=(PHIMAX-PHIMIN)/REAL(NPHI) DO J=1,NPHI PHI(J)=PHIMIN+DELTA*(REAL(J)-0.5) ENDDO C C** Specify weight factors WGTA, WGTB C (Note: weight function WGTA = 4*pi*P/(NTHETA*NPHI), where C P=(probability/solid angle) of orientation in direction THETA,PHI . C The orientational averaging program automatically samples uniformly C in cos(theta) to allow for d(solid angle)=sin(theta)d(theta)d(phi). C C Present version assumes random orientations. C C When NTHETA is >1 and even, we use trapezoidal integration C When NTHETA is >1 and odd, we use Simpson's rule integration. C IF(NTHETA.EQ.1.OR.NTHETA.EQ.2*(NTHETA/2))THEN C C either 1 or even number of theta values: midpoints of intervals C WG=1./REAL(NTHETA*NPHI) DO J=1,NPHI DO I=1,NTHETA WGTA(I,J)=WG ENDDO ENDDO ELSE C C odd number >1 of theta values: use Simpson's rule weighting C WG=1./REAL(3*(NTHETA-1)*NPHI) DO J=1,NPHI WGTA(1,J)=WG WGTA(NTHETA,J)=WG ENDDO WG=4./REAL(3*(NTHETA-1)*NPHI) DO I=2,NTHETA,2 DO J=1,NPHI WGTA(I,J)=WG ENDDO ENDDO IF(NTHETA.GE.5)THEN WG=2./REAL(3*(NTHETA-1)*NPHI) DO I=3,NTHETA-1,2 DO J=1,NPHI WGTA(I,J)=WG ENDDO ENDDO ENDIF ENDIF C C Assign weights for rotation angle beta: C WG=1./REAL(NBETA) DO I=1,NBETA WGTB(I)=WG ENDDO C RETURN END SUBROUTINE REAPAR(cbinflag, cnetflag, AEFFA, & BETAMI,BETAMX, & CALPHA,CDIEL,CFLEPS,CFLPAR,CMDSOL,CMDFFT,CSHAPE, & CMDTRQ,CXE01,CXE02, & DEGRAD, & ICTHM,IDVERR,IDVOUT,INIT,IOPAR,IORTH,IPHM, & IWRKSC, & MXBETA,MXCOMP,MXMEM,MXPHI,MXRAD,MXSCA,MXTHET, & MXWAV, & NBETA,NCOMP,NPHI,NRAD,NSCA,NTHETA,NWAV, & PHIMIN,PHIMAX,PHIN, & THETAN,THETMI,THETMX,TOL, & WAVEA,SHPAR) C*********************************************************************** C Subroutine REAPAR handles the reading of input parameters from the C "ddscat.par" file, as well as elementary processing with those input C parameters to generate arrays. C Original version created by P.J.Flatau, Colorado State Univ. C Modified by B.T.Draine, Princeton Univ. Obs. C History: C 90.11.09 (BTD): Modified so as to remove sensitivity to errors C in sign of entered value of DTHETA. C 90.11.30 (BTD): Added test to compare CMDFFT and MXMEM C 90.12.14 (BTD): Changed option MKRECT to RCTNGL C 91.01.03 (BTD): Added MXBETA,MXPHI,MXTHET to argument list C and added checking of NBETA,NPHI,NTHETA C 91.05.08 (BTD): Added CALPHA to allow choice of prescription C for polarizabilities C 91.05.08 (BTD): Some rewriting done. C 91.05.14 (BTD): Added separate options LDRXYZ and LDRAVG C 91.05.15 (BTD): Added new option LDR000 C 91.05.23 (BTD): Added variable IWRKSC to argument list and C added code to read IWRKSC C 91.05.30 (BTD): Changed LDRXYZ->LDR100 C Added option LDR111 C 91.11.13 (BTD): Now three options: DRAI88, GOBR88, LATTDR C (subroutine ALPHA now allows for directional C dependence for LATTDR option) C 92.05.14 (BTD): Allow for experimental option LDRISO C 92.09.21 (BTD): Allow for option ANIELL (ellipsoid of anisotropic C material); change option GRPHIT -> UNICYL (cylinder C of uniaxial material) C 93.01.07 (BTD): Allow for options TWOELL and TWOAEL (two touching C ellipsoids, each either isotropic or anisotropic) C 93.03.11 (BTD): Deleted all code associated with variable CMACHN C which previously was used to identify machine/OS C This information was not used or needed. C 93.12.15 (BTD): Added code to check for compatibility of NCOMP C when CSHAPE=TWOELL,THRELL,TWOAEL,or THRAEL C 94.01.27 (BTD): Replaced SHPAR1,SHPAR2,SHPAR3 by SHPAR(1-6) C Added CONELL (two concentric ellipsoids) as target C option. C 94.06.20 (PJF): Added code to handle NEWTMP C 94.12.02 (BTD): Added options C LB1991 (Dungey-Bohren 1991 polarizability) C LR1994 (Lumme-Rahola 1994 polarizability) C 95.05.26 (BTD): Corrected construction of orthogonal polarization C state for case where E01 is complex (e.g., C circular polarization) C 95.06.14 (BTD): Added code to check that user is not requesting C f_ml when illuminating target with other then C linearly polarized light C 95.06.19 (BTD): Added argument CMDSOL to allow selection of C different CCG methods C Added argument CMDTRQ to allow torque calculations C to be skipped. C 96.01.05 (BTD): Added shape options BLOCKS and DW1996 C 96.01.25 (BTD): Added shape option TWOSPH C 96.01.25 (BTD): Disable code added 95.06.14 so that user can C calculate f_ml even when illuminating target with C elliptically polarized light. C 96.10.18 (BTD): Changed option NEWTMP to GPFAFT (Generalized Prime C Factor Algorithm for Fourier Transform). C 96.11.06 (BTD): Added provision to C (1) produce fatal error if CXE01(1).NE.(0.,0.) C (2) normalize CXE01 if necessary C 96.11.14 (PJF): cbinflag for binary file C 96.11.15 (PJF): netCDF netflag C 96.11.21 (BTD): modified handling of CBINFLAG and CNETFLAG to C be consistent with other option flags C 98.12.21 (BTD): changed dimension of CFLPAR from CHARACTER*40 to C CHARACTER*60 to allow longer file names C (also changed in DDSCAT.f and dielec.f) C end history 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*********************************************************************** C C Arguments: CHARACTER & CALPHA*(*),CDIEL*(*),CFLPAR*(*),CMDSOL*(*),CMDFFT*(*), & CMDTRQ*(*),CSHAPE*(*), cbinflag*(*), cnetflag*(*) INTEGER & ICTHM,IDVERR,IDVOUT,INIT,IOPAR,IORTH,IPHM,IWRKSC, & MXBETA,MXCOMP,MXMEM,MXPHI,MXRAD,MXSCA,MXTHET,MXWAV, & NBETA,NCOMP,NPHI,NRAD,NSCA,NWAV,NTHETA CHARACTER*60 CFLEPS(MXCOMP) REAL & BETAMI,BETAMX,DEGRAD,PHIMIN,PHIMAX,THETMI,THETMX,TOL REAL & AEFFA(MXRAD),PHIN(MXSCA),SHPAR(6),THETAN(MXSCA),WAVEA(MXWAV) COMPLEX CXE01(3),CXE02(3) C C Local variables: CHARACTER CDIVID*3,CLINE*70,CMSGNM*70 INTEGER J,NSCA0 REAL & AEFEND,AEFINI,DELTA,DTHETA,E1,PHI1,THETA1,THETA2, & WAVEND,WAVINI C C External routines: EXTERNAL DIVIDE,ERRMSG,WRIMSG C C*********************************************************************** OPEN (UNIT=IOPAR,FILE=CFLPAR,STATUS='OLD') C*********************************************************************** C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG('REAPAR',CLINE) READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) C C Specify whether torques are to be calculated: C C CMDTRQ*6 = 'DOTORQ' - calculate torque on grain C = 'NOTORQ' - do not calculate torque on grain C READ(IOPAR,FMT=*,ERR=40)CMDTRQ IF(CMDTRQ.EQ.'dotorq')CMDTRQ='DOTORQ' IF(CMDTRQ.EQ.'notorq')CMDTRQ='NOTORQ' IF(CMDTRQ.EQ.'DOTORQ')THEN WRITE(CMSGNM,FMT='(A,A)')CMDTRQ,' - compute torques ' CALL WRIMSG('REAPAR',CMSGNM) ELSEIF(CMDTRQ.EQ.'NOTORQ')THEN WRITE(CMSGNM,FMT='(A,A)')CMDTRQ,' - do not compute torques ' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMDTRQ') ENDIF C C Define method used for iterative solution of complex linear equations C C CMDSOL*6 = 'PETRKP' - Petravic & Kuo-Petravic method C = 'PBCGST' - PIM BiConjugate Gradient with Stabilization C READ(IOPAR,FMT=*,ERR=40)CMDSOL IF(CMDSOL.EQ.'petrkp')CMDSOL='PETRKP' IF(CMDSOL.EQ.'pbcgst')CMDSOL='PBCGST' IF(CMDSOL.EQ.'PETRKP'.OR. & CMDSOL.EQ.'PBCGST')THEN WRITE(CMSGNM,FMT='(A,A)')CMDSOL,' - CCG Method ' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMDSOL') ENDIF C C Define FFT method: C C CMDFFT*6 = 'BRENNR' - Use Brenner's FOURX FFT routine. C = 'TMPRTN' - Use Temperton's FFT routine. C = 'GPFAFT' - GPFA code of Temperton C = 'CONVEX' - Use native Convex veclib routine C3DFFT. C READ(IOPAR,FMT=*,ERR=40)CMDFFT IF(CMDFFT.EQ.'BRENNR'.OR. & CMDFFT.EQ.'TMPRTN'.OR. & CMDFFT.EQ.'CONVEX'.OR. & CMDFFT.EQ.'GPFAFT')THEN WRITE(CMSGNM,FMT='(A,A)')CMDFFT,' - Solution Method ' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMDFFT') ENDIF C*** Check whether compiled with enough memory for TMPRTN option: IF(CMDFFT.EQ.'TMPRTN'.AND.MXMEM.EQ.0)THEN CALL ERRMSG('FATAL','REAPAR', & ' Must recompile with MXMEM=1 to use TMPRTN option') ENDIF IF(CMDFFT.NE.'TMPRTN'.AND.MXMEM.EQ.1)THEN CALL ERRMSG('WARNING','REAPAR', & ' You do not need MXMEM=1 with BRENNR,CONVEX, or GPFAFT') ENDIF C*********************************************************************** C C Define prescription for computing polarizabilities: C C CALPHA*6 = 'DRAI88' - Draine (1988) prescription C = 'GOBR88' - Goedecke & Obrien (1988) prescription C = 'LATTDR' - Lattice Dispersion Relation C = 'LDRISO' - Isotropized Lattice Dispersion Relation C = 'DB1991' - Dungey & Bohren (1991) prescription C = 'LR1994' - Lumme & Rahola (1994) prescription C READ(IOPAR,FMT=*,ERR=40)CALPHA IF(CALPHA.EQ.'DRAI88')THEN CALL WRIMSG('REAPAR', & 'DRAI88 - Draine (1988) prescription for alpha') ELSEIF(CALPHA.EQ.'GOBR88')THEN CALL WRIMSG('REAPAR', & 'GOBR88 - Goedecke & OBrien (1988) prescription for alpha') ELSEIF(CALPHA.EQ.'LATTDR')THEN CALL WRIMSG('REAPAR', & 'LATTDR - Lattice Dispersion Relation for alpha') ELSEIF(CALPHA.EQ.'LDRISO')THEN CALL WRIMSG('REAPAR', & 'LATTDR - Isotropized Lattice Dispersion Relation for alpha') ELSEIF(CALPHA.EQ.'DB1991')THEN CALL WRIMSG('REAPAR', & 'DB1991 - Dungey & Bohren (1991) prescription for alpha') ELSEIF(CALPHA.EQ.'LR1994')THEN CALL WRIMSG('REAPAR', & 'LR1994 - Lumme & Rahola (1994) prescription for alpha') ELSE CALL ERRMSG('FATAL','REAPAR',' wrong definition of CALPHA') ENDIF C C binary file flag READ(IOPAR,FMT=*,ERR=40)CBINFLAG IF(CBINFLAG.EQ.'allbin')CBINFLAG='ALLBIN' IF(CBINFLAG.EQ.'oribin')CBINFLAG='ORIBIN' IF(CBINFLAG.EQ.'notbin')CBINFLAG='NOTBIN' IF(CBINFLAG.EQ.'ALLBIN' & .OR.CBINFLAG.EQ.'ORIBIN' & .OR.CBINFLAG.EQ.'NOTBIN')THEN WRITE(CMSGNM,FMT='(A,A)')CBINFLAG, & ' - Unformatted binary dump option' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' Wrong definition of CBINFLAG') ENDIF C C binary file flag READ(IOPAR,FMT=*,ERR=40)CNETFLAG IF(CNETFLAG.EQ.'allcdf')CNETFLAG='ALLCDF' IF(CNETFLAG.EQ.'oricdf')CNETFLAG='ORICDF' IF(CNETFLAG.EQ.'notcdf')CNETFLAG='NOTCDF' IF(CNETFLAG.EQ.'ALLCDF' & .OR.CNETFLAG.EQ.'ORICDF' & .OR.CNETFLAG.EQ.'NOTCDF')THEN WRITE(CMSGNM,FMT='(A,A)')CNETFLAG, & ' - netCDF file option' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' Wrong definition of CNETFLAG') ENDIF C C*********************************************************************** C C Define shape: C C CSHAPE*6 = 'FRMFIL' read shape data from file C = 'RCTNGL' make rectangular solid C = 'ELLIPS' make ellipsoidal solid C = 'CYLNDR' make cylindrical solid C = 'HEXGON' make hexagonal prism C = 'TETRAH' make regular tetrahedron C = 'UNICYL' make cylinder of unixaxial material C = 'ANIELL' make ellipsoid of anisotropic material C = 'TWOELL' make two touching ellipsoids of materials 1,2 C = 'TWOAEL' make two touching anisotropic ellipsoids of C materials 1-6 C = 'THRELL' make three touching isotropic ellipsoids of C materials 1,2,3 C = 'THRAEL' make three touching anisotropic ellipsoids of C materials 1-9 C = 'CONELL' make two concentric ellipsoids of C materials 1,2 C = 'BLOCKS' make target out of cubic blocks defined by C data in file 'blocks.par' C = 'DW1996' make 13-cube target used by Draine & C Weingartner 1996 C = 'TWOSPH' make two touching spheroids with symmetry axes at C specified angle C READ(IOPAR,FMT=*,ERR=40)CSHAPE IF( CSHAPE.EQ.'FRMFIL' & .OR.CSHAPE.EQ.'RCTNGL' & .OR.CSHAPE.EQ.'ELLIPS' & .OR.CSHAPE.EQ.'CYLNDR' & .OR.CSHAPE.EQ.'HEXGON' & .OR.CSHAPE.EQ.'TETRAH' & .OR.CSHAPE.EQ.'UNICYL' & .OR.CSHAPE.EQ.'ANIELL' & .OR.CSHAPE.EQ.'TWOELL' & .OR.CSHAPE.EQ.'TWOAEL' & .OR.CSHAPE.EQ.'THRELL' & .OR.CSHAPE.EQ.'THRAEL' & .OR.CSHAPE.EQ.'CONELL' & .OR.CSHAPE.EQ.'BLOCKS' & .OR.CSHAPE.EQ.'DW1996' & .OR.CSHAPE.EQ.'TWOSPH')THEN WRITE(CMSGNM,FMT='(A,A)')CSHAPE,' - Shape definition ' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' Wrong shape directive ') ENDIF C C Read shape parameters C Need 3 shape parameters for all options except C FRMFIL (here skip a line, later read target from file) C TETRAH, DW1996 (1 shape parameter) C CYLNDR, HEXGON, UNICYL (2 shape parameters) C TWOSPH (5 shape parameters) C CONELL (6 shape parameters) C IF(CSHAPE.EQ.'FRMFIL')THEN READ(IOPAR,FMT=*,ERR=40) ELSEIF( & CSHAPE.EQ.'TETRAH'.OR. & CSHAPE.EQ.'DW1996')THEN READ(IOPAR,FMT=*,ERR=40)SHPAR(1) ELSEIF( & CSHAPE.EQ.'CYLNDR'.OR. & CSHAPE.EQ.'HEXGON'.OR. & CSHAPE.EQ.'UNICYL')THEN READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2) ELSEIF(CSHAPE.EQ.'TWOSPH')THEN READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2),SHPAR(3), & SHPAR(4),SHPAR(5) ELSEIF(CSHAPE.EQ.'CONELL')THEN READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2),SHPAR(3), & SHPAR(4),SHPAR(5),SHPAR(6) ELSE READ(IOPAR,FMT=*,ERR=40)SHPAR(1),SHPAR(2),SHPAR(3) ENDIF C*********************************************************************** C C Obtain names of file(s) containing dielectric function(s) C NCOMP = number of different dielectric functions C CFLEPS(1-NCOMP) = names of files containing dielectric functions C READ(IOPAR,FMT=*,ERR=40)NCOMP C C*** Check that MXCOMP is large enough ********************************* IF(NCOMP.GT.MXCOMP) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP > MXCOMP: must increase MXCOMP in DDSCAT') C C*** Check that NCOMP=2 if CSHAPE=UNICYL ******************************* C 3 ANIELL C 2 TWOELL C 3 THRELL C 6 TWOAEL C 9 THRAEL C 2 CONELL C 2 TWOSPH IF(CSHAPE.EQ.'UNICYL'.AND.NCOMP.NE.2) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 2 for option UNICYL') IF(CSHAPE.EQ.'ANIELL'.AND.NCOMP.NE.3) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 3 for option ANIELL') IF(CSHAPE.EQ.'TWOELL'.AND.NCOMP.NE.2) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 2 for option TWOELL') IF(CSHAPE.EQ.'THRELL'.AND.NCOMP.NE.3) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 3 for option THRELL') IF(CSHAPE.EQ.'TWOAEL'.AND.NCOMP.NE.6) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 6 for option TWOAEL') IF(CSHAPE.EQ.'THRAEL'.AND.NCOMP.NE.9) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 9 for option THRAEL') IF(CSHAPE.EQ.'CONELL'.AND.NCOMP.NE.2) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 2 for option CONELL') IF(CSHAPE.EQ.'TWOSPH'.AND.NCOMP.NE.2) & CALL ERRMSG('FATAL','REAPAR', & 'NCOMP must be 2 for option TWOSPH') C*********************************************************************** READ(IOPAR,FMT=*,ERR=40)CDIEL IF(CDIEL.EQ.'TABLES')THEN DO 1000 J=1,NCOMP READ(IOPAR,FMT=*,ERR=40)CFLEPS(J) 1000 CONTINUE ELSEIF(CDIEL.EQ.'H2OICE'.OR.CDIEL.EQ.'H2OLIQ')THEN IF(NCOMP.GT.1)CALL ERRMSG('FATAL','REAPAR', & 'H2OICE/LIQ with NCOMP>1?') ELSE CALL ERRMSG('FATAL','REAPAR','Unrecognized option for CDIEL') ENDIF C*********************************************************************** C C Define INITIALIZATION: C C INIT = 0 to start with |X0>=0 C 1 to obtain |X0> from 4th-order expansion in polarizability C 2 to read |X0> from file solvp.in C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) READ(IOPAR,FMT=*,ERR=40)INIT IF(INIT.EQ.0)THEN CALL WRIMSG('REAPAR', & 'INIT=0 to start with |X0>=0 (CCG method)') ELSEIF(INIT.EQ.1)THEN CALL WRIMSG('REAPAR', & 'INIT=1 to obtain |X0> from 4th-order expansion (CCG)') ELSEIF(INIT.EQ.2)THEN CALL WRIMSG('REAPAR', & 'INIT=2 to read |X0> from file solvp.in') ELSE CALL ERRMSG('FATAL','REAPAR',' Wrong value of INIT') ENDIF C*********************************************************************** C C Define error tolerance: C C TOL= maximum acceptable value of |Ax-E|/|E| C READ(IOPAR,FMT=*,ERR=40)TOL IF(TOL.GT.1.E-15.AND.TOL.LT.1.E-1)THEN WRITE(CMSGNM,FMT='(1P,E10.3,A)')TOL, & ' = TOL = max. acceptable normalized residue |Ax-E|/|E|' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('WARNING','REAPAR',' Strange value of TOL ') ENDIF C*********************************************************************** C C Define NUMBER OF THETA AND PHI VALUES FOR EVALUATION OF CSCA AND G C C ICTHM = number of thetan values to use in evaluation of CSCA, G C (must be odd; recommend ICTHM=33) C IPHM = number of phin values to use in evaluation of CSCA, G C (even or odd; recommend IPHM=12 or more) C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) READ(IOPAR,FMT=*)ICTHM IF((ICTHM+1)/2*2 .NE. (ICTHM+1))THEN CALL ERRMSG('FATAL','REAPAR',' ICTHM SHOULD BE ODD ') ELSE WRITE(CMSGNM,FMT='(I4,A)')ICTHM, &' - ICTHM; number of THETAN values for calc.of g and Qsca by SCAT' CALL WRIMSG('REAPAR',CMSGNM) ENDIF READ(IOPAR,FMT=*)IPHM WRITE(CMSGNM,FMT='(I4,A)')IPHM, &' - IPHM; number of PHIN values for calc. of g and Qsca by SCAT' CALL WRIMSG('REAPAR',CMSGNM) C*********************************************************************** C C Define WAVELENGTH : C C WAVINI = first wavelength (physical units) C WAVEND = last wavelength (physical units) C NWAV = number of wavelengths C CDIVID = 'LIN', 'LOG', or 'INV' for equal increments in C wavelength, log(wavelength), or frequency C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) READ(IOPAR,FMT=*,ERR=40)WAVINI,WAVEND,NWAV,CDIVID WRITE(CMSGNM,FMT='(I3,A,F7.4,A,F7.4)') & NWAV,' wavelengths from ',WAVINI,' to ',WAVEND CALL WRIMSG('REAPAR',CMSGNM) CALL DIVIDE(CDIVID,WAVINI,WAVEND,NWAV,MXWAV,WAVEA) C*********************************************************************** C C Define AEFF = radius of equal-volume sphere (physical units) C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) READ(IOPAR,FMT=*)AEFINI,AEFEND,NRAD,CDIVID WRITE(CMSGNM,FMT='(I3,A,F7.4,A,F7.4)') & NRAD,' eff. radii from ',AEFINI,' to ',AEFEND CALL WRIMSG('REAPAR',CMSGNM) CALL DIVIDE(CDIVID,AEFINI,AEFEND,NRAD,MXRAD,AEFFA) C*********************************************************************** C C Define incident polarizations (in Lab Frame) C It is assumed that incident radiation is along x axis in Lab Frame C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) C Read complex polarization vector CXE01=e01 (normalize if necessary). READ(IOPAR,FMT=*,ERR=40)CXE01 E1=SQRT(CXE01(1)*CONJG(CXE01(1))) IF(E1.NE.0.)THEN CALL ERRMSG('FATAL','REAPAR',' CXE01(1) must be zero!') ENDIF C Normalize: E1=SQRT(CXE01(2)*CONJG(CXE01(2))+CXE01(3)*CONJG(CXE01(3))) CXE01(2)=CXE01(2)/E1 CXE01(3)=CXE01(3)/E1 C Construct orthogonal normalized polarization vector CXE02=e02 C using xhat cross e01* = e02 C CXE02(1)=(0.,0.) CXE02(2)=-CONJG(CXE01(3)) CXE02(3)=CONJG(CXE01(2)) C C IORTH=1 to calculate only for single polarization C =2 to also calculate for orthogonal polarization READ(IOPAR,FMT=*,ERR=40)IORTH C C*********************************************************************** C C Specify whether or not to write ".sca" files C C IWRKSC=0 to NOT write ".sca" file for each target orientation C =1 to write ".sca" file for each target orientation READ(IOPAR,FMT=*,ERR=40)IWRKSC C C*********************************************************************** C C Read information determining target rotations C READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) READ(IOPAR,FMT=*)BETAMI,BETAMX,NBETA READ(IOPAR,FMT=*)THETMI,THETMX,NTHETA READ(IOPAR,FMT=*)PHIMIN,PHIMAX,NPHI C C*** Check that values of NBETA,NTHETA,NPHI are OK: IF(NBETA.GT.MXBETA) & CALL ERRMSG('FATAL','REAPAR',' NBETA>MXBETA: must recompile') IF(NPHI.GT.MXPHI) & CALL ERRMSG('FATAL','REAPAR',' NPHI>MXPHI: must recompile') IF(NTHETA.GT.MXTHET) & CALL ERRMSG('FATAL','REAPAR',' NTHETA>MXTHET: must recompile') C C If NPHI>1, then set IORTH=2 regardless of value input. IF(IORTH.EQ.1.AND.NPHI.GT.1)THEN IORTH=2 CALL WRIMSG('REAPAR', & 'Set IORTH=2 since NPHI>1 ') ENDIF C IF(IORTH.EQ.1)THEN CALL WRIMSG('REAPAR', & 'Calculate only for single polarization ') ELSEIF(IORTH.EQ.2)THEN CALL WRIMSG('REAPAR', & 'Do orthogonal polarization for each case ') ELSE CALL ERRMSG('FATAL','REAPAR',' WRONG VALUE OF IORTH ') ENDIF WRITE(CMSGNM,FMT='(2F7.2,A,I4)')BETAMI,BETAMX, & ' Range of BETA values ; NBETA =',NBETA CALL WRIMSG('REAPAR',CMSGNM) WRITE(CMSGNM,FMT='(2F7.2,A,I4)')THETMI,THETMX, & ' Range of THETA values; NTHETA=',NTHETA CALL WRIMSG('REAPAR',CMSGNM) WRITE(CMSGNM,FMT='(2F7.2,A,I4)')PHIMIN,PHIMAX, & ' Range of PHI values ; NPHI =',NPHI CALL WRIMSG('REAPAR',CMSGNM) C Convert from degrees to radians BETAMI=BETAMI/DEGRAD BETAMX=BETAMX/DEGRAD THETMI=THETMI/DEGRAD THETMX=THETMX/DEGRAD PHIMIN=PHIMIN/DEGRAD PHIMAX=PHIMAX/DEGRAD C*********************************************************************** C C Specify scattering directions to be calculated C THETAN,PHIN = Direction of n from n0 in calculation of the C scattering matrix; cos(THETAN) is n0 \dot n C and PHIN is azimuthal angle of n from Lab xy plane C Note: THETA1, THETA2, and PHI for each scattering C plane are entered in degrees, and immediately C converted to radians. Arbitrary number C of scattering planes may be considered. C NSCA=0 READ(IOPAR,FMT=*,ERR=40)CLINE CALL WRIMSG(' ',CLINE) 22 READ(IOPAR,FMT=*,ERR=40,END=25)PHI1,THETA1,THETA2,DTHETA WRITE(CMSGNM,FMT='(3F7.1,A)')PHI1,THETA1,THETA2, & ' = phi, theta_min, theta_max for scattering plane' CALL WRIMSG('REAPAR',CMSGNM) C Convert to radians: PHI1=PHI1/DEGRAD THETA1=THETA1/DEGRAD THETA2=THETA2/DEGRAD DTHETA=DTHETA/DEGRAD C Allow for possibility that user did not enter correct sign for C DTHETA (if I did it, others will too...) IF(THETA2.LT.THETA1)DTHETA=-ABS(DTHETA) C compute theta values for this scattering plane NSCA0=1 DELTA=THETA2-THETA1 IF(DELTA.NE.0)NSCA0=1+INT(DELTA/DTHETA+0.5) C Check that there is enough storage space: IF((NSCA+NSCA0).GT.MXSCA)THEN WRITE(IDVOUT,FMT='(I6,I6,A)')MXSCA,NSCA0, & ' =MXSCA,NSCA0 (dirs in this scatt. plane)' CALL ERRMSG('FATAL','REAPAR', & 'MXSCA too small in DDSCAT') ENDIF THETAN(NSCA+1)=THETA1 IF(NSCA0.GT.2)THEN DO 23 J=2,NSCA0-1 THETAN(NSCA+J)=THETA1+(J-1)*DTHETA 23 CONTINUE ENDIF THETAN(NSCA+NSCA0)=THETA2 DO 24 J=1,NSCA0 PHIN(NSCA+J)=PHI1 24 CONTINUE NSCA=NSCA+NSCA0 WRITE(CMSGNM,FMT='(I4,A)')NSCA0, & ' = number of scattering angles in this scattering plane' CALL WRIMSG('REAPAR',CMSGNM) GOTO 22 25 CONTINUE C C Check that if user has requested computation of f_ml, user has C also requested linearly-polarized incident radiation. C 96.01.25 (BTD) We have disabled this "protection". So long as user C is aware that when using elliptical incident polarization the C f_ml are not the "usual" ones with index l referring to linearly C polarized incident modes, there is not reason not to allow user C to calculate f_ml for the the scattered radiation. C C IF(NSCA.GT.0)THEN C DO 26 J=1,3 C IF(AIMAG(CXE01(J)).NE.0.)CALL ERRMSG('FATAL','REAPAR', C & 'pol. vector E01 is complex: DDSCAT cannot compute f_ml, so C &do not ask it to!') C 26 CONTINUE C ENDIF C CLOSE (IOPAR) RETURN 40 CONTINUE CALL ERRMSG('FATAL','REAPAR',' Error reading ddscat.par file') RETURN END SUBROUTINE REASHP(CFLSHP,IDVOUT,IDVSHP,A1,A2,CDESCR,MXNAT,NAT, & IXYZ,ICOMP,MXN3,NAT3) C Arguments: INTEGER IDVOUT,IDVSHP,MXNAT,NAT,MXN3,NAT3 CHARACTER CDESCR*67,CFLSHP*(*) INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3) REAL A1(3),A2(3) C Local variables: INTEGER J,JJ C*********************************************************************** C Given: C CFLSHP = name of (unformatted) file containing target description C IDVSHP = device number to use in reading this file C IDVOUT = device number to use for "running" output C MXNAT,MXN3 = dimensioning information C Returns: C A1(3) = target axis 1 C A2(3) = target axis 2 C CDESCR = string describing target C NAT = number of dipoles in target C NAT3 = 3*NAT C IXYZ(1-NAT,1-3) = lattice coordinates of dipoles C ICOMP(1-NAT,1-3)= composition for directions x,y,z (in Target C Frame) of dipoles 1-NAT. C C Note: It is expected that many users will wish to modify this routine C to be consistent with whatever target generation software they C use. However, as an example (and for test purposes) this routine C is supplied with code compatible with the "target.out" files which C are generated by the "calltarget" program when called for C homogenous isotropic targets (shape.dat files do not include C information on target composition). C Also provided below (commented out) are C * a module of code appropriate for reading formatted target.out C file generated by the routine tarcel.f (concentric ellipsoids) C as an example of an inhomogeneous target. C * a module of code appropriate for reading an unformatted C "shape.dat" file for a general target, assuming user's C target generation routine writes an unformatted file with C the same structure. C C C*********************************************************************** C History: C 91.05.02 (BTD): Added IDVOUT to argument list. C 91.05.23 (BTD): Added A1, A2 to argument list. C 91.11.11 (BTD): Modified for compatibility with WRITE statements in C program CALLTARGET; Changed dimensioning of ICOMP from C ICOMP(MXN3) to ICOMP(MXNAT,3) C 93.09.24 (BTD): Modified to include examples of READ statements C for both unformatted and formatted shape.dat files. C Leave the formatted example (compatible with the C "target.out" files generated by "calltarget") in effect. C 95.04.07 (BTD): Changed CDESCR*60 -> CDESCR*67 for compatibility with C calling routine TARGET C 95.05.15 (BTD): Changed READ(DVSHP,FMT='(A60)') to ..(A67).. C Added module for inhomogeneous/anisotropic targets 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*********************************************************************** C C Enable following module of code if "shape.dat" file has been written C in formatted form consistent with the "target.out" files generated by C routines tar2el.f, tar3el.f, tarcyl.f, tarell.f, tarhex.f, tarrec.f, C and tartet.f (in particular, this assumes homogeneous target): C OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='FORMATTED') READ(IDVSHP,FMT='(A67)')CDESCR READ(IDVSHP,*)NAT READ(IDVSHP,*)A1 READ(IDVSHP,*)A2 READ(IDVSHP,*) READ(IDVSHP,*)(JJ,IXYZ(J,1),IXYZ(J,2),IXYZ(J,3),J=1,NAT) DO 0200 JJ=1,3 DO 0100 J=1,NAT ICOMP(J,JJ)=1 0100 CONTINUE 0200 CONTINUE CLOSE(IDVSHP) C C*********************************************************************** C C Enable following module of code if "shape.dat" file has been written C in formatted form consistent with the "target.out" files generated by C routine tarcel.f (inhomogenous target): C C OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='FORMATTED') C READ(IDVSHP,FMT='(A67)')CDESCR C READ(IDVSHP,*)NAT C READ(IDVSHP,*)A1 C READ(IDVSHP,*)A2 C READ(IDVSHP,*) C READ(IDVSHP,*)(JJ,IXYZ(J,1),IXYZ(J,2),IXYZ(J,3), C & ICOMP(J,1),ICOMP(J,2),ICOMP(J,3),J=1,NAT) C CLOSE(IDVSHP) C C*********************************************************************** C C Enable following module of code if "shape.dat" file has been written C in unformatted form consistent with these read statements: C C INTEGER JX,JY,JZ C OPEN(UNIT=IDVSHP,FILE=CFLSHP,STATUS='OLD',FORM='UNFORMATTED') C READ(IDVSHP)CDESCR,NAT,A1,A2 C READ(IDVSHP)((IXYZ(JX,JY),JX=1,NAT),JY=1,3), C & ((ICOMP(J,JJ),J=1,NAT),JJ=1,3) C CLOSE(IDVSHP) C C*********************************************************************** NAT3=3*NAT RETURN END SUBROUTINE REDUCE(CXV,IOCC,MXN3,MXNAT,NAT,NAT0) C*********************************************************************** C Given: C CXV(1-3*NAT) defined for NAT lattice sites C IOCC(1-NAT) = 0 or 1 depending on whether lattice site C is vacant or occupied C MXN3,MXNAT = dimensioning information C NAT = number of lattice sites C NAT0 = number of occupied lattice sites C Returns C CXV(1-3*NAT0) defined for NAT0 occupied lattice sites C C Arguments: INTEGER MXN3,MXNAT,NAT,NAT0 INTEGER*2 IOCC(MXNAT) COMPLEX CXV(MXN3) C Local variables: INTEGER J,JOC,M C C B.T.Draine, Princeton Univ. Observatory, 90/11/29 C History: C 90/12/03 (BTD): Modified for new ordering of vectors. C 90/12/05 (BTD): Corrected errors C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DO 20 M=0,2 JOC=0 DO 10 J=1,NAT IF(IOCC(J).EQ.1)THEN JOC=JOC+1 CXV(JOC+M*NAT0)=CXV(J+M*NAT) ENDIF 10 CONTINUE 20 CONTINUE RETURN END