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,CMDSOL,CMDFFT,CMDTRQ, & CXADIA,CXALPH,CXE,CXE01R,CXE02R,CXEPS, & CXF11,CXF12,CXF21,CXF22,CXPOL,CXSC,CXSCR1,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,QTRQAB,QTRQSC, & SCRRS1,SCRRS2,THETAN,TOL,TIMERS,MXTIMERS,NTIMERS) 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 [in Target Frame] C AK3 =(kd)**3 C AKSR(1-3,1-NSCAT)=(k_x,k_y,k_z)*d in TF for NSCAT scattering directions C CALPHA =descriptor of method used for assigning polarizabilities C =LATTDR or DRAI88 or GOBR88 C CMDSOL =descriptor of method used for iterative solution C CMDFFT =descriptor of method used for FFTs C CXE01R(1-3)=incident polarization state 1 at origin [in TF] C CXE02R(1-3)=incident polarization state 2 at origin [in TF] C CXEPS(1-NCOMP)=dielectric constant for compositions 1-NCOMP C EM1R(1-3,1-NSCAT)=unit scat. polarization vectors 1 in TF C EM2R(1-3,1-NSCAT)=unit scat. polarization vectors 2 in TF 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 =which phi value 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 (in TF) C with sites 1-NAT0 in list being the occupied sites. 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 CXSCR1(1-MXN3)=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 SCRRS2(1-MXNAT)=scratch space C C Returns: C CXADIA(1-NAT3)=diagonal elements of "A matrix" C CXALPH(1-NAT3)=x,y,z dipole polarizabilities (in TF) C CXE(1-NAT3)=incident x,y,z E field at dipoles 1-NAT (in TF) 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 in TF C QABS(1-2)=Q_abs=C_abs/(PIA2*d**2) for incident pols 1,2 C QBKSCA(1-2)=diff.scatt.cross section/(PIA2*d^2) for backscat C 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-3,1-2)=*Q_sca for incident pols 1,2 C *Q_sca for incident pols 1,2 C *Q_sca for incident pols 1,2 C [QSCAG(1-3,1-2) is in Lab Frame] C QTRQAB(1-3,1-2)=vector torque cross section/(PIA2*d^2) for C torque due to absorption, along axes 1-3, for C incident pols 1,2 C [QTRQAB(1-3,1-2) is in Lab Frame] C QTRQSC(1-3,1-2)=vector torque cross section/(PIA2*d^2) for C torque due to scattering, along axes 1-3, for C incident pols 1,2 C [QTRQSC(1-3,1-2) is in Lab Frame] C Requires: C MATVEC = external routine to compute A*x C CMATVEC = external routine to compute conjg(A')*x C where A is matrix such that A*pol=E C and x is arbitrary vector 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 94.06.20 (PJF): Added DTIME to argument list for TIMEIT C 95.06.14 (BTD): Changed CSCAG from scalar to CSCAG(3) vector C where CSCAG(1)=C_sca* C Changed QSCAG(1-2) to QSCAG(1-3,1-2) C Added CTRQSC(1-3) to argument list of SCAT C Added QTRQSC(1-3,1-2) to argument list of GETFML C 95.06.15 (BTD): Added variable CXE to argument list of SCAT C Added CTRQAB(1-3) to argument list of SCAT C Added QTRQAB(1-3,1-2) to argument list of GETFML C 95.06.16 (PJF+: Modified to eliminate call to DDACCG and replace C BTD): with code to call PETR (a modularized implementation C of the Petravic-Kuo CCG algorithm). With this change C it becomes easier to substitute other implementations C of CCG (or other iterative methods, such as QMR). C 95.06.19 (BTD): Added CMDSOL to argument list to select CCG method C : Added CMDTRQ to argument list to allow torque C calculation to be skipped C 95.07.10 (BTD): Changed calls to TIMEIT (was not doing timing C properly) C 95.07.17 (BTD): Relocated call to EVALQ. C Call REDUCE to "reduce" CXE as well as CXP, prior C to calling EVALQ and SCAT. C 95.07.20 (BTD): Added scratch space SCRRS2 for use by SCAT C Added scratch array CXZW(1-MXNAT,10-12) to SCAT C argument list C 95.07.27 (BTD): Revised way in which CXALPH is "reduced" prior to C calls to EVALQ and SCAT. C Revised to compute CXALPH for each incident E field C [Only reason this is needed is because of slight C dependence of polarizabilities on incident pol state] C Added scratch space CXSCR1 to argument list, to store C first polarization solution while solving for second C 95.08.11 (BTD): Relocated CALL CINIT to precede IF statements for C selection of solution method. C 95.08.14 (BTD): Added COMMON/NORMERR/ERRSCAL to support desired C normalization of error printed out by subroutine C PROGRESS (part of cgcommon.f) C 96.11.14 (PJF): ADD TIMERS, MXTIMERS, NTIMERS to formal parameters C modify code to return timer results C TIMERS contains information needed to time the C code (for example CPU time or number of iterations) C timers(1) --- CG CPU time for all iterations C timers(2) --- number of iterations needed C timers(3) --- not used C timers(4) --- scat CPU time C ... see code for timers(5-12) C 96.11.21 (BTD): declared MXTIMERS,NTIMERS, C 97.02.20 (BTD): removed SNORM2 from EXTERNAL statement C 98.11.18 (BTD): corrected misleading comments apropos EM1R,EM2R C 98.12.04 (BTD): located error which resulted in subsequent incorrect C evaluation of Mueller scattering matrix. C Problem arose because at some point it was decided to C convert to usual convention for incident polarization C states 1 and 2 for computation of f_ml (CXF11,...) C Unfortunately, this was forgotten, and subroutine C GETMUELLER was written assuming that incident pol. C states l=1 and 2 referred to by f_ml were the ones C specified by the user in ddscat.par. Because this in C fact does seem preferable, we conform to this C prescription and have disabled the portion of GETFML C which converted to the "usual" convention. Note C that since f_ml is now entirely for internal use, C there is no motivation to adopt usual convention. C 98.12.07 (BTD): Experiment to verify claimed accuracy from PIM C New module has been inserted into code below but C commented out after verifying that reported errors C appear to be close to "true" errors. C However, module can be reactivated by uncommenting C if such verification is desired in future. Search C below for string "98.12.07" to find the 2 modules. C Copyright (C) 1993,1994,1995,1996,1997,1998 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,CMDSOL*6,CMDFFT*6,CMDTRQ*6 C INTEGER ICTHM,IDVOUT,INIT,IORTH,IPHI,IPHIM,ITASK, & LACE,LAXI,LCLM,LGI,LPI,LSC0,LQI, & MXCOMP,MXCXSC,MXMEM,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXSCA, & MXTIMERS,NAT,NAT0,NAT3,NSCAT,NTIMERS,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(3,2), & QTRQAB(3,2),QTRQSC(3,2),SCRRS1(MXN3),SCRRS2(MXNAT), & THETAN(MXSCA), TIMERS(MXTIMERS) C COMPLEX CXADIA(MXN3),CXALPH(MXN3), & CXE(MXN3),CXE01R(3),CXE02R(3),CXEPS(MXCOMP), & CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),CXF22(MXSCA), & CXPOL(MXN3),CXSC(MXCXSC),CXSCR1(MXN3), & CXZC(MXNX+1,MXNY+1,MXNZ+1,6), & CXZW(MXNAT,8*(3+MXMEM)) C C Common: C INTEGER IDVOUT2 REAL ERRSCAL C----------------------------------------------------------------------- COMMON/NORMERR/ERRSCAL,IDVOUT2 C----------------------------------------------------------------------- C C Local variables C COMPLEX CXE0R(3) C INTEGER I,JO,NAT03 C INTEGER IPAR(13) C REAL & CABS,CBKSCA,CEXT,CPHA,CSCA,COSPHI, & DTIME,DUMMY,E02,FALB,SINPHI C REAL CSCAG(3),CTRQAB(3),CTRQSC(3),SPAR(6) C SAVE E02 C EXTERNAL CMATVEC,DIAGL,MATVEC,PCSUM,PROGRESS,PSCNRM2 C IDVOUT2=IDVOUT 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 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 I=1,3 E02=E02+CXE0R(I)*CONJG(CXE0R(I)) ENDDO C C Compute quantity used for "normalizing" error: C ERRSCAL=SQRT(NAT0*E02) 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 (note: this has to be within the loop over directions and C polarizations because Lattice Dispersion Relation polarizabilities C depend on direction of propagation and on polarization state) C CALL ALPHA(CALPHA,CXALPH,CXE0R,CXEPS,ICOMP,AKR,NAT0,NAT3, & MXCOMP,MXN3) C C C*** call EVALA to prepare A matrix elements C CALL EVALA(CXADIA,CXALPH,MXN3,NAT) C C*** first call to timeit: C CALL TIMEIT(' SOLVER',DTIME) C C Version 4c called DDACCG at this point. This call to DDACCG is C replaced with code to C (1) define parameters IPAR(1),IPAR(4),IPAR(10),SPAR(1) C (2) call nuller to "clean" CXE (i.e., set to zero at unoccupied sites). C (3) initialize first guess for cxpol=0 C C Notation: C ipar(1) LDA LDA (Leading dimension of a) C ipar(2) N N (Number of rows/columns of a) C ipar(3) BLKSZ N (Size of block of data; used when data is C partitioned using cyclic mode) C ipar(4) LOCLEN N (Number of elements stored locally; for C sequential=n) C ipar(5) BASISDIM C=basis=10 (Dimension of orthogonal basis, used in C GMRES) C ipar(6) NPROCS -1 (Number of processors) C ipar(7) PROCID -1 (Processor identification) C ipar(8) PRECONTYPE PRET (=1 Type of preconditioning C 0-no preconditioning, 1-left, 2-right, 3-symmetric C ipar(9) STOPTYPE STOPT (Type of stopping criteria used) C ipar(10) MAXIT MAXIT (=int(n/2) Maximum number of iterations allowed) C IPAR(11) = -1 itno (Number of iterations executed) C IPAR(12) = -1 status (On exit of iterative method; 0-converged) C IPAR(13) = -1 steperr C IPAR(1)=MXN3 IPAR(2)=NAT3 IPAR(3)=NAT3 IPAR(4)=NAT3 IPAR(6)=-1 IPAR(7)=-1 IPAR(10)=MIN(300,NAT03) SPAR(1)=TOL C C Clean CXE: C IF(NAT0.LT.NAT)CALL NULLER(CXE,IOCC,MXNAT,MXN3,NAT) C C Initialize CXPOL (first guess) C CALL CINIT(NAT3,CMPLX(0.,0.),CXPOL,1) C C Iterate to improve CXPOL C IF(CMDSOL.EQ.'PETRKP')THEN IPAR(5)=-1 IPAR(8)=1 C C Note: PETR cannot be used with STOPTYPE=2 -- code crashes when C calling STOPCRIT. Stick with STOPTYPE=5 C IPAR(9)=5 CALL PETR(CXPOL,CXE,CXSC,MXN3,IPAR,SPAR,MATVEC,CMATVEC) C C----------------------------------------------------------------------- C 98.12.07 BTD Experiment to verify claimed accuracy from PETR. C This module is superfluous and could be excised without C affecting computed results, thereby eliminating C unnecessary computations. C c CALL MATVEC(CXPOL,CXSC,IPAR) c DUMMY=0. c DO I=1,NAT c IF(IOCC(I).EQ.1)THEN c DO J=3*I+1,3*I+3 c DUMMY=DUMMY+ c & CONJG(CXSC(J)-CXE(J))*(CXSC(J)-CXE(J)) c ENDDO c ENDIF c ENDDO c DUMMY=SQRT(DUMMY/(NAT0*E02)) c WRITE(IDVOUT,9200)DUMMY C----------------------------------------------------------------------- C C second call to timeit: C CALL TIMEIT('PETRKP',DTIME) TIMERS(1)=DTIME TIMERS(2)=FLOAT(IPAR(11)) ELSEIF(CMDSOL.EQ.'PBCGST')THEN IPAR(5)=-1 IPAR(8)=1 C C PIMCBICGSTAB appears to be about 25% slower (per iteration) when C called with STOPTYPE=2 rather than 5. Therefore use STOPTYPE=5 C even though actual stopping criterion used is complicated by C preconditioning. C IPAR(9)=5 CALL PIMCBICGSTAB(CXPOL,CXE,CXSC,IPAR,SPAR,MATVEC, & DIAGL, & DUMMY,PCSUM,PSCNRM2,PROGRESS) C C---------------------------------------------------------------------- C 98.12.07 BTD Experiment to verify claimed accuracy from PIM. C This module is superfluous and could be excised without C affecting computed results, thereby eliminating C unnecessary computations. C c CALL MATVEC(CXPOL,CXSC,IPAR) c DUMMY=0. c DO I=1,NAT c IF(IOCC(I).EQ.1)THEN c DO J=3*I+1,3*I+3 c DUMMY=DUMMY+ c & CONJG(CXSC(J)-CXE(J))*(CXSC(J)-CXE(J)) c ENDDO c ENDIF c ENDDO c DUMMY=SQRT(DUMMY/(NAT0*E02)) c WRITE(0,9200)DUMMY C---------------------------------------------------------------------- C C second call to timeit: C CALL TIMEIT('PBCGST',DTIME) TIMERS(1)=DTIME TIMERS(2)=FLOAT(IPAR(11)) ELSE STOP'ERROR -- INVALID CMDSOL IN GETFML' ENDIF C C*** REDUCE THE VECTORS CXALPH, CXPOL AND CXE TO RETAIN ONLY C OCCUPIED SITES, TO ELIMINATE UNNECESSARY CALCULATIONS FOR C UNOCCUPIED SITES WHEN CALLING EVALQ AND SCAT C CALL REDUCE(CXALPH,IOCC,MXN3,MXNAT,NAT,NAT0) CALL REDUCE(CXE,IOCC,MXN3,MXNAT,NAT,NAT0) CALL REDUCE(CXPOL,IOCC,MXN3,MXNAT,NAT,NAT0) C C*** IF IPHIM.GT.1: STORE REDUCED POLARIZATION VECTOR P C FOR USE IN SUBSEQUENT COMPUTATIONS FOR OTHER PHI VALUES C CXSCR1 FOR JO=1 C CXSC(LACE) FOR JO=2 (THIS SCRATCH SPACE NO LONGER NEEDED BY SOLVER) C IF(IPHIM.GT.1)THEN IF(JO.EQ.1)THEN CALL COPYIT(CXPOL,CXSCR1,NAT03) ELSEIF(JO.EQ.2)THEN CALL COPYIT(CXPOL,CXSC(LACE),NAT03) ENDIF ENDIF C CALL EVALQ(CXALPH,AKR,NAT03,E02,CXE,CXPOL,CABS,CEXT,CPHA, & MXN3,1) QABS(JO)=CABS/PIA2 QEXT(JO)=CEXT/PIA2 QPHA(JO)=CPHA/PIA2 WRITE(IDVOUT,9010)QABS(JO),QEXT(JO) C C*** first call to timeit: C CALL TIMEIT(' SCAT ',DTIME) C IF(JO.EQ.1)THEN CALL SCAT(AKR,AKSR,EM1R,EM2R,E02,CMDTRQ, & CBKSCA,CSCA,CSCAG,CTRQAB,CTRQSC, & CXE,CXE01R,CXF11,CXF21,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7),CXZW(1,10), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,SCRRS2,THETAN,IXYZ0) ELSEIF(JO.EQ.2)THEN CALL SCAT(AKR,AKSR,EM1R,EM2R,E02,CMDTRQ, & CBKSCA,CSCA,CSCAG,CTRQAB,CTRQSC, & CXE,CXE01R,CXF12,CXF22,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7),CXZW(1,10), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,SCRRS2,THETAN,IXYZ0) ENDIF C C*** second call to timeit: C CALL TIMEIT(' SCAT ',DTIME) TIMERS(4)=DTIME 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 Compute QSCAG(1-3,JO) and QTRQSC(1-3,JO) C DO I=1,3 QSCAG(I,JO)=QSCA(JO)*CSCAG(I)/CSCA QTRQAB(I,JO)=CTRQAB(I)/PIA2 QTRQSC(I,JO)=CTRQSC(I)/PIA2 ENDDO ENDDO 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 CXSCR1 contains polarization vector for IPHI=1 and JO=1 C CXSC(LACE) contains polarization vector for IPHI=1 and JO=2 C C Incident polarization state 1 = CXE01R C Call EVALE to obtain E at all lattice sites C first call to timing routine C CALL TIMEIT(' EVALE',DTIME) C C*** Call EVALE for NAT0 occupied sites C CALL EVALE(CXE01R,AKR,IXYZ0,MXNAT,MXN3,NAT0,NAT0,NX,NY,NZ,CXE) C C*** second call to timing routine C CALL TIMEIT(' EVALE',DTIME) TIMERS(5)=DTIME C C Compute polarizabilities at NAT sites, then reduce to NAT0 occupied C sites: C CALL TIMEIT(' ALPHA',DTIME) CALL ALPHA(CALPHA,CXALPH,CXE01R,CXEPS,ICOMP,AKR,NAT0,NAT3, & MXCOMP,MXN3) CALL REDUCE(CXALPH,IOCC,MXN3,MXNAT,NAT,NAT0) CALL TIMEIT(' ALPHA',DTIME) TIMERS(6)=DTIME C C*** Construct polarization for incident pol. state 1 C COSPHI=COS(PHI(IPHI)-PHI(1)) SINPHI=SIN(PHI(IPHI)-PHI(1)) DO I=1,NAT03 CXPOL(I)=COSPHI*CXSCR1(I)-SINPHI*CXSC(LACE+I-1) ENDDO C C*** first call to timing routine C CALL TIMEIT(' EVALQ',DTIME) C CALL EVALQ(CXALPH,AKR,NAT03,E02,CXE,CXPOL,CABS,CEXT,CPHA, & MXN3,1) C C*** second call to timing routine C CALL TIMEIT(' EVALQ',DTIME) TIMERS(7)=DTIME C QABS(1)=CABS/PIA2 QEXT(1)=CEXT/PIA2 QPHA(1)=CPHA/PIA2 C C*** first call to timing routine C CALL TIMEIT(' SCAT',DTIME) C C*** Now call SCAT C CXPOL has been reduced to NAT03 elements C first NAT03 elements of IXYZ0 correspond to physical sites C CALL SCAT(AKR,AKSR,EM1R,EM2R,E02,CMDTRQ, & CBKSCA,CSCA,CSCAG,CTRQAB,CTRQSC, & CXE,CXE01R,CXF11,CXF21,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7),CXZW(1,10), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,SCRRS2,THETAN,IXYZ0) C C*** second call to timing routine C CALL TIMEIT(' SCAT',DTIME) TIMERS(8)=DTIME 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) DO I=1,3 QSCAG(I,1)=QSCA(1)*CSCAG(I)/CSCA QTRQAB(I,1)=CTRQAB(I)/PIA2 QTRQSC(I,1)=CTRQSC(I)/PIA2 ENDDO C C*** Polarization state 2: C C*** first call to timing routine C CALL TIMEIT(' EVALE',DTIME) C C*** Call EVALE for NAT0 occupied sites C to obtain appropriate E vector for pol. state 2 C CALL EVALE(CXE02R,AKR,IXYZ0,MXNAT,MXN3,NAT0,NAT0,NX,NY,NZ,CXE) C C*** second call to timing routine C CALL TIMEIT(' EVALE',DTIME) TIMERS(9)=DTIME C C Compute polarizabilities at NAT sites, then reduce to NAT0 occupied C sites: C CALL TIMEIT(' ALPHA',DTIME) CALL ALPHA(CALPHA,CXALPH,CXE02R,CXEPS,ICOMP,AKR,NAT0,NAT3, & MXCOMP,MXN3) CALL REDUCE(CXALPH,IOCC,MXN3,MXNAT,NAT,NAT0) CALL TIMEIT(' ALPHA',DTIME) TIMERS(10)=DTIME C C*** Construct incident polarization vector for pol. state 2 C DO I=1,NAT03 CXPOL(I)=COSPHI*CXSC(LACE+I-1)+SINPHI*CXSCR1(I) ENDDO C C*** first call to timing routine C CALL TIMEIT(' EVALQ ',DTIME) C CALL EVALQ(CXALPH,AKR,NAT03,E02,CXE,CXPOL,CABS,CEXT,CPHA, & MXN3,1) C C*** second call to timing routine C CALL TIMEIT(' EVALQ ',DTIME) TIMERS(11)=DTIME C QABS(2)=CABS/PIA2 QEXT(2)=CEXT/PIA2 QPHA(2)=CPHA/PIA2 C C*** first call to timing routine C CALL TIMEIT(' SCAT',DTIME) C C*** Call SCAT C CXPOL has NAT03 elements C IXYZ0 was initially given NAT03 elements C CALL SCAT(AKR,AKSR,EM1R,EM2R,E02,CMDTRQ, & CBKSCA,CSCA,CSCAG,CTRQAB,CTRQSC, & CXE,CXE01R,CXF12,CXF22,CXPOL, & CXZW(1,1),CXZW(1,4),CXZW(1,7),CXZW(1,10), & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT0,NAT03,NSCAT, & PHIN,SCRRS1,SCRRS2,THETAN,IXYZ0) C C*** second call to timing routine C CALL TIMEIT(' SCAT',DTIME) TIMERS(12)=DTIME C QBKSCA(2)=CBKSCA/PIA2 QSCA(2)=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(2)/(.03*QEXT(2)))**2 QSCA(2)=(QSCA(2)+(QEXT(2)-QABS(2))*FALB)/(1.+FALB) DO I=1,3 QSCAG(I,2)=QSCA(2)*CSCAG(I)/CSCA QTRQAB(I,2)=CTRQAB(I)/PIA2 QTRQSC(I,2)=CTRQSC(I)/PIA2 ENDDO ENDIF C RETURN 9010 FORMAT(1X,'Q_abs =',1PE10.3,' Q_ext= ',E10.3) 9200 FORMAT(' final true frac.err=',F12.7) END SUBROUTINE GETMUELLER(IBETA,IORTH,IPHI,ITHETA,MXBETA, & MXSCA,MXPHI,MXTHET,NSCAT, & CMDTRQ,PHIN,SM,SMORI, & S1111,S2121, & CX1121,CXE01,CXE02, & CXF11,CXF12,CXF21,CXF22, & CXS1,CXS2,CXS3,CXS4, & QABS,QABSUM,QBKSCA,QBKSUM,QEXSUM, & QEXT,QPHA,QPHSUM,QSCAG,QSCAT, & QSCGSUM,QSCSUM,QTRQAB,QTRQABSUM, & QTRQSC,QTRQSCSUM,WGTA,WGTB) C Arguments: INTEGER & IBETA,IORTH,IPHI,ITHETA, & MXBETA,MXSCA,MXPHI,MXTHET, & NSCAT CHARACTER & CMDTRQ*6 REAL WG REAL & PHIN(MXSCA),QABS(2),QABSUM(2),QBKSCA(2),QBKSUM(2),QEXSUM(2), & QEXT(2),QPHA(2),QPHSUM(2),QSCAG(3,2),QSCAT(2),QSCGSUM(3,2), & QSCSUM(2),QTRQAB(3,2),QTRQABSUM(3,2),QTRQSC(3,2), & QTRQSCSUM(3,2),S1111(MXSCA),S2121(MXSCA), & SM(MXSCA,4,4),SMORI(MXSCA,4,4), & WGTA(MXTHET,MXPHI),WGTB(MXBETA) COMPLEX & CX1121(MXSCA),CXE01(3),CXE02(3), & CXF11(MXSCA),CXF12(MXSCA),CXF21(MXSCA),CXF22(MXSCA), & CXS1(MXSCA),CXS2(MXSCA),CXS3(MXSCA),CXS4(MXSCA) C C Local variables: INTEGER J,K,JO,ND REAL COSPHI,SINPHI COMPLEX CXA,CXB,CXC,CXD,CXI DATA CXI/(0.,1.)/ SAVE CXI C*********************************************************************** C Given: C IBETA =index specifying rotation of target around axis a1 C IORTH =number of incident polarization states being calculated for C IPHI =index specifying phi value for orientation of target axis a1 C ITHETA=index specifying theta value for orientation of target axsis C a1 C MXBETA=dimensioning information C MXSCA =dimensioning information C MXPHI =dimensioning information C MXTHET=dimensioning information C NSCAT =number of scattering directions for evaluation of Mueller C matrix C CMDTRQ='DOTORQ' or 'NOTORQ' depending on whether or not torque C calculation is to be carried out C PHIN(1-NSCAT)=values of PHI (rad) for scattering directions 1-NSCAT C S1111(1-NSCAT)=weighted sum of |f_11|^2 over previous orientations C S2121(1-NSCAT)=weighted sum of |f_21|^2 over previous orientations C CXE01(1-3) =incident pol state 1 in Lab Frame C CXE02(1-3) = 2 C CXF11(1-NSCAT)=f_11 for current orientation C CXF12(1-NSCAT)=f_12 for current orientation C CXF21(1-NSCAT)=f_21 for current orientation C CXF22(1-NSCAT)=f_22 for current orientation C QABS(2) =Q_abs for two incident polarizations C QABSUM(2) =weighted sum of Qabs over previous orientations C QBKSCA(2) =Q_bk for two incident polarizations C QBKSUM(2) =weighted sum of Q_bk over previous orientations C QEXSUM(2) =weighted sum of Q_ext over previous orientations C QEXT(2) =Q_ext for two incident polarizations C QPHA(2) =Q_pha for two incident polarizations C QPHSUM(2) =weighted sum of Q_ph over previous orientations C QSCAG(3,2) =Q_pr vector for two incident polarizations C QSCAT(2) =Q_sca for two incident polarizations C QSCGSUM(3,2) =weighted sum of Q_pr vector over prev. orients. C QSCSUM(2) =weighted sum of Q_sca over previous orientations C QTRQAB(3,2) =absorptive part of Q_gamma for two inc. pols. C QTRQABSUM(3,2)=weighted sum of Q_gamma over previous orientations C QTRQSC(3,2) =scattering part of Q_gamma for two inc. pols. C QTRQSCSUM(3,2)=weighted sum of Q_gamma over previous orientations C WGTA(ITHETA,IPHI)=weight factor for theta,phi orientation C WGTB(IBETA) =weight factor for beta orientation C SMORI(1-NSCAT,4,4)=weighted sum of Mueller matrix over previous C orientations C If IORTH=1, returns: C S1111(1-NSCAT)=updated weighted sum of |f_11|^2 over NSCAT C scattering directions C S2121(1-NSCAT)=updated weighted sum of |f_21|^2 C CX1121(1-NSCAT)=updated weighted summ of f_11*conjg(f_21) C C If IORTH=2, returns: C SM(1-NSCAT,4,4)=Mueller matrix for current orientation, and NSCAT C scattering directions C CXS1(1-NSCAT)=amplitude scattering matrix element S_1 for current C orientation C CXS2(1-NSCAT)=amplitude scattering matrix element S_2 for current C orientation C CXS3(1-NSCAT)=amplitude scattering matrix element S_3 for current C orientation C CXS4(1-NSCAT)=amplitude scattering matrix element S_4 for current C orientations C QABSUM(2) =updated weighted sum of Q_abs C QBKSUM(2) =updated weighted sum of Q_bk C QEXSUM(2) =updated weighted sum of Q_ext C QSCGSUM(3,2) =updated weighted sum of Q_pr vector C QSCSUM(2) =updated weighted sum of Q_sca C QTRQABSUM(3,2)=updated weighted sum of absorptive contribution to C Q_gamma C QTRQSCSUM(3,2)=updated weighted sum of scattering contribution to C Q_gamma C C SMORI(1-NSCAT,4,4)=updated weighted sum of Mueller matrix elements C C History: C 96.11.11 (BTD) ready for use C 97.07.24 (BTD) corrected error in computation of amplitude scattering C matrix elements (CXS1,CXS2,CXS3,CXS4) and Mueller matrix C elements SM. Until now code mistakenly used PHI C defining target orientation rather than PHIN C defining scattering direction. Correction of error C required replacing PHI by PHIN in argument list, C thereby requiring corresponding change in DDSCAT. C All calculated Mueller matrix elements until now were C therefore incorrect except for special cases where C PHI and PHIN were coincidentally the same (e.g. C target orientation PHI=0 and scattering plane PHIN=0). C This error was pointed out by Henriette Lemke (Inst. C for Atmospheric Physics, GKSS - Research Center, C Geesthacht. C 97.07.24 (BTD) also corrected sign mistake in evaluation of S_14 C and typo in evaluation of S_31. Note that since S_31 C had been evaluated incorrectly, PPOL (evaluated by C subroutine WRITESCA) was then numerically incorrect. C 98.11.01 (BTD) Timo Nousiainen (tpnousia@lumi.meteo.helsinki.fi) C reported a problem with S_12 computed by DDSCAT for C a pseudosphere target when the scattering plane has C e.g. phi=90. Problem is real and must be tracked down C and corrected. C 98.11.16 (BTD) Modified code calculating complex amplitude scattering C matrix elements S1,S2,S3,S4 from f_ij. C end history C C********* Sum scattering properties over orientations **************** C C Note: the equations used to compute the scattering amplitude C matrix elements S_1,S_2,S_3,S_4 from the f_ml are derived and C discussed in section 22.1 of the "UserGuide for the Discrete Dipole C Approximation Code DDSCAT (Version 5a9)", by B.T. Draine and P.J. Flatau. C C Correspondence to notation of Bohren & Huffman (1983): C Our f_jk correspond to notation of Draine (1988). C For the special case where incident polarization modes 1 and 2 are C parallel and perpendicular to the scattering plane, we have the C relationship between the f_jk and the elements S_j of the amplitude C scattering matrix as defined by Bohren & Huffman (1983) is as follows: C S_1 = -i*f_22 C S_2 = -i*f_11 C S_3 = i*f_12 C S_4 = i*f_21 C In this case (incident pol. states parallel,perpendicular to the C scattering plane), the elements of the 4x4 Mueller matrix (see Bohren & C Huffman 1983) are related to our f_jk as follows: C S_11 = [|f_11|^2 + |f_22|^2 + |f_12|^2 + |f_21|^2]/2 C S_12 = [|f_11|^2 + |f_21|^2 - |f_22|^2 - |f_12|^2]/2 C S_13 = -Re[f_11*conjg(f_12) + f_22*conjg(f_21)] C S_14 = Im[f_22*conjg(f_21) - f_11*conjg(f_12)] C S_21 = [|f_11|^2 + |f_12|^2 - |f_22|^2 - |f_21|^2]/2 C S_22 = [|f_11|^2 + |f_22|^2 - |f_21|^2 - |f_12|^2]/2 C S_23 = Re[f_22*conjg(f_21) - f_11*conjg(f_12)] C S_24 = -Im[f_11*conjg(f_12) + f_22*conjg(f_21)] C S_31 = -Re[f_11*conjg(f_21) + f_22*conjg(f_22)] C S_32 = Re[f_22*conjg(f_12) - f_22*conjg(f_21)] C S_33 = Re[f_22*conjg(f_11) + f_12*conjg(f_21)] C S_34 = Im[f_11*conjg(f_22) + f_21*conjg(f_22)] C S_41 = -Im[f_21*conjg(f_11) + f_22*conjg(f_12)] C S_42 = Im[f_22*conjg(f_12) - f_22*conjg(f_11)] C S_43 = Im[f_22*conjg(f_11) - f_12*conjg(f_21)] C S_44 = Re[f_22*conjg(f_11) - f_12*conjg(f_21)] C C Notation internal to this program: C ij.ne.kl: CX_ijkl = \f_ij* \times f_kl (summed over orientations) C S_ijij = \f_ij* \times f_ij (summed over orientations) C 4 pure real elements: S1111,S1212,S2121,S2222 C 8 complex elements: CX1112,CX1121,CX1122,CX1221,CX1222,CX2122 C 4 redundant elts.: CX1211 = Conjg(CX1112) C CX2111 = Conjg(CX1121) C CX2112 = Conjg(CX1221) C CX2211 = Conjg(CX1122) C C In this routine we will carry out fully general calculation of the C 4x4 Mueller matrix elements, valid for arbitrary incident polarization C states j=1,2 for which we have the scattering matrix f_ij C C First, we compute the amplitude scattering matrix elements S_1,S_2, C S_3,S_4 as defined by Bohren & Huffman: C C WG=weighting factor for this orientation C WG=WGTA(ITHETA,IPHI)*WGTB(IBETA) IF(IORTH.EQ.2)THEN C C CXE01,CXE02 are unit incident polarization vectors in Lab Frame. C We assume CXE01 and CXE02 are normalized (REAPAR makes sure of this). C CXA,CXB,CXC,CXD are complex coefficients corresponding to dot product C of complex polarization vectors CXE01,CXE02 with y and z unit vectors C in Lab Frame. C CXA=CONJG(CXE01(2)) CXB=CONJG(CXE01(3)) CXC=CONJG(CXE02(2)) CXD=CONJG(CXE02(3)) C C Compute complex scattering amplitudes S_1,S_2,S_3,S_4 for this C particular target orientation and NSCAT scattering directions: C DO ND=1,NSCAT SINPHI=SIN(PHIN(ND)) COSPHI=COS(PHIN(ND)) CXS1(ND)= CXI*(CXF21(ND)*(CXA*SINPHI-CXB*COSPHI)+ & CXF22(ND)*(CXC*SINPHI-CXD*COSPHI)) CXS2(ND)=-CXI*(CXF11(ND)*(CXB*SINPHI+CXA*COSPHI)+ & CXF12(ND)*(CXD*SINPHI+CXC*COSPHI)) CXS3(ND)=-CXI*(CXF11(ND)*(CXA*SINPHI-CXB*COSPHI)+ & CXF12(ND)*(CXC*SINPHI-CXD*COSPHI)) CXS4(ND)= CXI*(CXF21(ND)*(CXB*SINPHI+CXA*COSPHI)+ & CXF22(ND)*(CXD*SINPHI+CXC*COSPHI)) ENDDO c*** diagnostic c write(0,*)'in getmueller: scattering direction 1' c write(0,*)'S_1=',cxs1(1) c write(0,*)'S_2=',cxs2(1) c write(0,*)'S_3=',cxs3(1) c write(0,*)'S_4=',cxs4(1) c*** C C Now compute Mueller matrix elements for this particular target C orientation: C DO ND=1,NSCAT SM(ND,1,1)=0.5*(CXS1(ND)*CONJG(CXS1(ND))+ & CXS2(ND)*CONJG(CXS2(ND))+ & CXS3(ND)*CONJG(CXS3(ND))+ & CXS4(ND)*CONJG(CXS4(ND))) SM(ND,1,2)=0.5*(CXS2(ND)*CONJG(CXS2(ND))- & CXS1(ND)*CONJG(CXS1(ND))+ & CXS4(ND)*CONJG(CXS4(ND))- & CXS3(ND)*CONJG(CXS3(ND))) SM(ND,1,3)=REAL(CXS2(ND)*CONJG(CXS3(ND))+ & CXS1(ND)*CONJG(CXS4(ND))) SM(ND,1,4)=AIMAG(CXS2(ND)*CONJG(CXS3(ND))- & CXS1(ND)*CONJG(CXS4(ND))) SM(ND,2,1)=0.5*(CXS2(ND)*CONJG(CXS2(ND))- & CXS1(ND)*CONJG(CXS1(ND))- & CXS4(ND)*CONJG(CXS4(ND))+ & CXS3(ND)*CONJG(CXS3(ND))) SM(ND,2,2)=0.5*(CXS2(ND)*CONJG(CXS2(ND))+ & CXS1(ND)*CONJG(CXS1(ND))- & CXS4(ND)*CONJG(CXS4(ND))- & CXS3(ND)*CONJG(CXS3(ND))) SM(ND,2,3)=REAL(CXS2(ND)*CONJG(CXS3(ND))- & CXS1(ND)*CONJG(CXS4(ND))) SM(ND,2,4)=AIMAG(CXS2(ND)*CONJG(CXS3(ND))+ & CXS1(ND)*CONJG(CXS4(ND))) SM(ND,3,1)=REAL(CXS2(ND)*CONJG(CXS4(ND))+ & CXS1(ND)*CONJG(CXS3(ND))) SM(ND,3,2)=REAL(CXS2(ND)*CONJG(CXS4(ND))- & CXS1(ND)*CONJG(CXS3(ND))) SM(ND,3,3)=REAL(CXS1(ND)*CONJG(CXS2(ND))+ & CXS3(ND)*CONJG(CXS4(ND))) SM(ND,3,4)=AIMAG(CXS2(ND)*CONJG(CXS1(ND))+ & CXS4(ND)*CONJG(CXS3(ND))) SM(ND,4,1)=AIMAG(CXS4(ND)*CONJG(CXS2(ND))+ & CXS1(ND)*CONJG(CXS3(ND))) SM(ND,4,2)=AIMAG(CXS4(ND)*CONJG(CXS2(ND))- & CXS1(ND)*CONJG(CXS3(ND))) SM(ND,4,3)=AIMAG(CXS1(ND)*CONJG(CXS2(ND))- & CXS3(ND)*CONJG(CXS4(ND))) SM(ND,4,4)=REAL(CXS1(ND)*CONJG(CXS2(ND))- & CXS3(ND)*CONJG(CXS4(ND))) ENDDO C C Augment Mueller matrix for orientational average C DO J=1,4 DO K=1,4 DO ND=1,NSCAT SMORI(ND,J,K)=SMORI(ND,J,K)+WG*SM(ND,J,K) ENDDO ENDDO ENDDO C*********************************************************************** C ELSE C C When IORTH=1, cannot compute full Mueller matrix. C Compute three scattering properties: C DO ND=1,NSCAT S1111(ND)=S1111(ND)+WG*CONJG(CXF11(ND))*CXF11(ND) S2121(ND)=S2121(ND)+WG*CONJG(CXF21(ND))*CXF21(ND) CX1121(ND)=CX1121(ND)+WG*CONJG(CXF11(ND))*CXF21(ND) ENDDO ENDIF C C**** Now augment sums of QABS,QBKSCA,QEXT,QSCA,G*QSCA over incident C directions and polarizations. C DO JO=1,IORTH QEXSUM(JO)=QEXSUM(JO)+QEXT(JO)*WG QABSUM(JO)=QABSUM(JO)+QABS(JO)*WG QBKSUM(JO)=QBKSUM(JO)+QBKSCA(JO)*WG QPHSUM(JO)=QPHSUM(JO)+QPHA(JO)*WG QSCSUM(JO)=QSCSUM(JO)+QSCAT(JO)*WG DO J=1,3 QSCGSUM(J,JO)=QSCGSUM(J,JO)+QSCAG(J,JO)*WG ENDDO IF(CMDTRQ.EQ.'DOTORQ')THEN DO J=1,3 QTRQABSUM(J,JO)=QTRQABSUM(J,JO)+QTRQAB(J,JO)*WG QTRQSCSUM(J,JO)=QTRQSCSUM(J,JO)+QTRQSC(J,JO)*WG ENDDO ENDIF ENDDO RETURN END