SUBROUTINE ALPHA(CALPHA,CXALPH,CXE0R,CXEPS,ICOMP,AKR,NAT0,NAT3, & MXCOMP,MXN3) C*** Arguments: CHARACTER*6 CALPHA INTEGER MXN3,MXCOMP,NAT0,NAT3 COMPLEX CXALPH(MXN3),CXE0R(3),CXEPS(MXCOMP) INTEGER*2 ICOMP(MXN3) REAL AKR(3) C C*** Local Variables: CHARACTER CMSGNM*70 INTEGER IC,L COMPLEX CXI,CXRR REAL AK1,AK2,B1,B2,EMKD,PI,SUM DATA CXI/(0.,1.)/ C*********************************************************************** C Given: C CALPHA = 'LATTDR' or 'DRAI88' or 'GOBR88' or 'LDRISO' C to specify method ('LATTDR' is recommended) C CXE0R(1-3) = polarization vector in lattice coordinates C (assumed to be normalized) C CXEPS(1-NCOMP)=distinct values of dielectric constant C ICOMP(1-NAT3)=composition identifier for each lattice site and C direction (ICOMP=0 if lattice site is unoccupied) C MXCOMP = dimensioning information C MXN3 = dimensioning information C NAT0 = number of sites in original target (i.e., nonvacuum sites) C NAT3 = 3*number of sites in extended target (including vacuum C sites) C Returns: C CXALPHA(1-NAT3)=(polarizability/d^3) for each of 3 directions at C each lattice site C*** C If CALPHA = LATTDR: C Compute dipole polarizability using "Lattice Dispersion Relation" C of Draine & Goodman (1993,ApJ,March 10). It is required that C polarizability be such that an infinite lattice of such dipoles C reproduce the continuum dispersion relation for radiation C propagating with direction and polarization of radiation incident C on the DDA target. This is the recommended option. C C If CALPHA = DRAI88: C Compute dipole polarizability using C (1) eq. (2.03) of Draine (1988) (Clausius-Mossotti plus C radiative reaction correction) C C If CALPHA = GOBR88: C Compute dipole polarizability using C (1) radiative reaction correction from Draine (1988) C (2) O[(kd)^2] correction term used by C Goedecke & OBrien (1988) and Hage & Greenberg (1990) C C If CALPHA = LDRISO: C Compute dipole polarizability using "Isotropized Lattice Dispersion C Relation" of Draine & Goodman (1993,ApJ,March 10). This replaces C direction- and polarization-dependent term in Lattice Dispersion C Relation with the average value for random direction and polarization. C C*** C Note: CXALPH = polarizability/d^3 C In the event that ICOMP=0, then we set CXALPH=1. C*********************************************************************** C B.T.Draine, Princeton Univ. Obs. C History: C 90. 9.13 (BTD): Corrected error in DO loop limit. C 90.11. 1 (BTD): Special treatment for "vacuum" sites in order C to allow FFT treatment. C 90.11. 2 (BTD): Set CXALPH=(1.,0.) at vacuum sites. C 91. 4.30 (BTD): Modified to include both O[(kd)^2] correction term C from Goedecke & Obrien (1988) in addition to C radiative reaction correction. C 91. 5. 7 (BTD): Modified to allow easy choice among three methods C for computing polarizability. C 91. 5. 7 (BTD): Added call to WRIMSG to record which method in use C 91.05.08 (BTD): Added CALPHA to argument list. C 91.05.09 (BTD): Added printing of kd and magnitude of correction C terms. C 91.05.13 (BTD): Experiment with use of approximate directional C average for coefficient of CXEPS in Lattice C Disperision Relation theory C 91.05.13 (BTD): Corrected error in radiative reaction correction C (error presumably introduced during last week) C 91.05.14 (BTD): Added separate options LDRXYZ and LDRAVG C 91.05.15 (BTD): Added option LDR000 to omit CXEPS-dependent C correction term C 91.05.30 (BTD): Changed LDRXYZ -> LDR100 C Added option LDR111 C 91.09.12 (BTD): Corrected error in numerical coefficient for C LDRAVG case C 91.09.17 (BTD): Eliminate options LDR000,LDRXYZ,LDRAVG C Introduce option LATTDR (LATTice Dispersion Relation) C which includes explicit dependence C on direction and polarization state C remove AK1 from argument list C add AKR, CXE0R to argument list C 92.05.14 (BTD): Introduce option LDRISO to experiment with C using average directional correction term C rather than using correction for incident direction C and polarization C 92.05.16 (BTD): Correct error in LDRISO option: average SUM should be C 3/15 rather than 4/15 C C Copyright (C) 1993,1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** PI=4.*ATAN(1.) AK2=AKR(1)*AKR(1)+AKR(2)*AKR(2)+AKR(3)*AKR(3) AK1=SQRT(AK2) CXRR=-(AK1*AK2/1.5)*CXI C*** EMKD = |m|*k_0*d , where |m|=refractive index C k_0 = wave vector in vacuo C d = lattice spacing EMKD=SQRT(CABS(CXEPS(1)))*AK1 C C*********************************************************************** C*** Draine 1988: C IF(CALPHA.EQ.'DRAI88')THEN WRITE(CMSGNM,FMT='(A,F6.4)') & 'Clausius-Mossotti plus radiative reaction for |m|kd=', & EMKD CALL WRIMSG('ALPHA ',CMSGNM) DO 1000 L=1,NAT3 IC=ICOMP(L) IF(IC.GT.0)THEN C*** First compute Clausius-Mossotti polarizability: CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.) C*** Now correct this for C (1) radiative-reaction CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR) ELSEIF(IC.EQ.0)THEN C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites. CXALPH(L)=1. ENDIF 1000 CONTINUE RETURN C C*********************************************************************** C*** Goedecke & OBrien 1988: C ELSEIF(CALPHA.EQ.'GOBR88')THEN B1=-(4.*PI/3.)**(1./3.)*AK2 WRITE(CMSGNM,FMT='(A,F6.4)') & 'Goedecke & OBrien (1988) correction for |m|kd=',EMKD CALL WRIMSG('ALPHA ',CMSGNM) DO 2000 L=1,NAT3 IC=ICOMP(L) IF(IC.GT.0)THEN C*** First compute Clausius-Mossotti polarizability: CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.) C*** Now correct this for C (1) O[(kd)**2] correction from Goedecke and Obrien (1988) CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*B1) C (2) radiative-reaction CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR) ELSEIF(IC.EQ.0)THEN C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites. CXALPH(L)=1. ENDIF 2000 CONTINUE RETURN C C*********************************************************************** C*** Lattice dispersion relation C ELSEIF(CALPHA.EQ.'LATTDR')THEN WRITE(CMSGNM,FMT='(A,F6.4)') & 'Lattice dispersion relation for |m|k_0d=',EMKD CALL WRIMSG('ALPHA ',CMSGNM) C*** Compute sum (a_j*e_j)^2 , where a_j=unit propagation vector C e_j=unit polarization vector SUM=0. DO 2700 L=1,3 SUM=SUM+(AKR(L)*CABS(CXE0R(L)))**2 2700 CONTINUE SUM=SUM/AK2 B1=-1.8915316*AK2 B2=(0.1648469-1.7700004*SUM)*AK2 DO 3000 L=1,NAT3 IC=ICOMP(L) IF(IC.GT.0)THEN C*** First compute Clausius-Mossotti polarizability: CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.) C*** Determine polarizability by requiring that infinite lattice of C dipoles have dipersion relation of continuum. CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*(B1+CXEPS(IC)*B2)) C*** Radiative-reaction correction: CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR) ELSEIF(IC.EQ.0)THEN C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites. CXALPH(L)=1. ENDIF 3000 CONTINUE RETURN ELSEIF(CALPHA.EQ.'LDRISO')THEN WRITE(CMSGNM,FMT='(A,F6.4)') & 'Isotropic Lattice disp. rel. for |m|k_0d=',EMKD CALL WRIMSG('ALPHA ',CMSGNM) C*** Compute sum (a_j*e_j)^2 , where a_j=unit propagation vector C e_j=unit polarization vector C Here experiment by taking average value of sum (a_j*e_j)^2 = 1/5 SUM=0.2D0 B1=-1.8915316*AK2 B2=(0.1648469-1.7700004*SUM)*AK2 DO 4000 L=1,NAT3 IC=ICOMP(L) IF(IC.GT.0)THEN C*** First compute Clausius-Mossotti polarizability: CXALPH(L)=(.75/PI)*(CXEPS(IC)-1.)/(CXEPS(IC)+2.) C*** Determine polarizability by requiring that infinite lattice of C dipoles have dipersion relation of continuum. CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*(B1+CXEPS(IC)*B2)) C*** Radiative-reaction correction: CXALPH(L)=CXALPH(L)/(1.+CXALPH(L)*CXRR) ELSEIF(IC.EQ.0)THEN C To avoid divisions by zero, etc., set CXALPH=1 for vacuum sites. CXALPH(L)=1. ENDIF 4000 CONTINUE RETURN ENDIF END SUBROUTINE PETR(X,B,WRK,LDA,IPAR,SPAR,MATVEC,CMATVEC) C solves Ax=b C C Input: C B(LOCLEN) - RHS (doesn't change on output) C WRK(LDA, LOCLEN) - scratch array C IPAR - defines integer paramteres (see documentation of PIM) C SPAR - define real parameters C LDA - leading dimension of arrays C MATVEC -- name of procedure for computing y=Ax C CMATVEC -- name of procedure for computing y= conj(A') x C C Output: C X(LOCLEN) C C Reference: C Petravic,M and Kuo-Petravic,G. 1979, JCP, 32,263 C History: C Late 80' (Coded by Draine) C Changes to include "cprod" (PJF+BTD) C 95.06.01 (PJF) re-written to conform to PIM C 95.08.11 (BTD) minor editing (comments, etc.) C C Parameters: INTEGER IACE,IAXI,IGI,IPI,IQI,IR PARAMETER (IACE=1,IGI=2,IPI=3,IQI=4,IAXI=5,IR=6) C Arguments: INTEGER LDA INTEGER IPAR(*) REAL SPAR(*) COMPLEX B(*), WRK(LDA,*), X(*) EXTERNAL CMATVEC,MATVEC C C Local variables COMPLEX CDUMMY INTEGER IDUMMY,ITNO,L,LOCLEN,MAXIT,STATUS REAL & ALPHAI,BETAI1,BNORM,DUMMY,EPSILON,EXITNORM, & GIGI,GI1GI1,QIQI,RHSSTOP REAL SCSETRHSSTOP EXTERNAL PSCNRM2 C C*********************************************************************** LOCLEN = IPAR(4) MAXIT = IPAR(10) EPSILON=SPAR(1) C C Compute RHSSTOP = EPSILON*|B| C RHSSTOP = SCSETRHSSTOP(B,CDUMMY,EPSILON,IPAR,CDUMMY,PSCNRM2) C C Compute conjg(A')*B C CALL CMATVEC(B,WRK(1,IACE),IDUMMY) BNORM = 0. DO 170 L = 1,LOCLEN BNORM = BNORM + B(L)*CONJG(B(L)) WRK(L,IGI) = WRK(L,IACE) WRK(L,IPI) = WRK(L,IGI) 170 CONTINUE C C Compute |QI>=A|PI> C CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY) QIQI = 0. GIGI = 0. DO 190 L = 1,LOCLEN GIGI = GIGI + CONJG(WRK(L,IGI))*WRK(l,IGI) QIQI = QIQI + CONJG(WRK(L,IQI))*WRK(L,IQI) 190 CONTINUE ALPHAI = GIGI/QIQI DO 200 L = 1,LOCLEN X(L) = X(L) + ALPHAI*WRK(L,IPI) 200 CONTINUE C C compute |AX1>: C CALL MATVEC(X,WRK(1,IAXI),IDUMMY) C DO 340 ITNO = 2,MAXIT IPAR(11) = ITNO C C Transfer -> : C GI1GI1 = GIGI C C compute |GI>=AC|E>-AC|AXI>: C CALL CMATVEC(WRK(1,IAXI), WRK(1,IGI), IDUMMY) C C Compute GIGI=: GIGI = 0. DO 230 L = 1,LOCLEN WRK(L,IGI) = WRK(L,IACE) - WRK(L,IGI) 230 CONTINUE DO 235 L = 1,LOCLEN GIGI = GIGI + CONJG(WRK(L,IGI))*WRK(L,IGI) 235 CONTINUE C C Compute BETAI-1=/: C BETAI1 = GIGI/GI1GI1 C C Compute |PI>=|GI>+BETAI-1|PI-1>: C DO 240 L = 1,LOCLEN WRK(L,IPI) = WRK(L,IGI) + BETAI1*WRK(L,IPI) 240 CONTINUE C C Compute |QI>=A|PI>: C CALL MATVEC(WRK(1,IPI), WRK(1,IQI), IDUMMY) C C Compute : C QIQI = 0. DO 250 L = 1,LOCLEN QIQI = QIQI + CONJG(WRK(L,IQI))*WRK(L,IQI) 250 CONTINUE C C Compute ALPHAI = /: C ALPHAI = GIGI/QIQI C C Compute |XI+1>=|XI>+ALPHAI*|PI>: C DO 260 L = 1,LOCLEN X(L) = X(L) + ALPHAI*WRK(L,IPI) 260 CONTINUE C C Except every 10TH iteration, compute |AXI+1>=|AXI>+ALPHAI*|QI> C Draine: warning this part perhaps not checked C IF (ITNO.NE.10* (ITNO/10)) THEN DO 270 L = 1,LOCLEN WRK(L,IAXI) = WRK(L,IAXI) + ALPHAI*WRK(L,IQI) 270 CONTINUE ELSE CALL MATVEC(X,WRK(1,IAXI),IDUMMY) ENDIF C C Compute residual vector |RI>=A|XI>-|B> C DO 300 L=1, LOCLEN WRK(L,IR) = WRK(L,IAXI)-B(L) 300 CONTINUE C C Call STOPCRIT to check whether to terminate iteration. C Criterion: terminate if < TOL**2 * C Return with STATUS=0 if this condition is satisfied. C CALL STOPCRIT(B,WRK(1,IR),CDUMMY,CDUMMY,CDUMMY,WRK, + RHSSTOP,IDUMMY,EXITNORM, + STATUS,IPAR,DUMMY,DUMMY,DUMMY,DUMMY,PSCNRM2) C C Call PROGRESS to report "error" = sqrt()/sqrt() C CALL PROGRESS(LOCLEN,ITNO,EXITNORM,CDUMMY,CDUMMY,CDUMMY) C IF(STATUS.EQ.0)GOTO 500 340 CONTINUE 500 CONTINUE IPAR(12) = 0 RETURN END SUBROUTINE COPYIT(CXA,CXB,NPTS) C*** Arguments: INTEGER NPTS COMPLEX CXA(*), CXB(*) C C*** Local variables: INTEGER IPTS C*********************************************************************** C Purpose: To copy a complex vector from one memory location to another. C*********************************************************************** DO 10 IPTS=1,NPTS CXB(IPTS)=CXA(IPTS) 10 CONTINUE RETURN END SUBROUTINE DIAGL(U,V,IPAR) C cg package C Left preconditioning subroutine - division by diagonal C Arguments: COMPLEX U(*), V(*) INTEGER IPAR(*) C C Common variables: CHARACTER*6 CMETHD INTEGER & IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ COMPLEX CXADIA(1) C----------------------------------------------------------------------- COMMON /M2/ CXADIA COMMON /M6/ MXNATF, MXNXF, MXNYF, MXNZF, NAT, NAT3, NAT0, $ NX, NY, NZ, MXMEMF, MXN3F, IDVOUT, CMETHD C----------------------------------------------------------------------- CALL CDIV(U,V,CXADIA,NAT3) RETURN END SUBROUTINE CDIV(U,V,CXA,N) C Arguments: INTEGER N COMPLEX U(*), V(*), CXA(*) C C Local variables: INTEGER I C C cg package C needed by left preconditioning subroutine DO 10 I=1,N V(I)=U(I)/CXA(I) 10 CONTINUE RETURN END SUBROUTINE MATVEC(CXX, CXY, IPAR) C Arguments: COMPLEX CXX(*), CXY(*), IPAR(*) C Common: INTEGER*2 iocc(1) INTEGER & IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ CHARACTER*6 CMETHD REAL AKR(3) COMPLEX CXADIA(1),CXZC(1), CXZW(1) C----------------------------------------------------------------------- COMMON /M1/ AKR, /M2/ CXADIA, /M3/ CXZC, /M4/CXZW, /M5/ IOCC COMMON /M6/ MXNATF, MXNXF, MXNYF, MXNZF, NAT, NAT3, NAT0, & NX, NY, NZ, MXMEMF, MXN3F, IDVOUT, CMETHD C----------------------------------------------------------------------- C C*********************************************************************** C matvec --> cxy = A cxx 'N' C tmatvec --> cxy = A' cxx 'T' C cmatvec --> cxy = conjg(A') cxx 'C' C C Notice that DDSCAT has cases 'N' and 'C' implemented in CPROD C Notice, however, that in DDSCAT 'T'='N'. C Some CG methods use cases N, C (like Petravic). C Some CG methods use cases N, T (like PIM) C Some CG use all three products N, C, T (parts of CCGPAK) C C Information about matrix A is transfered via commons /m1/-/m6/ with C the main program. Arrays cxadia, cxzc, cxzwm, iocc are properly C dimensioned inside "cprod". You may get warning during the compilation C about "misaligned common". Ignore these. Also notice mxnat-f, C mxn3-f,mxnx-f, mxny-f, mxnz-f. These are defined by "parameter C statement" in main code and this is the reason why they have "f" suffix. C C ipar - information array, not used, defined for compatibility C*********************************************************************** CALL CPROD(AKR,CMETHD,'N',CXADIA, & CXX,CXY,CXZC,CXZW, & IDVOUT,IOCC,MXN3f,MXNATf,MXNXf,MXNYf,MXNZf, & NAT,NAT0,NAT3,NX,NY,NZ) RETURN END SUBROUTINE CMATVEC(CXX, CXY, IPAR) C Arguments: COMPLEX CXX(*), CXY(*), IPAR(*) C C Common: INTEGER*2 iocc(1) INTEGER & IDVOUT,MXMEMF,MXN3F,MXNATF,MXNXF,MXNYF,MXNZF, & NAT,NAT0,NAT3,NX,NY,NZ CHARACTER*6 CMETHD REAL akr(3) COMPLEX CXADIA(1),CXZC(1), CXZW(1) C----------------------------------------------------------------------- COMMON /M1/ AKR, /M2/ CXADIA, /M3/ CXZC, /M4/CXZW, /M5/ IOCC COMMON /M6/ MXNATF, MXNXF, MXNYF, MXNZF, NAT, NAT3, NAT0, & NX, NY, NZ, MXMEMF, MXN3F, IDVOUT, CMETHD C----------------------------------------------------------------------- C Local variables: C*********************************************************************** C matvec --> cxy = A cxx C tmatvec --> cxy = A' cxx C cmatvec --> cxy = conjg(A') cxx C*********************************************************************** CALL CPROD(AKR,CMETHD,'C',CXADIA, & CXX,CXY,CXZC,CXZW, & IDVOUT,IOCC,MXN3f,MXNATf,MXNXf,MXNYf,MXNZf, & NAT,NAT0,NAT3,NX,NY,NZ) RETURN END SUBROUTINE CPROD(AK,CMETHD,CWHAT,CXADIA, & CXX,CXY,CXZC,CXZW, & IDVOUT,IOCC,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ) C C If CWHAT = 'N', this returns CXY(j)=A(j,k)*CXX(k) (sum over k) C = 'C', this returns CXY(j)=conjg(A(k,j))*CXX(k) (sum over k) C C It is assumed that matrix A(j,k) is symmetric C Diagonal elements of A(j,k) are in vector ADIAG(j) C Off-diagonal elements of A(j,k) are produced internally by C subroutine ESELF, which uses 3d-FFT to accomplish convolution C of CXX with A-matrix. C C Original version of CPROD created by P.J.Flatau, Colorado State Univ. C Subsequently modified by B.T.Draine, Princeton Univ. Obs. C History: C 90/11/ 4 (BTD): Major modification to incorporate FFT method (FFT C computations are handled by ESELF). C 90/11/ 8 (BTD): Remove DO...ENDDO structures; renumber DO...CONTINUEs C 90/11/21 (BTD): Remove "CCGMVM" option: only FFT stuff retained. C 90/12/03 (BTD): No longer necessary to change ordering of CXX and CXY C 90/12/05 (BTD): Changed call to ESELF to eliminate MXNX,MXNY,MXNZ C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C*** Arguments: INTEGER IDVOUT,MXN3,MXNAT,MXNX,MXNY,MXNZ, & NAT,NAT0,NAT3,NX,NY,NZ CHARACTER CMETHD*6,CWHAT*1 INTEGER*2 IOCC(MXNAT) REAL AK(3) COMPLEX CXADIA(MXN3), & CXX(MXN3),CXY(MXN3),CXZC(NX+1,NY+1,NZ+1,6), & CXZW(MXNX*MXNY*MXNZ,*) C C*** Local variables: INTEGER J REAL AKD C C Intrinsic functions: INTRINSIC CONJG C C*********************************************************************** AKD=SQRT(AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3)) C We assume that input vector CXX has zeroes for elements corresponding C to vacuum sites. IF(CWHAT.EQ.'N')THEN C********************************************************************** C CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,CXZC,CXZW,CXY) C C************************************************************************* C ESELF omitted contribution from diagonal terms: add these. C Also remember that ESELF computes -A_jk*x_k DO 160 J=1,NAT3 CXY(J)=CXADIA(J)*CXX(J)-CXY(J) 160 CONTINUE C C*** C ELSEIF(CWHAT.EQ.'C')THEN C Need to temporarily replace CXX by conjg(CXX) DO 170 J=1,NAT3 CXX(J)=CONJG(CXX(J)) 170 CONTINUE C********************************************************************** C CALL ESELF(CMETHD,CXX,NX,NY,NZ,AKD,CXZC,CXZW,CXY) C C************************************************************************ C C ESELF omitted contribution from diagonal terms: add these. C Also remember that ESELF computed -A_jk*conjg(x_k), while C we seek to compute conjg(A_kj)* x_k. Since A_kj=A_jk, we have C conjg(A_kj)*x_k=conjg(A_jk*conjg(x_k)) DO 230 J=1,NAT3 CXY(J)=CXADIA(J)*CXX(J)-CXY(J) CXY(J)=CONJG(CXY(J)) 230 CONTINUE C Restore CXX: DO 240 J=1,NAT3 CXX(J)=CONJG(CXX(J)) 240 CONTINUE ELSE WRITE(IDVOUT,FMT=*)' Error in CPROD: CWHAT= ',CWHAT ENDIF C C*** Now must "clean" product vector CXY: C (zero out elements corresponding to vacuum sites) C IF(NAT0.LT.NAT)CALL NULLER(CXY,IOCC,MXNAT,MXN3,NAT) RETURN END SUBROUTINE CFFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) C C PURPOSE PERFORMS MULTIPLE FAST FOURIER TRANSFORMS. THIS PACKAGE C WILL PERFORM A NUMBER OF SIMULTANEOUS COMPLEX PERIODIC C FOURIER TRANSFORMS OR CORRESPONDING INVERSE TRANSFORMS. C THAT IS, GIVEN A SET OF COMPLEX GRIDPOINT VECTORS, THE C PACKAGE RETURNS A SET OF COMPLEX FOURIER C COEFFICIENT VECTORS, OR VICE VERSA. THE LENGTH OF THE C TRANSFORMS MUST BE A NUMBER GREATER THAN 1 THAT HAS C NO PRIME FACTORS OTHER THAN 2, 3, AND 5. C C THE PACKAGE CFFT99 CONTAINS SEVERAL USER-LEVEL ROUTINES: C C SUBROUTINE CFTFAX C AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE C BEFORE A SEQUENCE OF CALLS TO CFFT99 C (PROVIDED THAT N IS NOT CHANGED). C C SUBROUTINE CFFT99 C THE ACTUAL TRANSFORM ROUTINE ROUTINE, CABABLE OF C PERFORMING BOTH THE TRANSFORM AND ITS INVERSE. C HOWEVER, AS THE TRANSFORMS ARE NOT NORMALIZED, C THE APPLICATION OF A TRANSFORM FOLLOWED BY ITS C INVERSE WILL YIELD THE ORIGINAL VALUES MULTIPLIED C BY N. C C C ACCESS *FORTRAN,P=XLIB,SN=CFFT99 C C C USAGE LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 0, C Q .GE. 0, AND R .GE. 0. THEN A TYPICAL SEQUENCE OF C CALLS TO TRANSFORM A GIVEN SET OF COMPLEX VECTORS OF C LENGTH N TO A SET OF (UNSCALED) COMPLEX FOURIER C COEFFICIENT VECTORS OF LENGTH N IS C C DIMENSION IFAX(13),TRIGS(2*N) C COMPLEX A(...), WORK(...) C C CALL CFTFAX (N, IFAX, TRIGS) C CALL CFFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) C C THE OUTPUT VECTORS OVERWRITE THE INPUT VECTORS, AND C THESE ARE STORED IN A. WITH APPROPRIATE CHOICES FOR C THE OTHER ARGUMENTS, THESE VECTORS MAY BE CONSIDERED C EITHER THE ROWS OR THE COLUMNS OF THE ARRAY A. C SEE THE INDIVIDUAL WRITE-UPS FOR CFTFAX AND C CFFT99 BELOW, FOR A DETAILED DESCRIPTION OF THE C ARGUMENTS. C C HISTORY THE PACKAGE WAS WRITTEN BY 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----------------------------------------------------------------------- C C SUBROUTINE CFTFAX (N,IFAX,TRIGS) C C PURPOSE A SET-UP ROUTINE FOR CFFT99. IT NEED ONLY BE C CALLED ONCE BEFORE A SEQUENCE OF CALLS TO CFFT99, C PROVIDED THAT N IS NOT CHANGED. C C ARGUMENT IFAX(13),TRIGS(2*N) C DIMENSIONS C C ARGUMENTS C C ON INPUT N C AN EVEN NUMBER GREATER THAN 1 THAT HAS NO PRIME FACTOR C GREATER THAN 5. N IS THE LENGTH OF THE TRANSFORMS (SEE C THE DOCUMENTATION FOR CFFT99 FOR THE DEFINITION OF C THE TRANSFORMS). C C IFAX C AN INTEGER ARRAY. THE NUMBER OF ELEMENTS ACTUALLY USED C WILL DEPEND ON THE FACTORIZATION OF N. DIMENSIONING C IFAX FOR 13 SUFFICES FOR ALL N LESS THAN 1 MILLION. C C TRIGS C A REAL ARRAY OF DIMENSION 2*N C C ON OUTPUT IFAX C CONTAINS THE FACTORIZATION OF N. IFAX(1) IS THE C NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED C IN IFAX(2),IFAX(3),... IF N HAS ANY PRIME FACTORS C GREATER THAN 5, IFAX(1) IS SET TO -99. C C TRIGS C AN ARRAY OF TRIGONOMETRIC FUNCTION VALUES SUBSEQUENTLY C USED BY THE CFT ROUTINES. C C----------------------------------------------------------------------- C C SUBROUTINE CFFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) C C PURPOSE PERFORM A NUMBER OF SIMULTANEOUS (UNNORMALIZED) COMPLEX C PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE C TRANSFORMS. GIVEN A SET OF COMPLEX GRIDPOINT C VECTORS, THE PACKAGE RETURNS A SET OF C COMPLEX FOURIER COEFFICIENT VECTORS, OR VICE C VERSA. THE LENGTH OF THE TRANSFORMS MUST BE A C NUMBER HAVING NO PRIME FACTORS OTHER THAN C 2, 3, AND 5. THIS ROUTINE IS C OPTIMIZED FOR USE ON THE CRAY-1. C C ARGUMENT COMPLEX A(N*INC+(LOT-1)*JUMP), WORK(N*LOT) C DIMENSIONS REAL TRIGS(2*N), INTEGER IFAX(13) C C ARGUMENTS C C ON INPUT A C A COMPLEX ARRAY OF LENGTH N*INC+(LOT-1)*JUMP CONTAINING C THE INPUT GRIDPOINT OR COEFFICIENT VECTORS. THIS ARRAY C OVERWRITTEN BY THE RESULTS. C C N.B. ALTHOUGH THE ARRAY A IS USUALLY CONSIDERED TO BE OF C TYPE COMPLEX IN THE CALLING PROGRAM, IT IS TREATED AS C REAL WITHIN THE TRANSFORM PACKAGE. THIS REQUIRES THAT C SUCH TYPE CONFLICTS ARE PERMITTED IN THE USER"S C ENVIRONMENT, AND THAT THE STORAGE OF COMPLEX NUMBERS C MATCHES THE ASSUMPTIONS OF THIS ROUTINE. THIS ROUTINE C ASSUMES THAT THE REAL AND IMAGINARY PORTIONS OF A C COMPLEX NUMBER OCCUPY ADJACENT ELEMENTS OF MEMORY. IF C THESE CONDITIONS ARE NOT MET, THE USER MUST TREAT THE C ARRAY A AS REAL (AND OF TWICE THE ABOVE LENGTH), AND C WRITE THE CALLING PROGRAM TO TREAT THE REAL AND C IMAGINARY PORTIONS EXPLICITLY. C C WORK C A COMPLEX WORK ARRAY OF LENGTH N*LOT OR A REAL ARRAY C OF LENGTH 2*N*LOT. SEE N.B. ABOVE. C C TRIGS C AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST. C C IFAX C AN ARRAY SET UP BY CFTFAX, WHICH MUST BE CALLED FIRST. C C C N.B. IN THE FOLLOWING ARGUMENTS, INCREMENTS ARE MEASURED C IN WORD PAIRS, BECAUSE EACH COMPLEX ELEMENT IS ASSUMED C TO OCCUPY AN ADJACENT PAIR OF WORDS IN MEMORY. C C INC C THE INCREMENT (IN WORD PAIRS) BETWEEN SUCCESSIVE ELEMENT C OF EACH (COMPLEX) GRIDPOINT OR COEFFICIENT VECTOR C (E.G. INC=1 FOR CONSECUTIVELY STORED DATA). C C JUMP C THE INCREMENT (IN WORD PAIRS) BETWEEN THE FIRST ELEMENTS C OF SUCCESSIVE DATA OR COEFFICIENT VECTORS. ON THE CRAY- C TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8 C (TO AVOID MEMORY BANK CONFLICTS). FOR CLARIFICATION OF C INC AND JUMP, SEE THE EXAMPLES BELOW. C C N C THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF C TRANSFORMS, BELOW). C C LOT C THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY. C C ISIGN C = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER C COEFFICIENTS. C = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO C GRIDPOINT VALUES. C C ON OUTPUT A C IF ISIGN = -1, AND LOT GRIDPOINT VECTORS ARE SUPPLIED, C EACH CONTAINING THE COMPLEX SEQUENCE: C C G(0),G(1), ... ,G(N-1) (N COMPLEX VALUES) C C THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH C CONTAINING THE CORRESPONDING N COEFFICIENT VALUES: C C C(0),C(1), ... ,C(N-1) (N COMPLEX VALUES) C C DEFINED BY: C C(K) = SUM(J=0,...,N-1)( G(J)*EXP(-2*I*J*K*PI/N) ) C WHERE I = SQRT(-1) C C C IF ISIGN = +1, AND LOT COEFFICIENT VECTORS ARE SUPPLIED, C EACH CONTAINING THE COMPLEX SEQUENCE: C C C(0),C(1), ... ,C(N-1) (N COMPLEX VALUES) C C THEN THE RESULT CONSISTS OF LOT COMPLEX VECTORS EACH C CONTAINING THE CORRESPONDING N GRIDPOINT VALUES: C C G(0),G(1), ... ,G(N-1) (N COMPLEX VALUES) C C DEFINED BY: C G(J) = SUM(K=0,...,N-1)( G(K)*EXP(+2*I*J*K*PI/N) ) C WHERE I = SQRT(-1) C C C A CALL WITH ISIGN=-1 FOLLOWED BY A CALL WITH ISIGN=+1 C (OR VICE VERSA) RETURNS THE ORIGINAL DATA, MULTIPLIED C BY THE FACTOR N. C C C EXAMPLE GIVEN A 64 BY 9 GRID OF COMPLEX VALUES, STORED IN C A 66 BY 9 COMPLEX ARRAY, A, COMPUTE THE TWO DIMENSIONAL C FOURIER TRANSFORM OF THE GRID. FROM TRANSFORM THEORY, C IT IS KNOWN THAT A TWO DIMENSIONAL TRANSFORM CAN BE C OBTAINED BY FIRST TRANSFORMING THE GRID ALONG ONE C DIRECTION, THEN TRANSFORMING THESE RESULTS ALONG THE C ORTHOGONAL DIRECTION. C C COMPLEX A(66,9), WORK(64,9) C REAL TRIGS1(128), TRIGS2(18) C INTEGER IFAX1(13), IFAX2(13) C C SET UP THE IFAX AND TRIGS ARRAYS FOR EACH DIRECTION: C C CALL CFTFAX(64, IFAX1, TRIGS1) C CALL CFTFAX( 9, IFAX2, TRIGS2) C C IN THIS CASE, THE COMPLEX VALUES OF THE GRID ARE C STORED IN MEMORY AS FOLLOWS (USING U AND V TO C DENOTE THE REAL AND IMAGINARY COMPONENTS, AND C ASSUMING CONVENTIONAL FORTRAN STORAGE): C C U(1,1), V(1,1), U(2,1), V(2,1), ... U(64,1), V(64,1), 4 NULLS, C C U(1,2), V(1,2), U(2,2), V(2,2), ... U(64,2), V(64,2), 4 NULLS, C C . . . . . . . . C . . . . . . . . C . . . . . . . . C C U(1,9), V(1,9), U(2,9), V(2,9), ... U(64,9), V(64,9), 4 NULLS. C C WE CHOOSE (ARBITRARILY) TO TRANSORM FIRST ALONG THE C DIRECTION OF THE FIRST SUBSCRIPT. THUS WE DEFINE C THE LENGTH OF THE TRANSFORMS, N, TO BE 64, THE C NUMBER OF TRANSFORMS, LOT, TO BE 9, THE INCREMENT C BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 1, C AND THE INCREMENT BETWEEN THE STARTING POINTS C FOR EACH TRANSFORM, JUMP, TO BE 66 (THE FIRST C DIMENSION OF A). C C CALL CFFT99( A, WORK, TRIGS1, IFAX1, 1, 66, 64, 9, -1) C C TO TRANSFORM ALONG THE DIRECTION OF THE SECOND SUBSCRIPT C THE ROLES OF THE INCREMENTS ARE REVERSED. THUS WE DEFIN C THE LENGTH OF THE TRANSFORMS, N, TO BE 9, THE C NUMBER OF TRANSFORMS, LOT, TO BE 64, THE INCREMENT C BETWEEN ELEMENTS OF EACH TRANSFORM, INC, TO BE 66, C AND THE INCREMENT BETWEEN THE STARTING POINTS C FOR EACH TRANSFORM, JUMP, TO BE 1 C C CALL CFFT99( A, WORK, TRIGS2, IFAX2, 66, 1, 9, 64, -1) C C THESE TWO SEQUENTIAL STEPS RESULTS IN THE TWO-DIMENSIONA C FOURIER COEFFICIENT ARRAY OVERWRITING THE INPUT C GRIDPOINT ARRAY, A. THE SAME TWO STEPS APPLIED AGAIN C WITH ISIGN = +1 WOULD RESULT IN THE RECONSTRUCTION OF C THE GRIDPOINT ARRAY (MULTIPLIED BY A FACTOR OF 64*9). C C C----------------------------------------------------------------------- C*** Modification 90/11/29 (BTD) C DIMENSION A(1),WORK(1),TRIGS(1),IFAX(1) DIMENSION A(*),WORK(*),TRIGS(*),IFAX(*) C C SUBROUTINE "CFFT99" - MULTIPLE FAST COMPLEX FOURIER TRANSFORM C C A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA C WORK IS AN AREA OF SIZE N*LOT C TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES C IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N C INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' C (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) C JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR C N IS THE LENGTH OF THE DATA VECTORS C LOT IS THE NUMBER OF DATA VECTORS C ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT C = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL C C C VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN C PARALLEL. C C C C C NN = N+N INK=INC+INC JUM = JUMP+JUMP NFAX=IFAX(1) JNK = 2 JST = 2 IF (ISIGN.GE.0) GO TO 30 C C THE INNERMOST TEMPERTON ROUTINES HAVE NO FACILITY FOR THE C FORWARD (ISIGN = -1) TRANSFORM. THEREFORE, THE INPUT MUST BE C REARRANGED AS FOLLOWS: C C THE ORDER OF EACH INPUT VECTOR, C C G(0), G(1), G(2), ... , G(N-2), G(N-1) C C IS REVERSED (EXCLUDING G(0)) TO YIELD C C G(0), G(N-1), G(N-2), ... , G(2), G(1). C C WITHIN THE TRANSFORM, THE CORRESPONDING EXPONENTIAL MULTIPLIER C IS THEN PRECISELY THE CONJUGATE OF THAT FOR THE NORMAL C ORDERING. THUS THE FORWARD (ISIGN = -1) TRANSFORM IS C ACCOMPLISHED C C FOR NFAX ODD, THE INPUT MUST BE TRANSFERRED TO THE WORK ARRAY, C AND THE REARRANGEMENT CAN BE DONE DURING THE MOVE. C JNK = -2 JST = NN-2 IF (MOD(NFAX,2).EQ.1) GOTO 40 C C FOR NFAX EVEN, THE REARRANGEMENT MUST BE APPLIED DIRECTLY TO C THE INPUT ARRAY. THIS CAN BE DONE BY SWAPPING ELEMENTS. C IBASE = 1 ILAST = (N-1)*INK NH = N/2 DO 20 L=1,LOT I1 = IBASE+INK I2 = IBASE+ILAST CDIR$ IVDEP DO 10 M=1,NH C SWAP REAL AND IMAGINARY PORTIONS HREAL = A(I1) HIMAG = A(I1+1) A(I1) = A(I2) A(I1+1) = A(I2+1) A(I2) = HREAL A(I2+1) = HIMAG I1 = I1+INK I2 = I2-INK 10 CONTINUE IBASE = IBASE+JUM 20 CONTINUE GOTO 100 C 30 CONTINUE IF (MOD(NFAX,2).EQ.0) GOTO 100 C 40 CONTINUE C C DURING THE TRANSFORM PROCESS, NFAX STEPS ARE TAKEN, AND THE C RESULTS ARE STORED ALTERNATELY IN WORK AND IN A. IF NFAX IS C ODD, THE INPUT DATA ARE FIRST MOVED TO WORK SO THAT THE FINAL C RESULT (AFTER NFAX STEPS) IS STORED IN ARRAY A. C IBASE=1 JBASE=1 DO 60 L=1,LOT C MOVE REAL AND IMAGINARY PORTIONS OF ELEMENT ZERO WORK(JBASE) = A(IBASE) WORK(JBASE+1) = A(IBASE+1) I=IBASE+INK J=JBASE+JST CDIR$ IVDEP DO 50 M=2,N C MOVE REAL AND IMAGINARY PORTIONS OF OTHER ELEMENTS (POSSIBLY IN C REVERSE ORDER, DEPENDING ON JST AND JNK) WORK(J) = A(I) WORK(J+1) = A(I+1) I=I+INK J=J+JNK 50 CONTINUE IBASE=IBASE+JUM JBASE=JBASE+NN 60 CONTINUE C 100 CONTINUE C C PERFORM THE TRANSFORM PASSES, ONE PASS FOR EACH FACTOR. DURING C EACH PASS THE DATA ARE MOVED FROM A TO WORK OR FROM WORK TO A. C C FOR NFAX EVEN, THE FIRST PASS MOVES FROM A TO WORK IGO = 110 C FOR NFAX ODD, THE FIRST PASS MOVES FROM WORK TO A IF (MOD(NFAX,2).EQ.1) IGO = 120 LA=1 DO 140 K=1,NFAX IF (IGO.EQ.120) GO TO 120 110 CONTINUE CALL VPASSM(A(1),A(2),WORK(1),WORK(2),TRIGS, * INK,2,JUM,NN,LOT,N,IFAX(K+1),LA) IGO=120 GO TO 130 120 CONTINUE CALL VPASSM(WORK(1),WORK(2),A(1),A(2),TRIGS, * 2,INK,NN,JUM,LOT,N,IFAX(K+1),LA) IGO=110 130 CONTINUE LA=LA*IFAX(K+1) 140 CONTINUE C C AT THIS POINT THE FINAL TRANSFORM RESULT IS STORED IN A. C RETURN END SUBROUTINE CFTFAX(N,IFAX,TRIGS) C Modification 90/11/29 (BTD): C DIMENSION IFAX(13),TRIGS(1) DIMENSION IFAX(13),TRIGS(*) C C THIS ROUTINE WAS MODIFIED FROM TEMPERTON"S ORIGINAL C BY DAVE FULKER. IT NO LONGER PRODUCES FACTORS IN ASCENDING C ORDER, AND THERE ARE NONE OF THE ORIGINAL 'MODE' OPTIONS. C C ON INPUT N C THE LENGTH OF EACH COMPLEX TRANSFORM TO BE PERFORMED C C N MUST BE GREATER THAN 1 AND CONTAIN NO PRIME C FACTORS GREATER THAN 5. C C ON OUTPUT IFAX C IFAX(1) C THE NUMBER OF FACTORS CHOSEN OR -99 IN CASE OF ERROR C IFAX(2) THRU IFAX( IFAX(1)+1 ) C THE FACTORS OF N IN THE FOLLOWIN ORDER: APPEARING C FIRST ARE AS MANY FACTORS OF 4 AS CAN BE OBTAINED. C SUBSEQUENT FACTORS ARE PRIMES, AND APPEAR IN C ASCENDING ORDER, EXCEPT FOR MULTIPLE FACTORS. C C TRIGS C 2N SIN AND COS VALUES FOR USE BY THE TRANSFORM ROUTINE C CALL FACT(N,IFAX) K = IFAX(1) IF (K .LT. 1 .OR. IFAX(K+1) .GT. 5) IFAX(1) = -99 IF (IFAX(1) .LE. 0 ) CALL ERRMSG('FATAL','CFFT99F', $ ' FFTFAX - INVALID N') CALL CFTRIG (N, TRIGS) RETURN END SUBROUTINE FACT(N,IFAX) C FACTORIZATION ROUTINE THAT FIRST EXTRACTS ALL FACTORS OF 4 DIMENSION IFAX(13) IF (N.GT.1) GO TO 10 IFAX(1) = 0 IF (N.LT.1) IFAX(1) = -99 RETURN 10 NN=N K=1 C TEST FOR FACTORS OF 4 20 IF (MOD(NN,4).NE.0) GO TO 30 K=K+1 IFAX(K)=4 NN=NN/4 IF (NN.EQ.1) GO TO 80 GO TO 20 C TEST FOR EXTRA FACTOR OF 2 30 IF (MOD(NN,2).NE.0) GO TO 40 K=K+1 IFAX(K)=2 NN=NN/2 IF (NN.EQ.1) GO TO 80 C TEST FOR FACTORS OF 3 40 IF (MOD(NN,3).NE.0) GO TO 50 K=K+1 IFAX(K)=3 NN=NN/3 IF (NN.EQ.1) GO TO 80 GO TO 40 C NOW FIND REMAINING FACTORS 50 L=5 MAX = SQRT(FLOAT(NN)) INC=2 C INC ALTERNATELY TAKES ON VALUES 2 AND 4 60 IF (MOD(NN,L).NE.0) GO TO 70 K=K+1 IFAX(K)=L NN=NN/L IF (NN.EQ.1) GO TO 80 GO TO 60 70 IF (L.GT.MAX) GO TO 75 L=L+INC INC=6-INC GO TO 60 75 K = K+1 IFAX(K) = NN 80 IFAX(1)=K-1 C IFAX(1) NOW CONTAINS NUMBER OF FACTORS RETURN END SUBROUTINE CFTRIG(N,TRIGS) C Modification 90/11/29 (BTD): C DIMENSION TRIGS(1) DIMENSION TRIGS(*) PI=2.0*ASIN(1.0) DEL=(PI+PI)/FLOAT(N) L=N+N DO 10 I=1,L,2 ANGLE=0.5*FLOAT(I-1)*DEL TRIGS(I)=COS(ANGLE) TRIGS(I+1)=SIN(ANGLE) 10 CONTINUE RETURN END SUBROUTINE VPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA) DIMENSION A(N),B(N),C(N),D(N),TRIGS(N) C C SUBROUTINE "VPASSM" - MULTIPLE VERSION OF "VPASSA" C PERFORMS ONE PASS THROUGH DATA C AS PART OF MULTIPLE COMPLEX (INVERSE) FFT ROUTINE C A IS FIRST REAL INPUT VECTOR C B IS FIRST IMAGINARY INPUT VECTOR C C IS FIRST REAL OUTPUT VECTOR C D IS FIRST IMAGINARY OUTPUT VECTOR C TRIGS IS PRECALCULATED TABLE OF SINES & COSINES C INC1 IS ADDRESSING INCREMENT FOR A AND B C INC2 IS ADDRESSING INCREMENT FOR C AND D C INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S C INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S C LOT IS THE NUMBER OF VECTORS C N IS LENGTH OF VECTORS C IFAC IS CURRENT FACTOR OF N C LA IS PRODUCT OF PREVIOUS FACTORS C DATA SIN36/0.587785252292473/,COS36/0.809016994374947/, * SIN72/0.951056516295154/,COS72/0.309016994374947/, * SIN60/0.866025403784437/ C c*** diagnostic c write(0,*)'in vpassm, n=',n c*** M=N/IFAC IINK=M*INC1 JINK=LA*INC2 JUMP=(IFAC-1)*JINK IBASE=0 JBASE=0 IGO=IFAC-1 IF (IGO.GT.4) RETURN GO TO (10,50,90,130),IGO C C CODING FOR FACTOR 2 C 10 IA=1 JA=1 IB=IA+IINK JB=JA+JINK DO 20 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 15 IJK=1,LOT C(JA+J)=A(IA+I)+A(IB+I) D(JA+J)=B(IA+I)+B(IB+I) C(JB+J)=A(IA+I)-A(IB+I) D(JB+J)=B(IA+I)-B(IB+I) I=I+INC3 J=J+INC4 15 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 20 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 40 K=LA1,M,LA KB=K+K-2 C1=TRIGS(KB+1) S1=TRIGS(KB+2) DO 30 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 25 IJK=1,LOT C(JA+J)=A(IA+I)+A(IB+I) D(JA+J)=B(IA+I)+B(IB+I) C(JB+J)=C1*(A(IA+I)-A(IB+I))-S1*(B(IA+I)-B(IB+I)) D(JB+J)=S1*(A(IA+I)-A(IB+I))+C1*(B(IA+I)-B(IB+I)) I=I+INC3 J=J+INC4 25 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 30 CONTINUE JBASE=JBASE+JUMP 40 CONTINUE RETURN C C CODING FOR FACTOR 3 C 50 IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK DO 60 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 55 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I)) C(JB+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I))) C(JC+J)=(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I))) D(JB+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I))) D(JC+J)=(B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I))) I=I+INC3 J=J+INC4 55 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 60 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 80 K=LA1,M,LA KB=K+K-2 KC=KB+KB C1=TRIGS(KB+1) S1=TRIGS(KB+2) C2=TRIGS(KC+1) S2=TRIGS(KC+2) DO 70 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 65 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IC+I)) C(JB+J)= * C1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))) * -S1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))) D(JB+J)= * S1*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))-(SIN60*(B(IB+I)-B(IC+I)))) * +C1*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))+(SIN60*(A(IB+I)-A(IC+I)))) C(JC+J)= * C2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))) * -S2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))) D(JC+J)= * S2*((A(IA+I)-0.5*(A(IB+I)+A(IC+I)))+(SIN60*(B(IB+I)-B(IC+I)))) * +C2*((B(IA+I)-0.5*(B(IB+I)+B(IC+I)))-(SIN60*(A(IB+I)-A(IC+I)))) I=I+INC3 J=J+INC4 65 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 70 CONTINUE JBASE=JBASE+JUMP 80 CONTINUE RETURN C C CODING FOR FACTOR 4 C 90 IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK ID=IC+IINK JD=JC+JINK DO 100 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 95 IJK=1,LOT C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)) D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I)) D(JC+J)=(B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I)) C(JB+J)=(A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I)) C(JD+J)=(A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I)) D(JB+J)=(B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I)) D(JD+J)=(B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I)) I=I+INC3 J=J+INC4 95 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 100 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 120 K=LA1,M,LA KB=K+K-2 KC=KB+KB KD=KC+KB C1=TRIGS(KB+1) S1=TRIGS(KB+2) C2=TRIGS(KC+1) S2=TRIGS(KC+2) C3=TRIGS(KD+1) S3=TRIGS(KD+2) DO 110 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 105 IJK=1,LOT C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) D(JA+J)=(B(IA+I)+B(IC+I))+(B(IB+I)+B(ID+I)) C(JC+J)= * C2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) * -S2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))) D(JC+J)= * S2*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) * +C2*((B(IA+I)+B(IC+I))-(B(IB+I)+B(ID+I))) C(JB+J)= * C1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))) * -S1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))) D(JB+J)= * S1*((A(IA+I)-A(IC+I))-(B(IB+I)-B(ID+I))) * +C1*((B(IA+I)-B(IC+I))+(A(IB+I)-A(ID+I))) C(JD+J)= * C3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))) * -S3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))) D(JD+J)= * S3*((A(IA+I)-A(IC+I))+(B(IB+I)-B(ID+I))) * +C3*((B(IA+I)-B(IC+I))-(A(IB+I)-A(ID+I))) I=I+INC3 J=J+INC4 105 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 110 CONTINUE JBASE=JBASE+JUMP 120 CONTINUE RETURN C C CODING FOR FACTOR 5 C 130 IA=1 JA=1 IB=IA+IINK JB=JA+JINK IC=IB+IINK JC=JB+JINK ID=IC+IINK JD=JC+JINK IE=ID+IINK JE=JD+JINK DO 140 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 135 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I)) C(JB+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))) C(JE+J)=(A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I))) D(JB+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))) D(JE+J)=(B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I))) C(JC+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))) C(JD+J)=(A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I))) D(JC+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))) D(JD+J)=(B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I))) I=I+INC3 J=J+INC4 135 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 140 CONTINUE IF (LA.EQ.M) RETURN LA1=LA+1 JBASE=JBASE+JUMP DO 160 K=LA1,M,LA KB=K+K-2 KC=KB+KB KD=KC+KB KE=KD+KB C1=TRIGS(KB+1) S1=TRIGS(KB+2) C2=TRIGS(KC+1) S2=TRIGS(KC+2) C3=TRIGS(KD+1) S3=TRIGS(KD+2) C4=TRIGS(KE+1) S4=TRIGS(KE+2) DO 150 L=1,LA I=IBASE J=JBASE CDIR$ IVDEP DO 145 IJK=1,LOT C(JA+J)=A(IA+I)+(A(IB+I)+A(IE+I))+(A(IC+I)+A(ID+I)) D(JA+J)=B(IA+I)+(B(IB+I)+B(IE+I))+(B(IC+I)+B(ID+I)) C(JB+J)= * C1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * -S1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) D(JB+J)= * S1*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * -(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * +C1*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * +(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) C(JE+J)= * C4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * -S4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) D(JE+J)= * S4*((A(IA+I)+COS72*(A(IB+I)+A(IE+I))-COS36*(A(IC+I)+A(ID+I))) * +(SIN72*(B(IB+I)-B(IE+I))+SIN36*(B(IC+I)-B(ID+I)))) * +C4*((B(IA+I)+COS72*(B(IB+I)+B(IE+I))-COS36*(B(IC+I)+B(ID+I))) * -(SIN72*(A(IB+I)-A(IE+I))+SIN36*(A(IC+I)-A(ID+I)))) C(JC+J)= * C2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * -S2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) D(JC+J)= * S2*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * -(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * +C2*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * +(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) C(JD+J)= * C3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * -S3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) D(JD+J)= * S3*((A(IA+I)-COS36*(A(IB+I)+A(IE+I))+COS72*(A(IC+I)+A(ID+I))) * +(SIN36*(B(IB+I)-B(IE+I))-SIN72*(B(IC+I)-B(ID+I)))) * +C3*((B(IA+I)-COS36*(B(IB+I)+B(IE+I))+COS72*(B(IC+I)+B(ID+I))) * -(SIN36*(A(IB+I)-A(IE+I))-SIN72*(A(IC+I)-A(ID+I)))) I=I+INC3 J=J+INC4 145 CONTINUE IBASE=IBASE+INC1 JBASE=JBASE+INC2 150 CONTINUE JBASE=JBASE+JUMP 160 CONTINUE RETURN END SUBROUTINE C3DFFT(CZ,I1,I2,I3,I4,I5,I6,I7) C*********************************************************************** C Note: This module should NOT be linked into DDSCAT when C running on a Convex with the vector library -- in that case C link with the -lveclib flag and use the native Convex C3DFFT C routine. C On other systems, however, you need to link in this routine C to avoid an error message from the compiler. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** COMPLEX CZ(1) INTEGER I1,I2,I3,I4,I5,I6,I7 CALL ERRMSG('FATAL','C3DFFT',' CONVEX option but dummy routine!') RETURN END SUBROUTINE WRITENET(CNETFILE,IDNC,CNETFLAG,CSTAMP,CDESCR,CMDFFT, & CALPHA,CSHAPE,CMDTRQ,MXWAV,MXRAD,MXTHET, & MXBETA,MXPHI,MXN3,MXNAT,MXSCA,MXTIMERS,IORTH, & NX,NY,NZ,NAT0,NAT,NAT3,NSCAT,NCOMP,NORI,NWAV, & NRAD,NTHETA,NPHI,NBETA,NTIMERS,MXNX,MXNY,MXNZ, & MXCOMP,BETMID,BETMXD,THTMID,THTMXD, & PHIMID,PHIMXD,ICTHM,IPHIM,TOL,WAVEA,AEFFA, & THETA,BETA,PHI,A1,A2,ICOMP,IXYZ0,IWAV,IRAD, & ITHETA,IBETA,IPHI,IORI,AEFF,WAVE,XX,AK1,BETAD, & THETAD,PHID,AKR,CXE01R,CXE02R,CXRFR,CXEPS, & QEXT,QABS,QSCAT,G,QBKSCA,QPHA,QAV,QTRQAB, & QTRQSC,SM,SMORI,PHIN,THETAN,TIMERS) C*********************************************************************** C Dummy version of subroutine WRITENET C (genuine subroutine WRITENET writes NetCDF portable binary file C to file cnetfile 'dd.cdf') C C This dummy version is provided in order that use of DDSCAT is C possible even if NetCDF software is not installed on system. C C History: C 96.11.15 (PJF): First version of NetCDF (Wolff's request) C 5a6 ---> 5a7 Converted to subroutine. C 96.11.20 (BTD): dummy version created from real WRITENET C end history C Copyright (C) 1996 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C .. Scalar Arguments .. REAL AEFF,AK1,BETAD,BETMID,BETMXD,PHID,PHIMID,PHIMXD,THETAD, & THTMID,THTMXD,TOL,WAVE,XX INTEGER IBETA,ICTHM,IDNC,IORI,IORTH,IPHI,IPHIM,IRAD,ITHETA,IWAV, & MXBETA,MXCOMP,MXN3,MXNAT,MXNX,MXNY,MXNZ,MXPHI,MXRAD,MXSCA, & MXTHET,MXTIMERS,MXWAV,NAT,NAT0,NAT3,NBETA,NCOMP,NORI,NPHI, & NRAD,NSCAT,NTHETA,NTIMERS,NWAV,NX,NY,NZ CHARACTER CALPHA*6,CMDFFT*6,CMDTRQ*6,CNETFLAG*6,CSHAPE*6, & CNETFILE*14,CSTAMP*22,CDESCR*67 C .. Array Arguments .. COMPLEX CXE01R(3),CXE02R(3),CXEPS(MXCOMP),CXRFR(MXCOMP) REAL A1(3),A2(3),AEFFA(MXRAD),AKR(3),BETA(MXBETA),G(2),PHI(MXPHI), & PHIN(MXSCA),QABS(2),QAV(17),QBKSCA(2),QEXT(2),QPHA(2), & QSCAT(2),QTRQAB(3,2),QTRQSC(3,2),SM(MXSCA,4,4), & SMORI(MXSCA,4,4),THETA(MXTHET),THETAN(MXSCA), & TIMERS(MXTIMERS),WAVEA(MXWAV) INTEGER*2 ICOMP(MXN3),IXYZ0(MXNAT,3) C C This subroutine should never be called, so if it is, bail out... C CALL ERRMSG('FATAL','DUMMYWRITENET', & ' NetCDF is not installed. Change to NOTCDF in ddscat.par!') RETURN END