SUBROUTINE CXFFT3(CXDATA,MX,MY,MZ,ISIGN,CXWORK) PARAMETER (MXTRIG=1024) C*** Arguments: INTEGER ISIGN,MX,MY,MZ COMPLEX CXDATA(MX,MY,MZ),CXWORK(MX,MY,MZ) C C*** Local variables: INTEGER ICALL,INC,JUMP,K,LOT CHARACTER CMSGNM*70 INTEGER IFAX(13),IFAY(13),IFAZ(13) REAL TRIG(MXTRIG) C C*** External Subroutines: EXTERNAL CFFT99,CFTFAX C C*** SAVE statements: SAVE ICALL C C*** DATA statements: DATA ICALL/1/ C*********************************************************************** C C Purpose: C Calculates direct and inverse 3-dimensional FFT (in place) C for complex matrix CXDATA. C C Input: C C CXDATA(NX, NY, NZ) --- data array (complex). Will be overwritten. C MX, MY, MZ --- dimensions of CXDATA; have to be of the form C 2**a 3**b 5**c C ISIGN = (-1,1) --- (1 for inverse, -1 direct transform). C Read documentation of "cfft99f.f" C code for definition C TRIGX(2*MX) C TRIGY(2*MY) C TRIGZ(2*MZ) --- storage arrays. CAN NOT BE OVERWRITTEN C between calls C CXWORK(MX,MY,MZ) --- complex work array. Can be overwritten. C C Output: C C CXDATA --- 3D Fourier transform C TRIGX, TRIGY, TRIGZ C C Local variables: C C IFAX(13) C IFAY(13) C IFAZ(13) --- internal integer info array. Dimension 13 C is sufficient for very large nx, ny, nz C ICALL --- switch. On first entry=1, ofterwards=2 C INC C JUMP C LOT --- define how many FFTs C (see "cfft99.f" documentation) C C External documentation for "cfft99f" C C THE fft99 PACKAGE WAS WRITTEN BY C CLIVE TEMPERTON AT ECMWF IN C NOVEMBER, 1978. IT WAS MODIFIED, DOCUMENTED, AND TESTED C FOR NCAR BY RUSS REW IN SEPTEMBER, 1980. IT WAS C FURTHER MODIFIED FOR THE FULLY COMPLEX CASE BY DAVE C FULKER IN NOVEMBER, 1980. C C Efficiency: C C The code vectorizes on CRAY-YMP and is more efficient than the C FOURN routine from "Numerical Recipes" by W. H. Press, B. P. Flannery C S. A. Teukolsky, W. T. Vettering. For 64x64x32 case this code C is >6 faster than FOURN. C C If speed is important for even values C of n=nx=ny=nz one should insert dummy C rows and columns to make dimension cxdata(n+1, n+1, n+1) to C avoid bank conflicts (on CRAY). C History: C Written by P.J.Flatau, 11/24/1990 with the help of C. Temperton 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 initialize trigs (only once) IF(ICALL.EQ.1)THEN NTRIG=2*(MX+MY+MZ) IF(NTRIG.GT.MXTRIG)THEN WRITE(CMSGNM,FMT='(2I5,A)')NTRIG,MXTRIG, & ' 2*NX+2*NY+2*NZ .GT. MXTRIG. Increase MXTRIG! ' CALL ERRMSG('FATAL','CXFFT3',CMSGNM) ENDIF LTRIGX=1 LTRIGY=1+2*MX LTRIGZ=1+2*MX+2*MY CALL CFTFAX(MX,IFAX,TRIG(LTRIGX)) CALL CFTFAX(MY,IFAY,TRIG(LTRIGY)) CALL CFTFAX(MZ,IFAZ,TRIG(LTRIGZ)) ICALL=2 ENDIF C --- first dimension INC=1 JUMP=MX LOT=MY*MZ CALL CFFT99(CXDATA,CXWORK,TRIG(LTRIGX),IFAX,INC,JUMP,MX,LOT,ISIGN) C --- second dimension INC=MX JUMP=1 LOT=MX C --- one plane at a time DO 10 K=1,MZ CALL CFFT99(CXDATA(1,1,K),CXWORK,TRIG(LTRIGY),IFAY,INC,JUMP, & MY,LOT,ISIGN) 10 CONTINUE C --- third dimension INC=MX*MY JUMP=1 LOT=MX*MY CALL CFFT99(CXDATA,CXWORK,TRIG(LTRIGZ),IFAZ,INC,JUMP,MZ,LOT,ISIGN) RETURN END SUBROUTINE DDACCG(AK,AK3,CMETHD, & CXACE,CXADIA,CXALPH,CXAXI, & CXCLMN,CXE,CXE00,CXGI,CXQI, & CXPI,CXXI, & CXZC,CXZW, & IDVOUT,INIT,IOCC, & MXMEM,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ, & PIA2,QABS,QEXT,QPHA,QSCA, & SCRRS1,TOL) C*********************************************************************** C Given: C NAT=number of point dipoles C NAT3=3*NAT C NX,NY,NZ=x,y,z length of extended target C MXN3, MXNAT, MXNX, MXNY, MXNZ =dimensioning information C CXALPH(1-NAT3)=dipole polarizabilities C CXE00(1-3)=incident E field at origin C CXE(1-NAT3)=incident E field at dipole locations C AK(1-3)=d*(incident k vector), where d=lattice spacing C AK3=(kd)**3 C IDVOUT=device number for running output C INIT=parameter to select initial guess |X0> for solution C TOL=fractional error tolerance C PIA2=pi*(a_{eff}/d)**2, where d=lattice spacing C CXADIA(1-NAT3)=diagonal elements of A matrix C Returns: C CXXI(1-NAT3)=solution polarization vector C QEXT=Cext/pi*aeff**2 C QABS=Cabs/pi*aeff**2 C QSCA=Csca/pi*aeff**2 (computed as (QEXT-QABS) ) C Scratch space: C CXACE,CXAXI,CXCLMN,CXGI,CXQI,CXPI,SCRRS1 C C Subroutine DDACCG proceeds by iteration to solve the self-consistent C grain polarization for an array of point dipoles, using the a complex C conjugate gradient algorithm. C C B. T. Draine, Princeton Univ. Obs., 87/1/7 C History: C 88/ 6/30 (BTD): Modified C 89/11/27 (BTD): Modified C 90/ 9/13 (BTD): Modified to experiment with use of FFT method C proposed by J.J.Goodman. C 90/11/ 2 (BTD): Modified to treat nonrectangular targets with FFT C (necessary to deal with "vacuum" sites) C 90/11/ 4 (BTD): Major modification: transfer both FFT and MVM C computations to subroutine CPROD. Also transfer C requisite dimensioning information for use by CPROD. C 90/11/21 (BTD): Remove CXAOFF and CXCLMN from calls to CPROD. C 90/11/30 (BTD): Changes to reduce storage requirements. C 90/12/03 (BTD): Reorder CXZW and use CXZW to replace CXSCR0 C Remove CXSCR0 from calling sequence C 90/12/13 (BTD): Add new argument to calls to EVALQ C 93/05/17 (BTD): Moved call to 'NULLER' so that code works with option C INIT=1 as well as INIT=0 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 We seek to solve the complex matrix equation A*X=E C where X(1-3*NAT) is the unknown polarization vector for the dipoles C E(1-3*NAT) is the incident electric field vector at dipoles C and the complex matrix A describes the coupling of the dipoles to C the incident field and to one another. C This subroutine implements the complex conjugate gradient method C described by Petravic and Kuo-Petravic (1979; J. Comp. Phys. 32, 263) C with the particular choice H=K=identity matrix C in which case N=AC*A, where AC=conjugate(A) C Notation follows Petravic and Kuo-Petravic, except C that we define an additional vector |Q>=A|P>: C |GI>=AC|E>-ACA|XI> C BETAI-1=/ C |PI>=|GI>+BETAI-1*|PI-1> C |QI>=A|PI> C ALPHAI=/ C |XI+1>=|XI>+ALPHAI*|PI> C C Initialization: C C If INIT=0, take |X0>=0 C If INIT=1, compute |X0> to fourth order in polarizability alpha C C*********************************************************************** C*** Arguments: CHARACTER*6 CMETHD REAL AK3,TOL,PIA2,QABS,QEXT,QPHA,QSCA INTEGER IDVOUT,INIT,MXMEM,MXNAT,MXN3,MXNX,MXNY,MXNZ,NAT, & NAT3,NAT0 INTEGER NX,NY,NZ INTEGER*2 IOCC(MXNAT) COMPLEX CXACE(MXN3),CXALPH(MXN3),CXAXI(MXN3),CXCLMN(MXN3), & CXE(MXN3),CXE00(3),CXGI(MXN3),CXPI(MXN3), & CXQI(MXN3),CXXI(MXN3), & CXADIA(MXN3), & CXZC(MXNX+1,MXNY+1,MXNZ+1,6),CXZW(MXNAT,8*(3+MXMEM)) REAL AK(3),SCRRS1(MXN3) C C*** Local variables: COMPLEX CXALP0 REAL ALPHAI,BETAI1,CABS,CEXT,CPHA,E02,EMAG2,FRAC,FRACI1, & FRACI2,GI1GI1,GIGI,QABI1,QABI2,QEXEST,QEXI1,QEXI2, & QIQI,QPHI1,QPHI2 INTEGER IDVMSG,IORD,IT,ITER,ITMAX,J C C*** External Subroutines: EXTERNAL CPROD,EVALQ C C*** Intrinsic Functions: INTRINSIC ABS,SQRT,CONJG C C*** Define output device number for messages: DATA IDVMSG/0/ C C*********************************************************************** C*** Compute magnitude of incident electric field E02=0. DO 10 J=1,3 E02=E02+CXE00(J)*CONJG(CXE00(J)) 10 CONTINUE C C Clean CXE: IF(NAT0.LT.NAT)CALL NULLER(CXE,IOCC,MXNAT,MXN3,NAT) C ITER=0 IF(INIT.EQ.0)GOTO 160 C C*** Compute |X0> to fourth order in polarizability alpha C WRITE(IDVOUT,FMT=8991) C C*** Compute GI(J)=CXALPHA(J)*EINC(J) DO 40 J=1,NAT3 CXGI(J)=CXALPH(J)*CXE(J) C To first order, XI=GI CXXI(J)=CXGI(J) 40 CONTINUE CALL EVALQ(CXALPH,AK,NAT3,E02,CXE,CXXI,CABS,CEXT,CPHA,MXN3,1) QABS=CABS/PIA2 QEXT=CEXT/PIA2 QPHA=CPHA/PIA2 IORD=1 WRITE(IDVOUT,FMT=9000)IORD,QEXT,QABS,QPHA C C*** Compute PI(J)=ALPHA(J)*A(J,K)*GI(K) (sum over k) C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXGI,CXPI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) DO 50 J=1,NAT3 CXPI(J)=CXALPH(J)*CXPI(J) CXXI(J)=2.*CXGI(J)-CXPI(J) 50 CONTINUE C CALL EVALQ(CXALPH,AK,NAT3,E02,CXE,CXXI,CABS,CEXT,CPHA,MXN3,1) QABS=CABS/PIA2 QEXT=CEXT/PIA2 QPHA=CPHA/PIA2 IORD=2 WRITE (IDVOUT,FMT=9000) IORD,QEXT,QABS,QPHA C C*** Compute QI(J)=ALPHA(J)*A(J,K)*PI(K) (sum over k) C*** and XI(J)=3*G(J)-3*P(J)+ALPHA(J)*Q(J) + order(alpha**4) C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXPI,CXQI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) DO 60 J=1,NAT3 CXQI(J)=CXALPH(J)*CXQI(J) CXXI(J)=3.* (CXGI(J)-CXPI(J)) + CXQI(J) 60 CONTINUE C CALL EVALQ(CXALPH,AK,NAT3,E02,CXE,CXXI,CABS,CEXT,CPHA,MXN3,1) QABS=CABS/PIA2 QEXT=CEXT/PIA2 QPHA=CPHA/PIA2 IORD=3 WRITE (IDVOUT,FMT=9000) IORD,QEXT,QABS,QPHA C C*** Compute ACE(J)=ALPHA(J)*A(J,K)*QI(K) (sum over k) C*** and XI(J)=4.D0*GI(J)-6.D0*PI(J)+4.D0*QI(J)-ACE(J) + order(alpha**5) C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXQI,CXACE,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) DO 70 J=1,NAT3 CXXI(J)=4.*CXGI(J) - 6.*CXPI(J) + 4.*CXQI(J) - & CXALPH(J)*CXACE(J) 70 CONTINUE C C Given guess |X0>, compute: C |G0>=AC|E>-ACA|X0> C |P0>=|G0> C |Q0>=A|P0> C C nonzero |X0>: compute A|X0>, EMAG2, and FRAC=fractional error 80 CONTINUE ITER=0 EMAG2=0. DO 90 J=1,NAT3 EMAG2=EMAG2 + CXE(J)*CONJG(CXE(J)) 90 CONTINUE CALL CPROD(AK,CMETHD,'N',CXADIA, & CXXI,CXAXI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) DO 100 J=1,NAT3 CXZW(J,1)=CXAXI(J) - CXE(J) SCRRS1(J)=CXZW(J,1)*CONJG(CXZW(J,1)) 100 CONTINUE FRAC=0. DO 110 J=1,NAT3 FRAC=FRAC + SCRRS1(J) 110 CONTINUE C FRAC=SQRT(FRAC/EMAG2) CALL EVALQ(CXALPH,AK,NAT3,E02,CXE,CXXI,CABS,CEXT,CPHA,MXN3,1) QABS=CABS/PIA2 QEXT=CEXT/PIA2 QPHA=CPHA/PIA2 IORD=4 WRITE(IDVOUT,FMT=9010)IORD,FRAC,QEXT,QABS,QPHA GOTO 130 C 120 CONTINUE WRITE(IDVOUT,FMT=9020)ITER,FRAC,QEXT,QABS,QPHA C C Check whether result satisfies error criterion: C 130 CONTINUE IF(FRAC.GT.TOL)GOTO 140 QSCA=QEXT-QABS WRITE(IDVOUT,FMT=9050)QEXT,QABS,QPHA,QSCA GOTO 350 C C Need to iterate to reduce error C C Compute AC|E> C |G0>=AC|E>-ACA|X0> C |P0>=|G0> 140 CONTINUE CALL CPROD(AK,CMETHD,'C',CXADIA, & CXE,CXACE,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) CALL CPROD(AK,CMETHD,'C',CXADIA, & CXAXI,CXGI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) DO 150 J=1,NAT3 CXGI(J)=-CXGI(J) + CXACE(J) CXPI(J)=CXGI(J) 150 CONTINUE C GOTO 180 C C |X0>=0: 160 CONTINUE C*** Note: all the calls to timeit in this routine are now disabled. C If they are desired (in order to cleanly get single-iteration C timing) it is necessary to disable the calls to timeit in getfml. C*** first call to timeit: C CALL TIMEIT(' A') C WRITE (IDVOUT,FMT=8990) C C*** Set |X0>=0 DO 165 J=1,NAT3 CXXI(J)=0. 165 CONTINUE C C Compute AC|E> C EMAG2 C |G0>=AC|E> C |P0>=|G0> C CALL CPROD(AK,CMETHD,'C',CXADIA, & CXE,CXACE,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C EMAG2=0. DO 170 J=1,NAT3 EMAG2=EMAG2 + CXE(J)*CONJG(CXE(J)) CXGI(J)=CXACE(J) CXPI(J)=CXGI(J) 170 CONTINUE C QABS=0. QEXT=0. QPHA=0. FRAC=1. C*** second call to timeit C CALL TIMEIT('DDACCG.0') C*** WRITE(IDVOUT,FMT=9020)ITER,FRAC,QEXT,QABS,QPHA 180 CONTINUE C*** first call to timeit C CALL TIMEIT('b') C*** C C Compute |Q0>=A|P0> C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXPI,CXQI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C C NOW COMPUTE ,,AND ALPHA0=/ QIQI=0. GIGI=0. DO 190 J=1,NAT3 GIGI=GIGI + CONJG(CXGI(J))*CXGI(J) QIQI=QIQI + CONJG(CXQI(J))*CXQI(J) 190 CONTINUE ALPHAI=GIGI/QIQI C C Now compute |X1>=|X0>+ALPHA0*|P0> DO 200 J=1,NAT3 CXXI(J)=CXXI(J) + ALPHAI*CXPI(J) 200 CONTINUE ITER=1 C C Save results for extrapolation machinery: QEXI1=QEXT QABI1=QABS QPHI1=QPHA FRACI1=FRAC C Now compute |AX1> and FRAC=fractional error FRAC=0. C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXXI,CXAXI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C DO 210 J=1,NAT3 CXALP0=CXAXI(J) - CXE(J) FRAC=FRAC + CXALP0*CONJG(CXALP0) 210 CONTINUE C FRAC=SQRT(FRAC/EMAG2) C C NOW COMPUTE CROSS SECTIONS CALL EVALQ(CXALPH,AK,NAT3,E02,CXE,CXXI,CABS,CEXT,CPHA,MXN3,1) QABS=CABS/PIA2 QEXT=CEXT/PIA2 QPHA=CPHA/PIA2 C*** second call to timeit: C CALL TIMEIT(' DDACCG.1') C*** WRITE(IDVOUT,FMT=9020)ITER,FRAC,QEXT,QABS,QPHA C C CHECK WHETHER WISH TO TERMINATE ITERATION AT THIS POINT IF(FRAC.GT.TOL)GOTO 220 C EXTRAPOLATE FRAC=FRAC/(FRACI1-FRAC) QEXT=QEXT+FRAC*(QEXT-QEXI1) QABS=QABS+FRAC*(QABS-QABI1) QPHA=QPHA+FRAC*(QPHA-QPHI1) QSCA=QEXT-QABS WRITE(IDVOUT,FMT=9030)QEXT,QABS,QPHA,QSCA GOTO 350 C C Begin iteration ITER=2,3,... to compute |X2>,|X3>,... C*********************************************************** 220 CONTINUE ITMAX=NAT3 DO 340 IT=2,ITMAX C*** first call to timeit: C CALL TIMEIT(' c') C*** ITER=IT C C Transfer -> GI1GI1=GIGI C C Compute vector |GI>=AC|E>-AC|AXI> C Compute GIGI= C CALL CPROD(AK,CMETHD,'C',CXADIA, & CXAXI,CXGI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C GIGI=0. DO 230 J=1,NAT3 CXGI(J)=CXACE(J)-CXGI(J) 230 CONTINUE DO 235 J=1,NAT3 GIGI=GIGI+CONJG(CXGI(J))*CXGI(J) 235 CONTINUE C C Compute BETAI-1=/ BETAI1=GIGI/GI1GI1 C C Compute |PI>=|GI>+BETAI-1|PI-1> DO 240 J=1,NAT3 CXPI(J)=CXGI(J) + BETAI1*CXPI(J) 240 CONTINUE C C Compute |QI>=A|PI> C Compute C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXPI,CXQI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C QIQI=0. DO 250 J=1,NAT3 QIQI=QIQI+CONJG(CXQI(J))*CXQI(J) 250 CONTINUE C C Compute ALPHAI ALPHAI=GIGI/QIQI C C Compute |XI+1>=|XI>+ALPHAI*|PI> DO 260 J=1,NAT3 CXXI(J)=CXXI(J) + ALPHAI*CXPI(J) 260 CONTINUE C Except every 10TH iteration, compute |AXI+1>=|AXI>+ALPHAI*|QI> IF(IT.EQ.10*(IT/10))GOTO 280 DO 270 J=1,NAT3 CXAXI(J)=CXAXI(J)+ALPHAI*CXQI(J) 270 CONTINUE GOTO 290 C C In order to avoid accumulation of roundoff error, every 10th C iteration compute |AXI+1> directly from |XI+1> 280 CONTINUE C CALL CPROD(AK,CMETHD,'N',CXADIA, & CXXI,CXAXI,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C 290 CONTINUE C C Update extrapolation machinery QEXI2=QEXI1 QEXI1=QEXT QABI2=QABI1 QABI1=QABS QPHI2=QPHI1 QPHI1=QPHA FRACI2=FRACI1 FRACI1=FRAC C Compute new fractional error: FRAC=0. DO 300 J=1,NAT3 CXZW(J,1)=CXAXI(J)-CXE(J) SCRRS1(J)=CXZW(J,1)*CONJG(CXZW(J,1)) 300 CONTINUE DO 310 J=1,NAT3 FRAC=FRAC+SCRRS1(J) 310 CONTINUE FRAC=SQRT(FRAC/EMAG2) C C Compute new Qabs and Qext C CALL EVALQ(CXALPH,AK,NAT3,E02,CXE,CXXI,CABS,CEXT,CPHA,MXN3,1) QABS=CABS/PIA2 QEXT=CEXT/PIA2 QPHA=CPHA/PIA2 C*** second call to timeit C CALL TIMEIT(' DDACCG.2') C*** WRITE(IDVOUT,FMT=9020)ITER,FRAC,QEXT,QABS,QPHA C C********** Extrapolation machinery ************ C Attempt to extrapolate to zero error, assuming Q values C to be linear functions of fractional error C Do not attempt extrapolation unless C fractional error < 100.*TOL C and fractional error is decreasing monotonically. C IF(FRAC.GT.1.E2*TOL)GOTO 340 IF(FRAC.GE.FRACI1)GOTO 340 IF(FRACI1.GE.FRACI2)GOTO 330 C If FRAC.LE.TOL, skip extrapolation test IF(FRAC.LE.TOL)GOTO 320 C Do not accept extrapolation unless linear function fits C three consecutive Qext values C Test success of extrapolation: try to estimate present result for Qest QEXEST=QEXI1+(FRAC-FRACI1)*(QEXI1-QEXI2)/(FRACI1-FRACI2) IF(ABS(QEXT-QEXEST).GT.TOL*QEXT)GOTO 340 C Proceed to extrapolate to final Qext 320 CONTINUE QEXT=QEXT-FRAC*(QEXT-QEXI2)/(FRAC-FRACI2) QABS=QABS-(QABS-QABI2)*FRAC/(FRAC-FRACI2) QPHA=QPHA-(QPHA-QPHI2)*FRAC/(FRAC-FRACI2) QSCA=QEXT-QABS WRITE(IDVOUT,FMT=9030)QEXT,QABS,QPHA,QSCA GOTO 350 C 330 CONTINUE IF(FRAC.GE.TOL)GOTO 340 QSCA=QEXT-QABS WRITE(IDVOUT,FMT=9050)QEXT,QABS,QPHA,QSCA GOTO 350 C 340 CONTINUE 350 CONTINUE RETURN 8990 FORMAT(' >DDACCG Start with |X0>=0') 8991 FORMAT(' >DDACCG Obtain |X0> from expansion in powers of alpha') 8992 FORMAT(' >DDACCG Obtain |X0> from file DDACCG.IN') 9000 FORMAT (' >CG ORDER =',I2,' NO FRAC.ERR.',5X,'QEXT=',1P,E10.3, & ' QABS=',E10.3,' QPHA=',E10.3) 9010 FORMAT (' >CG ORDER =',I2,' F.ERR=',1P,E10.3,' QEXT=',E10.3,' Q', & 'ABS=',E10.3,' QPHA=',E10.3) 9020 FORMAT (' >CG ITER=',I4,' F.ERR=',1P,E10.3,' QEXT=',E10.3,' QAB', & 'S=',E10.3,' QPHA=',E10.3) 9030 FORMAT (' >CG EXTRAP.:',18X,' QEXT=',1P,E10.3,' QABS=',E10.3,' Q', & 'PHA=',E10.3,/,' ',11X,'QSCA=QEXT-QABS=',E10.3) 9040 FORMAT (' >DDACCG extrap: QEXT=',1P,E10.3,' QABS=',E10.3,' QPHA=', & E10.3) 9050 FORMAT (' >CG NO EXTRAP.:',15X,' QEXT=',1P,E10.3,' QABS=',E10.3, & ' QPHA=',E10.3,/,' ',13X,'QSCA=QEXT-QABS=',E10.3) 9060 FORMAT (' >CG NO EXTRAP.: QEXT=',1P,E10.3,' QABS=',E10.3, & ' QPHA=',E10.3) END SUBROUTINE DIELEC(CDIEL,WAVE,IDVOUT,CFLEPS,CXEPS, & MXCOMP,MXWAVT,NCOMP,E1A,E2A,WVA) C Arguments: INTEGER MXCOMP,MXWAVT CHARACTER*6 CDIEL CHARACTER*40 CFLEPS(MXCOMP) COMPLEX CXEPS(MXCOMP) INTEGER IDVOUT,NCOMP REAL WAVE REAL E1A(MXWAVT),E2A(MXWAVT), & WVA(MXWAVT) C Local Variables: LOGICAL INIT CHARACTER*80 CDESCR COMPLEX CXI INTEGER I,J,JJ,NWAVT,IWV,IREN,IIMN,IEPS1,IEPS2,ICOL REAL DUM1,DUM2,E1,E2,TEMP REAL XIN(10) EXTERNAL ERRMSG SAVE CXI,IEPS1,INIT,NWAVT *********************************************************************** C Given: C CDIEL = 'TABLES' to read from input tables C 'H2OICE' for H2O ice C 'H2OLIQ' for H2O liquid C WAVE = wavelength (micron) C IDVOUT=output unit number C CFLEPS(1-NCOMP) = names of files containing dielectric data C MXCOMP,MXWAVT = dimensioning information C NCOMP = number of components C NWAV = number of wavelengths C E1A, E2A, WVA = scratch arrays C Returns: C CXEPS(1-NCOMP) = dielectric constant for materials C 1-NCOMP C NOTE: C If CDIEL='TABLES', then it is assumed that file(s) containing C the table(s) will have following format: C line1 = description (CHARACTER*80) C line2 = IWV, IREN, IIMN, IEPS1, IEPS2 C where IWV = column in which wavelengths are tabulated C IREN = col. in which Re(N) is tabulated (0 if not) C IIMN = col. in which Im(N) is tabulated (0 if not) C IEPS1= col. in which Re(EPS) is tabulated (0 if not) C IEPS2= col. in which Im(EPS) is tabulated (0 if not) C line3 = header line (will not be read) C line4... = data, according to columns specified on line2 C wavelengths must be monotonic, either increasing or decreasing C C B.T.Draine, Princeton Univ. Obs. C History: C 90.12.04 (BTD): Rewritten to allow H2OICE and H2OLIQ options. C 90.12.21 (BTD): Correct handling of 'TABLES' option when more than C one wavelength is considered. C 91.05.02 (BTD): Added IDVOUT to argument list for subr. INTERP C 91.09.12 (BTD): Moved declaration of MXCOMP and MWAVT ahead of C declaration of CFLEPS and CXEPS C 93.12.06 (BTD): Added IEPS1,INIT, and NWAVT to SAVE statement C (without this did not run properly on SGI Indigo) 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 CXI/(0.,1.)/,INIT/.TRUE./ IF(CDIEL.EQ.'TABLES')THEN IF(INIT.OR.NCOMP.GT.1)THEN DO 1000 J=1,NCOMP OPEN(UNIT=3,FILE=CFLEPS(J),STATUS='OLD') C Read header line: READ(3,9000)CDESCR WRITE(IDVOUT,9000)CDESCR C Read line specifying columns for wavelength,Re(n),Im(n),Re(eps),Im(eps): READ(3,*)IWV,IREN,IIMN,IEPS1,IEPS2 ICOL=IWV IF(IREN.GT.ICOL)ICOL=IREN IF(IIMN.GT.ICOL)ICOL=IIMN IF(IEPS1.GT.ICOL)ICOL=IEPS1 IF(IEPS2.GT.ICOL)ICOL=IEPS2 C Skip header line: READ(3,*) DO 0500 I=1,MXWAVT READ(3,*,END=0600)(XIN(JJ),JJ=1,ICOL) WVA(I)=XIN(IWV) IF(IEPS1.GT.0)THEN E1A(I)=XIN(IEPS1) E2A(I)=XIN(IEPS2) ELSE E1A(I)=XIN(IREN) E2A(I)=XIN(IIMN) ENDIF NWAVT=I 0500 CONTINUE C Check whether there is unread data remaining in file READ(3,*,END=0600)(XIN(JJ),JJ=1,ICOL) C If this point is reached, apparently unread data remains, so C issue warning: CALL ERRMSG('WARNING','DIELEC', & 'parameter MXWAVT not large enough to read full dielec file') 0600 CLOSE(3) C C*** Have completed reading in table for this composition. C Now interpolate CALL INTERP(WAVE,E1,E2,WVA,E1A,E2A, & IDVOUT,MXWAVT,NWAVT) IF(IEPS1.GT.0)THEN CXEPS(J)=E1+CXI*E2 ELSE CXEPS(J)=(E1**2-E2**2)+2.*CXI*E1*E2 ENDIF 1000 CONTINUE INIT=.FALSE. ELSE C*** Perform this only if NCOMP=1 and previously initialized: CALL INTERP(WAVE,E1,E2,WVA,E1A,E2A,IDVOUT,MXWAVT,NWAVT) IF(IEPS1.GT.0)THEN CXEPS(1)=E1+CXI*E2 ELSE CXEPS(1)=(E1**2-E2**2)+2.*CXI*E1*E2 ENDIF ENDIF ELSEIF(CDIEL.EQ.'H2OICE')THEN WRITE(CDESCR,FMT='(A)')' H2O ice at T=270K' WRITE(IDVOUT,9000)CDESCR IF(NCOMP.GT.1)CALL ERRMSG('FATAL','DIELE', & 'H2OICE with NCOMP>1??') C For the moment set T=250K TEMP=250. CALL REFICE(0,WAVE,TEMP,E1,E2,DUM1,DUM2) CXEPS(1)=E1**2-E2**2+2.*CXI*E1*E2 RETURN ELSEIF(CDIEL.EQ.'H2OLIQ')THEN WRITE(CDESCR,FMT='(A)')' H2O liquid at T=280K' WRITE(IDVOUT,9000)CDESCR IF(NCOMP.GT.1)CALL ERRMSG('FATAL','DIELE', & 'H2OLIQ with NCOMP>1??') C For the moment set T=280K TEMP=280. CALL REFWAT(0,WAVE,TEMP,E1,E2,DUM1,DUM2) CXEPS(1)=E1**2-E2**2+2.*CXI*E1*E2 RETURN ELSE CALL ERRMSG('FATAL','DIELE','unrecognized value of CDIEL') ENDIF RETURN 9000 FORMAT(A80) END SUBROUTINE DIVIDE(CDIVID,X1,X2,NX,MXARY,ARY) C Arguments: REAL X1,X2 INTEGER MXARY,NX CHARACTER CDIVID*3 REAL ARY(MXARY) C C Local variables: INTEGER I REAL DELTA C C External Subroutines: EXTERNAL ERRMSG,WRIMSG C C Intrinsic Functions: INTRINSIC REAL C*********************************************************************** C Given: C CDIVID='LIN','INV', or 'LOG' C X1 = lower limit to interval C X2 = upper limit to interval C NX = number of elements C MXARY=dimensioning information for array ARY C C Returns: C ARY(1-NX)=vector of points with C ARY(1)=X1 C ARY(NX)=X2 C ARY(J) values spaced either C linearly (if CDIVID.EQ.'LIN') C linearly in 1/X (if CDIVID.EQ.'INV') C logarithmically (if CDIVID.EQ.'LOG') 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(NX.GT.MXARY)THEN CALL ERRMSG('FATAL','DIVIDE',' NX .GT. MXARY ') ENDIF C IF(X1.GT.X2)THEN C CALL ERRMSG('FATAL','DIVIDE',' X1 .GT. X2 ') C ENDIF IF(NX.LT.1)THEN CALL ERRMSG('FATAL','DIVIDE',' NX .LT. 1 ') ENDIF IF(NX.EQ.1)THEN CALL WRIMSG('DIVIDE',' Only one element initialized ') ARY(1)=X1 ELSE IF(CDIVID.EQ.'LIN')THEN DELTA=(X2-X1)/REAL(NX-1) DO 10 I=1,NX ARY(I)=X1+DELTA*REAL(I-1) 10 CONTINUE ELSEIF(CDIVID.EQ.'INV')THEN DELTA=(1./X2-1./X1)/REAL(NX-1) DO 20 I=1,NX ARY(I)=1./(1./X1+DELTA*REAL(I-1)) 20 CONTINUE ELSEIF(CDIVID.EQ.'LOG')THEN DELTA=ALOG(X2/X1)/REAL(NX-1) DO 30 I=1,NX ARY(I)=EXP(ALOG(X1)+DELTA*REAL(I-1)) 30 CONTINUE ENDIF ENDIF RETURN END SUBROUTINE ERRMSG(CSTATS,CSUBRT,CMSGNM) C Arguments: CHARACTER CMSGNM*(*),CSTATS*(*),CSUBRT*(*) C C Local variables: INTEGER IOERR C C External Subroutines: EXTERNAL GETSET C*********************************************************************** C Given: C CSTATS = 'WARNING' or 'FATAL' C CSUBRT = name of subroutine C CMESGN = message C Prints a warning message in a standardized way, C and STOPs if CSTATS='FATAL' C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** CALL GETSET('GET','IOERR',IOERR) IF(CSTATS.EQ.'FATAL')THEN WRITE (IOERR,FMT=9000) CSUBRT 9000 FORMAT(/,' >>>>> FATAL ERROR IN PROCEDURE: ',A) WRITE(IOERR,FMT=9010)CMSGNM 9010 FORMAT(1X,A) WRITE(IOERR,FMT=9020) 9020 FORMAT(' >>>>> EXECUTION ABORTED ') STOP ELSEIF(CSTATS.EQ.'WARNING')THEN WRITE(IOERR,FMT=9030)CSUBRT 9030 FORMAT(/,' >>>>> WARNING IN PROCEDURE: ',A) WRITE(IOERR,FMT=9010)CMSGNM 9040 FORMAT(1X,A) ENDIF RETURN END