SUBROUTINE FOURX(DATA,NN,NDIM,ISIGN) C Modified by PJF 11/17/1990 C to conform to DDSCAT needs: c c the cooley-tukey fast fourier transform in usasi basic fortran c c transform(j1,j2,,,,) = sum(data(i1,i2,,,,)*w1**((i2-1)*(j2-1)) c *w2**((i2-1)*(j2-1))*,,,), c where i1 and j1 run from 1 to nn(1) and w1=exp(isign*2*pi= c sqrt(-1)/nn(1)), etc. there is no limit on the dimensionality c (number of subscripts) of the data array. if an inverse c transform (isign=+1) is performed upon an array of transformed c (isign=-1) data, the original data will reappear. c multiplied by nn(1)*nn(2)*,,, the array of input data must be c in complex format. however, if all imaginary parts are zero (i.e. c the data are disguised real) running time is cut up to forty per- c cent. (for fastest transform of real data, nn(1) should be even.) c the transform values are always complex and are returned in the c original array of data, replacing the input data. the length c of each dimension of the data array may be any integer. the c program runs faster on composite integers than on primes, and is c particularly fast on numbers rich in factors of two. c c timing is in fact given by the following formula. let ntot be the c total number of points (real or complex) in the data array, that c is, ntot=nn(1)*nn(2)*... decompose ntot into its prime factors, c such as 2**k2 * 3**k3 * 5**k5 * ... let sum2 be the sum of all c the factors of two in ntot, that is, sum2 = 2*k2. let sumf be c the sum of all other factors of ntot, that is, sumf = 3*k3*5*k5*.. c the time taken by a multidimensional transform on these ntot data c is t = t0 + ntot*(t1+t2*sum2+t3*sumf). on the cdc 3300 (floating c point add time = six microseconds), t = 3000 + ntot*(600+40*sum2+ c 175*sumf) microseconds on complex data. c c implementation of the definition by summation will run in a time c proportional to ntot*(nn(1)+nn(2)+...). for highly composite ntot c the savings offered by this program can be dramatic. a one-dimen- c sional array 4000 in length will be transformed in 4000*(600+ c 40*(2+2+2+2+2)+175*(5+5+5)) = 14.5 seconds versus about 4000* c 4000*175 = 2800 seconds for the straightforward technique. c c the fast fourier transform places three restrictions upon the c data. c 1. the number of input data and the number of transform values c must be the same. c 2. both the input data and the transform values must represent c equispaced points in their respective domains of time and c frequency. calling these spacings deltat and deltaf, it must be c true that deltaf=2*pi/(nn(i)*deltat). of course, deltat need not c be the same for every dimension. c 3. conceptually at least, the input data and the transform output c represent single cycles of periodic functions. c c the calling sequence is-- c call fourt(data,nn,ndim,isign,iform,work) c c data is the array used to hold the real and imaginary parts c of the data on input and the transform values on output. it c is a multidimensional floating point array, with the real and c imaginary parts of a datum stored immediately adjacent in storage c (such as fortran iv places them). normal fortran ordering is c expected, the first subscript changing fastest. the dimensions c are given in the integer array nn, of length ndim. isign is -1 c to indicate a forward transform (exponential sign is -) and +1 c for an inverse transform (sign is +). iform is +1 if the data are c complex, 0 if the data are real. if it is 0, the imaginary c parts of the data must be set to zero. as explained above, the c transform values are always complex and are stored in array data. c work is an array used for working storage. it is floating point c real, one dimensional of length equal to twice the largest array c dimension nn(i) that is not a power of two. if all nn(i) are c powers of two, it is not needed and may be replaced by zero in the c calling sequence. thus, for a one-dimensional array, nn(1) odd, c work occupies as many storage locations as data. if supplied, c work must not be the same array as data. all subscripts of all c arrays begin at one. c c example 1. three-dimensional forward fourier transform of a c complex array dimensioned 32 by 25 by 13 in fortran iv. c dimension data(32,25,13),work(50),nn(3) c complex data c data nn/32,25,13/ c do 1 i=1,32 c do 1 j=1,25 c do 1 k=1,13 c 1 data(i,j,k)=complex value c call fourt(data,nn,3,-1,1,work) c c example 2. one-dimensional forward transform of a real array of c length 64 in fortran ii, c dimension data(2,64) c do 2 i=1,64 c data(1,i)=real part c 2 data(2,i)=0. c call fourt(data,64,1,-1,0,0) c c there are no error messages or error halts in this program. the c program returns immediately if ndim or any nn(i) is less than one. c c program by norman brenner from the basic program by charles c rader, june 1967. the idea for the digit reversal was c suggested by ralph alter. c c this is the fastest and most versatile version of the fft known c to the author. a program called four2 is available that also c performs the fast fourier transform and is written in usasi basic c fortran. it is about one third as long and restricts the c dimensions of the input array (which must be complex) to be powers c of two. another program, called four1, is one tenth as long and c runs two thirds as fast on a one-dimensional complex array whose c length is a power of two. c c reference-- c ieee audio transactions (june 1967), special issue on the fft. c C Parameters: INTEGER MXWORK PARAMETER (MXWORK=1024) C Arguments: INTEGER IFORM,ISIGN,NDIM C*** Modification 90/11/29 (BTD): C REAL DATA(1) C INTEGER NN(1) REAL DATA(*) INTEGER NN(*) C C Local variables: REAL DIFI,DIFR,OLDSI,OLDSR,RTHLF,SUMI,SUMR,T2I,T2R,T3I,T3R,T4I, + T4R,TEMPI,TEMPR,THETA,TWOPI,TWOWR,U1I,U1R,U2I,U2R,U3I,U3R, + U4I,U4R,W2I,W2R,W3I,W3R,WI,WR,WSTPI,WSTPR INTEGER I,I1,I1MAX,I1RNG,I2,I2MAX,I3,ICALL,ICASE,ICONJ,IDIM, + IDIV,IF,IFMIN,IFP1,IFP2,IMAX,IMIN,INON2,IPAR,IQUOT,IREM, + J,J1,J1MAX,J1MIN,J2,J2MAX,J2MIN,J2RNG,J3,J3MAX,JMAX,JMIN, + K1,K2,K3,K4,KDIF,KMIN,KSTEP,L,LMAX,M,MMAX, + N,NHALF,NMAX,NP0,NP1,NP1HF,NP1TW,NP2,NP2HF,NPREV,NTOT, + NTWO,NWORK INTEGER IFACT(32) REAL WORK(MXWORK) C C Intrinsic functions: INTRINSIC COS,MAX0,REAL,SIN C C Data statements: DATA NP0/0/,NPREV/0/ DATA TWOPI/6.2831853071796/,RTHLF/0.70710678118655/ Cflatau change: Hardwire that the data is complex. DATA IFORM, ICALL/1, 1/ C .. Cflatau addition: Error check IF(ICALL.EQ.1) THEN ICALL=2 NMAX = MAX( 2*NN(1), 2*NN(2), 2*NN(3) ) IF(NMAX .GT. MXWORK) THEN CALL ERRMSG('FATAL','FOURX', $ ' Increase size of working array MXWORK ') END IF IF(NDIM.NE.3) THEN CALL ERRMSG('FATAL','FOURX', $ ' This has to be 3D transform') END IF END IF IF (NDIM-1) 232,101,101 101 NTOT = 2 DO 103 IDIM = 1,NDIM IF (NN(IDIM)) 232,232,102 102 NTOT = NTOT*NN(IDIM) 103 CONTINUE c c main loop for each dimension c NP1 = 2 DO 231 IDIM = 1,NDIM N = NN(IDIM) NP2 = NP1*N IF (N-1) 232,230,104 c c is n a power of two and if not, what are its factors c 104 M = N NTWO = NP1 IF = 1 IDIV = 2 105 IQUOT = M/IDIV IREM = M - IDIV*IQUOT IF (IQUOT-IDIV) 113,106,106 106 IF (IREM) 108,107,108 107 NTWO = NTWO + NTWO IFACT(IF) = IDIV IF = IF + 1 M = IQUOT GO TO 105 108 IDIV = 3 INON2 = IF 109 IQUOT = M/IDIV IREM = M - IDIV*IQUOT IF (IQUOT-IDIV) 115,110,110 110 IF (IREM) 112,111,112 111 IFACT(IF) = IDIV IF = IF + 1 M = IQUOT GO TO 109 112 IDIV = IDIV + 2 GO TO 109 113 INON2 = IF IF (IREM) 115,114,115 114 NTWO = NTWO + NTWO GO TO 116 115 IFACT(IF) = M c c separate four cases-- c 1. complex transform or real transform for the 4th, 9th,etc. c dimensions. c 2. real transform for the 2nd or 3rd dimension. method-- c transform half the data, supplying the other half by con- c jugate symmetry. c 3. real transform for the 1st dimension, n odd. method-- c set the imaginary parts to zero. c 4. real transform for the 1st dimension, n even. method-- c transform a complex array of length n/2 whose real parts c are the even numbered real values and whose imaginary parts c are the odd numbered real values. separate and supply c the second half by conjugate symmetry. c 116 ICASE = 1 IFMIN = 1 I1RNG = NP1 IF (IDIM-4) 117,122,122 117 IF (IFORM) 118,118,122 118 ICASE = 2 I1RNG = NP0* (1+NPREV/2) IF (IDIM-1) 119,119,122 119 ICASE = 3 I1RNG = NP1 IF (NTWO-NP1) 122,122,120 120 ICASE = 4 IFMIN = 2 NTWO = NTWO/2 N = N/2 NP2 = NP2/2 NTOT = NTOT/2 I = 1 DO 121 J = 1,NTOT DATA(J) = DATA(I) I = I + 2 121 CONTINUE c c shuffle data by bit reversal, since n=2**k. as the shuffling c can be done by simple interchange, no working array is needed c 122 IF (NTWO-NP2) 132,123,123 123 NP2HF = NP2/2 J = 1 DO 131 I2 = 1,NP2,NP1 IF (J-I2) 124,127,127 124 I1MAX = I2 + NP1 - 2 DO 126 I1 = I2,I1MAX,2 DO 125 I3 = I1,NTOT,NP2 J3 = J + I3 - I2 TEMPR = DATA(I3) TEMPI = DATA(I3+1) DATA(I3) = DATA(J3) DATA(I3+1) = DATA(J3+1) DATA(J3) = TEMPR DATA(J3+1) = TEMPI 125 CONTINUE 126 CONTINUE 127 M = NP2HF 128 IF (J-M) 130,130,129 129 J = J - M M = M/2 IF (M-NP1) 130,128,128 130 J = J + M 131 CONTINUE GO TO 142 c c shuffle data by digit reversal for general n c 132 NWORK = 2*N DO 141 I1 = 1,NP1,2 DO 140 I3 = I1,NTOT,NP2 J = I3 DO 138 I = 1,NWORK,2 IF (ICASE-3) 133,134,133 133 WORK(I) = DATA(J) WORK(I+1) = DATA(J+1) GO TO 135 134 WORK(I) = DATA(J) WORK(I+1) = 0. 135 IFP2 = NP2 IF = IFMIN 136 IFP1 = IFP2/IFACT(IF) J = J + IFP1 IF (J-I3-IFP2) 138,137,137 137 J = J - IFP2 IFP2 = IFP1 IF = IF + 1 IF (IFP2-NP1) 138,138,136 138 CONTINUE I2MAX = I3 + NP2 - NP1 I = 1 DO 139 I2 = I3,I2MAX,NP1 DATA(I2) = WORK(I) DATA(I2+1) = WORK(I+1) I = I + 2 139 CONTINUE 140 CONTINUE 141 CONTINUE c c main loop for factors of two. perform fourier transforms of c length four, with one of length two if needed. the twiddle factor c w=exp(isign*2*pi*sqrt(-1)*m/(4*mmax)). check for w=isign*sqrt(-1) c and repeat for w=w*(1+isign*sqrt(-1))/sqrt(2). c 142 IF (NTWO-NP1) 174,174,143 143 NP1TW = NP1 + NP1 IPAR = NTWO/NP1 144 IF (IPAR-2) 149,146,145 145 IPAR = IPAR/4 GO TO 144 146 DO 148 I1 = 1,I1RNG,2 DO 147 K1 = I1,NTOT,NP1TW K2 = K1 + NP1 TEMPR = DATA(K2) TEMPI = DATA(K2+1) DATA(K2) = DATA(K1) - TEMPR DATA(K2+1) = DATA(K1+1) - TEMPI DATA(K1) = DATA(K1) + TEMPR DATA(K1+1) = DATA(K1+1) + TEMPI 147 CONTINUE 148 CONTINUE 149 MMAX = NP1 150 IF (MMAX-NTWO/2) 151,174,174 151 LMAX = MAX0(NP1TW,MMAX/2) DO 173 L = NP1,LMAX,NP1TW M = L IF (MMAX-NP1) 156,156,152 152 THETA = -TWOPI*REAL(L)/REAL(4*MMAX) IF (ISIGN) 154,153,153 153 THETA = -THETA 154 WR = COS(THETA) WI = SIN(THETA) 155 W2R = WR*WR - WI*WI W2I = 2.*WR*WI W3R = W2R*WR - W2I*WI W3I = W2R*WI + W2I*WR 156 DO 169 I1 = 1,I1RNG,2 KMIN = I1 + IPAR*M IF (MMAX-NP1) 157,157,158 157 KMIN = I1 158 KDIF = IPAR*MMAX 159 KSTEP = 4*KDIF IF (KSTEP-NTWO) 160,160,169 160 DO 168 K1 = KMIN,NTOT,KSTEP K2 = K1 + KDIF K3 = K2 + KDIF K4 = K3 + KDIF IF (MMAX-NP1) 161,161,164 161 U1R = DATA(K1) + DATA(K2) U1I = DATA(K1+1) + DATA(K2+1) U2R = DATA(K3) + DATA(K4) U2I = DATA(K3+1) + DATA(K4+1) U3R = DATA(K1) - DATA(K2) U3I = DATA(K1+1) - DATA(K2+1) IF (ISIGN) 162,163,163 162 U4R = DATA(K3+1) - DATA(K4+1) U4I = DATA(K4) - DATA(K3) GO TO 167 163 U4R = DATA(K4+1) - DATA(K3+1) U4I = DATA(K3) - DATA(K4) GO TO 167 164 T2R = W2R*DATA(K2) - W2I*DATA(K2+1) T2I = W2R*DATA(K2+1) + W2I*DATA(K2) T3R = WR*DATA(K3) - WI*DATA(K3+1) T3I = WR*DATA(K3+1) + WI*DATA(K3) T4R = W3R*DATA(K4) - W3I*DATA(K4+1) T4I = W3R*DATA(K4+1) + W3I*DATA(K4) U1R = DATA(K1) + T2R U1I = DATA(K1+1) + T2I U2R = T3R + T4R U2I = T3I + T4I U3R = DATA(K1) - T2R U3I = DATA(K1+1) - T2I IF (ISIGN) 165,166,166 165 U4R = T3I - T4I U4I = T4R - T3R GO TO 167 166 U4R = T4I - T3I U4I = T3R - T4R 167 DATA(K1) = U1R + U2R DATA(K1+1) = U1I + U2I DATA(K2) = U3R + U4R DATA(K2+1) = U3I + U4I DATA(K3) = U1R - U2R DATA(K3+1) = U1I - U2I DATA(K4) = U3R - U4R DATA(K4+1) = U3I - U4I 168 CONTINUE KDIF = KSTEP KMIN = 4* (KMIN-I1) + I1 GO TO 159 169 CONTINUE M = M + LMAX IF (M-MMAX) 170,170,173 170 IF (ISIGN) 171,172,172 171 TEMPR = WR WR = (WR+WI)*RTHLF WI = (WI-TEMPR)*RTHLF GO TO 155 172 TEMPR = WR WR = (WR-WI)*RTHLF WI = (TEMPR+WI)*RTHLF GO TO 155 173 CONTINUE IPAR = 3 - IPAR MMAX = MMAX + MMAX GO TO 150 c c main loop for factors not equal to two. apply the twiddle factor c w=exp(isign*2*pi*sqrt(-1)*(j1-1)*(j2-j1)/(ifp1+ifp2)), then c perform a fourier transform of length ifact(if), making use of c conjugate symmetries. c 174 IF (NTWO-NP2) 175,201,201 175 IFP1 = NTWO IF = INON2 NP1HF = NP1/2 176 IFP2 = IFACT(IF)*IFP1 J1MIN = NP1 + 1 IF (J1MIN-IFP1) 177,177,184 177 DO 183 J1 = J1MIN,IFP1,NP1 THETA = -TWOPI*REAL(J1-1)/REAL(IFP2) IF (ISIGN) 179,178,178 178 THETA = -THETA 179 WSTPR = COS(THETA) WSTPI = SIN(THETA) WR = WSTPR WI = WSTPI J2MIN = J1 + IFP1 J2MAX = J1 + IFP2 - IFP1 DO 182 J2 = J2MIN,J2MAX,IFP1 I1MAX = J2 + I1RNG - 2 DO 181 I1 = J2,I1MAX,2 DO 180 J3 = I1,NTOT,IFP2 TEMPR = DATA(J3) DATA(J3) = DATA(J3)*WR - DATA(J3+1)*WI DATA(J3+1) = TEMPR*WI + DATA(J3+1)*WR 180 CONTINUE 181 CONTINUE TEMPR = WR WR = WR*WSTPR - WI*WSTPI WI = TEMPR*WSTPI + WI*WSTPR 182 CONTINUE 183 CONTINUE 184 THETA = -TWOPI/REAL(IFACT(IF)) IF (ISIGN) 186,185,185 185 THETA = -THETA 186 WSTPR = COS(THETA) WSTPI = SIN(THETA) J2RNG = IFP1* (1+IFACT(IF)/2) DO 200 I1 = 1,I1RNG,2 DO 199 I3 = I1,NTOT,NP2 J2MAX = I3 + J2RNG - IFP1 DO 197 J2 = I3,J2MAX,IFP1 J1MAX = J2 + IFP1 - NP1 DO 193 J1 = J2,J1MAX,NP1 J3MAX = J1 + NP2 - IFP2 DO 192 J3 = J1,J3MAX,IFP2 JMIN = J3 - J2 + I3 JMAX = JMIN + IFP2 - IFP1 I = 1 + (J3-I3)/NP1HF IF (J2-I3) 187,187,189 187 SUMR = 0. SUMI = 0. DO 188 J = JMIN,JMAX,IFP1 SUMR = SUMR + DATA(J) SUMI = SUMI + DATA(J+1) 188 CONTINUE WORK(I) = SUMR WORK(I+1) = SUMI GO TO 192 189 ICONJ = 1 + (IFP2-2*J2+I3+J3)/NP1HF J = JMAX SUMR = DATA(J) SUMI = DATA(J+1) OLDSR = 0. OLDSI = 0. J = J - IFP1 190 TEMPR = SUMR TEMPI = SUMI SUMR = TWOWR*SUMR - OLDSR + DATA(J) SUMI = TWOWR*SUMI - OLDSI + DATA(J+1) OLDSR = TEMPR OLDSI = TEMPI J = J - IFP1 IF (J-JMIN) 191,191,190 191 TEMPR = WR*SUMR - OLDSR + DATA(J) TEMPI = WI*SUMI WORK(I) = TEMPR - TEMPI WORK(ICONJ) = TEMPR + TEMPI TEMPR = WR*SUMI - OLDSI + DATA(J+1) TEMPI = WI*SUMR WORK(I+1) = TEMPR + TEMPI WORK(ICONJ+1) = TEMPR - TEMPI 192 CONTINUE 193 CONTINUE IF (J2-I3) 194,194,195 194 WR = WSTPR WI = WSTPI GO TO 196 195 TEMPR = WR WR = WR*WSTPR - WI*WSTPI WI = TEMPR*WSTPI + WI*WSTPR 196 TWOWR = WR + WR 197 CONTINUE I = 1 I2MAX = I3 + NP2 - NP1 DO 198 I2 = I3,I2MAX,NP1 DATA(I2) = WORK(I) DATA(I2+1) = WORK(I+1) I = I + 2 198 CONTINUE 199 CONTINUE 200 CONTINUE IF = IF + 1 IFP1 = IFP2 IF (IFP1-NP2) 176,201,201 c c complete a real transform in the 1st dimension, n even, by con- c jugate symmetries. c 201 GO TO (230,220,230,202) ICASE 202 NHALF = N N = N + N THETA = -TWOPI/REAL(N) IF (ISIGN) 204,203,203 203 THETA = -THETA 204 WSTPR = COS(THETA) WSTPI = SIN(THETA) WR = WSTPR WI = WSTPI IMIN = 3 JMIN = 2*NHALF - 1 GO TO 207 205 J = JMIN DO 206 I = IMIN,NTOT,NP2 SUMR = (DATA(I)+DATA(J))/2. SUMI = (DATA(I+1)+DATA(J+1))/2. DIFR = (DATA(I)-DATA(J))/2. DIFI = (DATA(I+1)-DATA(J+1))/2. TEMPR = WR*SUMI + WI*DIFR TEMPI = WI*SUMI - WR*DIFR DATA(I) = SUMR + TEMPR DATA(I+1) = DIFI + TEMPI DATA(J) = SUMR - TEMPR DATA(J+1) = -DIFI + TEMPI J = J + NP2 206 CONTINUE IMIN = IMIN + 2 JMIN = JMIN - 2 TEMPR = WR WR = WR*WSTPR - WI*WSTPI WI = TEMPR*WSTPI + WI*WSTPR 207 IF (IMIN-JMIN) 205,208,211 208 IF (ISIGN) 209,211,211 209 DO 210 I = IMIN,NTOT,NP2 DATA(I+1) = -DATA(I+1) 210 CONTINUE 211 NP2 = NP2 + NP2 NTOT = NTOT + NTOT J = NTOT + 1 IMAX = NTOT/2 + 1 212 IMIN = IMAX - 2*NHALF I = IMIN GO TO 214 213 DATA(J) = DATA(I) DATA(J+1) = -DATA(I+1) 214 I = I + 2 J = J - 2 IF (I-IMAX) 213,215,215 215 DATA(J) = DATA(IMIN) - DATA(IMIN+1) DATA(J+1) = 0. IF (I-J) 217,219,219 216 DATA(J) = DATA(I) DATA(J+1) = DATA(I+1) 217 I = I - 2 J = J - 2 IF (I-IMIN) 218,218,216 218 DATA(J) = DATA(IMIN) + DATA(IMIN+1) DATA(J+1) = 0. IMAX = IMIN GO TO 212 219 DATA(1) = DATA(1) + DATA(2) DATA(2) = 0. GO TO 230 c c complete a real transform for the 2nd or 3rd dimension by c conjugate symmetries. c 220 IF (I1RNG-NP1) 221,230,230 221 DO 229 I3 = 1,NTOT,NP2 I2MAX = I3 + NP2 - NP1 DO 228 I2 = I3,I2MAX,NP1 IMIN = I2 + I1RNG IMAX = I2 + NP1 - 2 JMAX = 2*I3 + NP1 - IMIN IF (I2-I3) 223,223,222 222 JMAX = JMAX + NP2 223 IF (IDIM-2) 226,226,224 224 J = JMAX + NP0 DO 225 I = IMIN,IMAX,2 DATA(I) = DATA(J) DATA(I+1) = -DATA(J+1) J = J - 2 225 CONTINUE 226 J = JMAX DO 227 I = IMIN,IMAX,NP0 DATA(I) = DATA(J) DATA(I+1) = -DATA(J+1) J = J - NP0 227 CONTINUE 228 CONTINUE 229 CONTINUE c c end of loop on each dimension c 230 NP0 = NP1 NP1 = NP2 NPREV = N 231 CONTINUE 232 RETURN c c revision history--- c c january 1978 deleted references to the *cosy cards and c added revision history c----------------------------------------------------------------------- END SUBROUTINE GETFML(AKR,AK3,AKSR, & CALPHA,CMETHD,CXADIA,CXALPH,CXE,CXE01R,CXE02R,CXEPS, & CXF11,CXF12,CXF21,CXF22,CXPOL,CXSC,CXZC,CXZW, & EM1R,EM2R, & ICOMP,ICTHM,IDVOUT,INIT,IOCC,IORTH,IPHI,IPHIM,ITASK,IXYZ0, & LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0, & MXCOMP,MXCXSC,MXMEM,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXSCA, & NAT,NAT0,NAT3,NSCAT,NX,NY,NZ, & PHI,PHIN,PIA2, & QABS,QBKSCA,QEXT,QPHA,QSCA,QSCAG, & SCRRS1,THETAN,TOL) C C*********************************************************************** C Function of GETFML is to return the scattering function fml in C preselected scattering directions. C C Given: C AKR(1-3)=(k_x,k_y,k_z)*d for incident wave C AK3 =(kd)**3 C AKSR(1-3,1-NSCAT)=(k_x,k_y,k_z)*d for NSCAT scattering directions C CALPHA =descriptor of method used for assigning polarizabilities C =LATTDR or DRAI88 or GOBR88 C CMETHD =descriptor of method used for FFTs C CXE01R(1-3)=incident polarization state 1 at origin C CXE02R(1-3)=incident polarization state 2 at origin C CXEPS(1-NCOMP)=dielectric constant for compositions 1-NCOMP C ICOMP(1-NAT3)=x,y,z "composition" for sites 1-NAT C ICTHM =number of theta values for computing scattering integrals C IDVOUT =device for running output C INIT =0,1,2 for choice of |x0> for CCG method C IOCC(1-NAT)=0,1 if site in extended target is vacant,occupied C IORTH =1,2 to do 1,2 orthogonal pol states C IPHI =number of phi values for target orientation C IPHIM =number of phi values for computing scatt. integrals C ITASK not used C IXYZ0(1-MXNAT,1-3)=lattice coordinates for sites 1-NAT C LACE,LAXI,LCLM,LGI,LPI,LQI,LSC0=integers used to assign locations C in scratch array C MXCOMP =dimensioning info: maximum number of allowed compositions C MXCXSC =dimensioning info for complex scratch space C MXMEM =dimensioning info (0 or 1) C MXN3 =dimensioning info (3*max number of sites) C MXNAT =dimensioning info (max number of sites) C MXNX =dimensioning info (max x extent of target) C MXNY =dimensioning info (max y extent of target) C MXNZ =dimensioning info (max z extent of target) C MXPHI =dimensioning info (max number of phi vals for targ. orient) C MXSCA =dimensioning info (max number of scat. dirs for f_ml) C NAT =number of sites in extended target C NAT0 =number of sites in original target C NAT3 =3*NAT C NSCAT =number of scat. dirs for f_ml C NX =x extent of extended target C NY =y extent of extended target C NZ =z extent of extended target C PHI(1-NPHI)=phi values target orientation C PHIN(1-NSCAT)=phi values for scat. dirs for f_ml comp. C PIA2 =\pi*(3*NAT/4\pi)^{2/3} (NAT=number of sites in original target) C THETAN(1-NSCAT)=theta values for scat. dirs for f_ml comp C TOL =tolerance for terminating conjugate gradient iteration C C and scratch space: C CXSC(1-MXCXSC)=scratch space C CXZC(1-MXNX+1,1-MXNY+1,1-MXNZ+1,1-6)=scratch space C CXZW(1-MXNAT,1-8*(3+MXMEM))=scratch space C SCRRS1(1-MXN3)=scratch space C C Returns: C CXADIA(1-NAT3)=diagnoal elements of "A matrix" C CXALPH(1-NAT3)=complex x,y,z polarizabilities of dipoles 1-NAT C CXE(1-NAT3)=incident x,y,z E field at dipoles 1-NAT C CXF11(1-NSCAT)=scattering matrix element f_11 for NSCAT dirs C CXF12(1-NSCAT)= f_12 C CXF21(1-NSCAT)= f_21 C CXF22(1-NSCAT)= f_22 C CXPOL(1-NAT3)=x,y,z polarization of dipoles 1-NAT C EM1R(1-3,1-NSCAT)=scattered wave for incident pol state 1 C EM2R(1-3,1-NSCAT)=scattered wave for incident pol state 2 C QABS(1-2)=Q_abs=C_abs/(PIA2*d**2) for incident pols 1,2 C QBKSCA(1-2)=diff.scatt.cross section for backscat for inc.pols.1,2 C QEXT(1-2)=Q_ext=C_ext/(PIA2*d**2) for incident pols 1,2 C QPHA(1-2)=phase shift cross section for incident pols 1,2 C QSCA(1-2)=Q_sca=C_sca/(PIA2*d**2) for incident pols 1,2 C QSCAG(1-2)=*Q_sca for incident pols 1,2 C C B.T.Draine, Princeton Univ. Observatory, 1990. C History: C 90.11.06 (BTD): modified to change argument list and pass array C dimensions. C 90.11.29 (BTD): further modifications, including arguments C changed calls to EVALE C 90.11.29 (BTD): changes to use REDUCE so that computations in C EVALQ and SCAT are restricted to occupied sites. C 90.11.30 (BTD): multiple changes to reduce storage requirements C 90.12.03 (BTD): change ordering of XYZ0 C 90.12.10 (BTD): remove XYZ0 and use IXYZ0 instead C 90.12.13 (BTD): add new argument to calls to EVALQ C 91.01.03 (BTD): modified so that when IORTH=1 no attempt is made C to compute f_{ml} for "standard" definition of C incident pol. states defined by scattering plane. C 91.05.09 (BTD): after call to SCAT, adjust QSCA and QSCAG to make C QSCA=QEXT-QABS C 91.05.09 (BTD): eliminate redundant write statement C 91.08.14 (BTD): add QBKSCA to argument list for GETFML C add CBKSCA to argument list for SCAT C compute QBKSCA from CBKSCA C 91.09.17 (BTD): move call to subroutine ALPHA from DDSCAT to GETFML C (since polarizability from Lattice Dipersion Relation C depends on direction of propagation and polarization) C add CALPHA, CXEPS, ICOMP, MXCOMP to argument list for C GETFML C 93.01.14 (BTD): change method for computing QSCA to use result from C SCAT when albedo << .01 , to use (QEXT-QABS) when C albedo >> .01, and to smoothly join these limits. C 93.03.11 (BTD): deleted all code associated with unused variable C CMACHN (originally included to identify machine/OS) C 93.07.09 (BTD): modify method for computing QSCA to use .03 as C the dividing point rather than .01 . C 93.11.23 (BTD): added SAVE E02 to preserve value of E02 between calls C [to eliminate problem encountered on Silicon Graphics C machines when IPHI>1] 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*6,CMETHD*6 C INTEGER ICTHM,IDVOUT,INIT,IORTH,IPHI,IPHIM,ITASK,JO, & LACE,LAXI,LCLM,LGI,LPI,LSC0,LQI, & MXCOMP,MXCXSC,MXMEM,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXSCA, & NAT,NAT0,NAT3,NSCAT,NX,NY,NZ C INTEGER*2 ICOMP(MXN3),IOCC(MXNAT),IXYZ0(MXNAT,3) C REAL AK3,PIA2,TOL C REAL AKR(3),AKSR(3,MXSCA),EM1R(3,MXSCA),EM2R(3,MXSCA), & PHI(MXPHI),PHIN(MXSCA), & QABS(2),QBKSCA(2),QEXT(2),QPHA(2),QSCA(2),QSCAG(2), & SCRRS1(MXN3),THETAN(MXSCA) C COMPLEX CXADIA(MXN3),CXALPH(MXN3), & CXE(MXN3),CXE01R(3),CXE02R(3),CXEPS(MXCOMP), & CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),CXF22(MXSCA), & CXSC(MXCXSC),CXPOL(MXN3),CXZC(MXNX+1,MXNY+1,MXNZ+1,6), & CXZW(MXNAT,8*(3+MXMEM)) C C Local variables COMPLEX CXE0R(3),CXVAR1,CXVAR2 INTEGER I,NAT03 REAL CABS,CBKSCA,CEXT,CPHA,CSCA,CSCAG,COSPHI,E02,FALB,SINPHI SAVE E02 C NAT03=3*NAT0 C C*** Compute fml directly only if IPHI=1 C Otherwise use previously computed fml values to obtain new fml C IF(IPHI.EQ.1)THEN DO 70 JO=1,IORTH IF(JO.EQ.1)THEN CXE0R(1)=CXE01R(1) CXE0R(2)=CXE01R(2) CXE0R(3)=CXE01R(3) ELSE CXE0R(1)=CXE02R(1) CXE0R(2)=CXE02R(2) CXE0R(3)=CXE02R(3) ENDIF E02=0. DO 10 I=1,3 E02=E02+CXE0R(I)*CONJG(CXE0R(I)) 10 CONTINUE C C*** call EVALE for NAT sites C CALL EVALE(CXE0R,AKR,IXYZ0,MXNAT,MXN3,NAT,NAT0, & NX,NY,NZ,CXE) C C*** call ALPHA to determine polarizabilities at NAT sites C CALL ALPHA(CALPHA,CXALPH,CXE0R,CXEPS,ICOMP,AKR,NAT3, & MXCOMP,MXN3) C C*** call EVALA to prepare A matrix elements C CALL EVALA(CXADIA,CXALPH,MXN3,NAT) C C*** first call to timeit: C (disable if putting calls into ddaccg itself) C CALL TIMEIT(' DDACCG ') C C*** If JO=2 and IPHIM.GT.1 then move polarization vector P previously C computed for JO=1 [stored in CXSC(LACE)] to storage in CXSC(LSC0) C IF(JO.EQ.2.AND.IPHIM.GT.1)THEN CALL COPYIT(CXSC(LACE),CXSC(LSC0),NAT3) ENDIF C CALL DDACCG(AKR,AK3,CMETHD, & CXSC(LACE),CXADIA,CXALPH,CXSC(LAXI), & CXSC(LCLM),CXE,CXE0R,CXSC(LGI),CXSC(LQI), & CXSC(LPI),CXPOL, & CXZC,CXZW, & IDVOUT,INIT,IOCC, & MXMEM,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ, & PIA2,QABS(JO),QEXT(JO),QPHA(JO),QSCA(JO), & SCRRS1,TOL) C C*** Now reduce the computed polarization vector CXPOL C to retain only real elements C CALL REDUCE(CXPOL,IOCC,MXN3,MXNAT,NAT,NAT0) C C*** second call to timeit: C disable if putting calls into ddaccg itself: CALL TIMEIT(' DDACCG ') C C*** If IPHIM.GT.1: Store reduced polarization vector P in CXSC(LACE) C for use in subsequent computations for other phi values C IF(IPHIM.GT.1)CALL COPYIT(CXPOL,CXSC(LACE),NAT03) C C*** first call to timeit: CALL TIMEIT(' SCAT ') C IF(JO.EQ.1)THEN CALL SCAT(AKR,AKSR,EM1R,EM2R,E02, & CBKSCA,CSCA,CSCAG,CXE01R,CXF11,CXF21,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,THETAN,IXYZ0) ELSEIF(JO.EQ.2)THEN CALL SCAT(AKR,AKSR,EM1R,EM2R,E02, & CBKSCA,CSCA,CSCAG,CXE01R,CXF12,CXF22,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,THETAN,IXYZ0) ENDIF C C*** second call to timeit CALL TIMEIT(' SCAT ') QBKSCA(JO)=CBKSCA/PIA2 QSCA(JO)=CSCA/PIA2 C C*** Since Q_ext and Q_abs are calculated "exactly", when particles C are "large" we obtain Q_sca from (Q_ext-Q_abs) rather than C scattering cross section computed by SCAT. This is particularly C advantageous when the particles are large with complex scattering C patterns, so that angular integration of the differential scattering C cross section would be either very time-consuming or inaccurate. C However, when particles are "small" the scattering cross section C is small compared to Q_ext and Q_abs, so that the difference C (Q_ext-Q_abs) may be inaccurate. Accordingly, when particles C are "small" we instead accept the scattering cross section C computed by SCAT (probably quite accurate when albedo is small because C the scattering pattern tends to be dipolar and therefore well-resolved C by a reasonable number of scattering directions ICTHM and IPHI). C We use an ad-hoc procedure to smoothly make the transitions C between these regimes, based on the albedo. C The "dividing" line of albedo=.03 is chosen on basis of C numerical experiments. If Q_ext and Q_abs are each calculated to C accuracy of 1 times 10^{-4} (say), then the C difference between Q_ext and Q_abs would be accurate to C 0.0001*Q_ext*sqrt(2) . When albedo=.03 (Q_ext-Q_abs)/Q_sca would C then be accurate to .0001*sqrt(2)/.03=.005 C FALB=(QSCA(JO)/(.03*QEXT(JO)))**2 QSCA(JO)=(QSCA(JO)+(QEXT(JO)-QABS(JO))*FALB)/(1.+FALB) C C Must now correct value of QSCAG C QSCAG(JO)=QSCA(JO)*CSCAG/CSCA C 70 CONTINUE C ELSEIF(IPHI.GE.2)THEN C C C*** Have previously computed dipole polarization vectors for C first phi value (and reduced them to NAT03 elements). C Use these two solutions to obtain Qabs, Qext, Qpha, Qsca, g*Qsca, C and fml for any additional values of target rotation phi. C CXSC(LSC0) contains polarization vector for IPHI=1 and JO=1 C CXSC(LACE) contains polarization vector for IPHI=1 and JO=2 C C*** Call EVALE to obtain appropriate E vector C*** first call to timing routine CALL TIMEIT(' EVALE') C C*** Call EVALE for NAT0 occupied sites CALL EVALE(CXE01R,AKR,IXYZ0,MXNAT,MXN3,NAT0,NAT0,NX,NY,NZ,CXE) C C*** second call to timing routine CALL TIMEIT(' EVALE') C C*** Construct polarization for incident pol. state 1 COSPHI=COS(PHI(IPHI)-PHI(1)) SINPHI=SIN(PHI(IPHI)-PHI(1)) DO 80 I=1,NAT03 CXPOL(I)=COSPHI*CXSC(LSC0+I-1)-SINPHI*CXSC(LACE+I-1) 80 CONTINUE C*** first call to timing routine CALL TIMEIT(' EVALQ') C C Need to construct reduced vector CXALPH for use by EVALQ C -- store in CXSC(LQI) CALL COPYIT(CXALPH,CXSC(LQI),NAT3) CALL REDUCE(CXSC(LQI),IOCC,MXN3,MXNAT,NAT,NAT0) CALL EVALQ(CXSC(LQI),AKR,NAT03,E02,CXE,CXPOL,CABS,CEXT,CPHA, & MXN3,1) C C*** second call to timing routine CALL TIMEIT(' EVALQ') C QABS(1)=CABS/PIA2 QEXT(1)=CEXT/PIA2 QPHA(1)=CPHA/PIA2 C C*** first call to timing routine CALL TIMEIT(' SCAT') C C*** Now call SCAT C CXPOL has been reduced to NAT03 elements C IXYZ0 was initially given NAT03 elements C CALL SCAT(AKR,AKSR,EM1R,EM2R,E02, & CBKSCA,CSCA,CSCAG,CXE01R,CXF11,CXF21,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,THETAN,IXYZ0) C C*** second call to timing routine CALL TIMEIT(' SCAT') C QBKSCA(1)=CBKSCA/PIA2 QSCA(1)=CSCA/PIA2 C C*** As above, compute Q_sca from either SCAT or (Q_ext-Q_abs) C depending on whether albedo << .03 or albedo >> .03 C FALB=(QSCA(1)/(.03*QEXT(1)))**2 QSCA(1)=(QSCA(1)+(QEXT(1)-QABS(1))*FALB)/(1.+FALB) QSCAG(1)=QSCA(1)*CSCAG/CSCA C C*** first call to timing routine CALL TIMEIT(' EVALE') C C*** Call EVALE for NAT0 occupied sites C to obtain appropriate E vector for pol. state 2 CALL EVALE(CXE02R,AKR,IXYZ0,MXNAT,MXN3,NAT0,NAT0,NX,NY,NZ,CXE) C C*** second call to timing routine CALL TIMEIT(' EVALE') C C*** Construct incident polarization vector for pol. state 2 COSPHI=COS(PHI(IPHI)-PHI(1)) SINPHI=SIN(PHI(IPHI)-PHI(1)) DO 90 I=1,NAT03 CXPOL(I)=COSPHI*CXSC(LACE+I-1)+SINPHI*CXSC(LSC0+I-1) 90 CONTINUE C*** first call to timing routine CALL TIMEIT(' EVALQ ') C C*** Recall that reduced vector CXALPH is still stored in CXSC(LQI) CALL EVALQ(CXSC(LQI),AKR,NAT03,E02,CXE,CXPOL,CABS,CEXT,CPHA, & MXN3,1) C*** second call to timing routine CALL TIMEIT(' EVALQ ') C QABS(2)=CABS/PIA2 QEXT(2)=CEXT/PIA2 QPHA(2)=CPHA/PIA2 C C*** first call to timing routine CALL TIMEIT(' SCAT') C C*** Call SCAT C CXPOL has NAT03 elements C IXYZ0 was initially given NAT03 elements C CALL SCAT(AKR,AKSR,EM1R,EM2R,E02, & CBKSCA,CSCA,CSCAG,CXE01R,CXF12,CXF22,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,THETAN,IXYZ0) C C*** second call to timing routine CALL TIMEIT(' SCAT') C QBKSCA(2)=CBKSCA/PIA2 QSCA(2)=CSCA/PIA2 C*** As above, compute Q_sca from either SCAT or (Q_ext-Q_abs) C depending on whether albedo << .03 or albedo >> .03 C FALB=(QSCA(2)/(.03*QEXT(2)))**2 QSCA(2)=(QSCA(2)+(QEXT(2)-QABS(2))*FALB)/(1.+FALB) QSCAG(2)=QSCA(2)*CSCAG/CSCA C ENDIF C C*** If IORTH=2, then we have results for two incident polarization C states CXE01R and CXE02R, and it is possible to compute f_{ml} for C standard definition of f_{ml} in terms of incident polarization C states 1 and 2 defined by the scattering plane. C If, however, computations were done for IORTH=1, then results C exist only for single incident polarization state CXE01R, and C it is not possible to construct results for incident polarization C states defined by the scattering plane. In this case we therefore C use results for f_{m1}, where polarization state 1 = CXE01R. C IF(IORTH.EQ.2)THEN C*** Values of CXF11, CXF12, CXF21, CXF22 are for incident polarization C states CXE01R and CXE02R. Now convert to standard definition of C f_{ml}: incident pol states 1 and 2 defined by scattering plane DO 100 I=1,NSCAT CXVAR1=CXF11(I) CXVAR2=CXF21(I) COSPHI=COS(PHIN(I)) SINPHI=SIN(PHIN(I)) CXF11(I)=COSPHI*CXVAR1+SINPHI*CXF12(I) CXF21(I)=COSPHI*CXVAR2+SINPHI*CXF22(I) CXF12(I)=COSPHI*CXF12(I)-SINPHI*CXVAR1 CXF22(I)=COSPHI*CXF22(I)-SINPHI*CXVAR2 100 CONTINUE ENDIF RETURN END