SUBROUTINE GETSET(CWHAT,CVAR,INT) C .. Parameters .. INTEGER ISET PARAMETER (ISET=3) C .. C .. Scalar Arguments .. INTEGER INT CHARACTER CVAR* (*), CWHAT* (*) C .. C .. Local Scalars .. INTEGER I, IOERR C .. C .. Local Arrays .. INTEGER IARY(ISET) CHARACTER CARY(ISET)*10 C .. C .. External Subroutines .. EXTERNAL ERRMSG C .. C .. Data statements .. DATA CARY/'MXNAT', 'IOOUT', 'IOERR'/ DATA IARY/0, 6, 6/ DATA IOERR/6/ C .. C PURPOSE: SETS AND GETS INTEGER VARIABLES ("BASKET PROCEDURE") C FOR USE IN ERROR CHECKING OF PARAMETER STATEMENTS, AND/OR C SETTING INPUT/OUTPUT INTEGERS. THUS THESE VARIABLES DON'T HAVE C TO BE TRANSFERED AS FORMAL PARAMETERS OF THE PROCEDURE. C EACH NEW CVAR HAS TO BE MANUALY SET IN THIS ROUTINE! C CWHAT --- 'SET' OR 'GET' C CVAR --- CURRENTLY SET C 'MXNAT' C 'IOOUT' C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** IF (CWHAT .EQ. 'SET') THEN DO 10 I = 1, ISET IF (CVAR .EQ. CARY(I)) THEN IARY(I) = INT RETURN END IF 10 CONTINUE CALL ERRMSG('FATAL','GETSET',' WRONG CVAR ') RETURN C ELSE IF (CWHAT .EQ. 'GET') THEN DO 20 I = 1, ISET IF (CVAR .EQ. CARY(I)) THEN INT = IARY(I) RETURN END IF 20 CONTINUE CALL ERRMSG('FATAL','GETSET',' WRONG CVAR ') RETURN ELSE CALL ERRMSG('FATAL','GETSET',' WRONG CWHAT ') END IF RETURN END 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 MULT3(A,B,C) C*********************************************************************** C Purpose: C Given: A(I,J)= 3x3 matrix C B(J,K)= 3x3 matrix C Returns:C(I,K)= matrix product of A and B C This is auxilary for "rotation" routines C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. *********************************************************************** C Arguments: REAL A(3,3),B(3,3),C(3,3) C Local variables: REAL S INTEGER I,J,K C DO 30 J = 1, 3 DO 20 I = 1, 3 S = 0. DO 10 K = 1, 3 S = S + A(I,K)*B(K,J) 10 CONTINUE C(I,J) = S 20 CONTINUE 30 CONTINUE 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 New feature 92/04/01: C If NPHI=1, first angle is THETMI, last angle is THETMX C If NPHI=1 and NTHETA is odd, use Simpson's rule to evaluate C angle-average C If NPHI=1 and NTHETA is even, use trapezoidal rule to evaluate C angle-average C If NPHI>1, then angles are midpoints of intervals of equal C range in cos(theta) subdividing range from THETMI to THETMX C and weight each angle equally to evaluate angle-average 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 C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DELTA=(BETAMX-BETAMI)/REAL(NBETA) BETA(1)=BETAMI+0.5*DELTA IF(NBETA.GT.1)THEN DO 1000 J=1,NBETA BETA(J)=BETA(1)+DELTA*REAL(J-1) 1000 CONTINUE ENDIF C IF(NPHI.EQ.1.AND.NTHETA.GT.1)THEN DELTA=(COS(THETMX)-COS(THETMI))/REAL(NTHETA-1) THETA(1)=THETMI ELSE DELTA=(COS(THETMX)-COS(THETMI))/REAL(NTHETA) THETA(1)=ACOS(COS(THETMI)+0.5*DELTA) ENDIF IF(NTHETA.GT.1)THEN DO 2000 J=1,NTHETA THETA(J)=ACOS(COS(THETA(1))+DELTA*REAL(J-1)) 2000 CONTINUE ENDIF C DELTA=(PHIMAX-PHIMIN)/REAL(NPHI) PHI(1)=PHIMIN+0.5*DELTA IF(NPHI.GT.1)THEN DO 3000 J=1,NPHI PHI(J)=PHI(1)+DELTA*REAL(J-1) 3000 CONTINUE ENDIF 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 uniform weighting (i.e., random orientation) C C** New version: in special case NPHI=1 and NTHETA.GT.1, we let C THETMI,THETMX be smallest and largest values of theta in list C When NTHETA is even, we use trapezoidal integration C When NTHETA is odd, we use Simpson's rule integration. IF(NPHI.EQ.1)THEN IF(NTHETA.EQ.1)THEN WGTA(1,1)=1. ELSEIF(NTHETA.EQ.2)THEN WGTA(1,1)=0.5 WGTA(2,1)=0.5 ELSE IF(2*(NTHETA/2).EQ.NTHETA)THEN WG=1./REAL(NTHETA-1) WGTA(1,1)=.5*WG WGTA(NTHETA,1)=.5*WG DO 4100 J=2,NTHETA-1 WGTA(J,1)=WG 4100 CONTINUE ELSE WG=1./REAL(3*(NTHETA-1)) WGTA(1,1)=WG WGTA(NTHETA,1)=WG DO 4110 J=2,NTHETA-1,2 WGTA(J,1)=4.*WG 4110 CONTINUE IF(NTHETA.GE.5)THEN DO 4120 J=3,NTHETA-2,2 WGTA(J,1)=2.*WG 4120 CONTINUE ENDIF ENDIF ENDIF C ELSE WG=1./REAL(NTHETA*NPHI) DO 4300 I=1,NPHI DO 4200 J=1,NTHETA WGTA(J,I)=WG 4200 CONTINUE 4300 CONTINUE ENDIF WG=1./REAL(NBETA) DO 4500 I=1,NBETA WGTB(I)=WG 4500 CONTINUE C RETURN END SUBROUTINE PROD3(RM,VIN,VOUT) C Given: RM(I,J)=3x3 matrix C VIN(J)=3-vector C Returns:VOUT(I)=RM*VIN C C Auxilary for "rotation" routines C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments; REAL RM(3,3), VIN(3), VOUT(3) C Local variables: INTEGER I,J C DO 20 I = 1, 3 VOUT(I) = 0. DO 10 J = 1, 3 VOUT(I) = VOUT(I) + RM(I,J)*VIN(J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE PROD3C(RM,CVIN,CVOUT) C Arguments: REAL RM(3,3) COMPLEX CVIN(3),CVOUT(3) C Local variables: INTEGER I,J C Given: C RM=3x3 rotation matrix C CVIN=complex 3-vector C Computes: C CVOUT=RM \dot CVIN 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 I=1,3 CVOUT(I)=(0.,0.) DO 10 J=1,3 CVOUT(I)=CVOUT(I)+RM(I,J)*CVIN(J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE PROD3V(RM,VIN,VOUT,NV,NVMAX) C Arguments: INTEGER NV,NVMAX REAL RM(3,3),VIN(3,NVMAX),VOUT(3,NVMAX) C Local variables: INTEGER I,J,N C Given: C RM = 3x3 rotation matrix C VIN = NV 3-vectors C Computes: C VOUT=RM \dot VIN for each of the NV 3-vectors 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 I=1,3 DO 10 N=1,NV 10 VOUT(I,N)=0. 20 CONTINUE C DO 50 I=1,3 DO 40 J=1,3 DO 30 N=1,NV VOUT(I,N)=VOUT(I,N)+RM(I,J)*VIN(J,N) 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE REAPAR(AEFFA, & BETAMI,BETAMX, & CALPHA,CDIEL,CFLEPS,CFLPAR,CMETHD,CSHAPE, & 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,SHPAR1,SHPAR2,SHPAR3) 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. 9 (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 CMETHD 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 C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C Arguments: CHARACTER & CALPHA*(*),CDIEL*(*),CFLPAR*(*),CMETHD*(*),CSHAPE*(*) 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*40 CFLEPS(MXCOMP) REAL & BETAMI,BETAMX,DEGRAD,PHIMIN,PHIMAX,SHPAR1,SHPAR2,SHPAR3, & THETMI,THETMX,TOL REAL AEFFA(MXRAD),PHIN(MXSCA),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,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 Define solution method: C C CMETHD*6 = 'BRENNR' - Use Brenner's FOURX FFT routine. C = 'TMPRTN' - Use Temperton's FFT routine. C = 'CONVEX' - Use native Convex veclib routine C3DFFT. C READ(IOPAR,FMT=*,ERR=40)CMETHD IF(CMETHD.EQ.'BRENNR'.OR. & CMETHD.EQ.'TMPRTN'.OR. & CMETHD.EQ.'CONVEX')THEN WRITE(CMSGNM,FMT='(A,A)')CMETHD,' - Solution Method ' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' wrong definition of CMETHD') ENDIF C*** Check whether compiled with enough memory for TMPRTN option: IF(CMETHD.EQ.'TMPRTN'.AND.MXMEM.EQ.0)THEN CALL ERRMSG('FATAL','REAPAR', & ' Must recompile with MXMEM=1 to use TMPRTN option') 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 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') ELSE CALL ERRMSG('FATAL','REAPAR',' wrong definition of CALPHA') ENDIF 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 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')THEN WRITE(CMSGNM,FMT='(A,A)')CSHAPE,' - Shape definition ' CALL WRIMSG('REAPAR',CMSGNM) ELSE CALL ERRMSG('FATAL','REAPAR',' Wrong shape directive ') ENDIF READ(IOPAR,FMT=*,ERR=40)SHPAR1,SHPAR2,SHPAR3 IF(CSHAPE.EQ.'FRMFIL')THEN SHPAR1=0. SHPAR2=0. SHPAR3=0. 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 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') 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 READ(IOPAR,FMT=*,ERR=40)CXE01 C Construct orthogonal polarization vector CXE02=e02 CXE02(1)=(0.,0.) CXE02(2)=-CXE01(3) CXE02(3)=CXE01(2) 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 e01 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 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, & IX,IY,IZ,ICOMP,MXN3,NAT3) C Arguments: INTEGER IDVOUT,IDVSHP,MXNAT,NAT,MXN3,NAT3 CHARACTER CDESCR*60,CFLSHP*(*) INTEGER*2 IX(MXNAT),IY(MXNAT),IZ(MXNAT),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 IX(1-NAT),IY(1-NAT),IZ(1-NAT) = lattice coordinates of dipoles C ICOMP(1-NAT,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 provided with code compatible with the "target.out" files which C are generated by the "calltarget" program. C Also provided below (commented out) is a module of code appropriate C for reading an unformatted "shape.dat" file. 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 C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C 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)(IX(JX),JX=1,NAT),(IY(JY),JY=1,NAT),(IZ(JZ),JZ=1,NAT), C & ((ICOMP(J,JJ),J=1,NAT),JJ=1,3) C 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 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='(A60)')CDESCR READ(IDVSHP,*)NAT READ(IDVSHP,*)A1 READ(IDVSHP,*)A2 READ(IDVSHP,*) READ(IDVSHP,*)(JJ,IX(J),IY(J),IZ(J),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*********************************************************************** 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