C************************************************************* C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C The development of the MTZ files and associated software mark 1 C is part of the masterplan of the ESF/EACBM Working Group 2.1 for C better Protein Crystallographic software for Europe. C C************************************************************** C C FFTBIG C C In-memory FFT for space group P1, using FFT package from C Andrew McLachlan C C Tue Sep 1 13:16:54 BST 1992 C Modifeid SETGRD & SETLIM C C Subroutines C C (a) bits exclusive to fftbig C fftbig fftb1 C C (b) routines from fft.for slightly modified (& renamed) for this program C cinffb: modified from cinfft, LSG = 1 (number of possible space-groups) C Expand call removed C sfipyb: modified from sfipyc to load up array differently C prnbyz: modified from prntyz to output sections differently C C (c) routines identical to those in fft, included here for completeness C appnds chksym factrz fndsmp gethkl C prmvec prmvet rcsfft setgrd setlim C C (d) routines forming FFT package from ADMcL. Note that although some C of these have the same names as routines in TenEyck's FFT package C on which they are based, they are NOT the same C cmplft diprp gridim hermft inv21 mdftkd C prnt3d r2cftk r3cftk r4cftk r5cftk r8cftk C realft rpcftk rsymft sdiad shrnk2 srfp C stftp1 undmft C PROGRAM FFTBIG C ============== C IMPLICIT NONE C CCC INTEGER ISIZE CCC PARAMETER (ISIZE=10000000) CCC REAL X(ISIZE) INTEGER LSYM PARAMETER (LSYM=192) INTEGER ICEN,IPROJ,ISQ,IW,KSYM,LIST,NABCH,NANOM,NF1,NF2,NH,NI1, + NNA,NNB,NP,NP2,NPH,NSECAX,NSSM,NSG,NS1,NS2,NSYMCH,NW,NW2, + NX,NY,NZ,T2,T3,T4,XMAX,X MIN,YMAX,YMIN,ZMAX,ZMIN,NUSED7, + NUSED8,NUSED9,NUSED10,IUVW,NUSED11,LSPGRP,NUSED6 REAL BIAS,CELL,CS,DEL,DEL1,F1LIM,F2LIM,F1MAX,F2MAX,IPHTR,RHMIN, + RHMAX,RHMEAN,SCLF1,SCLF2,SCMOD1,SCMOD2,SIG1,SIG2,SMAX,SMIN, + SN,T,TFF1,TFF2,CX,CY,CZ COMMON /IOFFT/T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /MPHDR/NX,NY,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,NUSED6, + NUSED7,NUSED8,NUSED9,NUSED10,IUVW(3),NUSED11,CELL(6), + LSPGRP,RHMIN,RHMAX,RHMEAN COMMON /SEL/NF1,NS1,NF2,NS2,NP,NW,IW,ISQ,NH,SMIN,SMAX,SCLF1,TFF1, + SCLF2,TFF2,BIAS,NSSM,NSG,IPROJ,LIST,ICEN,KSYM,NI1,CS(24), + SN(24),T(LSYM,3),NNA,NNB,NPH,NANOM,NW2,NP2,SIG1,SIG2,DEL, + DEL1,F1LIM,F2LIM,IPHTR,SCMOD1,SCMOD2,F1MAX,F2MAX,CX,CY,CZ C INTEGER LDUM,IFAIL EXTERNAL CCPFYP, CCPOPN, CCPRCS, CINFFB, CCPALC, LRCLOS, + CCPERR, MTZINI, FFTB1 INTEGER LENGTH (1) CHARACTER TYPE (1) DATA TYPE /'R'/ C CALL CCPFYP IFAIL = 0 CALL MTZINI C C *********************************** CALL CCPOPN(5,'DATA',5,1,LDUM,IFAIL) CALL CCPOPN(6,'PRINTER',6,1,LDUM,IFAIL) CALL CCPRCS(6,'FFTBIG','$Date: 2002/07/31 16:05:01 $') C *********************************** C C---- Read input including spacegroup number C CALL CINFFB C C---- P1 C LENGTH (1) = NX*NY*(NZ+2) CALL CCPALC (FFTB1, 1, TYPE, LENGTH) CCC CALL FFTB1(X,ISIZE) CALL LRCLOS(1) IF (LIST .NE. 0) CLOSE(UNIT=NABCH,STATUS='KEEP') CALL CCPERR(0,'Normal termination') C END C C C C ============================== SUBROUTINE APPNDS(LINE,STRING) C ============================== C C---- Append "string" to "line" C C .. Scalar Arguments .. CHARACTER LINE* (*),STRING* (*) C .. C .. External Functions .. INTEGER LENSTR EXTERNAL LENSTR C .. C LINE(LENSTR(LINE) + 1:) = STRING END C C C C ================= SUBROUTINE CHKSYM C ================= C C Check hand of permuted reciprocal space symmetry C N1 is first dimension of JS & JT C C---- Loop symmetry operations C C---- Check if the screw axes of permuted symmetry are compatible C with the FFT spacegroup. C It is too complicated to do much more. C C .. Parameters .. INTEGER LSYM PARAMETER (LSYM=192) C .. C .. Scalars in Common .. INTEGER NLAUE,NLAFFT,NSMFFT,NSPFFT,NSYM,NSYMP,NF1,NS1,NF2,NS2,NP, &NW,IW,ISQ,NH,NSSM,NSG,IPROJ,LIST,ICEN,KSYM,NI1,NNA,NNB,NPH,NANOM, &NW2,NP2 REAL SMIN,SMAX,SCLF1,TFF1,SCLF2,TFF2,BIAS,SIG1,SIG2,DEL,DEL1, &F1LIM,F2LIM,IPHTR,SCMOD1,SCMOD2,F1MAX,F2MAX,CX,CY,CZ C .. C .. Arrays in Common .. REAL PERM,RSMFFT,RSYM,RSYMT,CS,SN,T C .. C .. Local Scalars .. REAL D,DCHK INTEGER I,I1,I2,IDCHK,IFFTWN,ISF,ISGN,J,N C .. C .. Local Arrays .. REAL ROTTP(4,4),TR(3) C .. C .. External Subroutines .. EXTERNAL CCPERR C .. C .. Intrinsic Functions .. INTRINSIC ABS,NINT C .. C .. Common blocks .. COMMON /ATSYM/RSYM(4,4,LSYM),RSYMT(4,4,LSYM),NSYM,PERM(4,4),NSYMP, + NLAUE COMMON /FFTSPG/NSPFFT,NSMFFT,RSMFFT(4,4,LSYM),NLAFFT COMMON /SEL/NF1,NS1,NF2,NS2,NP,NW,IW,ISQ,NH,SMIN,SMAX,SCLF1,TFF1, + SCLF2,TFF2,BIAS,NSSM,NSG,IPROJ,LIST,ICEN,KSYM,NI1,CS(24), + SN(24),T(LSYM,3),NNA,NNB,NPH,NANOM,NW2,NP2,SIG1,SIG2,DEL, + DEL1,F1LIM,F2LIM,IPHTR,SCMOD1,SCMOD2,F1MAX,F2MAX,CX,CY,CZ C .. C .. Save statement .. SAVE /ATSYM/, /FFTSPG/, /SEL/ C .. C ISYSW = 6 DO 80 ISF = 1,NSMFFT C C---- Warning flag unset C IFFTWN = 0 C C write (isysw,fmt=8980) ((rsmfft(i,j,isf),j=1,4),i=1,3) C 6000 FORMAT (' Input FFT reciprocal symmetry matrix and realtransla', C + 'tion vector :',/3 (5X,3F8.3,10X,F8.3,/)) C DO 70 N = 1,NSYM C C---- This is being recalculated for each FFT symmetry operator C to reorder hkl multiply rsymt* perm(inverse) C ie perm(transpose) C to reorder real translation ( =rsym(j,4,n)) multiply by perm C DO 20 J = 1,3 TR(J) = PERM(J,1)*RSYM(1,4,N) + PERM(J,2)*RSYM(2,4,N) + + PERM(J,3)*RSYM(3,4,N) ROTTP(J,4) = TR(J) C DO 10 I = 1,3 ROTTP(I,J) = RSYMT(I,1,N)*PERM(J,1) + + RSYMT(I,2,N)*PERM(J,2) + + RSYMT(I,3,N)*PERM(J,3) 10 CONTINUE 20 CONTINUE C IF (ISF.EQ.1) THEN C C write (isysw,fmt=8990) ((rsymt(i,j,n),j=1,3), C + rsym(i,4,n),i=1,3) C 6002 FORMAT (' Input reciprocal symmetry matrix and real transl', C + 'ation vector :',/3 (5X,3F8.3,10X,F8.3,/)) C C write (isysw,fmt=9000) ((rottp(i,j),j=1,3),tr(i),i=1,3) C 6004 FORMAT (' Permuted reciprocal symmetry matrix and real tra', C + 'nslation vector :',/3 (5X,3F8.3,10X,F8.3,/)) C C---- Calculate determinant of permuted matrix C D = 0.0 C DO 30 I = 1,3 I1 = I + 1 IF (I.EQ.3) I1 = 1 I2 = 6 - I - I1 D = (ROTTP(2,I1)*ROTTP(3,I2)-ROTTP(2,I2)*ROTTP(3,I1))* + ROTTP(1,I) + D 30 CONTINUE C IF (ABS(ABS(D)-1.0).GT.0.001) WRITE (ISYSW,FMT=6006) 6006 FORMAT (' Determinant of symop matrix is not + or - 1') ISGN = 1 IF (D.LT.0.0) ISGN = -1 CCC IF (D.LT.0.0) WRITE (ISYSW,FMT=6008) CCC 6008 FORMAT (' ***** MATRIX REVERSES HAND *****') END IF C C---- Now Check if the permuted translations are compatible C with the FFT spacegroup. C IF (NSPFFT.NE.1 .AND. ISQ.EQ.0) THEN C C---- Check real translation elements. C There may be unit cell translations. C DO 40 I = 1,3 J = 4 C C write(isysw,'(a,4i3,2f8.2)') C 1 ' isf i j isgn rottp(i,4) rsmfft(i,4,isf)', C 1 isf,i,j,isgn,rottp(i,4),rsmfft(i,4,isf) C I think? I need to take account of hand of symmetry operator here C DCHK = ROTTP(I,4)*ISGN - RSMFFT(I,4,ISF) IDCHK = NINT(DCHK) DCHK = ABS(DCHK-IDCHK) IF (DCHK.LT.0.01) GO TO 40 C C---- Label 40 means check another hkl symmetry operator.. C GO TO 70 40 CONTINUE C C---- Found a good translation C IFFTWN = 1 C C---- Now check all symmetry matrix elements for this value of isf C DO 60 I = 1,3 DO 50 J = 1,3 C C write(isysw,'(a,4i3,2f8.2)') C 1 ' isf i j isgn rottp(i,j) rsmfft(i,j,isf)', C 1 isf,i,j,isgn,rottp(i,j),rsmfft(i,j,isf) C DCHK = ABS(ROTTP(I,J)*ISGN-RSMFFT(I,J,ISF)) IF (DCHK.LT.0.01) GO TO 50 50 CONTINUE 60 CONTINUE C ELSE IFFTWN = 1 END IF C C---- All tests passed. C 70 CONTINUE C C---- No tests passed. C IF (IFFTWN.EQ.0) CALL CCPERR(1, + ' *** FFT symmetry incompatible with space-group ***' + ) C 80 CONTINUE C END C C C C ================= SUBROUTINE CINFFB C ================= C C Card input for data selection and equivalent reflexion generation C for FFT. Special version for FFTBIG C Differences: C 1) LSG = 1 C 2) Expand bit removed C C Enter with JFLAG=0 for FFT C or .ne.0 for EXPAND (exits early) C C .. Parameters .. INTEGER LSG,LSYM PARAMETER (LSG=1,LSYM=192) INTEGER MCOLS PARAMETER (MCOLS=500) C .. C .. Scalars in Common .. REAL A,B,BIAS,DEL,DEL1,F000,F1LIM,F2LIM,HMIN, + RHOMAX,RHOMIN,RMAX,RMIN,SCLF1,SCLF2,SCMOD1,SCMOD2,SIG1,SIG2, + SMAX,SMIN,SSQ,TFF1,TFF2,V,F1MAX,F2MAX,GRDSMP,CELL,RHMIN, + RHMAX,RHMEAN,CX,CY,CZ INTEGER HMAX,ICEN,IDIG,ILIN,IPHTR,IPROJ,ISORT,ISP,ISQ,IW,J,K,KE, + KMAX,KMIN,KSYM,L,LCMAX,LCMIN,LIST,LMAX,LMIN,M, + NABCH,NANOM,NCOL,NF1,NF2,NH,NI1,NLAUE,NLINES,NNA,NNB,NP, + NP2,NPH,NS1,NS2,NSECAX,NSG,NSMFFT,NSPFFT,NSSM,NSYM,NSYMCH, + NSYMP,NUSED10,NUSED11,NUSED6,NUSED7,NUSED8,NUSED9,NW,NW2, + NX,NY,NZ,T2,T3,T4,TYPE,XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN, + NLAFFT,IUVW,LSPGRP CHARACTER NAMSPG*10,PGNAME*10 C .. C .. Arrays in Common .. REAL CS,PERM,RSMFFT,RSYM,RSYMT,SN,XYZLIM INTEGER KP,MP,T C .. C .. Local Scalars .. REAL XTMP,SAMPLE INTEGER I,I1,I2,I3,IFAIL,J1,LDUM,N, + NSCRCH,NSYMM CHARACTER PATNAM*10 C .. C .. Local Arrays .. REAL TR(3,LSYM) INTEGER JUNK(3),KGPERM(LSG),KSC(LSG), + MCEN(LSG),MSG(LSG),MSP(3,2) CHARACTER KXYZ(3)*1 C .. C .. External Subroutines .. EXTERNAL CCPDPN,CCPERR,CHKSYM,INVSYM,PRMVEC,RCSFFT,SETGRD, + SETLIM,LENSTR INTEGER LENSTR C .. C .. Intrinsic Functions .. INTRINSIC COS,MAX,MIN,NINT,SIN C .. C .. Common blocks .. COMMON /ATSYM/RSYM(4,4,LSYM),RSYMT(4,4,LSYM),NSYM,PERM(4,4),NSYMP, + NLAUE COMMON /FFTPRM/KP(3) COMMON /FFTSPG/NSPFFT,NSMFFT,RSMFFT(4,4,LSYM),NLAFFT COMMON /FOUT/LCMIN,LCMAX,ISP,IDIG,ILIN COMMON /HKLF/A,B,ISORT,J,K,L,M,SSQ COMMON /HKLLIM/HMAX,KMAX,LMAX,HMIN,KMIN,LMIN COMMON /IOFFT/T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /LFTOVR/MP(3),KE,RMAX,RMIN,XYZLIM(2,3),GRDSMP COMMON /MAPFMT/TYPE,NLINES,NCOL,V,F000,RHOMAX,RHOMIN COMMON /MPHDR/NX,NY,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,NUSED6, + NUSED7,NUSED8,NUSED9,NUSED10,IUVW(3),NUSED11,CELL(6), + LSPGRP,RHMIN,RHMAX,RHMEAN COMMON /SEL/NF1,NS1,NF2,NS2,NP,NW,IW,ISQ,NH,SMIN,SMAX,SCLF1,TFF1, + SCLF2,TFF2,BIAS,NSSM,NSG,IPROJ,LIST,ICEN,KSYM,NI1,CS(24), + SN(24),T(LSYM,3),NNA,NNB,NPH,NANOM,NW2,NP2,SIG1,SIG2,DEL, + DEL1,F1LIM,F2LIM,IPHTR,SCMOD1,SCMOD2,F1MAX,F2MAX,CX,CY,CZ COMMON /SPACGP/ NAMSPG, PGNAME CC .. C .. Save statement .. SAVE C .. C .. Data statements .. C DATA MSG/1/ + MCEN/0/ + KSC/1/ DATA KGPERM/1/ DATA MSP/3,1,2,2,1,3/ DATA KXYZ/'X','Y','Z'/ C .. C C---- Space group number MSG, centric flag MCEN C KSC =1 for y sections , =2 for z sections C KGPERM = 1 if permutation of axes allowed, = 0 not C C Number of space-groups recognized C C P1 P-1 P2 P21 P2/M P222 P2221 P21212 P21212 P212121 C2221 C222 C I222 PMMM P41212 P43212 P31 P3 P3121 P3221 C C---- Program axis order : fast,medium,slow axes for y, z sections C ISYSW = 6 NSYMCH = 4 NABCH = 22 NSCRCH = 23 T2 = NSCRCH C DO 10 I = 1,3 KP(I) = I 10 CONTINUE C BIAS = 0.0 F000 = 0.0 IDIG = 4 ILIN = 1 ISQ = 0 IPHTR = 0 SCMOD1 = 0.0 SCMOD2 = 0.0 ISORT = 0 IPROJ = 0 LIST = 0 LCMAX = 0 LCMIN = 0 NCOL = 0 DEL = 0.0 DEL1 = 0.0 NLINES = 0 SIG1 = 0.0 SIG2 = 0.0 F1LIM = 0.0 F2LIM = 0.0 F1MAX = 0.0 F2MAX = 0.0 NSYM = 0 RHOMAX = 999.0 RHOMIN = -999.0 SCLF1 = 1.0 SCLF2 = 1.0 TYPE = 0 T3 = 1 T4 = 0 V = 0.0 NNA = 0 NNB = 0 NF1 = 0 NI1 = 0 NS1 = 0 NF2 = 0 NS2 = 0 NP = 0 NH = 0 NPH = 0 NW = 0 NANOM = 0 NW2 = 0 NP2 = 0 C C---- Control data error counter C KE = 0 C C---- Set one symmetry matrix to identity C DO 30 I = 1,4 DO 20 J = 1,4 IF (I.EQ.J) THEN RSYM(J,I,1) = 1 RSYMT(J,I,1) = 1.0 PERM(J,I) = 1.0 ELSE RSYM(J,I,1) = 0 RSYMT(J,I,1) = 0.0 PERM(J,I) = 0.0 END IF 20 CONTINUE 30 CONTINUE C C---- Table of Sin and Cos in 15 deg steps C DO 40 I = 1,24 CS(I) = COS(0.261799167*I) SN(I) = SIN(0.261799167*I) 40 CONTINUE C C---- read data cards C C ********************* CALL RCSFFT(MSG,LSG) C ********************* C C---- Check & write out input commands C C---- Patterson C IF (ISQ.EQ.1) THEN WRITE (ISYSW,FMT=6000) 6000 FORMAT (' PATTERSON calculated - amplitudes used = F1 squared -' + /' ** All translational elements of symmetry set to zero', + ' **') C Write out harker vectors CALL HARKER ( NSYM,RSYM) C Get Patterson spacegroup CALL PATSGP(NAMSPG, PGNAME, PATNAM, I) IF (I .LE. 0) THEN WRITE (ISYSW, FMT=6001) NAMSPG 6001 FORMAT(//' **** WARNING - no Patterson spacegroup known as ', $ 'equivalent to spacegroup ',A, '****'/ $ ' Spacegroup for map left as Fourier spacegroup'//) ELSE LSPGRP = I WRITE (ISYSW, FMT=6003) PATNAM(1:LENSTR(PATNAM)), LSPGRP 6003 FORMAT(/' Patterson spacegroup is ',A, $ ', number ',I6/) ENDIF END IF C C---- Phased translation function C IF(IPHTR.NE.0) THEN IF(IPHTR.GT.0) THEN WRITE (ISYSW,FMT=6006) + 'Phased translation function: W * F1 * F2 exp (i (PHI - PHI2))' ELSE WRITE (ISYSW,FMT=6006) + 'Phased translation function on other hand: ', + 'W * F1 * F2 exp (i (PHI + PHI2mod))' 6006 FORMAT (//' Coefficients used for Fourier calculation:-',/' **', + '*****************************************',//1X,2A,/) END IF C IF (NF1.EQ.0 .OR. NF2.EQ.0 .OR. NP.EQ.0 .OR. NP2.EQ.0) THEN WRITE (ISYSW,FMT='(A)') + ' ** !!! Define F1 F2 PHI PHI2 for PHTRANSLATION option **' KE = KE + 1 END IF ENDIF IF ((NF1.GT.0.AND.NS1.LE.0) .OR. + (NF2.GT.0.AND.NS2.LE.0)) BIAS = 0.0 C WRITE (ISYSW,FMT=6004) RMIN,RMAX 6004 FORMAT (/' Resolution limits: ',2F9.2) C WRITE (ISYSW,FMT=6116) 'F1',SCLF1,TFF1 6116 FORMAT (' Scale & B for ',A,' : ',2F9.5) C IF (NF2.GT.0) WRITE (ISYSW,FMT=6116) 'F2',SCLF2,TFF2 IF (ABS(BIAS).GT.0.0001) WRITE (ISYSW,FMT=6008) BIAS 6008 FORMAT (' Anomalous bias : ',F5.1,' DO NOT USE unless you K', + 'NOW your Sig(F) are sensible') C WRITE (ISYSW,FMT=6010) 6010 FORMAT (/' * F used = Scale * exp(- B * (sin theta/lambda)**2);', + ' ** Scale and B applied to F BEFORE squaring for Patter', + 'son **',/) C C---- Check on exclusions C IF (NF2 .EQ. 0) THEN IF (DEL1.NE.0.0 .AND. F1MAX.NE.0.0) THEN CALL CCPERR(3, + ' WARNING: possible inconsistency in EXCL: F1MAX cutoff taken') DEL1 = 0.0 ENDIF ELSE IF (ISQ.NE.0 .AND. (DEL.NE.0.0 .AND. DEL1.NE.0.0)) THEN CALL CCPERR(3, + ' WARNING: possible inconsistency in EXCL: DIFF cutoff taken') DEL1 = 0.0 ENDIF ENDIF C---- Only write out those tests which are actually being applied C IF (SIG1.GT.0.0 .OR. SIG2.GT.0.0 .OR. ABS(F1LIM).GT.0.0 .OR. + F2LIM.GT.0.0 .OR. F1MAX.GT.0. .OR. F2MAX.GT.0.0 .OR. + DEL.GT.0.0 .OR. DEL1.GT.0.0) THEN C WRITE (ISYSW,FMT=6012) 6012 FORMAT (/,' *** reflections excluded ***',/) C IF (SIG1 .GT. 0.0) WRITE (6,FMT=6101) SIG1 6101 FORMAT (21X,'if |F1| < ',F6.1,'* |sd(F1)|') IF (SIG2 .GT. 0.0) WRITE (6,FMT=6102) SIG2 6102 FORMAT (21X,'if |F2| < ',F6.1,'* |sd(F2)|') IF (NF1.GT.0 .AND.F1LIM .GT. 0.0) WRITE (6,FMT=6103) F1LIM 6103 FORMAT (21X,'if |F1| < ',F10.1) IF (NI1.GT.0 .AND.ABS(F1LIM) .GT. 0.0) WRITE (6,FMT=6123) F1LIM 6123 FORMAT (21X,'if I < ',F10.1) IF (F2LIM .GT. 0.0) WRITE (6,FMT=6104) F2LIM 6104 FORMAT (21X,'if |F2| < ',F10.1) IF (F1MAX .GT. 0.0) THEN IF(NF1.GT.0) WRITE (6,FMT=6013) F1MAX IF(NI1.GT.0) WRITE (6,FMT=6023) F1MAX 6013 FORMAT(21X,'if F1 > ',F10.1) 6023 FORMAT(21X,'if I > ',F10.1) ENDIF IF (F2MAX .GT. 0.0) THEN WRITE (6,FMT=6014) F2MAX 6014 FORMAT(21X,'or F2 .gt. ',F10.1) ENDIF IF (DEL .GT. 0.0) WRITE (6,FMT=6011) DEL IF (DEL1 .GT. 0.0) WRITE (6,FMT=6105) DEL1 6011 FORMAT (21X,'if differences |F1-F2| > ',F12.2) 6105 FORMAT (21X,'if absolute magnitude of Sqrt(A*A +B*B) > ', + F12.2) WRITE (6,FMT=6106) 6106 FORMAT (/,/) C ELSE WRITE (6,FMT=6107) 6107 FORMAT (/,' *** No reflections excluded ***',/) ENDIF C IFAIL = 0 C IF (LIST.NE.0) THEN CALL CCPDPN(NABCH,'ABCOEFFS','UNKNOWN','F',LDUM, IFAIL) WRITE(NABCH,'(A,A)') + ' H K L k1*F1 PHI k2*F2 PHI2 ', + 'W1 W2 A B ATAN(B,A) SQRT(AA+BB)' ENDIF C IF (LSPGRP.EQ.0) LSPGRP = NSPFFT WRITE (ISYSW,FMT=6015) NSYM,NSPFFT,LSPGRP,LIST 6015 FORMAT (///' Number of symmetry operations =',I2,' Space-group f', + 'or FFT =',I4,' True space-group =',I4,' List terms gt', + I6) C C---- check error flags C IF (NF1.GT.0 .OR. NF2.LE.0) THEN IF (NH.LE.0 .OR. (NF1.GT.0.AND.NF2.GT.0.AND.NP.GT.0)) THEN IF (NW.GT.0 .OR. IW.EQ.0) THEN IF ((NW2.LE.0.AND.NP2.LE.0) .OR. NF2.NE.0) GO TO 50 END IF END IF END IF C WRITE (ISYSW,FMT=6016) 6016 FORMAT (' ** !!! Invalid combination of data items **') KE = KE + 1 C 50 NSYMM = NSYM IF (NSYMM.EQ.0) NSYMM = 1 C DO 70 N = 1,NSYMM DO 60 I = 1,3 C C--- Set translational part of symmetry operator = 0 for Patterson C IF (ABS(ISQ).EQ.1.OR.IPHTR.NE.0) RSYM(I,4,N) = 0.0 TR(I,N) = RSYM(I,4,N) T(N,I) = NINT(TR(I,N)*24.0) 60 CONTINUE C C ******************************** CALL INVSYM(RSYM(1,1,N),RSYMT(1,1,N)) C ******************************** C 70 CONTINUE C KSYM = NSYM C DO 80 NSG = 1,LSG IF (MSG(NSG).EQ.NSPFFT) GO TO 90 80 CONTINUE C WRITE (ISYSW,FMT=6018) 6018 FORMAT (' ** !!! Illegal space group **') NSG = 1 KE = KE + 1 90 CONTINUE C ICEN = MCEN(NSG) C C---- Get axis permutation, read real axes for fast, medium slow C IF (MP(1).EQ.0 .OR. MP(2).EQ.0 .OR. MP(3).EQ.0) THEN C C---- Blank permutation, use default according to whether program is C sectioning along y or z C J = KSC(NSG) C DO 100 J1 = 1,3 MP(J1) = MSP(J1,J) 100 CONTINUE ELSE C C---- Only permit permutation in certain spacegroups C IF (KGPERM(NSG).EQ.0) THEN WRITE (ISYSW,FMT='(A)') + ' *** !!Axis permutation not allowed in this spacegroup ***' KE = KE + 1 END IF END IF C I1 = MP(1) I2 = MP(2) I3 = MP(3) WRITE (ISYSW,FMT=6020) KXYZ(I1),KXYZ(I2),KXYZ(I3) 6020 FORMAT (/' Axis order in map: fast=',A1,', medium= ',A1,', slo', + 'w(section)= ',A1,/) C C---- Now generate program axis permutation C y-section spacegroups have fast medium slow axes z,x,y C z-section y,x,z C J = KSC(NSG) C C---- NSECAX = 2 for "printy" calls; 3 for "printz" calls C NSECAX = J + 1 C DO 110 I = 1,3 J1 = MSP(I,J) KP(J1) = MP(I) 110 CONTINUE C C Fudge Laue group for certain cases IF (NSPFFT .EQ. 47) NLAUE = 14 IF (NSPFFT .EQ. 10) NLAUE = 5 C C---- Set default grid, if not already set C IF (NX.LE.0) THEN C C SETGRD needs to work on true axes C SETGRD uses Laue group, returns sample grid C C NXMIN etc = 2hmax + 1 - set JUNK(.) = this. JUNK(1) = 2*HMAX + 1 JUNK(2) = 2*KMAX + 1 JUNK(3) = 2*LMAX + 1 C Default desired sample = approx 3 * maximum index c (for SETGRD, SAMPLE is 1/2 desired sampling) SAMPLE = GRDSMP/2.0 C C ***************************************************** CALL SETGRD(NLAUE,SAMPLE,JUNK(1),JUNK(2),JUNK(3),NX,NY,NZ) C ***************************************************** C IF (NX.LE.0) THEN WRITE (ISYSW,FMT='(A,I5)') + ' ** !!!Unrecognized Laue group ** ', + NLAUE KE = KE + 1 END IF END IF C C---- Set up map volume if not set C IF (XYZLIM(1,1).LT.-0.9) THEN C C---- No volume has been set by the user from the keywords C The default is to give the whole unit cell C XYZLIM(1,1) = 0.0 XYZLIM(2,1) = 1.0 XYZLIM(1,2) = 0.0 XYZLIM(2,2) = 1.0 XYZLIM(1,3) = 0.0 XYZLIM(2,3) = 1.0 C WRITE (ISYSW,FMT=7003) XYZLIM 7003 FORMAT (/' ** The limits of the map have been set to those of', + /' the unit cell:', + //' Min X Max X Min Y Max Y Min Z Max Z', + /,1X,6F8.3, + //' (expressed in fractional coordinates and with x,y,z', + /' referring to the unpermutated axes.)',/) C END IF C C---- Set up limits in grid units C XTMP = (XYZLIM(2,1)-XYZLIM(1,1))* (XYZLIM(2,2)-XYZLIM(1,2))* + (XYZLIM(2,3)-XYZLIM(1,3)) C IF ((XYZLIM(2,1)-XYZLIM(1,1)).LT.0.0000001) THEN XTMP = (XYZLIM(2,2)-XYZLIM(1,2))* (XYZLIM(2,3)-XYZLIM(1,3)) ELSE IF ((XYZLIM(2,2)-XYZLIM(1,2)).LT.0.0000001) THEN XTMP = (XYZLIM(2,1)-XYZLIM(1,1))* (XYZLIM(2,3)-XYZLIM(1,3)) ELSE IF ((XYZLIM(2,3)-XYZLIM(1,3)).LT.0.0000001) THEN XTMP = (XYZLIM(2,2)-XYZLIM(1,2))* (XYZLIM(2,1)-XYZLIM(1,1)) END IF C IF (XTMP.GT.5.0) THEN C C---- more than 5 "gridpoints" in map, limits were input in grid units C XMIN = NINT(XYZLIM(1,1)) XMAX = NINT(XYZLIM(2,1)) YMIN = NINT(XYZLIM(1,2)) YMAX = NINT(XYZLIM(2,2)) ZMIN = NINT(XYZLIM(1,3)) ZMAX = NINT(XYZLIM(2,3)) ELSE C C---- limits were input as fractional limits C XMIN = NINT(XYZLIM(1,1)*NX) XMAX = INT(XYZLIM(2,1)*NX) YMIN = NINT(XYZLIM(1,2)*NY) YMAX = INT(XYZLIM(2,2)*NY) ZMIN = NINT(XYZLIM(1,3)*NZ) ZMAX = INT(XYZLIM(2,3)*NZ) END IF C C---- Force limits into range 0 -> NX-1 etc C XMIN = MAX(XMIN,0) XMAX = MIN(XMAX,NX-1) YMIN = MAX(YMIN,0) YMAX = MIN(YMAX,NY-1) ZMIN = MAX(ZMIN,0) ZMAX = MIN(ZMAX,NZ-1) C C---- If projection flag, fix up value for section axis C IPROJ = unpermuted projection axis C IF (IPROJ.NE.0) THEN IPROJ = MP(3) WRITE (ISYSW,FMT=6022) KXYZ(IPROJ) 6022 FORMAT (/' Map will be a projection along ',A1,/ + ' SECTION LIMITS set to zero') C C---- Get allowed minimum sample intervals C CALL SETGRD(NLAUE,1.0,0,0,0,JUNK(1),JUNK(2),JUNK(3)) C IF (IPROJ.EQ.1) THEN XMIN = 0 XMAX = 0 NX = JUNK(KP(1)) HMAX = 0 END IF IF (IPROJ.EQ.2) THEN YMIN = 0 YMAX = 0 NY = JUNK(KP(2)) KMAX = 0 END IF IF (IPROJ.EQ.3) THEN ZMIN = 0 ZMAX = 0 NZ = JUNK(KP(3)) LMAX = 0 END IF END IF C I1 = KP(1) I2 = KP(2) I3 = KP(3) WRITE (ISYSW,FMT=6024) KXYZ(I1),KXYZ(I2),KXYZ(I3) 6024 FORMAT (/' Permutation of axes: (program=input) x= ',A1,', y= ', + A1,', z= ',A1,/) C C---- Permute symmetry operations if symmetry operations are given C C---- reset permutation matrix C DO 130 I = 1,4 DO 120 J = 1,4 PERM(J,I) = 0.0 120 CONTINUE 130 CONTINUE C PERM(1,I1) = 1.0 PERM(2,I2) = 1.0 PERM(3,I3) = 1.0 PERM(4,4) = 1.0 C IF (NSYM.NE.0) THEN C DO 140 N = 1,NSYM C C ********************* CALL PRMVEC(PERM,T,N,LSYM) C ********************* C 140 CONTINUE C C ****** CALL CHKSYM C ****** C END IF C IF (ABS(BIAS).GT.0.0001) BIAS = BIAS**2 C C---- Exit for EXPAND C DL 15/11/95: This comment looks bogus; the file number (IFIL) C indicates whence GETHKL was called. C Use hmax=0 as flag that GETHKL has been called from expand C IF (NSYM.NE.0) ISORT = 1 IF (I1.NE.1 .OR. I2.NE.2 .OR. I3.NE.3) THEN ISORT = 1 IF (NSYM.EQ.0) THEN WRITE (ISYSW,FMT= + '('' ** !!! You need symmetry for axis permutation **'')') KE = KE + 1 END IF END IF C IF (LSPGRP.NE.NSPFFT) ISORT = 1 C C---- Check minimum & maximum limits C IF (XMIN.GE.0) THEN IF (YMIN.GE.0) THEN IF (ZMIN.GE.0) THEN IF (XMAX.LT.NX) THEN IF (YMAX.LT.NY) THEN IF (ZMAX.LT.NZ) GO TO 170 END IF END IF END IF END IF END IF C WRITE (ISYSW,FMT=6026) 6026 FORMAT (/' !!!!!! Illegal limits : XMIN,YMIN,ZMIN must be .ge. 0', + ', XMAX,YMAX,ZMAX must be .lt. NX,NY,NZ !!!!!!',/) KE = KE + 1 C C---- Check for sampling too coarse C 170 IF (NX.GE.2*HMAX+1) THEN IF (NY.GE.2*KMAX+1) THEN IF (NZ.GE.2*LMAX+1) GO TO 180 END IF END IF C WRITE (ISYSW,FMT=6028) 6028 FORMAT (/'!!!!!! Sampling too coarse : NX,NY,NZ must be .gt. 2*H', + 'MAX+1, 2*KMAX+1, 2*LMAX+1 !!!!!!',/) KE = KE + 1 C C---- Print input limits etc before permutation C 180 CONTINUE WRITE (ISYSW,FMT=6030) 6030 FORMAT (' Before permutation :') C WRITE (ISYSW,FMT=6032) HMAX,KMAX,LMAX,NX,NY,NZ,XMIN,XMAX, + YMIN,YMAX,ZMIN,ZMAX 6032 FORMAT (11X,'Maximum indices hkl ...................',3I5,/11X, + 'Sampling intervals on xyz .............',3I5,/11X, + 'Map limits in grid points on xyz ......',3 (2I5,3X),/) C C write (isysw,fmt=9140) type,nlines,ncol,v,f000,rhomin,rhomax,lcmin, C + lcmax,isp,idig,ilin CCC 6034 FORMAT (/' TYPE NLINES NCOL',7X,'V',8X,'F000',8X,'RHOMIN RHOMAX', CCC + ' MIN MAX ISP IDIG ILIN',/1X,I3,I7,I5,2X,2E12.5,2F9.3, CCC + 5I5,/) C IF (T3.NE.0) WRITE (ISYSW,FMT=6036) T3 6036 FORMAT (' Binary map will be written to stream',I4) IF (T4.NE.0) WRITE (ISYSW,FMT=6038) 6038 FORMAT (' Map will be printed to line-printer') C IF (ISORT.EQ.0) WRITE (ISYSW,FMT=6040) 6040 FORMAT (/' Data is assumed to be sorted') IF (ISORT.NE.0) WRITE (ISYSW,FMT=6042) 6042 FORMAT (/' Data is assumed to be unsorted') C IF (ISORT.EQ.0 .AND. PERM(1,1).NE.1) WRITE (ISYSW,FMT=6044) C 6044 FORMAT (//' ???? !!!!! Data is assumed to be sorted, but ', C + 'you are permuting indices. Are you sure you have resor', C + 'ted data???!!!!') C IF (ISORT.EQ.0 .AND. PERM(1,1).NE.1) NSYM = 1 IF (ISORT.EQ.0 .AND. PERM(1,1).NE.1) KSYM = 1 C C---- Store sampling intervals before permutation for output C C---- Now permute limits etc and print C C---- Patch included to conform to F77 Standards C It is not permitted to pass to a subroutine a C single variable (where an array is expected) and C rely on the program selecting the appropriate number C of variables from a labelled COMMON BLOCK C JUNK(1) = NX JUNK(2) = NY JUNK(3) = NZ CALL PRMVEC(PERM,JUNK,1,1) NX = JUNK(1) NY = JUNK(2) NZ = JUNK(3) JUNK(1) = HMAX JUNK(2) = KMAX JUNK(3) = LMAX CALL PRMVEC(PERM,JUNK,1,1) HMAX = JUNK(1) KMAX = JUNK(2) LMAX = JUNK(3) JUNK(1) = XMIN JUNK(2) = YMIN JUNK(3) = ZMIN CALL PRMVEC(PERM,JUNK,1,1) XMIN = JUNK(1) YMIN = JUNK(2) ZMIN = JUNK(3) JUNK(1) = XMAX JUNK(2) = YMAX JUNK(3) = ZMAX CALL PRMVEC(PERM,JUNK,1,1) XMAX = JUNK(1) YMAX = JUNK(2) ZMAX = JUNK(3) WRITE (ISYSW,FMT=6046) 6046 FORMAT (' After permutation :') WRITE (ISYSW,FMT=6032) HMAX,KMAX,LMAX,NX,NY,NZ,XMIN,XMAX, + YMIN,YMAX,ZMIN,ZMAX C IF (KE.NE.0) THEN WRITE (ISYSW,FMT=6048) KE 6048 FORMAT (' Note',I3,' fatal error messages above') CALL CCPERR(1,' ** Fatal error in input **') END IF C END C VAXPDS NAME=CMPLFT C*********************************************************************** SUBROUTINE CMPLFT(X,Y,N,IDIM,NWRITE) C----------------------------------------------------------------------- C COMPLEX FOURIER TRANSFORM OF 3-DIMENSIONAL ARRAY C LAST UPDATED A.D. MCLACHLAN 27 AUG 1987 C----------------------------------------------------------------------- INTEGER N,IDIM,NWRITE REAL X(*), Y(*) DIMENSION IDIM(10) C---------------------- INTEGER NDIM(5),LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW EQUIVALENCE (LSEPU,NDIM(1)),(LSEPV,NDIM(2)),(LSEPW,NDIM(3)) EQUIVALENCE (LGLIMV,NDIM(4)),(LGLIMW,NDIM(5)) C----- INTEGER MFACTR,MFACT1,M2GRP,MPORD,MPMAX PARAMETER(MFACTR=32,MFACT1=MFACTR+1) PARAMETER(M2GRP=8,MPORD=8,MPMAX=19) LOGICAL ERROR INTEGER NPMAX, NPSYM, N2GRP, NPORD INTEGER NFACTR(MFACT1), NSYM(MFACT1), NUNSYM(MFACT1) INTEGER NPRLST(MPORD) DATA NPRLST/2,3,5,7,11,13,17,19/ C----------------------------------------------------------------------- C IN/OUT -- X(*) REAL PART OF DATA/TRANSFORM C IN/OUT -- Y(*) IMAGINARY PART OF DATA/TRANSFORM C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- GRIDIM CALCULATE DIMENSIONING CONSTANTS C CALLS -- SRFP FACTORISATION OF N C CALLS -- MDFTKD FFT DRIVER C CALLS -- DIPRP REORDERING OF OUTPUT C----------------------------------------------------------------------- C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C------------------------------------------------------------------------------- C NOTE -- NEW CALLS AND INTEGER VARIABLES C NOTE -- TRANSFORMS ONE DIMENSION OF MULTI-DIMENSIONAL DATA C NOTE -- MODIFIED BY L. F. TEN EYCK FROM A ONE-DIMENSIONAL VERSION WRITTEN C NOTE -- BY G. T. SANDE, 1969. DIMENSIONING CHANGED BY A.D.MCLACHLAN, 1981. C NOTE -- C NOTE -- THIS PROGRAM CALCULATES THE TRANSFORM C NOTE -- (X(T) + I*Y(T))*(COS(2*PI*T/N) - I*SIN(2*PI*T/N)) C NOTE -- C----- C NOTE -- NPMAX IS THE LARGEST PRIME FACTOR THAT WILL BE TOLERATED BY THIS C NOTE -- PROGRAM. C NOTE -- NPMAX IS THE NPORD'TH PRIME NUMBER IN NPRLST LIST C NOTE -- N2GRP IS THE LARGEST POWER OF TWO THAT IS TREATED AS A SPECIAL C NOTE -- CASE. C NOTE -- NFACTR(MFACT1) ALLOWS UP TO MFACTR FACTORS FOR THE ARRAY C NOTE -- DIMENSION N C----------------------------------------------------------------------------- 20 FORMAT (1X,'**CMPLFT** ERROR: INVALID NUMBER OF POINTS N =', I10) 21 FORMAT(1X,'**CMPLFT** ERROR IN DIMENSIONS OF ARRAYS N=',I5, 1 /1X,' IDIM= ',3I5,1X,3I5,1X,3I5,1X,I5) C----------------------------------------------------------------------------- NPMAX = MPMAX NPORD=MPORD N2GRP= M2GRP ERROR=.FALSE. IF (N .LE. 1) GO TO 100 C --SET GRID DIMENSIONS CALL GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) IF(ERROR) GO TO 400 C --FACTORISE N, COLLECTING FACTORS IN PAIRS CALL SRFP(N,NPMAX,NPORD,N2GRP,NFACTR,NSYM,NPSYM,NUNSYM, 1 NPRLST, ERROR,NWRITE) IF (ERROR) GO TO 200 C --TRANSFORMS BY THE VARIOUS FACTORS IN TURN CALL MDFTKD(N,NFACTR,NDIM,X,Y,NWRITE) C --REODERING OF PERMUTED FOURIER COEFFICIENTS CALL DIPRP(N,NSYM,NPSYM,NUNSYM,NDIM,X,Y,NWRITE) 100 CONTINUE RETURN C---------------------------- C --FACTORS OF "N" IN ERROR 200 CONTINUE WRITE (NWRITE, 20) N CALL CCPERR (1, 'FFTBIG: Internal error') C --GRID DIMENSIONS IN ERROR 400 CONTINUE WRITE(NWRITE,21) N, IDIM CALL CCPERR (1, 'FFTBIG: Internal error') END C C VAXPDS NAME=DIPRP C******************************************************************** SUBROUTINE DIPRP(NPTS,NSYM,NPSYM,NUNSYM,NDIM,X,Y, 1 NWRITE) C-------------------------------------------------------------------- C DOUBLE IN PLACE REORDERING PROGRAMME C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987. C-------------------------------------------------------------------- INTEGER MFCMAX,MFCMX1 PARAMETER(MFCMAX=32,MFCMX1=MFCMAX+1) INTEGER NPTS,NPSYM REAL X(*), Y(*) INTEGER NSYM(MFCMX1), NUNSYM(MFCMX1),NDIM(5),NWRITE C------------------------------- REAL T LOGICAL ONEMOD INTEGER MODULO (MFCMX1) INTEGER NDK,JJ,KK,NLK,MODS,MULT,NFCMAX,NPNSYM,NTEST INTEGER DELTA, P, P1, P5, LDXV, LDXW, LEVEL INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW INTEGER LINDU,LINDV,NFCMX1,J,NDEEP,N,K,JL,MS,I,L,M INTEGER ISTEP(MFCMX1),IHIGH(MFCMX1) INTEGER ILOOP(MFCMX1),ILOW(MFCMX1) C-------------------------------------------------------------------- C INPUT -- NPTS LENGTH OF TRANSFORM C INPUT -- NSYM(MFCMX1) LIST OF PAIRED FACTORS C INPUT -- NPSYM PRODUCT OF PAIRED FACTORS ONCE EACH C INPUT -- NUNSYM(MFCMX1) LIST OF SINGLE FACTORS (GOING UP) C INPUT -- NDIM(5) DIMENSIONING CONSTANTS C IN/OUT -- X(*) REAL PART OF DATA C IN/OUT -- Y(*) IMAGINARY PART OF DATA C INPUT -- NWRITE PRINT UNIT (**NOT USED**) C-------------------------------------------------------------------- C CALLS -- *** C-------------------------------------------------------------------- C NOTE -- DIMENSIONS REVISED A.D.MCLACHLAN AUGUST 1981 C NOTE -- NESTED LOOPS REVISED A.D. MCLACHLAN DEC 1984 C-------------------------------------------------------------------- C NOTE -- NSYM(J) LISTS PAIRED FACTORS TWICE OVER: C NOTE -- (1) NP FACTORS DESCENDING C NOTE -- (2) MIDDLE TERM NPTS/(NPSYM**2) IF NQ.NE.0 C NOTE -- (3) NP FACTORS ASCENDING C NOTE -- NUNSYM(J) LISTS SINGLE FACTORS ASCENDING C NOTE -- EXCEPT IF ONLY ONE EXISTS C NOTE -- NPSYM IS PRODUCT OF PAIRED FACTORS ONCE EACH C-------------------------------------------------------- C --UP TO MFCMAX FACTORS ALLOWED NFCMAX=MFCMAX NFCMX1=NFCMAX+1 LSEPU=NDIM(1) LSEPV=NDIM(2) LSEPW=NDIM(3) LGLIMV=NDIM(4) LGLIMW=NDIM(5) C --DEAL WITH ALL REPEATED FACTORS (IF ANY) C --THERE ARE NDEEP OF THEM C --FIND NDEEP DO 100 J=1,NFCMX1 IF(NSYM(J).EQ.0) GO TO 110 100 CONTINUE J=NFCMX1 110 CONTINUE IF(J.EQ.1) GO TO 500 NDEEP=J-1 C --SET LIMITS FOR NESTED LOOPS C --IN REVERSE ORDER N=NPTS DO 200 J=1,NDEEP JJ=NDEEP+1-J IHIGH(JJ)=N N=N/NSYM(J) ISTEP(JJ)=N 200 CONTINUE C----- JJ=0 C%% WRITE(NWRITE,*)(NSYM(J),J=1,10) C%% WRITE(NWRITE,*)(NUNSYM(J),J=1,10) C%% WRITE(NWRITE,*) NPTS,NPSYM,NDIM C%% WRITE(NWRITE,*) IHIGH C%% WRITE(NWRITE,*) ISTEP C%% WRITE(NWRITE,*) NDEEP C--------------------------------------------------- C --NESTED LOOPS DONE AS A STACK WITH NDEEP LEVELS C --START ILOW(1)=1 ISTEP(1)=1 LEVEL=1 C --SET LOOP INDEX ON ENTRY BUT DO NOT ADVANCE 300 CONTINUE ILOOP(LEVEL)=ILOW(LEVEL) GO TO 315 C --ADVANCE INDEX 310 CONTINUE ILOOP(LEVEL)=ILOOP(LEVEL)+ISTEP(LEVEL) 315 CONTINUE C --TEST FOR LOOP FINISHED IF(ILOOP(LEVEL).GT.IHIGH(LEVEL)) GO TO 410 C --DEEPEN IF(LEVEL.LT.NDEEP) THEN LEVEL=LEVEL+1 C --SET LOWER LIMIT FOR CURRENT LEVEL ILOW(LEVEL)=ILOOP(LEVEL-1) GO TO 300 END IF C --INNERMOST OPERATION AT DEEPEST LEVEL N=ILOOP(NDEEP) JJ=JJ+1 C%% WRITE(NWRITE,*)(ILOOP(J),J=1,NDEEP),JJ IF (JJ.GE.N) GO TO 400 DELTA = (N-JJ)*LSEPU P1 = (JJ-1)*LSEPU + 1 LINDU=P1-LSEPV-LSEPW DO 350 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 350 LDXW=LSEPW,LGLIMW,LSEPW P=LINDV+LDXW P5 = P + DELTA T = X(P) X(P) = X(P5) X(P5) = T T = Y(P) Y(P) = Y(P5) Y(P5) = T 350 CONTINUE 400 CONTINUE C --LOOP REPEAT GO BACK GO TO 310 C --LOOP FINISHED 410 CONTINUE C --SET NEXT LEVEL UP LEVEL=LEVEL-1 C --TEST FOR OUTERMOST FINISH IF(LEVEL.EQ.0) GO TO 420 C --CONTINUE LOOPING ON CURRENT LEVEL GO TO 310 420 CONTINUE C--------------------------------------------- 500 CONTINUE IF (NUNSYM(1).EQ.0) GO TO 1900 C --DEAL WITH SINGLE FACTORS C --NPNSYM IS PRODUCT OF SINGLE FACTORS C --USE IHIGH(K) TO STORE MODULI NPNSYM=NPTS/(NPSYM**2) MULT=NPNSYM/NUNSYM(1) NTEST=(NUNSYM(1)*NUNSYM(2)-1)*MULT*NPSYM NLK=MULT NDK=MULT DO 600 K=2,NFCMAX IF (NUNSYM(K).EQ.0) GO TO 700 NLK=NLK*NUNSYM(K-1) NDK=NDK/NUNSYM(K) IHIGH(K)=(NLK-NDK)*NPSYM MODS=K 600 CONTINUE 700 CONTINUE ONEMOD=MODS.LT.3 IF (ONEMOD) GO TO 900 DO 800 J=3,MODS JJ=MODS+3-J MODULO(JJ)=IHIGH(J) 800 CONTINUE 900 CONTINUE MODULO(2)=IHIGH(2) JL=(NPNSYM-3)*NPSYM MS=NPNSYM*NPSYM C DO 1800 J=NPSYM,JL,NPSYM K=J C 1000 CONTINUE K=K*MULT IF (ONEMOD) GO TO 1200 C --TAKE REMAINDERS OF K IN TURN DO 1100 I=3,MODS K=K-(K/MODULO(I))*MODULO(I) 1100 CONTINUE 1200 CONTINUE IF (K.GE.NTEST) GO TO 1300 K=K-(K/MODULO(2))*MODULO(2) GO TO 1400 1300 CONTINUE K=K-(K/MODULO(2))*MODULO(2)+MODULO(2) 1400 CONTINUE IF (K.LT.J) GO TO 1000 C IF (K.EQ.J) GO TO 1700 DELTA = (K-J)*LSEPU DO 1600 L=1,NPSYM DO 1500 M=L,NPTS,MS P1 = (M+J-1)*LSEPU + 1 LINDU=P1-LSEPV-LSEPW DO 1500 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 1500 LDXW=LSEPW,LGLIMW,LSEPW JJ=LINDV+LDXW KK = JJ + DELTA T=X(JJ) X(JJ)=X(KK) X(KK)=T T=Y(JJ) Y(JJ)=Y(KK) Y(KK)=T 1500 CONTINUE 1600 CONTINUE 1700 CONTINUE 1800 CONTINUE C 1900 CONTINUE RETURN END C SUBROUTINE FFTB1(SIZE,X) C ======================= C C---- Fourier synthesis for space group P1, incore version C C Use Andrew McLachlan's in-core FFT routines C C .. Scalar Arguments .. INTEGER SIZE C .. C .. Array Arguments .. REAL X(SIZE) C .. C .. Scalars in Common .. INTEGER HMAX,HMIN,IXMAX,IXMIN,IYMAX,IYMIN,IZMAX,IZMIN,KMAX,KMIN, + LMAX,LMIN,NABCH,NSECAX,NSYMCH,NX,NY,NZ,T2,T3,T4 CHARACTER TITLE*80 C .. C .. Arrays in Common .. INTEGER JMPHDR C .. C .. Local Scalars .. REAL RHMAX,RHMEAN,RHMIN,RHRMS INTEGER LMODE,LSPGRP,M,NSEC,NU1,NU2,NV1,NV2,NW1,R,NNX,NNY,NNZ, + NWRITE,IPRINT,IERR C .. C .. Local Arrays .. REAL CELL(6) INTEGER IUVW(3),MXYZ(3) C .. C .. External Subroutines .. EXTERNAL MRDHDR,PRNBYZ C .. C .. Intrinsic Functions .. INTRINSIC MAX C .. C .. Common blocks .. COMMON /HKLLIM/HMAX,KMAX,LMAX,HMIN,KMIN,LMIN COMMON /IOFFT/T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /MPHDR/NX,NY,NZ,IXMIN,IXMAX,IYMIN,IYMAX,IZMIN,IZMAX, + JMPHDR(19) COMMON /MPHDRR/TITLE C .. SAVE /HKLLIM/, /IOFFT/, /MPHDR/, /MPHDRR/ C ISYSW = 6 WRITE (ISYSW,FMT=6000) SIZE C C---- Set index limits C HMIN = -HMAX KMIN = -KMAX LMIN = 0 C C---- Check sampling C IF (NZ.EQ.NZ/2*2) THEN C IF (NX*NY*NZ .GT. SIZE) THEN WRITE (ISYSW,FMT=6006) SIZE, NX*NY*NZ 6006 FORMAT (//' Buffer size too small, ',I8, $ ' needs to be .gt. ',I12/) CALL CCPERR(1,' ** Array too small **') ELSE C C---- transform on k C C---- read in data store as l h k C C ******************** CALL SFIPYB(X,NZ/2+1,NX,NY,M) C ******************** IF (M .LT. 0) THEN WRITE(ISYSW,'(A)') ' **** Indices outside limits ****' CALL CCPERR(1,' ** FFT failure, index limits **') ENDIF C NWRITE=6 IPRINT=0 CALL STFTP1(.TRUE.,.FALSE.,SIZE,NZ,NX,NY,IERR,NWRITE,IPRINT) IF (IERR .NE. 0) THEN WRITE(ISYSW,'(A)') ' ** Error from STFTP1 **' CALL CCPERR(1,' ** Error from STFTP1 **') ENDIF CALL FURINV(X,X) CALL SHRNK2(X,NZ,NX,NY) C C---- map calculated for p2 sections C R = IYMIN C *don't* use nx, ny, nz as args to avoid illegal aliasing NNX = NX NNY = NY NNZ= NZ CALL PRNBYZ(X,NNZ,NNX,NNY,R,T3,T4) C C ********************************************* IF (T3.GT.0) CALL MRDHDR(T3,'MAPOUT',TITLE,NSEC,IUVW,MXYZ,NW1, + NU1,NU2,NV1,NV2,CELL,LSPGRP,LMODE, + RHMIN,RHMAX,RHMEAN,RHRMS) C ********************************************* C END IF ELSE CALL CCPERR (1, ' !!!!!! NZ must be even !!!!!!') END IF C C---- Format statements C 6000 FORMAT (/////30X,'Fourier synthesis for space-group P1', $ //' Buffer size =',I12/) CCC 6004 FORMAT (' NZ,NX,P2,R,T3,T4,IYMAX',7I7) C END C C C C ================= SUBROUTINE GETHKL C ================= C C---- get reflexion record from file and generate equivalents C C .. Parameters .. REAL TWOPI PARAMETER (TWOPI=2.0*3.14159265359) INTEGER LSYM PARAMETER (LSYM=192) INTEGER MCOLS PARAMETER (MCOLS=500) C .. C .. Scalars in Common .. REAL A,B,BIAS,DEL,DEL1,F000,F1LIM,F2LIM,RHOMAX,RHOMIN,SCLF1,SCLF2, + SCMOD1,SCMOD2,SIG1,SIG2,SMAX,SMIN,SSQ,TFF1,TFF2,F1MAX,F2MAX, + V,CX,CY,CZ INTEGER ICEN,IFreeRexcludeVal,IH,IK,IL,IPHTR,IPROJ,ISORT,ISQ,IW, + KSYM,LIST,M,NABCH,NANOM,NCOL,NF1,NF2,NH,NI1,NLAUE,NLINES, + NLPRGI,NNA,NNB,NP,NP2,NPH,NS1,NS2,NSECAX,NSG,NSSM,NSYM, + NSYMCH,NSYMP,NUM_FREE_COLUMN,NW,NW2,TYPE,T2,T3,T4 LOGICAL USEFILL,USEFREE C .. C .. Arrays in Common .. REAL CS,PERM,RSYM,RSYMT,SN INTEGER LOOKUP,T C .. C .. Local Scalars .. REAL A1,A2,AA,B1,B2,BB,F,F1,F2,FLIST,FMAX,KF1,KF2,FSQ,PHI1,PHI2, + PHIH,PHIPH,PHIT,S,SC,VV,W1,W2,SUMFS INTEGER IH1,IHX,IK1,IKX,IL1,ILX,ISN,IT,JH,JK,JL,KF,NCOLS,NUMFS, + IFIL,ICHK,I LOGICAL EOF, LNAN C .. C .. Local Arrays .. REAL ADATA(MCOLS) INTEGER IN(3),KN(3) LOGICAL LOGMSS(MCOLS) C .. C .. External Subroutines .. EXTERNAL CCPERR,LRNCOL,LRREFL,LRREWD C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,COS,EXP,MAX,MOD,NINT,SIN,SQRT C .. C .. Common blocks .. COMMON /ATSYM/RSYM(4,4,LSYM),RSYMT(4,4,LSYM),NSYM,PERM(4,4),NSYMP, + NLAUE COMMON /HKLF/A,B,ISORT,IH,IK,IL,M,SSQ COMMON /IOFFT/T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /MAPFMT/TYPE,NLINES,NCOL,V,F000,RHOMAX,RHOMIN COMMON /SEL/NF1,NS1,NF2,NS2,NP,NW,IW,ISQ,NH,SMIN,SMAX,SCLF1,TFF1, + SCLF2,TFF2,BIAS,NSSM,NSG,IPROJ,LIST,ICEN,KSYM,NI1,CS(24), + SN(24),T(LSYM,3),NNA,NNB,NPH,NANOM,NW2,NP2,SIG1,SIG2,DEL, + DEL1,F1LIM,F2LIM,IPHTR,SCMOD1,SCMOD2,F1MAX,F2MAX,CX,CY,CZ COMMON /FREESTUFF/ USEFREE, NUM_FREE_COLUMN,IFreeRexcludeVal, + USEFILL COMMON /LOOKCOM/ NLPRGI, LOOKUP (MCOLS) C .. Save statement .. SAVE C .. DATA FMAX,SUMFS,NUMFS/0.0,0.0,0/ C .. ISYSW = 6 IFIL = 1 C C ****************** CALL LRNCOL(IFIL,NCOLS) C ****************** C M = 0 10 CONTINUE C C---- gen equiv or read new reflexion ? C IF (KSYM.GE.NSYMP .OR. NSYM.LE.0) THEN 20 CONTINUE C C---- Remember nf1 nna etc etc are now equal to lookup(jjj) C C *************************** CALL LRREFL(IFIL,SSQ,ADATA,EOF) C *************************** C IF (EOF) THEN C---- End of file M = -1 C C---- modify v for phased translation function to apply correct???? C scale C IF (IPHTR.NE.0) THEN TYPE = -1 F000 = 0 V = SQRT(SCMOD1*SCMOD2) END IF C IF (NUMFS.EQ.0) + CALL CCPERR (1, 'No reflexions pass acceptance criteria!' // + ' Check RESOLUTION, EXCLUDE, missing data.') C WRITE (ISYSW,FMT=6000) IHX,IKX,ILX,FMAX,SUMFS/NUMFS 6000 FORMAT (//' **** Largest F Value *** Reflection ',3I3, F15.2 + // 39X, ' Mean F ',F15.2/) RETURN END IF C S = SSQ M = M + 1 C C Check for missing data. CALL LRREFM(IFIL,LOGMSS) ICHK = 0 DO 21 I=1,NLPRGI IF (USEFILL .AND. + (LOOKUP(I).EQ.NF1 .OR.LOOKUP(I).EQ.NS1)) GO to 21 IF (LOOKUP (I) .GT. 0) THEN IF (LOGMSS(LOOKUP(I))) ICHK = ICHK + 1 END IF 21 CONTINUE IF (ICHK .GT. 0) GOTO 20 IN(1) = NINT(ADATA(1)) IN(2) = NINT(ADATA(2)) IN(3) = NINT(ADATA(3)) C C---- check resolution limits C IF (S.LT.SMIN .OR. S.GT.SMAX) GO TO 20 C C---- to calculate map with data left out for free R-factor C NB NUM_FREE_COLUMN=0 if .NOT.USEFREE C IF (USEFREE) THEN IF (NINT (ADATA(NUM_FREE_COLUMN)) .EQ. IFreeRexcludeVal) + GO TO 20 ENDIF C C---- Check Reflection cut offs C IF (NF1 .GT. 0) THEN IF (.NOT.LOGMSS(NF1)) THEN IF (F1LIM.GT.0.0) THEN IF(ABS(ADATA(NF1)).LT.F1LIM) GO TO 20 ENDIF IF (F1MAX.GT.0.0) THEN IF(ABS(ADATA(NF1)).GT.F1MAX) GO TO 20 ENDIF ENDIF ENDIF IF (NI1 .GT. 0) THEN IF (.NOT.LOGMSS(NI1)) THEN IF (ABS(F1LIM).GT.0.0) THEN IF(ADATA(NI1).LT.F1LIM) GO TO 20 ENDIF IF (F1MAX.GT.0.0) THEN IF(ABS(ADATA(NI1)).GT.F1MAX) GO TO 20 ENDIF ENDIF ENDIF IF (NF2 .GT. 0) THEN IF (F2LIM.GT.0.0) THEN IF(ABS(ADATA(NF2)).LT.F2LIM) GO TO 20 ENDIF IF (F2MAX.GT.0.0) THEN IF(ABS(ADATA(NF2)).GT.F2MAX) GO TO 20 ENDIF ENDIF C C---- Check standard deviations C C F1/SIG1 pair IF (NS1.GT.0 .AND. NF1.GT.0) THEN IF (.NOT. LOGMSS(NS1) .AND. .NOT. LOGMSS(NF1)) THEN IF (LNAN) THEN IF (LOGMSS(NS1) .NEQV. LOGMSS(NF1)) . WRITE(ISYSW,'(/,A,/,A,/)') . ' *** Warning: either F1/SIG1=MNF but other is not', . ' This must mean the colun pair is inconsistent ***' LNAN = .FALSE. ENDIF C old convention for missing datum: IF (ABS(ADATA(NS1)).LT.0.000001) GO TO 20 C apply sigma cutoff: IF (ABS(ADATA(NF1)).LT.SIG1*ABS(ADATA(NS1))) GO TO 20 END IF END IF C F2/SIG2 pair C N.B. SIG2 could be assigned to same column as SIG1 and so be a MNF IF (NS2.GT.0 .AND. NF2.GT.0) THEN IF (.NOT. LOGMSS(NS2) .AND. .NOT. LOGMSS(NF2)) THEN C old convention for missing datum: IF (ABS(ADATA(NS2)).LT.0.000001) GO TO 20 C apply sigma cutoff: IF (ABS(ADATA(NF2)).LT.SIG2*ABS(ADATA(NS2))) GO TO 20 END IF END IF C C I1/SIG1 pair IF (NS1.GT.0 .AND. NI1.GT.0) THEN IF (.NOT. LOGMSS(NS1) .AND. .NOT. LOGMSS(NI1)) THEN C old convention for missing datum: IF (ABS(ADATA(NS1)).LT.0.000001) GO TO 20 C apply sigma cutoff: IF (ABS(ADATA(NI1)).LT.SIG1*ABS(ADATA(NS1))) GO TO 20 END IF END IF C C---- Reset reflexion available flag C KF = 0 KSYM = 0 END IF KSYM = KSYM + 1 C C---- Generate equivalent reflexion C C [HKL'](ROW) = [HKL] *[RSYMT] * [PERM]**(-1)= ....*[PERM](TRANSPOSE) C IH1 = IN(1)*RSYMT(1,1,KSYM) + IN(2)*RSYMT(2,1,KSYM) + + IN(3)*RSYMT(3,1,KSYM) IK1 = IN(1)*RSYMT(1,2,KSYM) + IN(2)*RSYMT(2,2,KSYM) + + IN(3)*RSYMT(3,2,KSYM) IL1 = IN(1)*RSYMT(1,3,KSYM) + IN(2)*RSYMT(2,3,KSYM) + + IN(3)*RSYMT(3,3,KSYM) IH = PERM(1,1)*IH1 + PERM(1,2)*IK1 + PERM(1,3)*IL1 IK = PERM(2,1)*IH1 + PERM(2,2)*IK1 + PERM(2,3)*IL1 IL = PERM(3,1)*IH1 + PERM(3,2)*IK1 + PERM(3,3)*IL1 C---- store unpermuted indices for projection test KN(1) = IH1 KN(2) = IK1 KN(3) = IL1 C C---- Save indices for phase C JH = IH JK = IK JL = IL C C---- Check for projection C IF (IPROJ.NE.0) THEN C C---- NB KN equal to IH1,IK1,IL1 (unpermuted indices) C IF (KN(IPROJ).NE.0) GO TO 10 END IF GO TO (260,260,40,40,40,80,80,80,80, + 80,80,80,80,80,100,100,120,120, + 120,160,160,180,180) NSG C C 1 2 3 4 10 16 17 1018 18 19 20 C 21 23 47 92 96 144 145 146 152 154 169 170 C GO TO 250 C C----MONOCLINIC,Check K and L.GE.0, and H.GE.0 for HK0 C 40 ISN = 1 50 CONTINUE IF (IK.GE.0) THEN IF (IL) 70,60,280 60 IF (IH.GE.0) RETURN END IF 70 IF (ISN.EQ.2) THEN GO TO 10 ELSE IH = -IH IK = -IK IL = -IL ISN = 2 GO TO 50 END IF C C----ORTHORHOMBIC.Check H,K,L.ge.0 C 80 ISN = 1 90 CONTINUE IF (IH.GE.0) THEN IF (IK.GE.0) THEN IF (IL.GE.0) RETURN END IF END IF IF (ISN.EQ.2) THEN GO TO 10 ELSE IH = -IH IK = -IK IL = -IL ISN = 2 GO TO 90 END IF C C----TETRAGONALP41212/P43212. CHECK H,K,L.GE.0 H.GE.K C 100 ISN = 1 110 CONTINUE IF (IH.GE.0) THEN IF (IK.GE.0) THEN IF (IL.GE.0) THEN IF (IH-IK.GE.0) RETURN END IF END IF END IF IF (ISN.EQ.2) THEN GO TO 10 ELSE IH = -IH IK = -IK IL = -IL ISN = 2 GO TO 110 END IF C C----THREEFOLD CHECK H.GE.0 K.GT.0 TAKE ALL L AND 0 0 +L C 120 ISN = 1 130 CONTINUE IF (IH.LT.0) THEN GO TO 140 ELSE IF (IH.GT.0) THEN GO TO 150 END IF C IF (IK.LT.0) THEN GO TO 140 ELSE IF (IK.GT.0) THEN GO TO 280 END IF C IF (IL.GT.0) RETURN 140 IF (ISN.EQ.2) THEN GO TO 10 ELSE IH = -IH IK = -IK IL = -IL ISN = 2 GO TO 130 END IF C 150 IF (IK.LE.0) THEN GO TO 10 ELSE RETURN END IF C C---- P3221/P3121 CHECK H,K .GE. 0 H.GE.K C 160 ISN = 1 170 CONTINUE IF (IH.GE.0) THEN IF (IK.LT.0) THEN GO TO 10 ELSE IF (IH-IK.GE.0) THEN RETURN ELSE END IF END IF IF (ISN.EQ.2) THEN GO TO 10 ELSE IH = -IH IK = -IK IL = -IL ISN = 2 GO TO 170 END IF C C---- six fold check h.ge.0 k.gt.0 take l.GE.0 and 0 0 +l C 180 ISN = 1 190 CONTINUE C IF (IH.LT.0) THEN GO TO 220 ELSE IF (IH.GT.0) THEN GO TO 200 END IF C IF (IK.GE.0) THEN GO TO 210 ELSE GO TO 220 END IF C 200 IF (IK.LE.0) GO TO 10 210 IF (IL.GE.0) RETURN C 220 IF (ISN.EQ.2) THEN GO TO 10 ELSE IH = -IH IK = -IK IL = -IL ISN = 2 GO TO 190 END IF C 250 CONTINUE C C---- TRICLINIC, Check L.GE.0, or H.GE.0 for HK0, or K.GT.0 for 0K0 C 260 ISN = 1 IF (IL.LT.0) THEN GO TO 270 ELSE IF (IL.GT.0) THEN RETURN END IF C IF (IH.LT.0) THEN GO TO 270 ELSE IF (IH.GT.0) THEN RETURN END IF C IF (IK.GE.0) RETURN 270 IK = -IK IL = -IL IH = -IH ISN = 2 280 RETURN C C---- Entry point to calculate structure factors after FFT C has checked hkl. Restructured more logically while adding C variables NW2 NP2 13-2-88 EJD C C---- I want to be able to take vector differences between C two structure factors. C ENTRY GGETF() C ========== C C---- check reflexion available flag C M = 0 IF (KF.EQ.0) THEN IF (NF1.GT.0) THEN C C---- Set weights first C W1 = 1.0 IF (NW.GT.0) W1 = ADATA(NW) W2 = W1 IF (NW2.GT.0) W2 = ADATA(NW2) C C---- Set phases now - zero for Pattersons.. C PHI1 = 0.0 PHI2 = 0.0 C EJD comments about removing the following test: C Agreed - you cannnot have non zero patterson phases. C However in GGETF there are quite complicated vector C differences worked out, and it is important NOT to throw away C the phase information till you have completed the vector C summation, ans only then square the term, which involves C setting A = Magnitude **2, and B=0 C IF (ISQ.EQ.0) THEN C DL: Should we be allowing zero phase distinct from missing? IF (NP.GT.0) PHI1 = 0.0174532788*ADATA(NP) C C---- Anomalous difference Fourier, retard phase by 90 degrees C IF (NANOM.GT.0) PHI1 = PHI1 - 1.5707963 PHI2 = PHI1 IF (NP2.GT.0) PHI2 = 0.0174532788*ADATA(NP2) C END IF C C---- Set F1 now - and perhaps variance C SC = W1*SCLF1*EXP(-TFF1*S/4.0) F1 = 0.0 IF (.NOT.LOGMSS(NF1)) F1 = ADATA(NF1)*SC A1 = COS(PHI1)*F1 B1 = SIN(PHI1)*F1 C C---- Variance of first F C VV = 0.0 IF (ABS(BIAS).GT.0.0001 .AND. NS1.GT.0) THEN IF (.NOT.LOGMSS(NS1)) VV = (ADATA(NS1)*SC)**2 ENDIF C C---- Set F2 now if NF2 gt 0 - and perhaps variance C F2 = 0 A2 = 0 B2 = 0 IF (NF2.GT.0) THEN SC = W2*SCLF2*EXP(-TFF2*S/4.0) F2 = ADATA(NF2)*SC C---- CHECK del between f1 and f2 C C---- Check magnitude of DIFFERENCE C IF (DEL .GT. 0.0) THEN IF (.NOT.LOGMSS(NF1) .AND. ABS(F1-F2).GT.DEL) THEN WRITE (ISYSW,FMT=6004) IH,IK,IL,S,F1,F2,DEL 6004 FORMAT (' Difference too large - h k l f1 f2 ', + '"f" cutoff',/10X,3I4,4X,F7.3,3F12.4) F2 = F1 END IF ENDIF A2 = COS(PHI2)*F2 B2 = SIN(PHI2)*F2 C C---- Variance of second F C IF (BIAS.NE.0.0 .AND. NS2.NE.0) VV = VV + (ADATA(NS2)*SC)**2 C C---- Double difference map for heavy atoms?? C IF (NH.GT.0) THEN PHIH = ADATA(NPH)*0.0174532788 C---- Calculate Fph(calc) = Fp + Fh A2 = ADATA(NH)*SC*COS(PHIH) + A2 B2 = ADATA(NH)*SC*SIN(PHIH) + B2 C---- phi(ph) from Fph(calc) IF (ABS(A2).GT.0.00001 .OR. ABS(B2).GT.0.00001) THEN PHIPH = ATAN2(B2,A2) ELSE C---- default phi(ph) = phi(p) if undetermined PHIPH = PHI1 END IF C---- Reset first complex F to Fph(obs) exp(i phi(ph)) A1 = COS(PHIPH)*F1 B1 = SIN(PHIPH)*F1 C---- Then final F = (|Fph(obs)| - | Fph(calc) |) exp (i phi (ph)) END IF END IF C C---- Evaluate final A and B, F and Phase C A = A1 - A2 B = B1 - B2 C C---- If F1=MNF then replace with F2 if FILLIN keyword specified C IF (USEFILL .AND. LOGMSS(NF1)) THEN SC = (W1*SCLF1*EXP(-TFF1*S/4.0)-W2*SCLF2*EXP(-TFF2*S/4.0)) F = SC*ADATA(NF2) A = F*COS(PHI2) B = F*SIN(PHI2) END IF C FSQ = A*A + B*B F = 0.0 IF (FSQ.GT.0.000000001) F = SQRT(FSQ) C Ian Tickle says NO IF (BIAS.NE.0.0) F = F*F*F/ (F*F+BIAS*VV) C C---- Patterson C C IF (ISQ.NE.0) THEN IF(ISQ.GT.0) F = MAX(F*F - BIAS*VV, 0.0) A = F B = 0 END IF C C---- Extra bit for phased translation function C IF (IPHTR.NE.0) THEN SCMOD1 = F1*F1 + SCMOD1 SCMOD2 = F2*F2 + SCMOD2 F = F1*F2 C To change hand PHI2 becomes TWOPI(IH*CX+IK+CY+IL*CZ)) -PHI2 IF(IPHTR.GT.0) THEN PHIT = (PHI1-PHI2) ELSE PHIT = (PHI1 - TWOPI*(IH*CX+IK+CY+IL*CZ) +PHI2) ENDIF A = COS(PHIT)*F B = SIN(PHIT)*F END IF C C---- Check magnitude of terms C IF (DEL1.GT.0.0 .AND. ABS(F).GT.DEL1) THEN C IF (NF2.EQ.0) F2 = 0.0 6002 FORMAT (' Term too large - h k l f1 f2 "f" cutoff', + /10X,3I4,4X,F7.3,2F12.1,3X,2F12.1) C WRITE (ISYSW,FMT=6002) IH,IK,IL,S,F1,F2,F,DEL1 F = 0 END IF C C----- C IF (ABS(F).EQ.0.0) THEN KSYM = NSYM M = -1 RETURN END IF C ELSE IF (NI1.GT.0) THEN C C A = ADATA(NI1) F = ADATA(NI1) B = 0 C C---- Check magnitude of terms C IF (DEL1.GT.0.0 .AND. ABS(A).GT.DEL1) THEN C C WRITE (ISYSW,FMT=6002) IH,IK,IL,S,F1,F2,F,DEL1 F = 0 END IF C C----- C IF (ABS(F).EQ.0.0) THEN KSYM = NSYM M = -1 RETURN END IF C C C---- Previously calculated A and B parts C ELSE IF (NNA.GT.0) THEN A = ADATA(NNA) B = 0 IF (NNB.GT.0) B = ADATA(NNB) ELSE C C---- nf1 = 0 and nna = 0 - GARBAGE C CALL CCPERR(1,'***** Essential columns not assigned ****') END IF C C---- Deal with phase changes for different symmetry equivalents C IF (NSYM.LE.0) THEN GO TO 290 ELSE AA = A BB = B KF = 1 END IF END IF C C---- Calc phase shift C IT = T(KSYM,1)*JH + T(KSYM,2)*JK + T(KSYM,3)*JL IF (IT.NE.0) THEN IT = MOD(IT,24) IF (IT.LE.0) IT = IT + 24 C C---- Apply phase shift due to space group translations C A = CS(IT)*AA - SN(IT)*BB B = SN(IT)*AA + CS(IT)*BB ELSE A = AA B = BB END IF IF (ISN.EQ.2) B = -B C C---- Optional output of Fourier coefficients to file ABCOEFFS C 290 FLIST = SQRT(A*A+B*B) IF(ISQ.EQ.-1) FLIST = A NUMFS = NUMFS + 1 SUMFS = SUMFS + FLIST IF (FLIST.GT.FMAX) THEN IHX = IH1 IKX = IK1 ILX = IL1 FMAX = FLIST ENDIF IF (LIST.NE.0) THEN PHIT = ATAN2(B,A)/0.0174532788 IF (W1 .EQ. 0.0) THEN KF1 = 0.0 ELSE KF1 = F1/W1 ENDIF IF (W2 .EQ. 0.0) THEN KF2 = 0.0 ELSE KF2 = F2/W2 ENDIF WRITE (NABCH,FMT=6006) $ IH1,IK1,IL1,KF1,PHI1,KF2,PHI2,W1,W2,A,B,PHIT,FLIST 6006 FORMAT (3I4,2(F11.2,F8.2),2F8.3,2F11.2,F8.2,F11.2) END IF RETURN C ENTRY REWND() C =========== C C---- rewind hkl file for reread C C ****** CALL LRREWD(1) C ****** C END C C C C ======================================== SUBROUTINE HARKER(NSYM,RSYM) C ======================================== parameter (msym = 96) dimension rsym(4,4,*),rsympatt(4,4,msym),hsum(msym) character symchs(msym)*80,sect(3)*1 data sect/'X','Y','Z'/ ISYSW = 6 nsm = 1 call symtr3( nsm, rsym, symchs(1),0) C check harker vectors do 1 isym = 2,nsym do 2 i= 1,4 do 3 j= 1,4 rsympatt(i,j,isym) = rsym(i,j,1)- rsym(i,j,isym) 3 continue hsum(i) = rsympatt(i,1,isym)*rsympatt(i,1,isym) + + rsympatt(i,2,isym)*rsympatt(i,2,isym) + + rsympatt(i,3,isym)*rsympatt(i,3,isym) 2 continue rsympatt(4,4,4) = 1.0 do 4 i= 1,3 if(hsum(i).le.0.1) then C C call symtr3(nsm, rsym(1,1,isym),symchs(isym),0) C---- write a message C WRITE (ISYSW,'(//A,2x,I2,3A,I5,2A)') 'Symmetry operators:', + nsm,': ',symchs(1) (1:MIN(350,LENSTR(SYMCHS(1)))),':', + isym,':',SYMCHS(isym) (1:MIN(350,LENSTR(SYMCHS(isym)))) C C call symtr3( 1, rsympatt(1,1,isym),symchs(isym),0) C---- write a message rsympatt(i,4,isym) = MOD( (rsympatt(i,4,isym)+10.0),1.0) C WRITE (ISYSW,FMT='(A,I3,5X,A)') ' Harker Vector:',isym, + SYMCHS(isym) (1:MIN(350,LENSTR(SYMCHS(isym)))) C WRITE (ISYSW,FMT='(3A,F5.2)') + ' Harker Section is: ',SECT(I),' = ',rsympatt(i,4,isym) C WRITE (ISYSW,FMT='(/,A,4(/,8x,4F6.2))') + ' Harker vector Matrix:',((RSYMPATT(JDO50,J,isym),J=1,4), + JDO50 = 1,4) end if 4 continue 1 continue 99 continue end C C C VAXPDS NAME=GRIDIM C********************************************************************** SUBROUTINE GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) C---------------------------------------------------------------------- C CALCULATE ARRAY DIMENSIONING CONSTANTS C A.D. MCLACHLAN 1981. LAST UPDATED 26 AUG 1987. C---------------------------------------------------------------------- INTEGER IDIM,NDIM,NGRID,IUVW,LSEP,IU,IV,IW DIMENSION IDIM(10),NDIM(3),NGRID(3),IUVW(3),LSEP(3) EQUIVALENCE(IU,IUVW(1)) EQUIVALENCE(IV,IUVW(2)) EQUIVALENCE(IW,IUVW(3)) LOGICAL ERROR INTEGER N,LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW,NWRITE,I,ICMPLX, + LGRIDU,LGRIDV,LGRIDW C----------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C OUTPUT -- ERROR .L. .TRUE. IF ERROR IS DETECTED C OUTPUT -- LSEPU SEPARATION OF ELEMENTS ALONG IU AXIS C OUTPUT -- LSEPV SEPARATION OF ELEMENTS ALONG IV AXIS C OUTPUT -- LSEPW SEPARATION OF ELEMENTS ALONG IW AXIS C OUTPUT -- LGLIMV UPPER LIMIT FOR LINDV LOOP C OUTPUT -- LGLIMW UPPER LIMIT FOR LINDW LOOP C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- *** C----------------------------------------------------------------------- C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C------------------------------------------------------------------------------- 20 FORMAT(1X,' GRIDIM ERROR *** NGRID ',3I5,' GT NDIM ',3I5) 21 FORMAT(1X,' GRIDIM ERROR *** INVALID AXES ',3I5) 22 FORMAT(1X,' GRIDIM ERROR *** N.GT.LGRIDU ',2I5) C-------------------------------- ERROR=.FALSE. DO 100 I=1,3 NDIM(I)=IDIM(I+3) NGRID(I)=IDIM(I+6) IUVW(I)=IDIM(I) 100 CONTINUE ICMPLX=IDIM(10) DO 105 I=1,3 IF(NGRID(I).GT.NDIM(I)) GO TO 300 105 CONTINUE DO 110 I=1,3 IF((IUVW(I).LT.1).OR.(IUVW(I).GT.3)) GO TO 400 110 CONTINUE IF((IU.EQ.IV).OR.(IV.EQ.IW).OR.(IW.EQ.IU)) GO TO 400 LSEP(1)=ICMPLX LSEP(2)=ICMPLX*NDIM(1) LSEP(3)=ICMPLX*NDIM(1)*NDIM(2) LGRIDU=NGRID(IU) LGRIDV=NGRID(IV) LGRIDW=NGRID(IW) IF(LGRIDU.LT.N) GO TO 500 LSEPU=LSEP(IU) LSEPV=LSEP(IV) LSEPW=LSEP(IW) LGLIMV=LSEPV*LGRIDV LGLIMW=LSEPW*LGRIDW RETURN 300 CONTINUE C GRID SIZE ERROR WRITE(NWRITE,20) NGRID,NDIM ERROR=.TRUE. GO TO 600 400 CONTINUE C AXES ERROR WRITE(NWRITE,21) IUVW ERROR=.TRUE. GO TO 600 500 CONTINUE C FOURIER RANGE TOO BIG FOR ARRAY DIMENSION WRITE(NWRITE,22) N,LGRIDU ERROR= .TRUE. 600 CONTINUE RETURN END C C VAXPDS NAME=HERMFT C***************************************************************************** SUBROUTINE HERMFT (X, Y, N, IDIM,NWRITE) C----------------------------------------------------------------------------- C HERMITIAN SYMMETRIC FOURIER TRANSFORM C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C----------------------------------------------------------------------------- REAL X(*), Y(*) INTEGER N INTEGER IDIM(10) C----- INTEGER DIM(5) LOGICAL ERROR REAL A, B, C, D, E, F, CO, SI INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW EQUIVALENCE (LSEPU,DIM(1)),(LSEPV,DIM(2)),(LSEPW,DIM(3)) EQUIVALENCE (LGLIMV,DIM(4)),(LGLIMW,DIM(5)) INTEGER I, J, K, N OVER 2, I0, LINDU,LINDV,LDXV,LDXW INTEGER K1,NWRITE,L DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBTWON DATA DBPI2/6.28318530717958623D0/ C-- C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------- C IN/OUT -- X(*) REAL PART OF DATA/TRANSFORM C IN/OUT -- Y(*) IMAGINARY PART OF DATA/TRANSFORM C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- GRIDIM CALCULATE DIMENSIONING CONSTANTS C CALLS -- CMPLFT COMPLEX FOURIER TRANSFORM C----------------------------------------------------------------------- C NOTE -- INDEXING REVISED AND SIN COS RECALCULATED A.D.MCLACHLAN AUGUST 1981 C----------------------------------------------------------------------------- C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C------------------------------------------------------------------------------- C NOTE -- GIVEN THE UNIQUE TERMS OF A HERMITIAN SYMMETRIC SEQUENCE OF LENGTH C NOTE -- 2N THIS SUBROUTINE CALCULATES THE 2N REAL NUMBERS WHICH ARE ITS C NOTE -- FOURIER TRANSFORM. THE EVEN NUMBERED ELEMENTS OF THE TRANSFORM C NOTE -- (0, 2, 4, . . ., 2N-2) ARE RETURNED IN X AND THE ODD NUMBERED C NOTE -- ELEMENTS (1, 3, 5, . . ., 2N-1) IN Y. C NOTE -- C NOTE -- A FINITE HERMITIAN SEQUENCE OF LENGTH 2N CONTAINS N + 1 UNIQUE C NOTE -- REAL NUMBERS AND N - 1 UNIQUE IMAGINARY NUMBERS. FOR CONVENIENCE C NOTE -- THE REAL VALUE FOR X(N) IS NOW STORED AT X(N+1) INSTEAD OF Y(0). C NOTE -- THE ROUTINE RETURNS WITH X(N+1) SET TO 0.0 C------------------------------------------------------------------------------ DBTWON=DBLE(2*N) ERROR=.FALSE. CALL GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) IF(ERROR) GO TO 600 J=N*LSEPU LINDU=1-LSEPV-LSEPW C-- DO 100 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 100 LDXW=LSEPW,LGLIMW,LSEPW I=LINDV+LDXW L = I + J A = X(I) B = X(L) X(L)=0.0 X(I) = A + B Y(I) = A - B 100 CONTINUE C-- N OVER 2 = N/2 + 1 IF (N OVER 2 .LT. 2) GO TO 500 DBANGL=DBPI2/DBTWON DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 400 I0 = 2, N OVER 2 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 CO=DBC1 SI=DBS1 K = (N + 2 - 2*I0)*LSEPU K1 = (I0 - 1)*LSEPU + 1 LINDU=K1-LSEPV-LSEPW DO 300 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 300 LDXW=LSEPW,LGLIMW,LSEPW I=LINDV+LDXW J = I + K A = X(I) + X(J) B = X(I) - X(J) C = Y(I) + Y(J) D = Y(I) - Y(J) E = B*CO + C*SI F = B*SI - C*CO X(I) = A + F X(J) = A - F Y(I) = E + D Y(J) = E - D 300 CONTINUE 400 CONTINUE C----- CALL CMPLFT(X,Y,N,IDIM,NWRITE) C----- 500 CONTINUE RETURN 600 CONTINUE WRITE(NWRITE,700) N,IDIM CALL CCPERR (1, 'FFTBIG: internal error') 700 FORMAT(1X,'**HERMFT** ERROR IN DIMENSIONS OF ARRAYS N=,IDIM=', 1 I5,5X,3I5,1X,3I5,1X,3I5,1X,I5) C---------- END C C VAXPDS NAME=INV21 C************************************************************************ SUBROUTINE INV21(X,Y,N,IDIM,NWRITE) C------------------------------------------------------------------------ C INVERTS FOURIER TRANSFORM ALONG A SCREW DIAD. THE RESULT IS SCALED BY N C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C------------------------------------------------------------------------ REAL X(*),Y(*) INTEGER N,IDIM(10) C----- INTEGER DIM(5) LOGICAL ERROR INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW EQUIVALENCE(LSEPU,DIM(1)),(LSEPV,DIM(2)),(LSEPW,DIM(3)) EQUIVALENCE(LGLIMV,DIM(4)),(LGLIMW,DIM(5)) INTEGER LINDV,LDXV,LDXW,NWRITE,NOVER2,LL,KK,LINDU, + J,K,L,I,J1 REAL A,B,C,S DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBTWON DATA DBPI2/6.28318530717958623D0/ C C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------- C IN/OUT -- X(*) REAL PART OF DATA/TRANSFORM C IN/OUT -- Y(*) IMAGINARY PART OF DATA/TRANSFORM C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- GRIDIM CALCULATE DIMENSIONING CONSTANTS C CALLS -- CMPLFT COMPLEX FOURIER TRANSFORM C----------------------------------------------------------------------- C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C------------------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED AND SIN COS RECALC A.D.MCLACHLAN C NOTE -- AUGUST 1981 C NOTE -- CORRECTION BECAUSE RESULTS TOO SMALL BY FACTOR OF 2.0 C------------------------------------------------------------------------- DBTWON=DBLE(2*N) ERROR=.FALSE. CALL GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) IF(ERROR) GO TO 500 C-- NOVER2 = N/2 LL = N*LSEPU KK = NOVER2*LSEPU LINDU=1-LSEPV-LSEPW DO 100 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 100 LDXW=LSEPW,LGLIMW,LSEPW J=LINDV+LDXW L = LL+J K = KK+J X(L) = X(J)+X(K) X(K) = X(K)+Y(K) Y(L) = 0.0 100 Y(K) = 0.0 C-- DBANGL=DBPI2/DBTWON DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 200 I = 2, NOVER2 KK = (N+2-2*I)*LSEPU LL = (N+1-I)*LSEPU DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C=DBC1 S=DBS1 J1 = (I-1)*LSEPU+1 LINDU=J1-LSEPV-LSEPW DO 150 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 150 LDXW=LSEPW,LGLIMW,LSEPW J=LINDV+LDXW L = J+LL K = J+KK X(L) = X(L)+X(J)+X(K) X(J) = X(J)+S*Y(J) X(K) = X(K)+S*Y(K) Y(J) = C*Y(J) 150 Y(K) = -C*Y(K) 200 CONTINUE C-- MULTIPLY X AND Y BY 2.0 N1=N+1 DO 220 I=1,N1 J1=(I-1)*LSEPU + 1 LINDU=J1-LSEPV-LSEPW DO 220 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 220 LDXW=LSEPW,LGLIMW,LSEPW J=LINDV+LDXW X(J)=2.0*X(J) Y(J)=2.0*Y(J) 220 CONTINUE C------------------------------------- CALL CMPLFT(X,Y,N,IDIM,NWRITE) C------------------------------------- DO 300 I = 1,NOVER2 KK = (N+1-2*I)*LSEPU LL = KK+I*LSEPU J1 = (I-1)*LSEPU+1 LINDU=J1-LSEPV-LSEPW DO 250 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 250 LDXW=LSEPW,LGLIMW,LSEPW J=LINDV+LDXW K = J+KK L = J+LL A = X(J)-X(L) B = Y(J)+Y(L) X(J) = X(L) Y(J) = -Y(L) X(L) = X(K)+A Y(L) = Y(K)-B X(K) = A 250 Y(K) = B 300 CONTINUE C- M = N-2 DO 400 I = 1,M K = I 320 J = K K = J/2 IF(2*K.NE.J) K = N-1-K IF(K-I) 320,400,340 340 KK = (K-I)*LSEPU J1 = I*LSEPU+1 LINDU=J1-LSEPV-LSEPW DO 360 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 360 LDXW=LSEPW,LGLIMW,LSEPW J=LINDV+LDXW K = J+KK A = X(K) B = Y(K) X(K) =X(J) Y(K) = Y(J) X(J) = A 360 Y(J) = B 400 CONTINUE C----- RETURN C----- 500 CONTINUE WRITE(NWRITE,600) N,IDIM CALL CCPERR (1, 'FFTBIG: internal error') 600 FORMAT(1X,'**INV21** ERROR IN ARRAY DIMENSIONS N=,IDIM=', 1 I5,5X,3I5,1X,3I5,1X,3I5,1X,I5) END C C VAXPDS NAME=MDFTKD C************************************************************************* SUBROUTINE MDFTKD (N,NFACTR,NDIM,X,Y,NWRITE) C------------------------------------------------------------------------- C MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL DRIVER C SMALL REVISIONS DEC 1984. A.D. MCLACHLAN. LAST UPDATED 26 AUG 1987 C------------------------------------------------------------------------- INTEGER N INTEGER NFACTR (*),NDIM(5) REAL X(*), Y(*) C----- INTEGER NF, M, NP, NR, LSEPU C------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- NFACTR(*) LIST OF PRIME FACTORS IN SPECIAL ORDER C INPUT -- NDIM(5) ARRAY DIMENSIONING CONSTANTS C IN/OUT -- X(*) REAL PART OF DATA C IN/OUT -- Y(*) IMAGINARY PART OF DATA C INPUT -- NWRITE PRINT UNIT C------------------------------------------------------------------------- C CALLS -- R2CFTK RADIX 2 KERNEL C CALLS -- R3CFTK RADIX 3 KERNEL C CALLS -- R4CFTK RADIX 4 KERNEL C CALLS -- R5CFTK RADIX 5 KERNEL C CALLS -- R8CFTK RADIX 8 KERNEL C CALLS -- RPCFTK PRIME RADIX KERNEL C------------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED A.D. MCLACHLAN AUGUST 1981 C------------------------------------------------------------------------- 20 FORMAT (1X,'***MDFTKD*** ERROR- ILLEGAL FACTOR DETECTED =',I8) C------------------------------------------------------------------------- LSEPU=NDIM(1) NF = 0 M = N C --GO THROUGH ALL FACTORS, SKIPPING NP=1 100 CONTINUE NF = NF + 1 NP = NFACTR(NF) IF (NP.EQ.0) RETURN M = M/NP NR = M*LSEPU IF (NP.GT.8) GO TO 900 GO TO (100, 200, 300, 400, 500, 600, 700, 800), NP GO TO 1000 C --FACTOR OF 2 200 CONTINUE CALL R2CFTK(N, M, X(1), Y(1), X(NR+1), Y(NR+1), NDIM) GO TO 100 C --FACTOR OF 3 300 CONTINUE CALL R3CFTK(N, M, X(1), Y(1), X(NR+1), Y(NR+1), 1 X(2*NR+1), Y(2*NR+1), NDIM) GO TO 100 C --FACTOR OF 4 400 CONTINUE CALL R4CFTK(N, M, X(1), Y(1), X(NR+1), Y(NR+1), 1 X(2*NR+1), Y(2*NR+1), X(3*NR+1), Y(3*NR+1), NDIM) GO TO 100 C --FACTOR OF 5 500 CONTINUE CALL R5CFTK(N, M, X(1), Y(1), X(NR+1), Y(NR+1), 1 X(2*NR+1), Y(2*NR+1), X(3*NR+1), Y(3*NR+1), 2 X(4*NR+1), Y(4*NR+1), NDIM) GO TO 100 C --FACTOR 6 IS ILLEGAL 600 CONTINUE GO TO 1000 C --7 TREATED AS A GENERAL PRIME 700 CONTINUE GO TO 900 C --FACTOR OF 8 800 CONTINUE CALL R8CFTK (N, M, X(1), Y(1), X(NR+1), Y(NR+1), 1 X(2*NR+1), Y(2*NR+1), X(3*NR+1), Y(3*NR+1), 2 X(4*NR+1), Y(4*NR+1), X(5*NR+1), Y(5*NR+1), 3 X(6*NR+1), Y(6*NR+1), X(7*NR+1), Y(7*NR+1), 4 NDIM) GO TO 100 C --GENERAL PRIME FACTOR 900 CONTINUE CALL RPCFTK(N, M, NP, NR, X, Y, NDIM) GO TO 100 C --ERROR STOP 1000 CONTINUE WRITE(NWRITE,20) NP CALL CCPERR (1, 'FFTBIG: internal error') END C C C C =============================== SUBROUTINE PRMVEC(PERM,JV,N,N1) C =============================== C C Permute vector JV(N,3) by permutation matrix PERM C N1 is first dimension of JV C C---- Permute C C .. Scalar Arguments .. INTEGER N,N1 C .. C .. Array Arguments .. REAL PERM(4,4) INTEGER JV(N1,3) C .. C .. Local Scalars .. INTEGER I C .. C .. Local Arrays .. INTEGER LV(3) C .. DO 10 I = 1,3 LV(I) = PERM(I,1)*JV(N,1) + PERM(I,2)*JV(N,2) + + PERM(I,3)*JV(N,3) 10 CONTINUE C C---- Copy back C DO 20 I = 1,3 JV(N,I) = LV(I) 20 CONTINUE END C C C C =============================== SUBROUTINE PRMVET(PERM,JV,N,N1) C =============================== C C Permute vector JV(N,3) by transpose of permutation matrix PERM C N1 is first dimension of JV C C---- Permute C C .. Scalar Arguments .. INTEGER N,N1 C .. C .. Array Arguments .. REAL PERM(4,4) INTEGER JV(N1,3) C .. C .. Local Scalars .. INTEGER I C .. C .. Local Arrays .. INTEGER LV(3) C .. DO 10 I = 1,3 LV(I) = PERM(1,I)*JV(N,1) + PERM(2,I)*JV(N,2) + + PERM(3,I)*JV(N,3) 10 CONTINUE C C---- Copy back C DO 20 I = 1,3 JV(N,I) = LV(I) 20 CONTINUE END C C C C =========================================== SUBROUTINE PRNBYZ(X,NF,NM,NS,ISEC0,BIN,BCD) C =========================================== C C FFT map output subroutine for section defined by 'kp3' C This version writes the binary file with 'kp2' running fastest C C*C C *** This is a special version of PRNTYZ for FFTBIG *** C C Differences are flagged with C*C, and merely concern the fact that C here we enter with all sections from 0 always, so the definition C of KSEC1 & KSEC2 is different, also the calculation of SEC. C Note that ISEC0 still enters as minimum section number for output C*C C .. Parameters .. INTEGER LSYM PARAMETER (LSYM=192) C .. C .. Scalar Arguments .. INTEGER BCD,BIN,ISEC0,NF,NM,NS C .. C .. Array Arguments .. REAL X(NF,NM,NS) C .. C .. Scalars in Common .. REAL F000,RHMAX,RHMEAN,RHMIN,RHOMAX,RHOMIN,T2,T3,T4,V INTEGER IDIG,ILIN,ISP,KP1,KP2,KP3,LCMAX,LCMIN,LSPGRP,NABCH,NCOL, + NLINES,NSECAX,NSYM,NSYMCH,NUSED10,NUSED11,NUSED12,NUSED6, + NUSED7,NUSED9,TYPE,XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN CHARACTER TITLE*80 C .. C .. Arrays in Common .. REAL CELL,PERM,ROT,ROTT INTEGER IDUMMY,NXYZ C .. C .. Local Scalars .. REAL R,RMAX,RMAXOV,RMIN,RMINOV,RMS INTEGER I,J,K,KX1,KX1MAX,KX1MIN,KX1MNV,KX1MXV,KX2,KX2MAX,KX2MIN, + KX2MNV,KX2MXV,KX3,KX3MNV,KX3MXV,L,MODE,N,NFF,NMM,NPT,NSEC, + NSS,SEC,XL,XLLP,XU,XULP,YZL,YZLLP,YZU,YZULP,KSEC1,KSEC2 LOGICAL SPEW CHARACTER AST*1,SPC*1 C .. C .. Local Arrays .. INTEGER IUVW(3),JUNK(3),LINE(120),XYZMN(3),XYZMX(3) CHARACTER KXYZ(3)*1,IFMT(8)*4,JDIG(9)*4,JLIN(4)*4 C .. C .. External Subroutines .. EXTERNAL MSPEW,MSYPUT,MWCLOSE,MWRHDR,MWRSEC,PRMVEC,PRMVET C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN,NINT,REAL,SQRT C .. C .. Common blocks .. COMMON /ATSYM/ROT(4,4,LSYM),ROTT(4,4,LSYM),NSYM,PERM(4,4),NSYMP, + NLAUE COMMON /FFTPRM/KP1,KP2,KP3 COMMON /FOUT/LCMIN,LCMAX,ISP,IDIG,ILIN COMMON /IOFFT/T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /MAPFMT/TYPE,NLINES,NCOL,V,F000,RHOMAX,RHOMIN COMMON /MPHDR/NXYZ(3),XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,NUSED6,NUSED7, + NUSED9,NUSED10,NUSED11,IDUMMY(3),NUSED12,CELL(6),LSPGRP, + RHMIN,RHMAX,RHMEAN COMMON /MPHDRR/TITLE C .. C .. Save statement .. SAVE C .. C .. Data statements .. DATA KXYZ/'X','Y','Z'/ DATA AST/'*'/,SPC/' '/ DATA IFMT/'(',' ','1X,I','3,1H','*,','128I','4',')'/ DATA JDIG/'1','2','3','4','5','6','7','8','9'/ DATA JLIN/' ','/ ','// ','/// '/ C .. C C---- Take care: nf nm ns are " nz nx p2" for FFT1 ,etc C AND nxyz(3) are the permuted values of NX NY NZ C SO nf and nm can be reset by PRMVET and this can mess up the C definition of SPEW C ISYSW = 6 NFF = NF NMM = NM NSS = NS C C---- Now permute "nx" "ny" "nz" to get nxyz(.) related to A B C C---- PRMVET will return NXYZ AND XMIN. to original order C C ********************* CALL PRMVET(PERM,NXYZ,1,1) C ********************* C C---- Patch included to conform to F77 Standards C It is not permitted to pass to a subroutine a C single variable (where an array is expected) and C rely on the program selecting the appropriate number C of variables from a labelled COMMON BLOCK C JUNK(1) = XMIN JUNK(2) = YMIN JUNK(3) = ZMIN C C ********************* CALL PRMVET(PERM,JUNK,1,1) C ********************* C XMIN = JUNK(1) YMIN = JUNK(2) ZMIN = JUNK(3) C JUNK(1) = XMAX JUNK(2) = YMAX JUNK(3) = ZMAX C C ********************* CALL PRMVET(PERM,JUNK,1,1) C ********************* C XMAX = JUNK(1) YMAX = JUNK(2) ZMAX = JUNK(3) C C---- if this is the first time, derive the scale factors C XYZMN(1) = XMIN XYZMN(2) = YMIN XYZMN(3) = ZMIN C XYZMX(1) = XMAX XYZMX(2) = YMAX XYZMX(3) = ZMAX C C---- use nsecax to choose whether this is a "printy" look-alike C or a "printz" one. C C---- FOR "PRINTZ" PUT kx1=kp1 kx2=kp2 kx3=kp3 C KX1 = KP1 KX2 = KP2 KX3 = KP3 C C---- FOR "PRINTY" PUT kx1=kp1 kx2=kp3 kx3=kp2 C IF (NSECAX.EQ.2) THEN KX2 = KP3 KX3 = KP2 END IF C C*C C Section limits KSEC1 = XYZMN(KX3)+1 KSEC2 = XYZMX(KX3)+1 C*C IF (ISEC0.EQ.XYZMN(KX3)) THEN IF (NLINES.LE.0) NLINES = XYZMX(KX1) - XYZMN(KX1) + 1 I = 128/IDIG IF (NCOL.LE.0 .OR. NCOL.GT.I) NCOL = I C IF (BIN.GT.0) THEN C C---- Write out binary map file header C Title, NSEC (number of sections in map), C fast, medium and sMIN axis C numbers, NNX,NNY,NNZ (sampling intervals along real xyz) C NSEC = XYZMX(KX3) - XYZMN(KX3) + 1 C C---- Set up map header in common C Axis order C IUVW(1) = KX2 IUVW(2) = KX1 IUVW(3) = KX3 MODE = 2 C C *************************************************** CALL MWRHDR(BIN,TITLE,NSEC,IUVW,NXYZ,XYZMN(KX3),XYZMN(KX2), + XYZMX(KX2),XYZMN(KX1),XYZMX(KX1),CELL,LSPGRP,MODE) C *************************************************** C C---- Copy symmetry operations from library file stream 4 to output map C C ************************* CALL MSYPUT(NSYMCH,LSPGRP,BIN) C ************************* C END IF C C---- XL XU refer to column numbers C YZL YZU refer to row numbers C XL = XYZMN(KX1) + 1 XU = XYZMX(KX1) + 1 YZL = XYZMN(KX2) + 1 YZU = XYZMX(KX2) + 1 RMINOV = 10000.0 RMAXOV = -10000.0 SPEW = (XL.EQ.1) .AND. (XU.EQ.NMM) .AND. (YZL.EQ.1) .AND. + (YZU.EQ.NFF) .AND. (BIN.NE.0) IF (V.EQ.0.0) V = 1.0 C IF (RHOMIN.GT.RHOMAX) THEN R = RHOMIN RHOMIN = RHOMAX RHOMAX = R END IF C IF (TYPE.GE.0) THEN RMAX = X(YZL,XL,1) RMIN = RMAX C DO 30 I = KSEC1, KSEC2 DO 20 J = XL,XU DO 10 K = YZL,YZU R = X(K,J,I) IF (R.GT.RMAX) RMAX = R IF (R.LT.RMIN) RMIN = R 10 CONTINUE 20 CONTINUE 30 CONTINUE C IF (TYPE.GT.0) THEN V = (RMAX-RMIN)/ (RHOMAX-RHOMIN) F000 = V*RHOMAX - RMAX ELSE RMAX = MAX(ABS(RMIN),RMAX) V = RMAX/RHOMAX F000 = 0.0 END IF END IF C WRITE (ISYSW,FMT=6000) F000,V 6000 FORMAT (' F000 = ',1P,E12.3,10X,'V = ',E12.3) IF (TYPE.GE.0) WRITE (ISYSW,FMT=6002) NSS 6002 FORMAT (' Scaling Based on Values in the First ',I3, + ' Sections.') WRITE (ISYSW,FMT=6004) KXYZ(KX3) 6004 FORMAT (/' Section axis is ',A1) IF (BIN.NE.0) WRITE (ISYSW,FMT=6006) KXYZ(KX2),KXYZ(KX1) 6006 FORMAT (/' Binary output has ',A1,' running fastest, ',A1,' ne', + 'xt') IF (BCD.NE.0) WRITE (ISYSW,FMT=6008) KXYZ(KX2),KXYZ(KX1) 6008 FORMAT (/' Line-printer output has ',A1,' across, ',A1,' down ', + 'the page') C C---- Initialise max,min,mean C RHMEAN = 0.0 RHMIN = (X(YZL,XL,1)+F000)/V RHMAX = RHMIN RMS = 0.0 NPT = 0 END IF C C---- put out scaled "yz" sections C DO 110 K = KSEC1,KSEC2 C*C SEC = K - 1 C*C C---- scale this section C XL = XYZMN(KX1) + 1 XU = XYZMX(KX1) + 1 YZL = XYZMN(KX2) + 1 YZU = XYZMX(KX2) + 1 RMIN = (X(YZL,XL,K)+F000)/V RMAX = RMIN KX1MIN = XL - 1 KX2MIN = YZL - 1 KX1MAX = XL - 1 KX2MAX = YZL - 1 C DO 50 J = XL,XU DO 40 I = YZL,YZU R = (X(I,J,K)+F000)/V X(I,J,K) = R RHMEAN = RHMEAN + R RMS = R*R + RMS NPT = NPT + 1 C C---- find min and max values for this section : find position they occu C and also overall min and max. C IF (R.LT.RMIN) THEN RMIN = R KX1MIN = J - 1 KX2MIN = I - 1 END IF C IF (R.GT.RMAX) THEN RMAX = R KX1MAX = J - 1 KX2MAX = I - 1 END IF C 40 CONTINUE 50 CONTINUE C C---- Find overall max min C IF (RMAX.GT.RMAXOV) THEN RMAXOV = RMAX KX1MXV = KX1MAX KX2MXV = KX2MAX KX3MXV = SEC END IF C IF (RMIN.LT.RMINOV) THEN RMINOV = RMIN KX1MNV = KX1MIN KX2MNV = KX2MIN KX3MNV = SEC END IF C RHMIN = MIN(RHMIN,RMIN) RHMAX = MAX(RHMAX,RMAX) RHMIN = MIN(RHMIN,RMIN) RHMAX = MAX(RHMAX,RMAX) C C---- binary output C IF (BIN.NE.0) THEN IF (SPEW) THEN C C ******************* CALL MSPEW(BIN,X(1,1,K)) C ******************* C ELSE C C ****************************************** CALL MWRSEC(BIN,X(1,1,K),NFF,NMM,YZL,YZU,XL,XU) C ****************************************** C END IF END IF C C---- decimal output C IF (BCD.NE.0) THEN IF (BCD.EQ.0) THEN GO TO 120 ELSE C C---- construct format C IFMT(2) = JLIN(ILIN) IFMT(7) = JDIG(IDIG) YZULP = XYZMN(KX2) 60 CONTINUE C YZLLP = YZULP + 1 YZULP = YZULP + NCOL IF (YZULP.GT.XYZMX(KX2)) YZULP = XYZMX(KX2) + 1 XULP = XYZMN(KX1) 70 CONTINUE XLLP = XULP + 1 XULP = XULP + NLINES IF (XULP.GT.XYZMX(KX1)) XULP = XYZMX(KX1) + 1 I = YZLLP - 1 J = YZULP - 1 WRITE (BCD,FMT=6010) TITLE,KXYZ(KX3),SEC,KXYZ(KX2),I,J 6010 FORMAT (////10X,A80,5X,A1,' =',I3,', ',A1,' =',I3,' TO',I4, + /) L = 0 C DO 80 J = YZLLP,YZULP L = L + 1 LINE(L) = J - 1 80 CONTINUE C WRITE (BCD,FMT=IFMT) SEC, (LINE(J),J=1,L) N = IDIG - 1 WRITE (BCD,FMT=6012) ((SPC,I=1,N),AST,J=1,L) 6012 FORMAT (5X,128A1) C DO 100 J = XLLP,XULP LINE(1) = J - 1 L = 1 C DO 90 I = YZLLP,YZULP L = L + 1 LINE(L) = NINT(X(I,J,K)) IF (LINE(L).GT.LCMIN .AND. LINE(L).LT.LCMAX) LINE(L) = 0 90 CONTINUE C WRITE (BCD,FMT=IFMT) (LINE(I),I=1,L) 100 CONTINUE C IF (XULP.LE.XYZMX(KX1)) GO TO 70 IF (YZULP.LE.XYZMX(KX2)) GO TO 60 END IF END IF C IF (BCD.GT.0) WRITE (ISYSW,FMT=6014) 6014 FORMAT (1X) WRITE (ISYSW,FMT=6016) KXYZ(KX3),SEC,RMIN,KXYZ(KX1),KX1MIN, + KXYZ(KX2),KX2MIN,RMAX,KXYZ(KX1),KX1MAX,KXYZ(KX2),KX2MAX 6016 FORMAT (' Section ',A1,' = ',I3,' minimum rho =',F15.4,' at ', + A1,' = ',I3,' ',A1,' = ',I3,' maximum rho =',F15.4, + ' at ',A1,' = ',I3,' ',A1,' = ',I3) 110 CONTINUE C IF (ISEC0+NSS.GT.XYZMX(KX3)) THEN C C---- Close map file - MCLOSE prints out density limits C C ***************** IF (BIN.NE.0) CALL MWCLOSE(BIN) C ***************** C C---- Print max min information C WRITE (ISYSW,FMT=6020) RMAXOV,KXYZ(KX3),KX3MXV,KXYZ(KX1),KX1MXV, + KXYZ(KX2),KX2MXV 6020 FORMAT (//' Overall maximum rho = ',F15.4,' on section ',A1, + ' = ',I3,' at ',A1,' = ',I3,1X,A1,' = ',I3) RHMEAN = RHMEAN/REAL(NPT) WRITE (ISYSW,FMT=6022) RMINOV,KXYZ(KX3),KX3MNV,KXYZ(KX1),KX1MNV, + KXYZ(KX2),KX2MNV,RHMEAN 6022 FORMAT (//' Overall minimum rho = ',F15.4,' on section ',A1, + ' = ',I3,' at ',A1,' = ',I3,1X,A1,' = ',I3,//' Average', + ' density on this map is ',F19.8) END IF C 120 CONTINUE C C---- Restore permuted order of nxyz xmin etc C C ********************* CALL PRMVEC(PERM,NXYZ,1,1) C ********************* C C---- Patch included to conform to F77 Standards C It is not permitted to pass to a subroutine a C single variable (where an array is expected) and C rely on the program selecting the appropriate number C of variables from a labelled COMMON BLOCK C JUNK(1) = XMIN JUNK(2) = YMIN JUNK(3) = ZMIN C C ********************* CALL PRMVEC(PERM,JUNK,1,1) C ********************* C XMIN = JUNK(1) YMIN = JUNK(2) ZMIN = JUNK(3) C JUNK(1) = XMAX JUNK(2) = YMAX JUNK(3) = ZMAX C C ********************* CALL PRMVEC(PERM,JUNK,1,1) C ********************* C XMAX = JUNK(1) YMAX = JUNK(2) ZMAX = JUNK(3) C END C VAXPDS NAME=PRNT3D C********************************************************* SUBROUTINE PRNT3D(RHO,NU,NV,NW,HEAD,NWRITE) C--------------------------------------------------------- C --A.D. MCLACHLAN SEP 1980. LAST UPDATED 27 JAN 1988. C PRINTS A 3-D ARRAY IN SIMPLE FORMAT. SEVERAL ENTRIES C---------------------------------------------------------- C ENTRY -- PRNT3D(RHO,NU,NV,NW,HEAD) PRINT REAL ARRAY C ENTRY -- PRN3DI(IRHO,NU,NV,NW,HEAD) PRINT INTEGERS C next two commented-out for standardisation reasons C ENTRY -- PRN3DJ(JRHO,NU,NV,NW,HEAD) PRINT INTEGER*2 C ENTRY -- PRN3DK(KRHO,NU,NV,NW,HEAD) PRINT BYTE C ENTRY -- PRFTMR(FMTR,MROW) SET FORMAT FOR REAL PRNT3D C ENTRY -- PRFTMI(FMTI,MROW) SET FORMAT FOR INTEGER PRINTS C------------------------------------------------------------- C INPUT -- RHO(NU,NV,NW) REAL ARRAY C INPUT -- IRHO(NU,NV,NW) INTEGER ARRAY C INPUT -- JRHO(NU,NV,NW) INTEGER*2 ARRAY C INPUT -- KRHO(NU,NV,NW) BYTE ARRAY C INPUT -- NU X DIMENSION C INPUT -- NV Y DIMENSION C INPUT -- NW Z DIMENSION C INPUT -- HEAD CH*(*) HEADER LABEL C INPUT -- FMTR CH*(*) FORMAT FOR REAL NUMBER OUTPUT C INPUT -- FMTI CH*(*) FORMAT FOR INTEGER OUTPUT C INPUT -- MROW ITEMS PER ROW C INPUT -- NWRITE PRINT UNIT C------------------------------------------------------------- C NOTE -- DEFAULT MROW=16 FMTR = (1X,I4,' )',8(1X,2F7.4)) C NOTE -- DEFAULT FMTI = (1X,I4,' )',8(1X,2I5)) C%%C NOTE -- DEFAULT MROW=8 FMTR=(1X,I4,8(1PG12.4E2)) C%%C NOTE -- FMTI = (1X,I4,8I12) C NOTE -- IF THE ARRAY IS WIDER THAN MROW IT IS PRINTED IN PAGES OF MROW COLS C------------------------------------------------------------- REAL RHO(NU,NV,NW) INTEGER IRHO(NU,NV,NW) CCC INTEGER*2 JRHO(NU,NV,NW) C BYTE KRHO(NU,NV,NW) CHARACTER*(*) HEAD CHARACTER*(*) FMTR,FMTI CHARACTER*30 FORMR CHARACTER*30 FORMI C%% CHARACTER*30 FORMR/'(1X,I4,8(1PG12.4E2))'/ C%% CHARACTER*30 FORMI/'(1X,I4,8(1X,2I5))'/ DATA NROW/16/ C%% DATA NROW/8/ DATA FORMR/'(1X,I4,2H ),8(1X,2F7.4))'/ DATA FORMI/'(1X,I4,2H ),8(I12))'/ C----------------------------------------------------------- 20 FORMAT(1X) 22 FORMAT(1X,'<',A,'>',' : DIM=',3I5) 23 FORMAT(1X,10X,' EACH PLANE IN PAGES OF ',I5,' COLUMNS ') 24 FORMAT(1X,' PLANE IW=',I5,': - COLUMNS IX=',I5,' TO ',I5) 25 FORMAT(1X,4X,'-',I5,' -') C------------------------------- NVTYPE=1 GO TO 100 C++++++++++++++++++++++++++++++++++++++++++++++ ENTRY PRN3DI(IRHO,NU,NV,NW,HEAD,NWRITE) C++++++++++++++++++++++++++++++++++++++++++++++ NVTYPE=2 GO TO 100 C++++++++++++++++++++++++++++++++++++++++++++++ CCC ENTRY PRN3DJ(JRHO,NU,NV,NW,HEAD,NWRITE) CCCC++++++++++++++++++++++++++++++++++++++++++++++ CCC NVTYPE=3 CCC GO TO 100 C++++++++++++++++++++++++++++++++++++++++++++++ CCC ENTRY PRN3DK(KRHO,NU,NV,NW,HEAD,NWRITE) CCCC++++++++++++++++++++++++++++++++++++++++++++++ CCC NVTYPE=4 CCC GO TO 100 C------------------------------------- 100 CONTINUE WRITE(NWRITE,22) HEAD,NU,NV,NW IF(NU.GT.NROW) WRITE(NWRITE,23) NROW DO 150 IW=1,NW WRITE(NWRITE,20) IU1=1 130 CONTINUE IU2=IU1+NROW-1 IF(IU2.GT.NU) IU2=NU IF(NU.GT.NROW) THEN WRITE(NWRITE,24) IW,IU1,IU2 ELSE WRITE(NWRITE,25) IW END IF DO 140 IV=1,NV GO TO (131,132,133,134),NVTYPE 131 CONTINUE WRITE(NWRITE,FORMR) IV,(RHO(IU,IV,IW),IU=IU1,IU2) GO TO 140 132 CONTINUE WRITE(NWRITE,FORMI) IV,(IRHO(IU,IV,IW),IU=IU1,IU2) GO TO 140 133 CONTINUE CCC WRITE(NWRITE,FORMI) IV,(JRHO(IU,IV,IW),IU=IU1,IU2) GO TO 140 134 CONTINUE CCC WRITE(NWRITE,FORMI) IV,(KRHO(IU,IV,IW),IU=IU1,IU2) GO TO 140 140 CONTINUE IU1=IU1+NROW IF(IU1.LE.NU) GO TO 130 150 CONTINUE WRITE(NWRITE,20) RETURN C+++++++++++++++++++++++++++++++++++++ ENTRY PRFMTR(FMTR,MROW,NWRITE) C+++++++++++++++++++++++++++++++++++++ FORMR=FMTR NROW=MROW RETURN C+++++++++++++++++++++++++++++++++++++ ENTRY PRFMTI(FMTI,MROW,NWRITE) C+++++++++++++++++++++++++++++++++++++ FORMI=FMTI NROW=MROW RETURN C---------------- END C C VAXPDS NAME=R2CFTK C************************************************************************** SUBROUTINE R2CFTK (N, M, X0, Y0, X1, Y1, DIM) C-------------------------------------------------------------------------- C RADIX 2 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987. C-------------------------------------------------------------------------- INTEGER N, M, DIM(5) REAL X0(*),Y0(*),X1(*),Y1(*) C----- LOGICAL FOLD,ZERO INTEGER J,K,K0,M2,M OVER 2 INTEGER KK, MM2 INTEGER NS REAL C,IS,IU,RS,RU,S DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1 DOUBLE PRECISION DBANGL,DBFM2,DBCC DATA DBPI2/6.28318530717958623D0/ C C%% DBPI2=8.0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- M N DIVIDED BY CURRENT FACTORS C IN/OUT -- X0(*)...X1(*) REAL VALUES AT 2 EQUALLY SEPARATED DATA POINTS C IN/OUT -- Y0(*)...Y1(*) IMAGINARY VALUES AT 2 EQUALLY SEPARATED DATA POINTS C INPUT -- DIM(5) ARRAY DIMENSIONING CONSTANTS C----------------------------------------------------------------------------- C CALLS -- *** C----------------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED AND SIN COS RECALC A.D. MCLACHLAN AUG 1981 C-------------------------------------------------------------------------- LSEPU=DIM(1) LSEPV=DIM(2) LSEPW=DIM(3) LGLIMV=DIM(4) LGLIMW=DIM(5) NS=N*LSEPU M2=M*2 DBFM2=DBLE(M2) DBANGL=DBPI2/DBFM2 DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 M OVER 2=M/2+1 MM2 = LSEPU*M2 C-- DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*LSEPU + 1 ZERO= J.EQ.1 IF (ZERO) GO TO 200 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C=DBC1 S=DBS1 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*LSEPU + 1 C=-C 200 CONTINUE C-- DO 500 KK=K0,NS,MM2 LINDU=KK-LSEPV-LSEPW DO 440 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 420 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW RS=X0(K)+X1(K) IS=Y0(K)+Y1(K) RU=X0(K)-X1(K) IU=Y0(K)-Y1(K) X0(K)=RS Y0(K)=IS IF (ZERO) GO TO 300 X1(K)=RU*C+IU*S Y1(K)=IU*C-RU*S GO TO 400 300 CONTINUE X1(K)=RU Y1(K)=IU 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C-- RETURN END C C VAXPDS NAME=R3CFTK C************************************************************************** SUBROUTINE R3CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, DIM) C-------------------------------------------------------------------------- C RADIX 3 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C-------------------------------------------------------------------------- INTEGER N, M, DIM(5) REAL X0(*),Y0(*),X1(*),Y1(*),X2(*),Y2(*) C----- LOGICAL FOLD,ZERO INTEGER J,K,K0,M3,M OVER 2 INTEGER KK,MM3 INTEGER NS REAL A,B,C1,C2,S1,S2,T REAL I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBFM3 DATA DBPI2/6.28318530717958623D0/ DATA A/-0.5/, B/0.86602540/ C----------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- M N DIVIDED BY CURRENT FACTORS C IN/OUT -- X0(*)...X2(*) REAL VALUES AT 3 EQUALLY SEPARATED DATA POINTS C IN/OUT -- Y0(*)...Y2(*) IMAGINARY VALUES AT 3 EQUALLY SEPARATED DATA POINTS C INPUT -- DIM(5) ARRAY DIMENSIONING CONSTANTS C----------------------------------------------------------------------------- C CALLS -- *** C-------------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED AND SIN COS RECALC A.D. MCLACHLAN AUG 1981 C----------------------------------------------------------------------------- LSEPU=DIM(1) LSEPV=DIM(2) LSEPW=DIM(3) LGLIMV=DIM(4) LGLIMW=DIM(5) NS = N*LSEPU M3=M*3 DBFM3=DBLE(M3) DBANGL=DBPI2/DBFM3 DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 MM3 = LSEPU*M3 M OVER 2=M/2+1 C-- DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*LSEPU + 1 ZERO= J.EQ.1 IF (ZERO) GO TO 200 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C1=DBC1 S1=DBS1 C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*LSEPU + 1 T=C1*A+S1*B S1=C1*B-S1*A C1=T T=C2*A-S2*B S2=-C2*B-S2*A C2=T 200 CONTINUE C-- DO 500 KK = K0, NS, MM3 LINDU=KK-LSEPV-LSEPW DO 440 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 420 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW R0=X0(K) I0=Y0(K) RS=X1(K)+X2(K) IS=Y1(K)+Y2(K) X0(K)=R0+RS Y0(K)=I0+IS RA=R0+RS*A IA=I0+IS*A RB=(X1(K)-X2(K))*B IB=(Y1(K)-Y2(K))*B IF (ZERO) GO TO 300 R1=RA+IB I1=IA-RB R2=RA-IB I2=IA+RB X1(K)=R1*C1+I1*S1 Y1(K)=I1*C1-R1*S1 X2(K)=R2*C2+I2*S2 Y2(K)=I2*C2-R2*S2 GO TO 400 300 CONTINUE X1(K)=RA+IB Y1(K)=IA-RB X2(K)=RA-IB Y2(K)=IA+RB 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C-- RETURN END C C VAXPDS NAME=R4CFTK C************************************************************************** SUBROUTINE R4CFTK (N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,DIM) C-------------------------------------------------------------------------- C RADIX 4 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C-------------------------------------------------------------------------- INTEGER N, M, DIM(5) REAL X0(*),Y0(*),X1(*),Y1(*) REAL X2(*),Y2(*),X3(*),Y3(*) C----- LOGICAL FOLD,ZERO INTEGER J,K,K0,M4,M OVER 2 INTEGER KK,MM4 INTEGER NS REAL C1,C2,C3,S1,S2,S3,T REAL I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1 DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBFM4 DATA DBPI2/6.28318530717958623D0/ C C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- M N DIVIDED BY CURRENT FACTORS C IN/OUT -- X0(*)...X3(*) REAL VALUES AT 4 EQUALLY SEPARATED DATA POINTS C IN/OUT -- Y0(*)...Y3(*) IMAGINARY VALUES AT 4 EQUALLY SEPARATED DATA POINTS C INPUT -- DIM(5) ARRAY DIMENSIONING CONSTANTS C----------------------------------------------------------------------------- C CALLS -- *** C----------------------------------------------------------------------------- C NOTE -- INDEXING REVISED AND SIN COS RECALC A.D. MCLACHLAN AUG 1981 C-------------------------------------------------------------------------- LSEPU=DIM(1) LSEPV=DIM(2) LSEPW=DIM(3) LGLIMV=DIM(4) LGLIMW=DIM(5) NS = N*LSEPU M4=M*4 DBFM4=DBLE(M4) DBANGL=DBPI2/DBFM4 DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 MM4 = LSEPU*M4 M OVER 2=M/2+1 C-- DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*LSEPU + 1 ZERO= J.EQ.1 IF (ZERO) GO TO 200 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C1=DBC1 S1=DBS1 C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 C3=C2*C1-S2*S1 S3=S2*C1+C2*S1 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*LSEPU + 1 T=C1 C1=S1 S1=T C2=-C2 T=C3 C3=-S3 S3=-T 200 CONTINUE C-- DO 500 KK = K0, NS, MM4 LINDU=KK-LSEPV-LSEPW DO 440 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 420 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW RS0=X0(K)+X2(K) IS0=Y0(K)+Y2(K) RU0=X0(K)-X2(K) IU0=Y0(K)-Y2(K) RS1=X1(K)+X3(K) IS1=Y1(K)+Y3(K) RU1=X1(K)-X3(K) IU1=Y1(K)-Y3(K) X0(K)=RS0+RS1 Y0(K)=IS0+IS1 IF (ZERO) GO TO 300 R1=RU0+IU1 I1=IU0-RU1 R2=RS0-RS1 I2=IS0-IS1 R3=RU0-IU1 I3=IU0+RU1 X2(K)=R1*C1+I1*S1 Y2(K)=I1*C1-R1*S1 X1(K)=R2*C2+I2*S2 Y1(K)=I2*C2-R2*S2 X3(K)=R3*C3+I3*S3 Y3(K)=I3*C3-R3*S3 GO TO 400 300 CONTINUE X2(K)=RU0+IU1 Y2(K)=IU0-RU1 X1(K)=RS0-RS1 Y1(K)=IS0-IS1 X3(K)=RU0-IU1 Y3(K)=IU0+RU1 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C-- RETURN END C C VAXPDS NAME=R5CFTK C********************************************************************** SUBROUTINE R5CFTK (N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,DIM) C---------------------------------------------------------------------- C RADIX 5 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C---------------------------------------------------------------------- INTEGER N, M, DIM(5) REAL X0(*),Y0(*),X1(*),Y1(*),X2(*),Y2(*) REAL X3(*),Y3(*),X4(*),Y4(*) C----- LOGICAL FOLD,ZERO INTEGER J,K,K0,M5,M OVER 2 INTEGER KK,MM5 INTEGER NS REAL A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T REAL R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2 REAL I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2 DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBFM5 DATA A1/0.30901699/, B1/0.95105652/, 1 A2/-0.80901699/, B2/0.58778525/ DATA DBPI2/6.28318530717958623D0/ C C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- M N DIVIDED BY CURRENT FACTORS C IN/OUT -- X0(*)...X4(*) REAL VALUES AT 4 EQUALLY SEPARATED DATA POINTS C IN/OUT -- Y0(*)...Y4(*) IMAGINARY VALUES AT 4 EQUALLY SEPARATED DATA POINTS C INPUT -- DIM(5) ARRAY DIMENSIONING CONSTANTS C----------------------------------------------------------------------------- C CALLS -- *** C---------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED AND SIN COS RECALCULATED C NOTE -- A.D. MCLACHLAN. AUG 1981 C---------------------------------------------------------------------- LSEPU=DIM(1) LSEPV=DIM(2) LSEPW=DIM(3) LGLIMV=DIM(4) LGLIMW=DIM(5) NS = N*LSEPU M5=M*5 DBFM5=DBLE(M5) DBANGL=DBPI2/DBFM5 DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 MM5 = LSEPU*M5 M OVER 2=M/2+1 C-- DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*LSEPU + 1 ZERO= J.EQ.1 IF (ZERO) GO TO 200 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C1=DBC1 S1=DBS1 C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 C3=C2*C1-S2*S1 S3=S2*C1+C2*S1 C4=C2*C2-S2*S2 S4=S2*C2+C2*S2 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*LSEPU + 1 T=C1*A1+S1*B1 S1=C1*B1-S1*A1 C1=T T=C2*A2+S2*B2 S2=C2*B2-S2*A2 C2=T T=C3*A2-S3*B2 S3=-C3*B2-S3*A2 C3=T T=C4*A1-S4*B1 S4=-C4*B1-S4*A1 C4=T 200 CONTINUE C-- DO 500 KK = K0, NS, MM5 LINDU=KK-LSEPV-LSEPW DO 440 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 420 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW R0=X0(K) I0=Y0(K) RS1=X1(K)+X4(K) IS1=Y1(K)+Y4(K) RU1=X1(K)-X4(K) IU1=Y1(K)-Y4(K) RS2=X2(K)+X3(K) IS2=Y2(K)+Y3(K) RU2=X2(K)-X3(K) IU2=Y2(K)-Y3(K) X0(K)=R0+RS1+RS2 Y0(K)=I0+IS1+IS2 RA1=R0+RS1*A1+RS2*A2 IA1=I0+IS1*A1+IS2*A2 RA2=R0+RS1*A2+RS2*A1 IA2=I0+IS1*A2+IS2*A1 RB1=RU1*B1+RU2*B2 IB1=IU1*B1+IU2*B2 RB2=RU1*B2-RU2*B1 IB2=IU1*B2-IU2*B1 IF (ZERO) GO TO 300 R1=RA1+IB1 I1=IA1-RB1 R2=RA2+IB2 I2=IA2-RB2 R3=RA2-IB2 I3=IA2+RB2 R4=RA1-IB1 I4=IA1+RB1 X1(K)=R1*C1+I1*S1 Y1(K)=I1*C1-R1*S1 X2(K)=R2*C2+I2*S2 Y2(K)=I2*C2-R2*S2 X3(K)=R3*C3+I3*S3 Y3(K)=I3*C3-R3*S3 X4(K)=R4*C4+I4*S4 Y4(K)=I4*C4-R4*S4 GO TO 400 300 CONTINUE X1(K)=RA1+IB1 Y1(K)=IA1-RB1 X2(K)=RA2+IB2 Y2(K)=IA2-RB2 X3(K)=RA2-IB2 Y3(K)=IA2+RB2 X4(K)=RA1-IB1 Y4(K)=IA1+RB1 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C-- RETURN END C C VAXPDS NAME=R8CFTK C************************************************************************ SUBROUTINE R8CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4, 1 X5,Y5,X6,Y6,X7,Y7, DIM) C------------------------------------------------------------------------ C RADIX 8 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C------------------------------------------------------------------------ INTEGER N, M, DIM(5) REAL X0(*),Y0(*),X1(*),Y1(*),X2(*),Y2(*) REAL X3(*),Y3(*),X4(*),Y4(*),X5(*),Y5(*) REAL X6(*),Y6(*),X7(*),Y7(*) C----- LOGICAL FOLD,ZERO INTEGER J,K,K0,M8,M OVER 2 INTEGER KK,MM8 INTEGER NS REAL C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T REAL R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3 REAL I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3 REAL RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1 REAL ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1 DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1 DOUBLE PRECISION DBANGL,DBFM8 DATA DBPI2/6.28318530717958623D0/ DATA E/0.70710678/ C%% DBPI2=8.0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- M N DIVIDED BY CURRENT FACTORS C IN/OUT -- X0(*)...X7(*) REAL VALUES AT 8 EQUALLY SEPARATED DATA POINTS C IN/OUT -- Y0(*)...Y7(*) IMAGINARY VALUES AT 8 EQUALLY SEPARATED DATA POINTS C INPUT -- DIM(5) ARRAY DIMENSIONING CONSTANTS C----------------------------------------------------------------------------- C CALLS -- *** C----------------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED AND SIN COS RECALC A.D. MCLACHLAN AUG 1981 C----------------------------------------------------------------------------- LSEPU=DIM(1) LSEPV=DIM(2) LSEPW=DIM(3) LGLIMV=DIM(4) LGLIMW=DIM(5) NS = N*LSEPU M8=M*8 DBFM8=DBLE(M8) DBANGL=DBPI2/DBFM8 DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 MM8 = LSEPU*M8 M OVER 2=M/2+1 C-- DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*LSEPU + 1 ZERO= J.EQ.1 IF (ZERO) GO TO 200 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C1=DBC1 S1=DBS1 C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 C3=C2*C1-S2*S1 S3=S2*C1+C2*S1 C4=C2*C2-S2*S2 S4=S2*C2+C2*S2 C5=C4*C1-S4*S1 S5=S4*C1+C4*S1 C6=C4*C2-S4*S2 S6=S4*C2+C4*S2 C7=C4*C3-S4*S3 S7=S4*C3+C4*S3 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*LSEPU + 1 T=(C1+S1)*E S1=(C1-S1)*E C1=T T=S2 S2=C2 C2=T T=(-C3+S3)*E S3=(C3+S3)*E C3=T C4=-C4 T=-(C5+S5)*E S5=(-C5+S5)*E C5=T T=-S6 S6=-C6 C6=T T=(C7-S7)*E S7=-(C7+S7)*E C7=T 200 CONTINUE C-- DO 500 KK = K0, NS, MM8 LINDU=KK-LSEPV-LSEPW DO 440 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 420 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW RS0=X0(K)+X4(K) IS0=Y0(K)+Y4(K) RU0=X0(K)-X4(K) IU0=Y0(K)-Y4(K) RS1=X1(K)+X5(K) IS1=Y1(K)+Y5(K) RU1=X1(K)-X5(K) IU1=Y1(K)-Y5(K) RS2=X2(K)+X6(K) IS2=Y2(K)+Y6(K) RU2=X2(K)-X6(K) IU2=Y2(K)-Y6(K) RS3=X3(K)+X7(K) IS3=Y3(K)+Y7(K) RU3=X3(K)-X7(K) IU3=Y3(K)-Y7(K) RSS0=RS0+RS2 ISS0=IS0+IS2 RSU0=RS0-RS2 ISU0=IS0-IS2 RSS1=RS1+RS3 ISS1=IS1+IS3 RSU1=RS1-RS3 ISU1=IS1-IS3 RUS0=RU0-IU2 IUS0=IU0+RU2 RUU0=RU0+IU2 IUU0=IU0-RU2 RUS1=RU1-IU3 IUS1=IU1+RU3 RUU1=RU1+IU3 IUU1=IU1-RU3 T=(RUS1+IUS1)*E IUS1=(IUS1-RUS1)*E RUS1=T T=(RUU1+IUU1)*E IUU1=(IUU1-RUU1)*E RUU1=T X0(K)=RSS0+RSS1 Y0(K)=ISS0+ISS1 IF (ZERO) GO TO 300 R1=RUU0+RUU1 I1=IUU0+IUU1 R2=RSU0+ISU1 I2=ISU0-RSU1 R3=RUS0+IUS1 I3=IUS0-RUS1 R4=RSS0-RSS1 I4=ISS0-ISS1 R5=RUU0-RUU1 I5=IUU0-IUU1 R6=RSU0-ISU1 I6=ISU0+RSU1 R7=RUS0-IUS1 I7=IUS0+RUS1 X4(K)=R1*C1+I1*S1 Y4(K)=I1*C1-R1*S1 X2(K)=R2*C2+I2*S2 Y2(K)=I2*C2-R2*S2 X6(K)=R3*C3+I3*S3 Y6(K)=I3*C3-R3*S3 X1(K)=R4*C4+I4*S4 Y1(K)=I4*C4-R4*S4 X5(K)=R5*C5+I5*S5 Y5(K)=I5*C5-R5*S5 X3(K)=R6*C6+I6*S6 Y3(K)=I6*C6-R6*S6 X7(K)=R7*C7+I7*S7 Y7(K)=I7*C7-R7*S7 GO TO 400 300 CONTINUE X4(K)=RUU0+RUU1 Y4(K)=IUU0+IUU1 X2(K)=RSU0+ISU1 Y2(K)=ISU0-RSU1 X6(K)=RUS0+IUS1 Y6(K)=IUS0-RUS1 X1(K)=RSS0-RSS1 Y1(K)=ISS0-ISS1 X5(K)=RUU0-RUU1 Y5(K)=IUU0-IUU1 X3(K)=RSU0-ISU1 Y3(K)=ISU0+RSU1 X7(K)=RUS0-IUS1 Y7(K)=IUS0+RUS1 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C-- RETURN END C C ================================ SUBROUTINE RCSFFT(MSG,LSG) C ================================ C C This reads in the input cards, for fft, using keywords. C C MSG(LSG) list of allowed FFT spacegroups C C .. Parameters .. INTEGER NPAR PARAMETER (NPAR=200) INTEGER MCOLS,MXCOLS,NLI,ECOLS C NLI is number of input program labels PARAMETER (MCOLS=500,MXCOLS=500,NLI=19) PARAMETER (ECOLS=MCOLS-NLI) INTEGER LSYM PARAMETER (LSYM=192) INTEGER NKEY PARAMETER (NKEY=25) C .. C .. Scalar Arguments .. INTEGER LSG C .. C .. Array Arguments .. INTEGER MSG(LSG) C .. C .. Scalars in Common .. REAL A,B,BIAS,DEL,DEL1,F000,F1LIM,F2LIM,RHMAX,RHMEAN,RHMIN, + RHOMAX,RHOMIN,RMAX,RMIN,SCLF1,SCLF2,SCMOD1,SCMOD2,SIG1,SIG2, + SMAX,SMIN,SSQ,TFF1,TFF2,V,F1MAX,F2MAX,GRDSMP,CX,CY,CZ INTEGER HMAX,HMIN,ICEN,IDIG,IFreeRexcludeVal,ILIN,IPHTR,IPROJ, + ISORT,ISP,ISQ,IW,J,K,KE,KMAX,KMIN,KSYM,L,LCMAX,LCMIN, + LIST,LMAX,LMIN,LSPGRP,M,NABCH,NANOM,NCOL,NF1,NF2,NH,NI1, + NLAUE,NLINES,NLPRGI,NNA,NNB,NP,NP2,NPH,NS1,NS2,NSECAX,NSG, + NSMFFT,NSPFFT,NSSM,NSYM,NSYMCH,NSYMP,NUM_FREE_COLUMN, + NUSED10,NUSED11,NUSED6,NUSED7,NUSED8,NUSED9,NW,NW2,NX,NY, + NZ,T2,T3,T4,TYPE,XMAX,XMIN,YMAX,YMIN,ZMAX,ZMIN,NLAFFT CHARACTER TITLE*80,NAMSPG*10,PGNAME*10 LOGICAL USEFILL,USEFREE C .. C .. Arrays in Common .. REAL CELL,CS,PERM,RSMFFT,RSYM,RSYMT,SN,XYZLIM INTEGER IUVW,LOOKUP,MP,T C .. C .. Local Scalars .. REAL BB,RHOTEM,SCAL,VOL,XTMP INTEGER I,IERR,IFIL,II,ILPRGI,ISYSW,ITOK,JJ,JTOK,N, + NCHK,NCOLS,NREFLS,NSTART,NTOK,TMP,IPRINT,TEMPSPG LOGICAL LEND,LASU,LAXIS CHARACTER LTYPE*1,KEY*4,WORD*4,DUMMY*10,TEMPNAM*10, + LAUNAM*10,FILIN*70,LINE*400 C .. C .. Local Arrays .. REAL CRANGE(2,MXCOLS),FVALUE(NPAR),RR(3,3,6),RSMINV(4,4) INTEGER IBEG(NPAR),IDEC(NPAR),IEND(NPAR),ITYP(NPAR) CHARACTER CTPRGI(MXCOLS)*1,KXYZ(3)*1, + CVALUE(NPAR)*4,KEYWRD(NKEY)*4,CDUMMY(MXCOLS)*30, + LSPRGI(MCOLS)*30,LSUSRJ(MCOLS)*30 C .. C .. External Functions .. INTEGER LENSTR EXTERNAL LENSTR C .. C .. External Subroutines .. EXTERNAL APPNDS,CCPERR,CCPUPC,GTPINT,GTPREA,INVSYM,LKYIN,LKYSET, + LRASSN,LRCELL,LRCLAB,LRINFO,LROPEN,LRRSOL,LRSYMI,LRSYMM, + PARSER,RBFRO1,RDRESO,RDSCAL,RDSYMM C .. C .. Intrinsic Functions .. INTRINSIC INT,MAX,MOD,SQRT C .. C .. Common blocks .. COMMON /ATSYM/RSYM(4,4,LSYM),RSYMT(4,4,LSYM),NSYM,PERM(4,4),NSYMP, + NLAUE COMMON /FFTSPG/NSPFFT,NSMFFT,RSMFFT(4,4,LSYM),NLAFFT COMMON /FOUT/LCMIN,LCMAX,ISP,IDIG,ILIN COMMON /HKLF/A,B,ISORT,J,K,L,M,SSQ COMMON /HKLLIM/HMAX,KMAX,LMAX,HMIN,KMIN,LMIN COMMON /IOFFT/T2,T3,T4,NSYMCH,NABCH,NSECAX COMMON /LFTOVR/MP(3),KE,RMAX,RMIN,XYZLIM(2,3),GRDSMP COMMON /MAPFMT/TYPE,NLINES,NCOL,V,F000,RHOMAX,RHOMIN COMMON /MPHDR/NX,NY,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,NUSED6, + NUSED7,NUSED8,NUSED9,NUSED10,IUVW(3),NUSED11,CELL(6), + LSPGRP,RHMIN,RHMAX,RHMEAN COMMON /MPHDRR/TITLE COMMON /SEL/NF1,NS1,NF2,NS2,NP,NW,IW,ISQ,NH,SMIN,SMAX,SCLF1,TFF1, + SCLF2,TFF2,BIAS,NSSM,NSG,IPROJ,LIST,ICEN,KSYM,NI1,CS(24), + SN(24),T(LSYM,3),NNA,NNB,NPH,NANOM,NW2,NP2,SIG1,SIG2,DEL, + DEL1,F1LIM,F2LIM,IPHTR,SCMOD1,SCMOD2,F1MAX,F2MAX,CX,CY,CZ COMMON /SPACGP/ NAMSPG, PGNAME COMMON /FREESTUFF/ USEFREE, NUM_FREE_COLUMN,IFreeRexcludeVal, + USEFILL COMMON /LOOKCOM/ NLPRGI, LOOKUP (MCOLS) C .. C .. Save statement .. SAVE C .. C .. Data statements .. DATA LSPRGI/'H','K','L','A','B','F1','SIG1','PHI','F2', + 'SIG2','FH','PHIH','W','DANO','W2','PHI2','FREE','I','SIGI', + ECOLS*' '/ DATA CTPRGI/3*'H', 2*'R', ' ', 'Q', 'P', ' ', 'Q', ' ', 'P', 'W', + 'D', 'W', 'P', 'I', 'J','Q', ECOLS*' '/ DATA FILIN/'HKLIN'/ DATA KXYZ/'X','Y','Z'/ C---- Key Words DATA KEYWRD/'LIST','AXIS','BIAS','BINM','EXCL','SCAL','FFTS', + 'FILL','GRID','END ','SYMM','LPFM','LPMA','LABI','PROJ', + 'PATT','RESO','RHOL','NOCH','TITL','XYZL','VF00','PHTR', + 'HKLM','FREE'/ C .. C NLPRGI = number of input program labels NLPRGI = NLI DO 1, I=1,3 LOOKUP (I) = -1 1 CONTINUE DO 2, I=4,MCOLS LOOKUP (I) = 0 2 CONTINUE USEFREE= .FALSE. USEFILL=.FALSE. ISYSW = 6 NSYMP = -1 NSPFFT = 1 C Set NLAFFT = 3 - P1 value. NLAFFT = 3 NSMFFT = -1 IPROJ = 0 IFreeRexcludeVal = 0 RMIN = -1.0 NX = -1 NY = -1 NZ = -1 XYZLIM(1,1) = -1.0 C Set HKLMAX flags for explicit input HMAX = -1 KMAX = -1 LMAX = -1 C Set flag for AXIS keyword LAXIS = .FALSE. C Set flag for XYZLIM ASU option LASU = .FALSE. C Default grid sample = 3 ie 1/3 resolution GRDSMP = 3.0 C C---- default map scale uses true volume and f000 = 0.0 C TYPE = -1 F000 = 0.0 C C---- set vol = 0.0 for starters C V = 0.0 TITLE = ' Map from fft' C C---- Start of parse loop C 10 CONTINUE LINE = ' ' KEY = ' ' NTOK = NPAR C C *********************************************************** CALL PARSER(KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK,LEND, + .TRUE.) C *********************************************************** C C---- End of file? C IF (LEND) THEN GO TO 1000 ELSE C C---- extra check for comment line (!) C IF (KEY(1:1).EQ.'!' .OR. NTOK .EQ. 0) THEN GO TO 10 ELSE C C---- Loop over possible key-words C DO 20 I = 1,NKEY IF (KEY.EQ.KEYWRD(I)) GO TO 30 20 CONTINUE C C---- Keyword not found C WRITE (ISYSW,FMT=6000) LINE(IBEG(1) :IEND(1)) 6000 FORMAT (/' *** ERROR - unrecognized keyword: ',A) KE = KE + 1 GO TO 10 C 30 CONTINUE C DATA KEYWRD/'LIST','AXIS','BIAS','BINM','EXCL','SCAL','FFTS', C + 'FILL','GRID','END ','SYMM','LPFM','LPMA','LABI','PROJ', C + 'PATT','RESO','RHOL','NOCH','TITL','XYZL','VF00','PHTR', C + 'HKLM','FREE'/ GO TO ( 40, 50, 90, 100, 110, 130, 150, + 155, 160, 1000, 170, 180, 190, 200, 230, + 240, 250, 260, 270, 280, 290, 320, 330, + 340, 345) I C C---- Listing "LIST" - lists information about some reflections C ========================================================== C 40 CONTINUE C LIST = 1 C ****************************** CALL GTPINT(2,LIST,NTOK,ITYP,FVALUE) C ****************************** GO TO 10 C C---- Axis orientations (Z,X,Y) "AXIS" C ================================== C 50 IF (NTOK.EQ.4) THEN C DO 80 I = 2,4 C C ************ CALL CCPUPC(LINE) C ************ C DO 60 J = 1,3 IF (LINE(IBEG(I) :IEND(I)).EQ.KXYZ(J)) GO TO 70 60 CONTINUE C WRITE (ISYSW,FMT=6002) LINE(IBEG(2) :IEND(NTOK)) 6002 FORMAT (' *** There is an error in "AXIS"',A10,/) KE = KE + 1 GO TO 80 70 MP(I-1) = J 80 CONTINUE ELSE WRITE (ISYSW,FMT='(A)') ' ** Too few tokens on AXIS **' KE = KE + 1 END IF LAXIS = .TRUE. GO TO 10 C C---- Bias "BIAS" C ============ C 90 CONTINUE C C ****************************** CALL GTPREA(2,BIAS,NTOK,ITYP,FVALUE) C ****************************** C GO TO 10 C C---- Binary map flag "BINM" C ====================== C C Subsidiary keyword NONE to switch off C 100 T3 = 1 C IF (NTOK.GT.1) THEN KEY = CVALUE(2) C C *********** CALL CCPUPC(KEY) C *********** C IF (KEY.EQ.'NONE') THEN T3 = 0 END IF END IF GO TO 10 C C---- Exclude reflections "EXCLUDE" C ============================= C C Subsidiary key words SIG1 sig1 SIG2 f1lim f2lim term diff f1max f2max C 110 CONTINUE NSTART = 2 NCHK = MOD(NTOK,2) C IF (NCHK.NE.1) THEN WRITE (ISYSW,FMT='(A)') ' ** Error on EXCLUDE card **' KE = KE + 1 GO TO 10 END IF C DO 120 ITOK = NSTART,NTOK,2 KEY = CVALUE(ITOK) C C *********** CALL CCPUPC(KEY) C *********** C IF (KEY.EQ.'SIG1') THEN SIG1 = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'SIG2') THEN SIG2 = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'F1LI') THEN F1LIM = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'ILIM') THEN F1LIM = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'F2LI') THEN F2LIM = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'DIFF') THEN DEL = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'TERM') THEN DEL1 = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'F1MA') THEN F1MAX = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'IMAX') THEN F1MAX = FVALUE(ITOK+1) ELSEIF (KEY.EQ.'F2MA') THEN F2MAX = FVALUE(ITOK+1) ENDIF 120 CONTINUE GO TO 10 C C---- SCALE Label Scale B-factor C =========================== C C F1 or F2 Scale and temperature factor C 130 CONTINUE C ITOK = 2 140 IF (ITOK .LE. NTOK .AND. ITOK .GT. 0) THEN C C---- Subroutine RDSCAL extracts from line, beginning at field ITOK, C sets of 3 arguments:- Label (= F1 or F2); Scale; Bfactor C C Input: LSPRGI list of NLPRGI possible labels C Output: ILPRGI label number found C SCAL scale C BB temperature factor C ITOK next item to read, = 0 if finished, .lt. 0 if error C I = ITOK C *************************************************** CALL RDSCAL(ITOK,LINE,IBEG,IEND,ITYP,FVALUE,NTOK,NLPRGI, + LSPRGI,ILPRGI,SCAL,BB) C *************************************************** IF (ITOK .EQ. I) THEN C Allow for old version of RDSCAL which doesn't reset ITOK ITOK = ITOK + 3 ELSEIF (ITOK .EQ. -1) THEN WRITE (ISYSW, '(A)') $ ' **** No scale factor given ****' KE = KE + 1 ELSEIF (ITOK .EQ. -2) THEN WRITE (ISYSW, '(A)') $ ' **** Unrecognized program label ****' KE = KE + 1 ENDIF C IF (ITOK .GE. 0) THEN WRITE (ISYSW,FMT=6004) $ LSPRGI(ILPRGI)(1:LENSTR(LSPRGI(ILPRGI))),SCAL,BB 6004 FORMAT (// $ ' Scale factor and Temperature factor for',7X $ ,A,4X,2F9.4) WORD = LSPRGI(ILPRGI) C IF (WORD(1:2).EQ.'F1') THEN SCLF1 = SCAL TFF1 = BB ELSE IF (WORD(1:2).EQ.'F2') THEN SCLF2 = SCAL TFF2 = BB END IF GO TO 140 ENDIF ENDIF GO TO 10 C C---- FFT Space group (should be 2,10,47 for patterson) "FFTS" C =========================================================== C 150 CONTINUE WRITE (ISYSW, + '('' *** FFT SG ignored in this P1-only version'')') GO TO 10 C C---- FILLIN - F1 is substituted by F2 is F1 is missing C ================================================== C 155 CONTINUE USEFILL=.TRUE. GOTO 10 C C---- "GRID" - Number of divisions along each whole unit cell edge C ============================================================ C 160 CONTINUE KEY = CVALUE(2) CALL CCPUPC(KEY) IF (KEY .EQ. 'SAMP') THEN CALL GTPREA(3,GRDSMP,NTOK,ITYP,FVALUE) IF (GRDSMP .LT. 2.0) THEN WRITE (ISYSW,FMT='(A)') $ ' ** GRID sampling must be at least 2 **' KE = KE + 1 ENDIF ELSEIF (NTOK.EQ.4) THEN C C ****************************** CALL GTPINT(2,NX,NTOK,ITYP,FVALUE) CALL GTPINT(3,NY,NTOK,ITYP,FVALUE) CALL GTPINT(4,NZ,NTOK,ITYP,FVALUE) C ****************************** C ELSE WRITE (ISYSW,FMT='(A)') ' ** Too few tokens on GRID **' KE = KE + 1 END IF GO TO 10 C C---- SYMMETRY the true space group for the data C =========================================== C 170 CONTINUE JTOK = 2 C C ************************************************** CALL RDSYMM(JTOK,LINE,IBEG,IEND,ITYP,FVALUE,NTOK,NAMSPG, + LSPGRP,PGNAME,NSYM,NSYMP,RSYM) C ************************************************** C C---- Get NLAUE C C Default unset point group to P1 IF (PGNAME .EQ. ' ') THEN PGNAME = 'PG1' ENDIF CALL PGDEFN(PGNAME,NSYMP,NSYM,RSYM,.FALSE.) CALL PGNLAU(PGNAME,NLAUE,LAUNAM) C C For fft hkl generation we only need the primitive symmetry operators. C The output map reads the SYMOP.LIB again and copies all the symmetry C operators into the output file. C NSYM = NSYMP GO TO 10 C C---- Line printer MAP FORMAT "LPFM" C ============================== C 180 CONTINUE C C ****************************** CALL GTPINT(2,NLINES,NTOK,ITYP,FVALUE) CALL GTPINT(3,NCOL,NTOK,ITYP,FVALUE) CALL GTPINT(4,IDIG,NTOK,ITYP,FVALUE) CALL GTPINT(5,ILIN,NTOK,ITYP,FVALUE) CALL GTPINT(6,LCMIN,NTOK,ITYP,FVALUE) CALL GTPINT(7,LCMAX,NTOK,ITYP,FVALUE) C ****************************** C IF (IDIG.LT.1 .OR. IDIG.GT.4) IDIG = 4 IF (ILIN.GT.3 .OR. ILIN.LT.1) ILIN = 1 C IF (LCMAX.LT.LCMIN) THEN TMP = LCMAX LCMAX = LCMIN LCMIN = TMP END IF C GO TO 10 C C---- "LPMA" Line printer map flag C ============================ C 190 T4 = 6 GO TO 10 C C---- "LABI" i.e. set input labels C ================================ C 200 CONTINUE C ITOK = 2 C C ****************************************************** CALL LKYSET(LSPRGI,NLPRGI,LSUSRJ,LOOKUP,ITOK,NTOK,LINE, + IBEG,IEND) C ****************************************************** C C---- Set up Label Keys. Needed for LRASSN C C **************************************** CALL LKYIN(1,LSPRGI,NLPRGI,NTOK,LINE,IBEG,IEND) C **************************************** C GO TO 10 C C---- "PROJ" Map axis projection - down section axis C ================================================ C 230 CONTINUE C C---- Choose sectioning axis at end C IPROJ = 1 GO TO 10 C C---- "PATT" Patterson C ================= C 240 ISQ = 1 GO TO 10 C C---- "RESO"lution C ============== C 250 CONTINUE JTOK = 2 C C ************************************************* CALL RDRESO(JTOK,ITYP,FVALUE,NTOK,RMIN,RMAX,SMIN,SMAX) C ************************************************* C GO TO 10 C C---- "RHOL" Rho limits C ==================== C 260 CONTINUE RHOMAX = 0.0 TYPE = 0 RHOMIN = 0.0 C C ********************************* CALL GTPREA(2,RHOMAX,NTOK,ITYP,FVALUE) C ********************************* IF (NTOK.GT.2) THEN TYPE = 1 C ********************************* CALL GTPREA(3,RHOMIN,NTOK,ITYP,FVALUE) C ********************************* END IF C IF (RHOMAX.LT.RHOMIN) THEN RHOTEM = RHOMAX RHOMAX = RHOMIN RHOMIN = RHOTEM END IF C GO TO 10 C C---- "NOCHeck" C ========= C C - You do not wish to generate any symmetry C equivalent data; If you change axis order, C or do anything fancy it will be overridden C 270 CONTINUE ISORT = 0 NSYM = 0 GO TO 10 C C---- "TITLE" Title C =============== C 280 CONTINUE IF(NTOK.GE.2)TITLE = LINE(IBEG(2) :IEND(NTOK)) GO TO 10 C C---- "XYZLIM" X Y and Z Limits, grid or fractional C ============================================== C 290 CONTINUE C C 22/12/98 - XYZLIM ASU option will get map limits from a call to C s/r SETLIM, otherwise limits set as normal. C IF (NTOK.EQ.2) THEN C IF (CVALUE(2).EQ.'ASU '.OR. CVALUE(2).EQ.'asu ') THEN WRITE (ISYSW,FMT='(A,A)') + ' Map limits will be set to the CCP4 asymmetric', + ' unit' LASU = .TRUE. ELSE WRITE (ISYSW,FMT='(A,/,2A)') + ' ** Unrecognised argument for XYZLIM **', + ' Must be ASU or else 6 numbers specifying', + ' map limits' END IF C ELSE IF (NTOK.NE.7) THEN C KE = KE + 1 WRITE (ISYSW,FMT='(A)') ' ** Too few tokens on XYZLIM **' GO TO 10 C ELSE C I = 1 C DO 310 II = 1,3 DO 300 JJ = 1,2 I = I + 1 C C *************************************** CALL GTPREA(I,XYZLIM(JJ,II),NTOK,ITYP,FVALUE) C *************************************** C 300 CONTINUE C C---- Reorder if needed C IF (XYZLIM(2,II).LT.XYZLIM(2,II)) THEN XTMP = XYZLIM(1,II) XYZLIM(1,II) = XYZLIM(2,II) XYZLIM(2,II) = XTMP END IF C 310 CONTINUE C END IF C GO TO 10 C C---- "VF00" V and F000 Rho scaling C =============================== C 320 CONTINUE F000 = 0.0 C C ****************************** CALL GTPREA(2,V,NTOK,ITYP,FVALUE) CALL GTPREA(3,F000,NTOK,ITYP,FVALUE) C ****************************** C GO TO 10 C C---- "PHTRANSLATION" C ================== C C Phased translation calculation - terms used are: C w*F1*F2 exp(phi1 - phi2) or w*F1*F2 exp(phi1 + PHI2) C 330 IPHTR = 1 IF(NTOK.GT.2)THEN CALL CCPUPC(CVALUE(2)) IF(CVALUE(2).EQ.'HAND') + CALL GTPINT(3,IHAND,NTOK,ITYP,FVALUE) IF(IHAND.EQ.2) CALL HANDCHANGE(LSPGRP,CX,CY,CZ) IF(IHAND.EQ.2) IPHTR = -1 END IF GO TO 10 C C---- "HKLMAX" C ======== C C Set hkl maxima: only for anisotropic resolution C Ignore empty fields, to leave HMAX, KMAX, LMAX C to be set later from resolution limits C (note that HMAX, KMAX, LMAX are set = -1 at beginning C of this routine) 340 IF (ITYP(2) .EQ. 2) CALL GTPINT(2,HMAX,NTOK,ITYP,FVALUE) IF (ITYP(3) .EQ. 2) CALL GTPINT(3,KMAX,NTOK,ITYP,FVALUE) IF (ITYP(4) .EQ. 2) CALL GTPINT(4,LMAX,NTOK,ITYP,FVALUE) GO TO 10 C C---- Free R factor - Value in table to exclude. C ===== C 345 CONTINUE C C ********************************* CALL GTPINT(2,IFreeRexcludeVal,NTOK,ITYP,FVALUE) C ********************************* GO TO 10 C END IF END IF C C---- "END " Must be at end C ========================= C 1000 CONTINUE C C For fft hkl generation we only need the primitive symmetry operators. C The output map reads the SYMOP.LIB again and copies all the symmetry C operators into the output file. C IFIL = 1 IPRINT = 2 IERR = -1 C C ****************************** CALL LROPEN(IFIL,FILIN,IPRINT,IERR) C ****************************** C IF (IERR.EQ.-1) THEN WRITE (ISYSW,FMT='(A)') ' *** NO SUCH FILE !**' C C ************************************** CALL CCPERR(1,' ** File does not exist **') C ************************************** C END IF C C---- Read file info, mostly to get hkl range C C ************************************** CALL LRINFO(IFIL,DUMMY,NCOLS,NREFLS,CRANGE) C ************************************** C C---- Set up true symmetry C IF (NSYMP.LE.0) THEN C C No true symmetry read from control data, so read from HKLIN file C For rest of program, we need only C NSYM = NSYMP (primitive only) C RSYM symmetry operations (primitive) for expansion C LSPGRP space group number to write symmetry to map file C C ********************************************* CALL LRSYMI(IFIL,NSYMP,LTYPE,LSPGRP,NAMSPG,PGNAME) CALL LRSYMM(IFIL,NSYM,RSYM) C---- Get NLAUE CALL PGDEFN (PGNAME, NSYMP, NSYM, RSYM, .FALSE.) CALL PGNLAU(PGNAME,NLAUE,LAUNAM) C ********************************************* C If PATT set this may need changing. C NSYM = NSYMP END IF C C---- If LASU then get XYZLIM limits from SETLIM to correspond C to the CCP4 asymmetric unit for the spacegroup of the map C IF (LASU) THEN C C---- Check for Patterson first...want to set limits for map C IF (ISQ. EQ. 0) THEN TEMPSPG = LSPGRP ELSE CALL PATSGP(NAMSPG, PGNAME, TEMPNAM, TEMPSPG) END IF C CALL SETLIM (TEMPSPG, XYZLIM) WRITE (ISYSW,FMT=7002)TEMPSPG, XYZLIM 7002 FORMAT (/' ** The limits of the map have been set to those of', + /' the CCP4 asymmetric unit for appropriate spacegroup', + //' For map spacegroup ',I4,' these limits are:' + //' Min X Max X Min Y Max Y Min Z Max Z', + /,1X,6F8.3, + //' expressed in fractional coordinates and with x,y,z', + /' referring to the unpermutated axes.',/) END IF C C---- Set up FFT symmetry C If not explicitly given set = true symmetry C IF (NSMFFT.LE.0) THEN C C---- Check if true symmetry is one of allowed spacegroups C c DO 1350 N = 1,LSG c IF (MSG(N).EQ.LSPGRP) GO TO 1360 c 1350 CONTINUE C C---- No, so set spacegroup = P1 or P-1 C IF (ISQ.EQ.0 .OR. LSG.EQ.1) THEN C C--- P1 C NSMFFT = 1 NSPFFT = 1 NLAFFT = 3 ELSE C C-----P-1 C NSMFFT = 1 NSPFFT = 2 NLAFFT = 3 END IF C DO 1353 N = 1,NSMFFT DO 1352 II = 1,4 DO 1351 JJ = 1,4 RSMFFT(II,JJ,N) = RSYM(II,JJ,N) 1351 CONTINUE 1352 CONTINUE 1353 CONTINUE C GO TO 1370 C C---- True spacegroup OK, use it C 1360 CONTINUE NSMFFT = NSYM NLAFFT = NLAUE IF (ISQ.EQ.0) THEN C -- Not a patterson NSPFFT = LSPGRP C FFTSG matrices from input symmtery DO 1400 N = 1,NSMFFT DO 1390 II = 1,4 DO 1380 JJ = 1,4 RSMFFT(II,JJ,N) = RSYM(II,JJ,N) 1380 CONTINUE 1390 CONTINUE 1400 CONTINUE ELSE C These pattersons can be set. No others.. IF(NLAFFT.EQ.4) NSPFFT = 10 IF(NLAFFT.EQ.6) NSPFFT = 47 IF(NLAFFT.EQ.3) NSPFFT = 2 IF(NSPFFT.LE.0) + CALL CCPERR(1,' Please set FFTSPG for Patterson') C Get FFTSG symmetry matrices from reset SG C the translation parts will be set back to 0.0 later CALL MSYGET(24, NSPFFT, NSMFFT, RSMFFT) END IF 1370 CONTINUE END IF C C FFT/FFTBIG: the default axis ordering is Z,X,Y for P1, monoclinic C space groups, and space groups 16, 17 and 18, and Y,X,Z for all C higher symmetry space groups. Default map volume in FFTBIG is now C always the whole unit cell. C IF (NSPFFT.EQ.1 .AND. LSPGRP.GE.19 .AND. .NOT.LAXIS) THEN MP(1) = 2 MP(2) = 1 MP(3) = 3 END IF C C---- Get inverse for symmetry - leave RSMFFT(4,??) as the real vector C DO 1430 N = 1,NSMFFT C C ********************************* CALL INVSYM(RSMFFT(1,1,N),RSMINV(1,1)) C ********************************* C DO 1420 II = 1,3 DO 1410 JJ = 1,3 RSMFFT(II,JJ,N) = RSMINV(II,JJ) 1410 CONTINUE 1420 CONTINUE 1430 CONTINUE C C---- if vol not set AND TYPE = -1 work it out C C ******************** CALL LRCELL(IFIL,CELL) CALL RBFRO1(CELL,VOL,RR) C ******************** C IF (V.EQ.0.0) V = VOL C C---- Set up resolution limits if not set C IF (RMIN.LT.0.0) THEN C C ********************** CALL LRRSOL(IFIL,SMIN,SMAX) C ********************** C C RMIN = 1.0/SQRT(SMIN) C RMAX = 1.0/SQRT(SMAX) RMIN = 10000 RMAX = 0.5 IF (SMIN.GT.0.0) RMIN = 1.0/SQRT(SMIN) IF (SMAX.GT.0.0) RMAX = 1.0/SQRT(SMAX) C C---- If resolution limit taken from file, C set index limits == limits from file, unless explicitly set C IF (HMAX .LT. 0) $ HMAX = INT(MAX(-CRANGE(1,1),CRANGE(2,1),CELL(1)/RMAX)) IF (KMAX .LT. 0) $ KMAX = INT(MAX(-CRANGE(1,2),CRANGE(2,2),CELL(2)/RMAX)) IF (LMAX .LT. 0) $ LMAX = INT(MAX(-CRANGE(1,3),CRANGE(2,3),CELL(3)/RMAX)) ELSE C . . otherwise set them from resolution command IF (HMAX .LT. 0) $ HMAX = INT(CELL(1)/RMAX) IF (KMAX .LT. 0) $ KMAX = INT(CELL(2)/RMAX) IF (LMAX .LT. 0) $ LMAX = INT(CELL(3)/RMAX) END IF C IERR = -1 C C **************************************** CALL LRASSN(IFIL,LSPRGI,NLPRGI,LOOKUP,CTPRGI) C **************************************** C C---- Actual column types returned from LRASSN in program order C C---- Get actual column types (& labels) in file order C C ******************************** CALL LRCLAB(IFIL,CDUMMY,CTPRGI,NCOLS) C ******************************** C C------ Set up lookup pointers, & check column types C C Count assignments other than A, B C NCHK = 0 C DO 1440 I = 6,NLPRGI NCHK = LOOKUP(I) + NCHK 1440 CONTINUE C IF (LOOKUP(4).NE.0) NNA = LOOKUP(4) IF (LOOKUP(5).NE.0) NNB = LOOKUP(5) C IF (NNA.EQ.0) THEN IF (NNB.NE.0) THEN C C---- Can't have B without A C GO TO 1450 ELSE IF (NCHK.EQ.0) THEN C C---- If not A&B must have something else assigned C GO TO 1450 END IF ELSE C C---- A is assigned, must have B as well C IF (NNB.EQ.0) THEN GO TO 1450 ELSE C C---- Both A & B, check types, should both be real C IF (CTPRGI(NNA).NE.'R' .OR. CTPRGI(NNB).NE.'R') THEN WRITE (ISYSW,FMT='(A)') + ' ** Columns A & B must both be of type R **' KE = KE + 1 END IF END IF END IF C Pregenerated patterson terms IF (LOOKUP(18) .NE.0) THEN ISQ = -1 NI1 = LOOKUP(18) IF (CTPRGI(NI1).NE.'J') THEN WRITE (ISYSW,FMT='(A)') + ' ** Column I must be of type J **' KE = KE + 1 END IF IF (NI1.NE.0) THEN IF(LOOKUP(4).NE.0 .OR.LOOKUP(6) .NE.0 .OR.LOOKUP(14).NE.0)THEN WRITE (ISYSW,FMT='(A)') + ' ** Column I set; none of F or A or DANO allowed **' KE = KE + 1 END IF END IF IF(LOOKUP(19).NE.0) THEN NS1=LOOKUP(19) C---- SIG column can be of any type, but print warning if not Q, L or M C IF (CTPRGI(NS1).NE.'Q' .AND. CTPRGI(NS1).NE.'L'.AND. + CTPRGI(NS1).NE.'M') THEN WRITE (ISYSW,FMT='(A)') + ' ** WARNING ** column SIG1 is not of type Q, L or M' END IF END IF IF (SIG1.GT.0.0 .AND. LOOKUP(19).EQ.0) + CALL CCPERR(1,' EXCLUDE on SIG1 requested but SIG1 unassigned') END IF C C ---- F1 IF (LOOKUP(6).NE.0) THEN NF1 = LOOKUP(6) IF (CTPRGI(NF1).NE.'F' .AND. CTPRGI(NF1).NE.'J'.AND. + CTPRGI(NF1).NE.'D' .AND. CTPRGI(NF1).NE.'G'.AND. + CTPRGI(NF1).NE.'K') THEN WRITE (ISYSW,FMT='(A)') + ' ** Column F1 must be of type F or J or D or G or K **' KE = KE + 1 END IF IF (SIG1.GT.0.0 .AND. LOOKUP(7).EQ.0) + CALL CCPERR(1,' EXCLUDE on SIG1 requested but SIG1 unassigned') END IF C ---- SIG1 IF (LOOKUP(7).NE.0) THEN NS1 = LOOKUP(7) C C---- SIG column can be of any type, but print warning if not Q, L or M C IF (CTPRGI(NS1).NE.'Q' .AND. CTPRGI(NS1).NE.'L'.AND. + CTPRGI(NS1).NE.'M') THEN WRITE (ISYSW,FMT='(A)') + ' ** WARNING ** column SIG1 is not of type Q, L or M' END IF END IF C ---- PHI IF (LOOKUP(8).NE.0) THEN NP = LOOKUP(8) IF (CTPRGI(NP).NE.'P') THEN WRITE (ISYSW,FMT='(A)') ' ** Column PHI must be of type P **' KE = KE + 1 END IF END IF C ---- F2 IF (LOOKUP(9).NE.0) THEN NF2 = LOOKUP(9) IF (CTPRGI(NF2).NE.'F' .AND. CTPRGI(NF2).NE.'J' .AND. + CTPRGI(NF2).NE.'D' .AND. CTPRGI(NF2).NE.'G' .AND. + CTPRGI(NF2).NE.'K') THEN WRITE (ISYSW,FMT='(A)') + ' ** Column F2 must be of type F or J or D or G or K **' KE = KE + 1 END IF IF (SIG2.GT.0.0 .AND. LOOKUP(10).EQ.0) + CALL CCPERR(1,' EXCLUDE on SIG2 requested but SIG2 unassigned') END IF C IF (LOOKUP(10).NE.0) THEN NS2 = LOOKUP(10) C C---- SIG column can be of any type, but print warning if not Q, L or M C IF (CTPRGI(NS2).NE.'Q' .AND.CTPRGI(NS2).NE.'L' .AND. + CTPRGI(NS2).NE.'M') THEN WRITE (ISYSW,FMT='(A)') + ' ** WARNING ** column SIG2 is not of type Q, K or M' END IF END IF C IF (LOOKUP(11).NE.0) THEN NH = LOOKUP(11) IF (CTPRGI(NH).NE.'F') THEN WRITE (ISYSW,FMT='(A)') ' ** Column FH must be of type F **' KE = KE + 1 END IF END IF C IF (LOOKUP(12).NE.0) THEN NPH = LOOKUP(12) IF (CTPRGI(NPH).NE.'P') THEN WRITE (ISYSW,FMT='(A)') ' ** Column PHIH must be of type P **' KE = KE + 1 END IF END IF C IF (LOOKUP(13).NE.0) THEN NW = LOOKUP(13) IW = 1 IF (CTPRGI(NW).NE.'W') THEN WRITE (ISYSW,FMT='(A)') ' ** Column W must be of type W **' KE = KE + 1 END IF END IF C C Weight will be multiplied by FREE flag (1 or 0) to exclude reflections. NUM_FREE_COLUMN = 0 IF (LOOKUP(17).NE.0) THEN NUM_FREE_COLUMN = LOOKUP(17) USEFREE = .TRUE. END IF IF (NUM_FREE_COLUMN .eq. 0) + Write(ISYSW,*)' FREE column NOT assigned' C IF (LOOKUP(14).NE.0) THEN NANOM = LOOKUP(14) IF (CTPRGI(NANOM).NE.'D') THEN WRITE (ISYSW,FMT='(A)') ' ** Column DANO must be of type D **' KE = KE + 1 END IF END IF C IF (LOOKUP(15).NE.0) THEN NW2 = LOOKUP(15) IF (CTPRGI(NW2).NE.'W') THEN WRITE (ISYSW,FMT='(A)') ' ** Column W2 must be of type W **' KE = KE + 1 END IF END IF C IF (LOOKUP(16).NE.0) THEN NP2 = LOOKUP(16) IF (CTPRGI(NP2).NE.'P') THEN WRITE (ISYSW,FMT='(A)') ' ** Column PHI2 must be of type P **' KE = KE + 1 END IF END IF C IF (NANOM.NE.0) THEN C C---- Anomalous difference Fourier (DANO) C Set column number for F1 C Flag NANOM used later to retard phase by 90 degrees C NF1 = NANOM LSUSRJ(6) = LSUSRJ(14) END IF C C---- Set flag for previously prepared A & B parts C IF (NNA.GT.0) NF1 = -1 C C---- Check to see if FILLIN is sensible C IF(ABS(ISQ).EQ.1 .OR. NF2.EQ.0) USEFILL = .FALSE. IF (USEFILL) THEN WRITE(ISYSW,'(A,/,A,A,/,A,/,A)') + '*****************************************', + ' Coefficients used for Fourier calculation will have ', + 'missing data replaced:', + ' ie: If F1 is missing, then (SC1 -SC2)* F2 replaces it.', + ' NB - ONLY to be used for (n+1)*F1 - n*F2 type maps', + '*****************************************' IF(ABS(SCLF1-SCLF2-1.0).GT.0.1 .OR. + ABS(SCLF1-SCLF2).LT.0.9) WRITE (ISYSW,'(//,a,/,a,/,a,//)') + '*****************************************', + ' BEWARE - Scales seem strange for a FILLIN procedure ', + '*****************************************' END IF C C--- Print the type of Fourier coefficients to be used C C---- Rare cases first C IF (IPHTR.NE.0) THEN C Set search area whole unit cell XYZLIM(1,1) = 0.0 XYZLIM(2,1) = 1.0 XYZLIM(1,2) = 0.0 XYZLIM(2,2) = 1.0 XYZLIM(1,3) = 0.0 XYZLIM(2,3) = 1.0 LASU = .FALSE. C C---- Phased translation function C WRITE (ISYSW,FMT=6006) + 'Phased translation function: W * F1 * F2 exp (i (PHI - PHI2))' 6006 FORMAT (//' Coefficients used for Fourier calculation:-',/' **', + '*****************************************',//1X,A,/) ELSE IF (NNA.NE.0) THEN WRITE (ISYSW,FMT=6006) 'Fourier coefficients A & B' ELSE IF (NANOM.NE.0) THEN C C---- Anomalous Fourier C WRITE (ISYSW,FMT=6006) + 'Anomalous Fourier: DANO exp (i (PHI - 90))' ELSE IF (NH.NE.0) THEN C C---- Double difference Fourier C WRITE (ISYSW,FMT=6006) LINE = 'Double difference Fourier: ' IF (NW.NE.0) CALL APPNDS(LINE,' W * ') C C ****************************************************** CALL APPNDS(LINE, + ' (k1 * F1 - k2 * | F2exp(i PHI2) + FHexp(i PHIH)|)' + //' exp (i AlphaPH))') C ****************************************************** C WRITE (ISYSW,FMT=6008) LINE(1:LENSTR(LINE)) 6008 FORMAT (1X,A,/) ELSE IF (NP2.NE.0) THEN C C---- Vector difference C WRITE (ISYSW,FMT=6006) LINE = 'Vector difference: ' C IF (NW.NE.0 .AND. NW2.EQ.0) THEN CALL APPNDS(LINE,' W * (') ELSE IF (NW.NE.0) THEN CALL APPNDS(LINE,'( W * ') ELSE CALL APPNDS(LINE,' ( ') END IF CALL APPNDS(LINE,' k1 * F1 exp(i PHI) - ') IF (NW2.NE.0) CALL APPNDS(LINE,'W2 * ') CALL APPNDS(LINE,' k2 * F2 exp(i PHI2))') WRITE (ISYSW,FMT=6008) LINE(1:LENSTR(LINE)) ELSE C C---- Usual options C WRITE (ISYSW,FMT=6006) IF (ISQ.NE.0) THEN LINE = 'Patterson synthesis: ' ELSE LINE = 'Fourier synthesis: ' END IF C IF (NW.NE.0 .AND. NW2.EQ.0) THEN CALL APPNDS(LINE,' W * (') ELSE IF (NW.NE.0) THEN CALL APPNDS(LINE,' (W * ') ELSE CALL APPNDS(LINE,' ( ') END IF C IF(NI1.NE.0)CALL APPNDS(LINE,' k1 * I ') C IF(NF1.NE.0)CALL APPNDS(LINE,' k1 * F1 ') C IF (NF2.NE.0) THEN CALL APPNDS(LINE,' -') IF (NW2.NE.0) CALL APPNDS(LINE,' W2 * ') CALL APPNDS(LINE,' k2 * F2)') ELSE CALL APPNDS(LINE,')') END IF C IF (ISQ.NE.0) THEN CALL APPNDS(LINE,'**2') ELSE CALL APPNDS(LINE,' exp( i PHI)') END IF WRITE (ISYSW,FMT=6008) LINE(1:LENSTR(LINE)) END IF C C---- Echo column labels C IF (NF1.GT.0) THEN LINE = ' F1 = '//LSUSRJ(6) IF (NS1.NE.0) CALL APPNDS(LINE, + ' checked by SIG1 = '//LSUSRJ(7)) WRITE (ISYSW,FMT='(1x,a)') LINE(1:LENSTR(LINE)) END IF C IF (NF2.GT.0) THEN LINE = ' F2 = '//LSUSRJ(9) IF (NS2.NE.0) CALL APPNDS(LINE, + ' checked by SIG2 = '//LSUSRJ(10)) WRITE (ISYSW,FMT='(1x,a)') LINE(1:LENSTR(LINE)) END IF C IF (NP.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' PHI = '//LSUSRJ(8) IF (NW.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' W = '//LSUSRJ(13) IF (NW2.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' W2 = '//LSUSRJ(15) IF (NP2.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' PHI2 = '//LSUSRJ(16) IF (NH.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' FH = '//LSUSRJ(11) IF (NPH.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' PHIH = '//LSUSRJ(12) IF (NANOM.NE.0) WRITE (ISYSW,FMT='(1x,a)') ' DANO = '// + LSUSRJ(14) C WRITE (ISYSW,FMT='(/)') C RETURN C C---- Error in input C C ******************************************************* 1450 CALL CCPERR(1, $ 'A & B parts must either be both given,'// + 'or neither, & may not be combined with other items') C ******************************************************* C END C C VAXPDS NAME=REALFT C********************************************************************** SUBROUTINE REALFT (EVEN, ODD, N, IDIM,NWRITE) C---------------------------------------------------------------------- C REAL FOURIER TRANSFORM C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C---------------------------------------------------------------------- REAL EVEN(*), ODD(*) INTEGER N INTEGER IDIM(10) C------ INTEGER DIM(5) REAL A, B, C, D, E, F, CO, SI LOGICAL ERROR INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW EQUIVALENCE (LSEPU,DIM(1)),(LSEPV,DIM(2)),(LSEPW,DIM(3)) EQUIVALENCE (LGLIMV,DIM(4)),(LGLIMW,DIM(5)) INTEGER I, J, K, L, N OVER 2, I0, LINDU,LINDV,LDXV,LDXW DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBTWON DATA DBPI2/6.28318530717958623D0/ C C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------- C IN/OUT -- EVEN(*) EVEN PART OF DATA/TRANSFORM C IN/OUT -- ODD(*) ODD PART OF DATA/TRANSFORM C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- GRIDIM CALCULATE DIMENSIONING CONSTANTS C CALLS -- CMPLFT COMPLEX FOURIER TRANSFORM C---------------------------------------------------------------------- C NOTE -- INDEXING REVISED AND SIN COS RECALC A.D.MCLACHLAN AUGUST 1981 C NOTE -- C NOTE -- GIVEN A REAL SEQUENCE OF LENGTH 2N THIS SUBROUTINE CALCULATES THE C NOTE -- UNIQUE PART OF THE FOURIER TRANSFORM. THE FOURIER TRANSFORM HAS C NOTE -- N + 1 UNIQUE REAL PARTS AND N - 1 UNIQUE IMAGINARY PARTS. SINCE C NOTE -- THE REAL PART AT X(N) IS FREQUENTLY OF INTEREST, THIS SUBROUTINE C NOTE -- STORES IT AT X(N) RATHER THAN IN Y(0). THEREFORE X AND Y MUST BE C NOTE -- OF LENGTH N + 1 INSTEAD OF N. NOTE THAT THIS STORAGE ARRANGEMENT C NOTE -- IS NOW SAME AS THAT EMPLOYED BY THE HERMITIAN FOURIER TRANSFORM C NOTE -- SUBROUTINE. C NOTE -- C NOTE -- FOR CONVENIENCE THE DATA IS PRESENTED IN TWO PARTS, THE FIRST C NOTE -- CONTAINING THE EVEN NUMBERED REAL TERMS AND THE SECOND CONTAINING C NOTE -- THE ODD NUMBERED TERMS (NUMBERING STARTING AT 0). ON RETURN THE C NOTE -- REAL PART OF THE TRANSFORM REPLACES THE EVEN TERMS AND THE C NOTE -- IMAGINARY PART OF THE TRANSFORM REPLACES THE ODD TERMS. C------------------------------------------------------------------------------ C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C------------------------------------------------------------------------------- DBTWON=DBLE(2*N) ERROR=.FALSE. CALL GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) IF(ERROR) GO TO 700 C----------------------------------------- CALL CMPLFT(EVEN,ODD,N,IDIM,NWRITE) C----------------------------------------- N OVER 2 = N/2 + 1 C-- IF (N OVER 2 .LT. 2) GO TO 400 DBANGL=DBPI2/DBTWON DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 300 I = 2, N OVER 2 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 CO=DBC1 SI=DBS1 I0 = (I - 1)*LSEPU + 1 J = (N + 2 - 2*I)*LSEPU LINDU=I0-LSEPV-LSEPW DO 200 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 200 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + J A = (EVEN(L) + EVEN(K))/2.0 C = (EVEN(L) - EVEN(K))/2.0 B = (ODD(L) + ODD(K))/2.0 D = (ODD(L) - ODD(K))/2.0 E = C*SI + B*CO F = C*CO - B*SI EVEN(K) = A + E EVEN(L) = A - E ODD(K) = F - D ODD(L) = F + D CCC 100 CONTINUE 200 CONTINUE 300 CONTINUE C-- 400 CONTINUE IF (N .LT. 1) GO TO 600 J = N*LSEPU LINDU=1-LSEPV-LSEPW DO 500 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 500 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + J EVEN(L) = EVEN(K) - ODD(K) ODD(L) = 0.0 EVEN(K) = EVEN(K) + ODD(K) ODD(K) = 0.0 500 CONTINUE C-- 600 CONTINUE RETURN C-- 700 CONTINUE WRITE(NWRITE,800) N,IDIM CALL CCPERR (1, 'FFTBIG: internal error') 800 FORMAT(1X,'**REALFT** ERROR IN DIMENSIONS OF ARRAYS N=,IDIM=', 1 I5,5X,3I5,1X,3I5,1X,3I5,1X,I5) END C C VAXPDS NAME=RPCFTK C****************************************************************************** SUBROUTINE RPCFTK (N, M, P, R, X, Y, DIM) C------------------------------------------------------------------------------ C RADIX PRIME MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C------------------------------------------------------------------------------ INTEGER N, M, P, R, DIM(5) REAL X(R,P), Y(R,P) C----- LOGICAL FOLD,ZERO REAL IS,IU,RS,RU,T,XT,YT INTEGER J,JJ,K0,K,M OVER 2,MP,PM,PP,U,V INTEGER KK,MMP INTEGER NS,LDXV,LDXW INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW,LINDU,LINDV DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBFP,DBFMP C-- PARAMETER(MPPMAX=9,MPMAX2=2*MPPMAX) REAL AA(MPPMAX,MPPMAX),BB(MPPMAX,MPPMAX) REAL A(MPMAX2),B(MPMAX2),C(MPMAX2),S(MPMAX2) REAL IA(MPPMAX),IB(MPPMAX),RA(MPPMAX),RB(MPPMAX) DATA DBPI2/6.28318530717958623D0/ C-- C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------------- C INPUT -- N LENGTH OF TRANSFORM C INPUT -- M N DIVIDED BY CURRENT FACTORS C INPUT -- P PRIME NUMBER C INPUT -- R SEPARATION BETWEEN ADJACENT ELEMENTS OF X,Y C IN/OUT -- X(R,P) REAL VALUES AT P EQUALLY SEPARATED DATA POINTS C IN/OUT -- Y(R,P) IMAGINARY VALUES AT P EQUALLY SEPARATED DATA POINTS C INPUT -- DIM(5) ARRAY DIMENSIONING CONSTANTS C----------------------------------------------------------------------------- C CALLS -- *** C------------------------------------------------------------------------------ C NOTE -- DIMENSIONING REVISED AND SIN COS RECALC A.D. MCLACHLAN AUGUST 1981 C NOTE -- THE LARGEST PRIME ALLOWED IS 2*MPPMAX+1 (HERE 19) C------------------------------------------------------------------------------ LSEPU=DIM(1) LSEPV=DIM(2) LSEPW=DIM(3) LGLIMV=DIM(4) LGLIMW=DIM(5) NS = N*LSEPU M OVER 2=M/2+1 MP=M*P DBFMP=DBLE(MP) MMP = LSEPU*MP PP=P/2 PM=P-1 DBFP=DBLE(P) DBANGL=DBPI2/DBFP DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 100 U=1,PP DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 JJ=P-U A(U)=DBC1 B(U)=DBS1 A(JJ)=A(U) B(JJ)=-B(U) 100 CONTINUE DO 300 U=1,PP DO 200 V=1,PP JJ=U*V-U*V/P*P AA(V,U)=A(JJ) BB(V,U)=B(JJ) 200 CONTINUE 300 CONTINUE C-- DBANGL=DBPI2/DBFMP DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 1500 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*LSEPU + 1 ZERO= J.EQ.1 IF (ZERO) GO TO 700 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C(1)=DBC1 S(1)=DBS1 DO 400 U=2,PM C(U)=C(U-1)*C(1)-S(U-1)*S(1) S(U)=S(U-1)*C(1)+C(U-1)*S(1) 400 CONTINUE GO TO 700 500 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*LSEPU + 1 DO 600 U=1,PM T=C(U)*A(U)+S(U)*B(U) S(U)=-S(U)*A(U)+C(U)*B(U) C(U)=T 600 CONTINUE 700 CONTINUE C-- DO 1400 KK = K0, NS, MMP LINDU=KK-LSEPV-LSEPW DO 1340 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 1320 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW XT=X(K,1) YT=Y(K,1) RS=X(K,2)+X(K,P) IS=Y(K,2)+Y(K,P) RU=X(K,2)-X(K,P) IU=Y(K,2)-Y(K,P) DO 800 U=1,PP RA(U)=XT+RS*AA(U,1) IA(U)=YT+IS*AA(U,1) RB(U)=RU*BB(U,1) IB(U)=IU*BB(U,1) 800 CONTINUE XT=XT+RS YT=YT+IS DO 1000 U=2,PP JJ=P-U RS=X(K,U+1)+X(K,JJ+1) IS=Y(K,U+1)+Y(K,JJ+1) RU=X(K,U+1)-X(K,JJ+1) IU=Y(K,U+1)-Y(K,JJ+1) XT=XT+RS YT=YT+IS DO 900 V=1,PP RA(V)=RA(V)+RS*AA(V,U) IA(V)=IA(V)+IS*AA(V,U) RB(V)=RB(V)+RU*BB(V,U) IB(V)=IB(V)+IU*BB(V,U) 900 CONTINUE 1000 CONTINUE X(K,1)=XT Y(K,1)=YT DO 1300 U=1,PP JJ=P-U IF (ZERO) GO TO 1100 XT=RA(U)+IB(U) YT=IA(U)-RB(U) X(K,U+1)=XT*C(U)+YT*S(U) Y(K,U+1)=YT*C(U)-XT*S(U) XT=RA(U)-IB(U) YT=IA(U)+RB(U) X(K,JJ+1)=XT*C(JJ)+YT*S(JJ) Y(K,JJ+1)=YT*C(JJ)-XT*S(JJ) GO TO 1200 1100 CONTINUE X(K,U+1)=RA(U)+IB(U) Y(K,U+1)=IA(U)-RB(U) X(K,JJ+1)=RA(U)-IB(U) Y(K,JJ+1)=IA(U)+RB(U) 1200 CONTINUE 1300 CONTINUE 1320 CONTINUE 1340 CONTINUE 1400 CONTINUE IF (FOLD) GO TO 500 1500 CONTINUE C RETURN END C C VAXPDS NAME=RSYMFT C****************************************************************************** SUBROUTINE RSYMFT(X,N,IDIM,NWRITE) C------------------------------------------------------------------------------ C REAL SYMMETRIC MULTIDIMENSIONAL FOURIER TRANSFORM C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C------------------------------------------------------------------------------ REAL X(*) INTEGER N INTEGER IDIM(10) C----- INTEGER DIM(5) LOGICAL ERROR REAL A, B, CO, SI EQUIVALENCE (LSEPU,DIM(1)),(LSEPV,DIM(2)),(LSEPW,DIM(3)) EQUIVALENCE (LGLIMV,DIM(4)),(LGLIMW,DIM(5)) INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW INTEGER LINDU,LINDV,LDXV,LDXW INTEGER I, J0, J1, K, K0, L, NN,NWRITE,C,D INTEGER N OVER 2, N OVER 4 INTEGER II, I0, J, M, MJ, MK, ML, MM, TWOD2 DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBTWON DATA DBPI2/6.28318530717958623D0/ C C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------- C IN/OUT -- X(*) REAL DATA/TRANSFORM C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- GRIDIM CALCULATE DIMENSIONING CONSTANTS C CALLS -- HERMFT HERMITIAN FOURIER TRANSFORM C----------------------------------------------------------------------- C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C------------------------------------------------------------------------------- C NOTE -- INDEXING REVISED AND SIN COS RECALCULATED A.D.MCLACHLAN SEP 1981 C NOTE -- C NOTE -- N MUST BE A MULTIPLE OF 4. THE TWO UNIQUE ELEMENTS ARE STORED AT C NOTE -- X(1) AND X(N+1). C NOTE -- C NOTE -- DECIMATION IN FREQUENCY APPLIED TO A REAL SYMMETRIC SEQUENCE OF C NOTE -- LENGTH 2N GIVES A REAL SYMMETRIC SEQUENCE OF LENGTH N, THE C NOTE -- TRANSFORM OF WHICH GIVES THE EVEN NUMBERED FOURIER COEFFICIENTS, C NOTE -- AND A HERMITIAN SYMMETRIC SEQUENCE OF LENGTH N, THE TRANSFORM OF C NOTE -- WHICH GIVES THE ODD NUMBERED FOURIER COEFFICIENTS. THE SUM OF C NOTE -- THE TWO SEQUENCES IS A HERMITIAN SYMMETRIC SEQUENCE OF LENGTH N, C NOTE -- WHICH MAY BE STORED IN N/2 COMPLEX LOCATIONS. THE TRANSFORM OF C NOTE -- THIS SEQUENCE IS N REAL NUMBERS REPRESENTING THE TERM BY TERM SUM C NOTE -- OF THE EVEN AND ODD NUMBERED FOURIER COEFFICIENTS. THIS SYMMETRIC C NOTE -- SEQUENCE MAY BE SOLVED IF ANY OF THE FOURIER COEFFICIENTS ARE C NOTE -- KNOWN. FOR THIS PURPOSE X0, WHICH IS SIMPLY THE SUM OF THE C NOTE -- ORIGINAL SEQUENCE, IS COMPUTED AND SAVED IN X(N+1). C------------------------------------------------------------------------------- IF (N .EQ. 1) GO TO 1300 N OVER 2 = N/2 N OVER 4 = N/4 IF (4*N OVER 4 .NE. N) GO TO 1400 ERROR=.FALSE. CALL GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) IF(ERROR) GO TO 1600 DBTWON=DBLE(2*N) TWOD2 = 2*LSEPU C-- K0 = N*LSEPU + 1 LINDU=K0-LSEPV-LSEPW DO 100 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 100 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW X(K) = X(K)/2.0 100 CONTINUE C-- DBANGL=DBPI2/DBTWON DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 300 I = 2, N OVER 2 DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 CO=DBC1 SI=DBS1 K0 = (I - 1)*LSEPU + 1 J0 = (N + 2 - 2*I)*LSEPU J1 = (N + 1 - I)*LSEPU LINDU=K0-LSEPV-LSEPW DO 200 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 200 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + J0 NN = K + J1 A = X(L) + X(K) B = X(L) - X(K) X(K) = A - B*CO X(L) = B*SI X(NN) = X(NN) + A 200 CONTINUE 300 CONTINUE C-- IF (N OVER 4 .EQ. 1) GO TO 600 J0 = N OVER 4 - 1 DO 500 I = 1, J0 K0 = (N OVER 2 + I)*LSEPU + 1 J1 = (N OVER 2 - 2*I)*LSEPU LINDU=K0-LSEPV-LSEPW DO 400 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 400 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + J1 A = X(K) X(K) = X(L) X(L) = A 400 CONTINUE 500 CONTINUE C-- 600 CONTINUE J0 = N OVER 2*LSEPU J1 = N*LSEPU LINDU=1-LSEPV-LSEPW DO 700 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 700 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW I = K + J0 L = K + J1 X(I) = 2.0*X(I) X(L) = X(K) + X(I) + 2.0*X(L) X(K) = 2.0*X(K) 700 CONTINUE C-- K = NOVER2*LSEPU + 1 C---------------------------------------------- CALL HERMFT(X(1),X(K),NOVER2,IDIM,NWRITE) C---------------------------------------------- C-- SOLVE THE EQUATIONS FOR ALL OF THE SEQUENCES C-- I0 = 1 - LSEPU MK = N OVER 2*LSEPU MJ = MK + LSEPU ML = N*LSEPU + LSEPU MM = ML DO 800 II = 1, N OVER 4 I0 = I0 + LSEPU MJ = MJ - TWOD2 ML = ML - TWOD2 MM = MM - LSEPU LINDU=I0-LSEPV-LSEPW DO 800 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 800 LDXW=LSEPW,LGLIMW,LSEPW I=LINDV+LDXW J = I + MJ K = I + MK L = I + ML M = I + MM A = X(I) - X(M) B = X(L) - A C = X(K) - B D = X(J) - C X(I) = X(M) X(J) = A X(K) = B X(L) = C X(M) = D 800 CONTINUE C------------------------------------------------------------------------ C THE RESULTS ARE NOW IN A SCRAMBLED DIGIT REVERSED ORDER, I.E. C X(1), X(5), X(9), ..., X(10), X(6), X(2), ..., X(3), X(7), X(11), C ..., X(12), X(8), X(4). THE FOLLOWING SECTION OF PROGRAM FOLLOWS C THE PERMUTATION CYCLES AND DOES THE NECESSARY INTERCHANGES. C------------------------------------------------------------------------ IF (N OVER 4 .EQ. 1) GO TO 1300 NN = N - 2 DO 1200 I = 1, NN K = I C-- 1000 CONTINUE K0 = K/4 L = K - K0*4 IF (L .NE. (L/2)*2) K0 = N OVER 4 - 1 - K0 K = K0 + L*N OVER 4 IF (K .LT. I) GO TO 1000 IF (K .EQ. I) GO TO 1200 C-- K0 = I*LSEPU + 1 J0 = (K - I)*LSEPU LINDU=K0-LSEPV-LSEPW DO 1100 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 1100 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + J0 A = X(K) X(K) = X(L) X(L) = A 1100 CONTINUE 1200 CONTINUE C-- 1300 CONTINUE RETURN C 1400 CONTINUE WRITE(NWRITE,1500) N CALL CCPERR (1, 'FFTBIG: internal error') C-- 1500 FORMAT(1X,' **RSYMFT** ERROR: N NOT A MULTIPLE OF 4. N =',I10) C-- 1600 CONTINUE WRITE(NWRITE,1700) N,IDIM CALL CCPERR (1, 'FFTBIG: internal error') 1700 FORMAT(1X,'**RSYMFT** ERROR IN DIMENSIONS OF ARRAYS N=,IDIM=', 1 I5,5X,3I5,1X,3I5,1X,3I5,1X,I5) END C C VAXPDS NAME=SDIAD C*************************************************************************** SUBROUTINE SDIAD(X,Y,N,IDIM,NWRITE) C--------------------------------------------------------------------------- C FOURIER SYNTHESIS ALONG A SCREW DIAD C LAST UPDATED A.D. MCLACHLAN 26 AUG 1987 C--------------------------------------------------------------------------- REAL X(*), Y(*) INTEGER N,NWRITE INTEGER IDIM(10) C----- INTEGER DIM(5) LOGICAL ERROR INTEGER LSEPU,LSEPV,LSEPW,LGLIMV,LGLIMW EQUIVALENCE(LSEPU,DIM(1)),(LSEPV,DIM(2)),(LSEPW,DIM(3)) EQUIVALENCE(LGLIMV,DIM(4)),(LGLIMW,DIM(5)) INTEGER LINDU,LINDV,LDXV,LDXW,MN LOGICAL FOLD REAL A, C, S INTEGER I, J, K, K0, L, M, NN, N OVER 2 DOUBLE PRECISION DBPI2,DBCOS1,DBSIN1,DBC1,DBS1,DBCC DOUBLE PRECISION DBANGL,DBTWON DATA DBPI2/6.28318530717958623D0/ C-- C%% DBPI2=8.0D0*DATAN2(1.0D0,1.0D0) C----------------------------------------------------------------------- C IN/OUT -- X(*) REAL PART OF DATA/TRANSFORM C IN/OUT -- Y(*) IMAGINARY PART OF DATA/TRANSFORM C INPUT -- N LENGTH OF TRANSFORM AXIS C INPUT -- IDIM(10) DIMENSIONING CONSTANTS FOR 1,2 OR 3-D ARRAY C INPUT -- NWRITE PRINT UNIT C----------------------------------------------------------------------- C CALLS -- GRIDIM CALCULATE DIMENSIONING CONSTANTS C CALLS -- CMPLFT COMPLEX FOURIER TRANSFORM C----------------------------------------------------------------------- C NOTE -- INDEXING -- THE ARRANGEMENT OF THE MULTI-DIMENSIONAL DATA IS C NOTE -- SPECIFIED BY THE INTEGER ARRAY D, THE VALUES OF WHICH ARE USED AS C NOTE -- CONTROL PARAMETERS IN DO LOOPS. WHEN IT IS DESIRED TO COVER ALL C NOTE -- ELEMENTS OF THE DATA FOR WHICH THE SUBSCRIPT BEING TRANSFORMED HAS C NOTE -- THE VALUE I0, THE FOLLOWING IS USED. C NOTE -- C NOTE -- I1 = (I0 - 1)*LSEPU + 1 C NOTE -- LINDU=I1-LSEPV-LSEPW C NOTE -- DO 100 LDXV=LSEPV,LGLIMV,LSEPV C NOTE -- LINDV=LINDU+LDXV C NOTE -- DO 100 LDXW=LSEPW,LGLIMW,LSEPW C NOTE -- I=LINDV+LDXW C NOTE -- . . . C NOTE -- . . . C NOTE -- 100 CONTINUE C NOTE -- HERE LSEPU=D(1) IS SEPARATION BETWEEN ADJACENT ELEMENTS OF X C NOTE -- (OR OF Y) ALONG THE INDEX IU WHICH IS BEING TRANSFORMED. C NOTE -- LSEPV=D(2) AND LSEPW=D(3) ARE SEPARATION OF ADJACENT ELEMENTS C NOTE -- ALONG THE SECOND AND THIRD DIRECTIONS INDEXED BY IV,IW. C NOTE -- LGLIMV=D(4)=LSEPV*LGRIDV AND LGLIMW=D(5)=LSEPW*LGRIDW WITH C NOTE -- LGRIDV, LGRIDW BEING THE NUMBER OF GRID POINTS IN THE X (OR Y) C NOTE -- ARRAYS ALONG THE IV AND IW AXES C NOTE -- WITH THIS INDEXING IT IS POSSIBLE TO USE A NUMBER OF ARRANGEMENTS C NOTE -- OF THE DATA, INCLUDING NORMAL FORTRAN COMPLEX NUMBERS C NOTE -- (LSEPU=D(1) = 2) C NOTE -- THE SUBROUTINE GRIDIM CALCULATES THE ELEMENTS OF D ARRAY FROM THE C NOTE -- IDIM ARRAY, WHERE C NOTE -- IDIM(1)=IU AXIS ALONG WHICH TRANSFORM IS DONE C NOTE -- IDIM(2)=IV 2ND AXIS NORMAL TO IU C NOTE -- IDIM(3)=IW 3RD AXIS NORMAL TO IU AND IV C NOTE -- IDIM(4)=NDIMX ARRAY DIMENSION OF X(OR Y) ALONG A AXIS C NOTE -- IDIM(5)=NDIMY ARRAY DIMENSION ALONG B AXIS C NOTE -- IDIM(6)=NDIMZ ARRAY DIMENSION ALONG C AXIS C NOTE -- IDIM(7)=NGRIDX NUMBER OF GRID POINTS USED ALONG A AXIS C NOTE -- IDIM(8)=NGRIDY NUMBER OF GRID POINTS USED ALONG B AXIS C NOTE -- IDIM(9)=NGRIDZ NUMBER OF GRID POINTS USED ALONG C AXIS C NOTE -- IDIM(10)=ICMPLX VALUE 1 FOR SEPARATE X,Y ARRAYS,2 FOR C NOTE -- INTERLEAVED ARRAYS(X,I*Y) C NOTE -- NOTE THAT IU,IV,IW MUST BE A PERMUTATION OF 1,2,3 AND EACH C NOTE -- NGRID.LE.NDIM. C NOTE -- ALSO NGRID FOR AXIS IU IS SAME AS N, LENGTH OF TRANSFORM. C NOTE -- THIS ALLOWS TRANSFORM TO BE DONE FOR ARRAYS EMBEDDED IN PART OF A C NOTE -- LARGER ARRAY AND CHOICE OF AXES IN ANY ORDER C--------------------------------------------------------------------------- C NOTE -- DIMENSIONING REVISED AND SIN COS RECALCULATED C NOTE -- A.D.MCLACHLAN AUG 1981 C NOTE -- MINOR CORRECTIONS MARCH 1981 BY A.D. MCLACHLAN TO REMOVE LIMITATION C NOTE -- THAT F(N) IS ASSUMED TO BE ZERO IN THE ORIGINAL PROGRAM. C NOTE -- USE IMAGINARY PART OF F(N) -REALLY ZERO- AS A SCRATCH AREA C NOTE -- OLD STATEMENTS RETAINED WITH 'C%%' IN FRONT C NOTE -- THIS SUBROUTINE COMPUTES HALF THE FOURIER SYNTHESIS ALONG A SCREW C NOTE -- DIAD LYING ALONG A CRYSTALLOGRAPHIC AXIS GIVEN HALF THE FOURIER C NOTE -- COEFFICIENTS. THAT IS, IT ASSUMES THAT F(T) = CONJG(F(-T)) FOR T C NOTE -- EVEN AND F(T) = -CONJG(F(-T)) FOR T ODD. N IS THE LENGTH OF THE C NOTE -- DESIRED HALF OF THE TRANSFORM. THE LOCATION X(N+1) IS REQUIRED AS C NOTE -- A SCRATCH LOCATION AND THEREFORE A VALUE IS ALSO RETURNED IN C NOTE -- X(N+1) AND Y(N+1). THE VALUE OF THE SECOND HALF OF THE TRANSFORM C NOTE -- MAY BE GENERATED FROM THE FIRST HALF BY THE FORMULA X(N+T) = X(T), C NOTE -- Y(N+T) = -Y(T). IN OTHER WORDS, THE LAST HALF OF THE TRANSFORM IS C NOTE -- THE COMPLEX CONJUGATE OF THE FIRST HALF. C NOTE -- C NOTE -- THE TRANSFORM IS CALCULATED BY FORMING THE SUM OF THE EVEN TERMS C NOTE -- AND THE ODD TERMS IN PLACE, USING THE SYMMETRY RELATIONS TO C NOTE -- OBTAIN THE VALUES FOR NEGATIVE SUBSCRIPTS. THE TRANSFORM OF THE C NOTE -- RESULTING SEQUENCE MAY BE SEPARATED BY USING THE FACT THAT THE C NOTE -- TRANSFORM OF THE EVEN TERMS IS REAL, WHILE THE PRODCT OF THE C NOTE -- TRANSFORM OF THE ODD TERMS AND (COS(PI*T/N) - I*SIN(PI*T/N)) IS C NOTE -- IMAGINARY. THE SCRATCH LOCATION IS REQUIRED BECAUSE THE FORMULA C NOTE -- FOR SEPARATING THE TWO TRANSFORMS BREAKS DOWN WHEN T = N/2. C------------------------------------------------------------------------------- NOVER2 = N/2 IF (2*NOVER2 .NE. N) GO TO 700 DBTWON=DBLE(2*N) ERROR=.FALSE. CALL GRIDIM(N,IDIM,ERROR,LSEPU,LSEPV,LSEPW, 1 LGLIMV,LGLIMW,NWRITE) IF(ERROR) GO TO 1100 C-- C *** NEW SIGN CORRECTION OF ERROR IN ORIGINAL ROUTINE C-- STORE +/- X(N) IN Y(N+1) S= -1.0 IF(NOVER2.EQ.(2*(NOVER2/2))) S= -S K0 = (N - 1)*LSEPU + 1 LINDU=K0-LSEPV-LSEPW DO 100 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 100 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + LSEPU Y(L) = S*X(K) C%% X(L) = X(K) 100 CONTINUE S = 1.0 NN = N - 2 DO 200 I = 1, NN, 2 S = -S MN = (N + 1 - I)*LSEPU K0 = (I - 1)*LSEPU + 1 LINDU=K0-LSEPV-LSEPW DO 150 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 150 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW J = K + LSEPU L = K + 2*LSEPU M = K + MN C-- UPDATE TOTAL IN SCRATCH LOCATION Y(M) =Y(M) + S*X(J) C%% X(M) = X(M) + S*X(J) X(K) = X(K) + X(J) X(J) = X(L) - X(J) Y(K) = Y(K) + Y(J) Y(J) = Y(J) - Y(L) 150 CONTINUE 200 CONTINUE K0 = (N - 2)*LSEPU + 1 LINDU=K0-LSEPV-LSEPW DO 250 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 250 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + LSEPU C-- NEXT STATEMENT IS NEW *** J=L+LSEPU X(K) = X(K) + X(L) Y(K) = Y(K) + Y(L) C-- USE THE REAL PART OF X(N+1) TO FORM THE (N)TH FOURIER INPUT X(L) = X(J) - X(L) C%% X(L) = -X(L) 250 CONTINUE C --REORDER SCRAMBLED FOURIER COEFFICIENTS DO 400 I = 1, NN K = I 300 CONTINUE K = 2*K IF (K .GT. N - 1) K = 2*N - 1 - K IF (K .LT. I) GO TO 300 IF (K .EQ. I) GO TO 400 J = (K - I)*LSEPU K0 = I*LSEPU + 1 LINDU=K0-LSEPV-LSEPW DO 350 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 350 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW L = K + J A = X(K) X(K) = X(L) X(L) = A A = Y(K) Y(K) = Y(L) Y(L) = A 350 CONTINUE 400 CONTINUE C---------------------------------------- CALL CMPLFT(X,Y,N,IDIM,NWRITE) C---------------------------------------- M = NOVER2 - 1 DBANGL=DBPI2/DBTWON DBCOS1=DCOS(DBANGL) DBSIN1=DSIN(DBANGL) DBC1=1.0D0 DBS1=0.0D0 DO 600 I = 1, M DBCC=DBC1 DBC1=DBC1*DBCOS1-DBS1*DBSIN1 DBS1=DBS1*DBCOS1+DBCC*DBSIN1 C=DBC1 S=DBS1 K0 = I*LSEPU + 1 FOLD = .TRUE. GO TO 500 C-- 450 CONTINUE C = -C K0 = (N - I)*LSEPU + 1 FOLD = .FALSE. C-- 500 CONTINUE LINDU=K0-LSEPV-LSEPW DO 550 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 550 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW A = Y(K)/C X(K) = X(K) + S*A Y(K) = A 550 CONTINUE IF (FOLD) GO TO 450 600 CONTINUE C-- M = NOVER2*LSEPU K0 = M + 1 C --SPECIAL FOURIER COEFFICIENT IS X(N/2),Y(N/2) AND USES SCRATCH SUM 'A' C --WHICH IS STORED IN Y(N+1) LINDU=K0-LSEPV-LSEPW DO 650 LDXV=LSEPV,LGLIMV,LSEPV LINDV=LINDU+LDXV DO 650 LDXW=LSEPW,LGLIMW,LSEPW K=LINDV+LDXW J = K - M L = K + M A = 2.0*Y(L) C%% A = 2.0*X(L) X(K) = X(K) + A Y(K) = A C --COPY X(1),Y(1) INTO X(N+1),Y(N+1) X(L) = X(J) Y(L) = -Y(J) 650 CONTINUE C-- RETURN C-- 700 CONTINUE WRITE(NWRITE,1000) N CALL CCPERR (1, 'FFTBIG: internal error') C-- 1000 FORMAT (1X,'**SDIAD** ERROR: N ODD. N =', I10) C-- 1100 CONTINUE WRITE(NWRITE,1200) N,IDIM 1200 FORMAT(1X,'**SDIAD** ERROR IN ARRAY DIMENSIONS N=,IDIM=', 1 I5,5X,3I5,1X,3I5,1X,3I5,1X,I5) CALL CCPERR (1, 'FFTBIG: internal error') END C SUBROUTINE SFIPYB(X,N1,N2,N3,MM) C ================================ C C Input routine for complex data for in-core FFT C The data is stored in the array X as lhk C Note that l .ge. 0 C N1 = NZ/2+1 C N2 = NX C N3 = NY C MM is return = -1 if error, else = 0 C C IMPLICIT NONE C C .. Scalar Arguments .. INTEGER MM,N1,N2,N3 C .. C .. Array Arguments .. CCC COMPLEX X(N1,N2,N3) REAL X(2*N1,N2,N3) C .. C .. Scalars in Common .. COMPLEX F INTEGER H,HMAX,HMIN,ISORT,K,KMAX,KMIN,L, + LMAX,LMIN,M REAL A,B,SSQ C .. C .. Local Scalars .. INTEGER IH,IK,IL,TOTAL,USED C .. C .. External Subroutines .. EXTERNAL GETHKL,GGETF C .. C .. Intrinsic Functions .. INTRINSIC CMPLX,CONJG C .. C .. Common blocks .. COMMON /HKLF/A,B,ISORT,H,K,L,M,SSQ COMMON /HKLLIM/HMAX,KMAX,LMAX,HMIN,KMIN,LMIN C .. SAVE /HKLF/, /HKLLIM/ C .. C DO 30 IK = 1,N3 DO 20 IH = 1,N2 DO 10 IL = 1,2*N1 X(IL,IH,IK) = 0.0 10 CONTINUE 20 CONTINUE 30 CONTINUE C ISYSW = 6 TOTAL = 0 USED = 0 MM = 1 C 60 CONTINUE C C ****** CALL GETHKL C ****** C Check end of file IF (M.LT.0) GO TO 100 C TOTAL = TOTAL + M C Check index out of range IF (H.GT.HMAX .OR. K.GT.KMAX .OR. L.GT.LMAX .OR. $ H.LT.HMIN .OR. K.LT.KMIN .OR. L.LT.LMIN) THEN MM = -1 RETURN ENDIF C C ***** CALL GGETF F = CMPLX(A,B) C ***** C IF (M.LT.0) GO TO 60 C USED = USED + 1 IH = H + 1 IF (IH.LE.0) IH = IH + N2 IK = K + 1 IF (IK.LE.0) IK = IK + N3 L = L + 1 CCC X(L,IH,IK) = F X(2*L-1,IH,IK) = A X(2*L,IH,IK) = B C C---- For space group P1, fill in rest of hk0 plane C IF (L.EQ.1) THEN C For hk0 reflections, we should receive h.ge.0, with k .ge. 0 if h=0 IF (H.EQ.0) THEN C . . if h=0, generate, 0 -k 0 CCC X(1,IH,N3+1-K) = CONJG(F) X(1,IH,N3+1-K) = A X(2,IH,N3+1-K) = -B ELSE C . . if h .ne. 0 generate -h -k 0 IF (K.GT.0) THEN CCC X(1,N2+1-H,N3+1-K) = CONJG(F) X(1,N2+1-H,N3+1-K) = A X(2,N2+1-H,N3+1-K) = -B ELSE CCC X(1,N2+1-H,1-K) = CONJG(F) X(1,N2+1-H,1-K) = A X(2,N2+1-H,1-K) = -B ENDIF ENDIF ENDIF GO TO 60 C 100 WRITE (ISYSW,FMT=6004) USED,TOTAL TOTAL = 0 MM = 0 RETURN C---- Format statements C 6004 FORMAT(1X,I10,' Structure factors used out of',I6, $ ' on the input file') C END C C VAXPDS NAME=SHRNK2 C************************************************************ SUBROUTINE SHRNK2(X,NX,NY,NZ) C-------------------------------------------------- C REDUCE DIMENSION OF A 3-D ARRAY FROM (NX+2)*NY*NZ TO (NX)*NY*NZ C ASSUMING LOCATIONS IX=NX+1 AND NX+2 ARE EMPTY C USED FOR TESTING FFT ROUTINES. A.D. MCLACHLAN. LAST UPDATED 27 AUG 1987 C------------------------------------------------------------- DIMENSION X(*) INTEGER NX,NY,NZ,NX2,NYZ,NYZ1,NSTEP,NT1,JS,JS2,I, + NTOT1,NTOT2,JCOL REAL X C------------------------------------------------------------- C ENTRY -- SHRNK2 SHRINK C ENTRY -- EXPND2 EXPAND C------------------------------------------------------------- C IN/OUT -- X(*) ARRAY C INPUT -- NX X-DIMENSION C INPUT -- NY Y-DIMENSION C INPUT -- NZ Z-DIMENSION C------------------------------------------------------------- C CALLS -- *** C------------------------------------------------------------- NX2=NX+2 NYZ=NY*NZ IF(NYZ.EQ.1) RETURN NYZ1=NYZ-1 NSTEP=2 NT1=NYZ1*NX DO 110 JS=NX,NT1,NX JS2=JS+NSTEP NSTEP=NSTEP+2 DO 100 I=1,NX X(JS+I)=X(JS2+I) 100 CONTINUE 110 CONTINUE C --CLEAR THE UNUSED LOCATIONS NTOT1=NX*NYZ +1 NTOT2=NX2*NYZ DO 120 I=NTOT1,NTOT2 X(I)=0.0 120 CONTINUE RETURN C+++++++++++++++++++++++++++++++++ ENTRY EXPND2(X,NX,NY,NZ) C+++++++++++++++++++++++++++++++++ C --REVERSE OPERATION EXPAND NX2=NX+2 NYZ=NY*NZ NTOT1=NX*NYZ+1 NTOT2=NX2*NYZ IF(NYZ.EQ.1) RETURN NYZ1=NYZ-1 JS=NTOT1 JS2=NTOT2-1 DO 210 JCOL=1,NYZ1 C --IN REVERSE ORDER TO AVOID OVERWRITING DO 200 I=1,NX X(JS2-I)=X(JS-I) 200 CONTINUE JS=JS-NX JS2=JS2-NX2 210 CONTINUE C --CLEAR THE IX=NX+1 AND IX=NX+2 POSITIONS OF THE NEW ARRAY NTOT2=NX2*NYZ DO 220 I=1,NTOT2,NX2 J=I+NX X(J)=0.0 X(J+1)=0.0 220 CONTINUE RETURN END C C VAXPDS NAME=SRFP C*************************************************************************** SUBROUTINE SRFP(NPTS,NPMAX,NPORD,N2GRP,NFACTR,NSYM, 1 NPSYM,NUNSYM,NPRLST,ERROR,NWRITE) C--------------------------------------------------------------------------- C SYMMETRIZED REORDERING FACTORING PROGRAMME C LAST UPDATED A.D. MCLACHLAN 27 AUG 1987 C--------------------------------------------------------------------------- INTEGER MPMAX,MPMAX1,MPORD PARAMETER(MPMAX=32,MPMAX1=MPMAX+1,MPORD=8) LOGICAL ERROR INTEGER NPTS,NPMAX,N2GRP,NPSYM INTEGER NFACTR(MPMAX1), NSYM(MPMAX1), NUNSYM(MPMAX1) INTEGER NPRLST(MPORD) C----- INTEGER NPP(MPMAX1), NQQ(MPMAX1),NPORD,NWRITE INTEGER NF,J,JJ,N,NFCMAX,NP,NPOW2,NQ,NR C--------------------------------------------------------------------------- C INPUT -- NPTS NUMBER OF POINTS C INPUT -- NPMAX LARGEST ALLOWED FACTOR C INPUT -- NPORD NUMBER OF PRIMES USED C INPUT -- N2GRP HIGHEST VALUE OF 2**N TREATED AS A SINGLE C INPUT -- SPECIAL FACTOR C OUTPUT -- NFACTR(MPMAX1) DOUBLE LIST OF FACTORS C OUTPUT -- NSYM(MPMAX1) LIST OF PAIRED FACTORS TWICE (DOWN THEN UP) C OUTPUT -- NPSYM PRODUCT OF PAIRED FACTORS ONCE EACH C OUTPUT -- NUNSYM(MPMAX1) LIST OF SINGLE FACTORS, SMALLEST FIRST C INPUT -- NPRLST(MPORD) LIST OF FIRST NPORD PRIMES /2,3,5,7,.../ C OUTPUT -- ERROR .L. .TRUE. IF ERROR IS FOUND C INPUT -- NWRITE PRINT CONTROL C--------------------------------------------------------------------------- C NOTE -- UP TO 32 FACTORS DEC 1984 A.D. MCLACHLAN. UPDATED 17 DEC 1984. C--------------------------------------------------------------------------- C NOTE -- THE PARAMETER MPMAX (HERE 32) SPECIFIES THE LARGEST NUMBER OF C NOTE -- FACTORS (POWERS OF 2 ARE CLUSTERED AS 4 AND 8) C NOTE -- THERE ARE NP PAIRED FACTORS, NQ SINGLE FACTORS C NOTE -- NQQ HAS DIMENSION MPMAX1 (UP TO NPORD+1 USED) C NOTE -- NPP HAS DIMENSION MPMAX1 (UP TO MPMAX1 USED) C-------------------------------------------------------------------------- 20 FORMAT (1X,'***SRFP*** ERROR:', 1 ' LARGEST FACTOR EXCEEDS ',I3,' N = ',I6) 21 FORMAT (1X,'***SRFP*** ERROR:', 1 ' FACTOR COUNT EXCEEDS ',I3,' N = ',I6) C-------------------------------------------------------------------------- NFCMAX=MPMAX N=NPTS NPSYM=1 IPMIN=1 NP=0 NQ=0 100 CONTINUE C --DIVIDE BY POSSIBLE PRIME FACTORS IN TURN IF (N.LE.1) GO TO 500 DO 200 IP=IPMIN,NPORD J=NPRLST(IP) IF (N.EQ.(N/J)*J) GO TO 300 200 CONTINUE C --FACTOR TOO BIG GO TO 1400 300 CONTINUE C --CHECK FOR TOO MANY FACTORS IF ((2*NP+NQ).GE.NFCMAX) GO TO 1600 IPMIN=IP NF=J N=N/NF C --CHECK IF STILL DIVISIBLE BY NF IF (N.EQ.(N/NF)*NF) GO TO 400 C --NOW NF IS A SINGLE FACTOR OF NPTS (UP TO NPORD DIFERENT ONES MAY EXIST) NQ=NQ+1 NQQ(NQ)=NF GO TO 100 400 CONTINUE C --NOW NF IS A PAIRED FACTOR C --DIVIDE BY NF A SECOND TIME TO COLLECT PAIRS OF FACTORS C --NPSYM IS THE PRODUCT OF ALL THE PAIRED FACTORS, ONCE EACH N=N/NF NP=NP+1 NPP(NP)=NF NPSYM=NPSYM*NF GO TO 100 C------------------- 500 CONTINUE C --RESERVE SPACE FOR A MIDDLE VALUE IN NSYM ARRAY NR=1 IF (NQ.EQ.0) NR=0 IF (NP.LT.1) GO TO 700 C --INDEX ALL THE DOUBLE FACTORS C --NSYM(J) LISTS LARGEST FIRST (1...NP) C --NFACTR(J) LISTS LARGEST FIRST (1...NP) C --NFACTR(NP+NQ+J) LISTS AGAIN SMALLEST FIRST C --NSYM(NP+NR+J) LISTS AGAIN SMALLEST FIRST DO 600 J=1,NP JJ=NP+1-J NSYM(J)=NPP(JJ) NFACTR(J)=NPP(JJ) JJ=NP+NQ+J NFACTR(JJ)=NPP(J) JJ=NP+NR+J NSYM(JJ)=NPP(J) 600 CONTINUE 700 CONTINUE IF (NQ.LT.1) GO TO 900 C --INDEX ALL THE SINGLE FACTORS (1...NQ) IN THE MIDDLE SECTION OF NFACTR C --NUNSYM LISTS SMALLEST FIRST C --NSYM(NP+1) IS PRODUCT OF ALL SINGLE FACTORS (IF ANY) DO 800 J=1,NQ JJ=NP+J NUNSYM(J)=NQQ(J) NFACTR(JJ)=NQQ(J) 800 CONTINUE NSYM(NP+1)=NPTS/(NPSYM**2) 900 CONTINUE C --NFACS IS TOTAL NUMBER OF FACTORS NFACS=2*NP+NQ NFACTR(NFACS+1)=0 NPOW2=1 J=0 C --GO THROUGH FACTORS COLLECTING CONSECUTIVE POWERS OF 2 IN CLUSTERS UP C --TO N2GRP 1000 CONTINUE J=J+1 C --TEST FOR LAST FACTOR IF (NFACTR(J).EQ.0) GO TO 1200 C --SEARCH FOR FACTOR OF 2 IF (NFACTR(J).NE.2) GO TO 1000 C --2 FOUND NPOW2=NPOW2*2 C --REDUCE CURRENT FACTOR TO 1 NFACTR(J)=1 C --IS CLUSTER FULL? IF (NPOW2.GE.N2GRP) GO TO 1100 C --GO BACK FOR MORE IF (NFACTR(J+1).EQ.2) GO TO 1000 C --SAVE CURRENT CLUSTER 1100 CONTINUE NFACTR(J)=NPOW2 NPOW2=1 GO TO 1000 C----- 1200 CONTINUE C --IF ONLY ONE SINGLE FACTOR ERASE IT FROM NUNSYM LIST IF (NP.EQ.0) NR=0 JJ=2*NP+NR NSYM(JJ+1)=0 IF (NQ.LE.1) NQ=0 NUNSYM(NQ+1)=0 ERROR=.FALSE. 1300 CONTINUE RETURN C C --ERROR: FACTOR TOO LARGE 1400 CONTINUE WRITE (NWRITE,20) NPMAX,NPTS ERROR=.TRUE. GO TO 1300 C --ERROR: TOO MANY FACTORS 1600 CONTINUE WRITE (NWRITE,21) NFCMAX,NPTS ERROR=.TRUE. GO TO 1300 END C C VAXPDS NAME=STFTP1 C************************************************************************ SUBROUTINE STFTP1(LZEQZI,LSCALE,NMAXI,NXI,NYI,NZI,IERR, 1 MWRITE,JPRINT) C------------------------------------------------------------ C SETTING ROUTINE FOR P1 FOURIER TRANSFORM WITH REAL DATA C A.D. MCLACHLAN 1980. LAST UPDATED 27 AUG 1987 C------------------------------------------------------------ SAVE DIMENSION NDIM(3) INTEGER NMAXI,NXI,NYI,NZI,IERR,MWRITE,JPRINT,NWRITE,IPRINT, + NMAX,NX,NY,NZ,NTOT,NXYZ,NDIM,KAXIS REAL ANXYZ,RROOTN,SCAFAC LOGICAL ZEQZFT,SCALED,CONJBF,CONJAF,ERROR,LSCALE,LZEQZI, LZEQZF REAL ZFT(*), Z(*) C------------------------------------------------------------- C INPUT -- LZEQZI .TRUE. IF Z & ZFT WILL BE SAME ARRAY (IN-PLACE) C INPUT -- LSCALE .TRUE. IF SCALE=1/SQRT(NXYZ), ELSE SCALE=1.0 C INPUT -- NMAXI TOTAL MEMORY SPACE AVAILABLE (REALS) C INPUT -- NXI CELL DIMENSION NX C INPUT -- NYI CELL DIMENSION NY C INPUT -- NZI CELL DIMENSION NZ C OUTPUT -- IERR (0,1) (0) NO ERRORS (1) ERROR DETECTED C INPUT -- MWRITE PRINT UNIT C INPUT -- JPRINT (0,1,2...) PRINT CONTROL C-------------------------------------------------------------- C ENTRY -- STFTP1 SET UP C ENTRY -- FURTRN TRANSFORM FROM REAL TO RECIPROCAL SPACE C ENTRY -- FURINV TRANSFORM BACK TO REAL SPACE C--------------------------------------------------------------- C CALLS -- UNDMFT ONE-DIMENSIONAL TRANSFORM OF 3-D ARRAY C CALLS -- PRNT3D PRINT 3-D ARRAY C--------------------------------------------------------------- C NOTE -- DIMENSION (NX,NY,NZ). NX MUST BE EVEN C NOTE -- DENSITY RHO(X,Y,Z) STORED AS REAL ARRAY SIZE C NOTE -- (NX+2,NY,NZ) C NOTE -- TRANSFORM F(H,K,L) STORED AS COMPLEX ARRAY SIZE (NX/2+1,NY,NZ) C NOTE -- HALF H ALL K,L C--------------------------------------------------------------- 20 FORMAT(1X,' ARRAYS TOO BIG. SIZE ',3I5,' = ',I9,' NMAX=',I9) 21 FORMAT(1X,'***STFTP1 STOPPED IN ERROR***') 22 FORMAT(1X,' NX SHOULD BE EVEN ',3I5) C--------------------------------------------------------------- C --CHECK INPUT C --SAVE ARGUMENTS LZEQZF=LZEQZI NWRITE=MWRITE IPRINT=JPRINT IF(IPRINT.GE.1) WRITE(NWRITE,*) ' <> ' NMAX=NMAXI NX=NXI NY=NYI NZ=NZI IF(IPRINT.GE.2) THEN WRITE(NWRITE,*) ' NX,NY,NZ ',NX,NY,NZ WRITE(NWRITE,*) ' NMAX ',NMAX END IF IERR=0 IF(NX.EQ.2*(NX/2)) GO TO 110 IERR=1 WRITE(NWRITE,22) NX,NY,NZ 110 CONTINUE NTOT=(NX+2)*NY*NZ IF(NTOT.LE.NMAX) GO TO 120 IERR=1 WRITE(NWRITE,20) NX,NY,NZ,NTOT,NMAX 120 CONTINUE NXYZ=NX*NY*NZ ANXYZ=NXYZ IF (LSCALE) THEN RROOTN=1.0/SQRT(ANXYZ) ELSE RROOTN=1.0 ENDIF IF(IERR.EQ.0) GO TO 130 WRITE(NWRITE,21) 130 CONTINUE NDIM(1)=NX/2 + 1 NDIM(2)=NY NDIM(3)=NZ CONJBF=.FALSE. CONJAF=.FALSE. IF(IPRINT.GE.1) WRITE(NWRITE,*) ' ' RETURN C####################################### ENTRY FURTRN( Z , ZFT ) C####################################### C --FORWARD TRANSFORM RHO(X,Y,Z) IF(IPRINT.GE.3) WRITE(NWRITE,*) ' <> ' ZEQZFT=LZEQZF SCALED=.TRUE. SCAFAC=RROOTN ERROR=.FALSE. C --TRANSFORM REAL ALONG 1 (X) KAXIS=1 CALL UNDMFT(Z,ZFT,KAXIS,NX,'REALFT ',NDIM,NDIM,ZEQZFT, 1 SCALED,CONJBF,CONJAF,SCAFAC,ERROR, 2 NWRITE,IPRINT-2) IF(IPRINT.GE.4) THEN CALL PRNT3D(ZFT,NX+2,NY,NZ,'ZFT X ',NWRITE) END IF C --TRANSFORM COMPLEX ALONG 2 (Y) KAXIS=2 ZEQZFT=.TRUE. SCALED=.FALSE. CALL UNDMFT(ZFT,ZFT,KAXIS,NY,'CMPLFT ',NDIM,NDIM,ZEQZFT, 1 SCALED,CONJBF,CONJAF,1.0 , ERROR, 2 NWRITE,IPRINT-2) IF(IPRINT.GE.4) THEN CALL PRNT3D(ZFT,NX+2,NY,NZ,'ZFT Y ',NWRITE) END IF C --TRANSFORM COMPLEX ALONG 3 (Z) KAXIS=3 CONJAF=.TRUE. CALL UNDMFT(ZFT,ZFT,KAXIS,NZ,'CMPLFT ',NDIM,NDIM,ZEQZFT, 1 SCALED,CONJBF,CONJAF,1.0,ERROR, 2 NWRITE,IPRINT-2) CONJAF=.FALSE. IF(IPRINT.GE.3) WRITE(NWRITE,*) ' ' RETURN C########################################## ENTRY FURINV( ZFT , Z ) C########################################## C --REVERSE TRANSFORM F(H,K,L) IF(IPRINT.GE.3) WRITE(NWRITE,*) ' <> ' ZEQZFT=.FALSE. SCALED=.TRUE. SCAFAC=RROOTN ERROR=.FALSE. C --TRANSFORM COMPLEX ALONG 3 (Z) KAXIS=3 CALL UNDMFT(ZFT,Z,KAXIS,NZ,'CMPLFT ',NDIM,NDIM,ZEQZFT, 1 SCALED,CONJBF,CONJAF,SCAFAC,ERROR, 2 NWRITE,IPRINT-2) IF(IPRINT.GE.4) THEN CALL PRNT3D(Z,NX+2,NY,NZ,'Z AFTERZ',NWRITE) END IF C --TRANSFORM COMPLEX ALONG 2 (Y) KAXIS=2 ZEQZFT=.TRUE. SCALED=.FALSE. CALL UNDMFT(Z,Z,KAXIS,NY,'CMPLFT ',NDIM,NDIM,ZEQZFT, 1 SCALED,CONJBF,CONJAF,1.0,ERROR, 2 NWRITE,IPRINT-2) IF(IPRINT.GE.4) THEN CALL PRNT3D(Z,NX+2,NY,NZ,'Z AFTERY',NWRITE) END IF C --TRANSFORM COMPLEX ALONG 1 (X) KAXIS=1 CALL UNDMFT(Z,Z,KAXIS,NX,'HERMFT ',NDIM,NDIM,ZEQZFT, 1 SCALED,CONJBF,CONJAF,1.0,ERROR, 2 NWRITE,IPRINT-2) IF(IPRINT.GE.3) WRITE(NWRITE,*) ' ' C--------------- RETURN END C C VAXPDS NAME=UNDMFT C*************************************************************** SUBROUTINE UNDMFT( Z , ZFT ,KAXIS,NLENG,FTSUB,NDIM, 1 NGRID,ZEQZFT,SCALED,CONJBF,CONJAF, 2 SCAFAC,ERROR,NWRITE,IPRINT) C--------------------------------------------------------------- C A.D. MCLACHLAN 1981. LAST UPDATED 27 JUL 1987. C TRANSFORMS ALONG ONE DIMENSION (KAXIS) OF A 1- 2- OR C 3-DIMENSIONAL ARRAY USING LYNN TEN EYCK'S 'CMPLFT' ROUTINES C-------------------------------------------------------------------- DIMENSION Z(*),ZFT(*),NDIM(3),NGRID(3) LOGICAL ZEQZFT,SCALED,CONJBF,CONJAF,ERROR INTEGER NDIM,NGRID REAL Z,ZFT CHARACTER*8 FTSUB C---------------------------------------------------- INTEGER IDIM,IU,IV,IW,MDIMX,MDIMY,MDIMZ,MGRIDX,MGRIDY, + MGRIDZ,INTLXY,IUVW,MDIM,MGRID DIMENSION IDIM(10) EQUIVALENCE (IDIM(1),IU),(IDIM(2),IV),(IDIM(3),IW) EQUIVALENCE (IDIM(4),MDIMX),(IDIM(5),MDIMY),(IDIM(6),MDIMZ) EQUIVALENCE (IDIM(7),MGRIDX),(IDIM(8),MGRIDY),(IDIM(9),MGRIDZ) EQUIVALENCE (IDIM(10),INTLXY) DIMENSION IUVW(3),MDIM(3),MGRID(3) EQUIVALENCE(IDIM(1),IUVW(1)) EQUIVALENCE(IDIM(4),MDIM(1)) EQUIVALENCE(IDIM(7),MGRID(1)) C --ARRAY IDIM IS USED BY 'CMPLFT' ETC TO CONTROL ARRAY INDEXING INTEGER KAXIS,NLENG,NWRITE,IPRINT,I,NSIZE,IFTSUB REAL SCAFAC CHARACTER*8 SUBNAM(7) DATA SUBNAM /'CMPLFT','REALFT','HERMFT','SDIAD', 1 'INV21','RSYMFT','NULL'/ C------------------------------------------------------ C INPUT -- Z(*) REAL ARRAY (COMPLEX TREATED AS REAL) C OUTPUT -- ZFT(*) FOURIER TRANSFORM OF Z C INPUT -- KAXIS (1,2,3) AXIS OF TRANSFORM C INPUT -- NLENG FULL CELL LENGTH (NOT HALF LENGTH IN REALFT ETC) C INPUT -- FTSUB CH*8 NAME OF BASIC FOURIER ROUTINE USED. C INPUT -- NDIM(3) DIMENSIONS OF ARRAY C INPUT -- NGRID(3) DIMENSIONS OF GRID ACTUALLY USED (NGRID.LE.NDIM) C INPUT -- ZEQZFT .L. .TRUE. IF Z AND ZFT ARE SAME ARRAY. C INPUT -- OTHERWISE Z IS COPIED TO ZFT C INPUT -- SCALED .L. .TRUE. IF SCALE FACTOR SCAFAC IS TO BE APPLIED C INPUT -- CONJBF .L. .TRUE. IF COMPLEX CONJ BEFORE TRANSFORM C INPUT -- CONJAF .L. .TRUE. IF COMPLEX CONJ AFTER TRANSFORM C INPUT -- SCAFAC SCALE FACTOR TO APPLY TO ZFT C OUTPUT -- ERROR .L. .TRUE. IF ERROR DETECTED C INPUT -- NWRITE PRINT UNIT C INPUT -- IPRINT (0,1,2...) PRINT CONTROL C----------------------------------------------------------------------------- C CALLS -- CMPLFT COMPLEX FOURIER TRANSFORM C CALLS -- REALFT REAL DATA TRANSFORM C CALLS -- HERMFT HERMITIAN DATA TRANSFORM C CALLS -- SDIAD SCREW DIAD TRANSFORM C CALLS -- INV21 INVERSE OF SDIAD C CALLS -- RSYMFT REAL SYMMETRIC DATA TRANSFORM C----------------------------------------------------------------------------- C NOTE -- FTSUB MUST BE ONE OF THE NAMES C NOTE -- CMPLFT C NOTE -- REALFT C NOTE -- HERMFT C NOTE -- SDIAD C NOTE -- INV21 C NOTE -- RSYMFT C NOTE -- NULL C NOTE -- NDIM(1,2,3) ARE DIMENSIONS OF ARRAY, TAKEN AS COMPLEX EXCEPT C NOTE -- FOR RSYMFT WHERE DIMENSIONS ARE TREATED AS REAL C-------------------------------------------------------------------- 20 FORMAT(1X,'***UNDMFT ERROR *** FTSUB=(',A8,') NOT KNOWN ') C------------------------------------------------------- IF(IPRINT.GE.1) WRITE(NWRITE,*) ' <> ' ERROR=.FALSE. IU=KAXIS IV=KAXIS+1 IW=KAXIS+2 IF(IV.GT.3) IV=IV-3 IF(IW.GT.3) IW=IW-3 DO 100 I=1,3 MDIM(I)=NDIM(I) MGRID(I)=NGRID(I) 100 CONTINUE C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C SUBSECTION FOR IDENTIFYING SUBROUTINE NAMES C IDENTIFY SUBROUTINE 'FTSUB' DO 210 IFTSUB=1,7 IF(SUBNAM(IFTSUB).EQ.FTSUB) GO TO 220 210 CONTINUE ERROR=.TRUE. WRITE(NWRITE,20) FTSUB 220 CONTINUE INTLXY=2 IF(IFTSUB.EQ.6) INTLXY=1 C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ IF(ERROR) RETURN C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ C FOURIER TRANSFORM SUBSECTION ENTRY A C IF ARRAYS Z ZFT NOT SAME COPY Z INTO ZFT NSIZE=MDIMX*MDIMY*MDIMZ*INTLXY IF(ZEQZFT) GO TO 280 DO 270 I=1,NSIZE ZFT(I)=Z(I) 270 CONTINUE 280 CONTINUE IF(.NOT.SCALED) GO TO 290 DO 285 I=1,NSIZE ZFT(I)=ZFT(I)*SCAFAC 285 CONTINUE 290 CONTINUE IF(.NOT.CONJBF) GO TO 300 DO 295 I=2,NSIZE,2 ZFT(I)=-ZFT(I) 295 CONTINUE C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 300 CONTINUE C FOURIER TRANSFORM SUBSECTION ENTRY B C USE THE REQUIRED FFT ROUTINE GO TO (400,410,420,430,440,450,500),IFTSUB 400 CONTINUE CALL CMPLFT(ZFT(1),ZFT(2),NLENG,IDIM,NWRITE) GO TO 500 410 CONTINUE CALL REALFT(ZFT(1),ZFT(2),NLENG/2,IDIM,NWRITE) GO TO 500 420 CONTINUE CALL HERMFT(ZFT(1),ZFT(2),NLENG/2,IDIM,NWRITE) GO TO 500 430 CONTINUE CALL SDIAD(ZFT(1),ZFT(2),NLENG/2,IDIM,NWRITE) GO TO 500 440 CONTINUE CALL INV21(ZFT(1),ZFT(2),NLENG/2,IDIM,NWRITE) GO TO 500 450 CONTINUE CALL RSYMFT(ZFT(1),NLENG/2,IDIM,NWRITE) GO TO 500 500 CONTINUE IF(.NOT.CONJAF) GO TO 510 DO 505 I=2,NSIZE,2 ZFT(I)=-ZFT(I) 505 CONTINUE 510 CONTINUE IF(IPRINT.GE.1) WRITE(NWRITE,*) ' ' C---------------------- RETURN END