	PROGRAM FOUR
C  GENERALIZED FOURIER PROGRAM   MAY 1974   BY L. W. FINGER
C	CONVERTED TO PDP-11, OCT 78, JULY 79
C	TO RSX-11M ON 28-MAY-80
C	TO VAX/VMS ON 9-NOV-81
C	2*FO - FC OPTION ADDED 23-NOV-82
C
	PARAMETER MAXST=30000
      DIMENSION FBUF(MAXST),FRMT(50),SINHX(MAXST),COSHX(MAXST)
	DIMENSION SINKY(MAXST),SINLZ(MAXST),COSLZ(MAXST),DHK(MAXST)
	DIMENSION EHK(MAXST),UHK(MAXST),VHK(MAXST),WHK(MAXST)
	DIMENSION COSKY(MAXST),XHK(MAXST),YH(MAXST),ZH(MAXST)
      REAL*8 LABEL(22)
      DIMENSION TITLE(20),JH(3),IH(3),IP(43),RMAT(3,3,24),IRMAT(3,
     13,24),TRANS(3,24),IEXTR(6),FX(21),FMULT(12),IFBUF(500)
      INTEGER OUT,S1,S2,S3
      EQUIVALENCE (IEXTR(1),IBX),(IEXTR(2),IBY),(IEXTR(3),IBZ)
      EQUIVALENCE (IEXTR(4),IMAX),(IEXTR(5),KMAX),(IEXTR(6),LMAX)
      EQUIVALENCE (FX,IP,FMULT),(FBUF,IFBUF)
	EQUIVALENCE (FBUF,SINHX,COSHX,SINKY,COSKY,SINLZ,COSLZ,DHK
	1,EHK,UHK,VHK,WHK,XHK,YH,ZH)
C
      DATA FRMT/4H(1X,,4HF4.0,3H,1X,47*1H /
      DATA LABEL/1H ,6H  DIFF,3HXYZ,3HXZY,3HYXZ,3HYZX,3HZYX,3HZXY,
     1 3HH00,3H00L,3HHK0,3HH0L,3H0KL,3HHH0,3HH0H,3H0KK,3HHHH,3HHHL,3HHKK
     2,3HHKL,6H  BETA,6H2Fo-Fc/
      DATA IN,OUT,S1,ZRO,AZRO,ENDD,S2,S3/5,6,12,3H,4X,3H,I4,1H),11,13/
      DATA TWOPI,ISPC,IBNK/6.28318531,1H0,1H /
C
    2 FORMAT(20A4)
    3 FORMAT(I8,E16.6/(5E16.6))
    4 FORMAT(7I1,I5,8I4,3F8.0/12F6.0)
    5 FORMAT(4I4,4F8.3)
    6 FORMAT(1X,F4.0)
    7 FORMAT(2A3,I2,3F8.5)
    8 FORMAT(1H1,20A4,5X,7HSECTION,I4,1H/,I4//F10.0,14F8.0)
    9 FORMAT(' PROBLEM TOO LARGE, INCREASE DIMENSIONS AND RERUN.')
   12 FORMAT(/5X,24HALL REFLECTIONS INCLUDED )
   13 FORMAT(/5X,43HREFLECTIONS WITH REJ. FLAG SET NOT INCLUDED )
   14 FORMAT(1H1,20A4//2X,A6,8H FOURIER,10X,A3,9H SECTIONS
     1//5X,9HF(000) = F10.3
     2//5X,19HMAP EXAGGERATION = F10.3
     3//5X,31HMULTIPLICITY OF GENERAL REFL = F10.3)
   15 FORMAT(/40X,14HMULTIPLICITIES/5X,11(A3,5X),A3//12F8.2//
     2//9X,50HINITIAL   FINAL    INTERVAL   MIN INDEX  MAX INDEX/)
   16 FORMAT(10H POINTS   I4,I8,I12,I10,I11/8H LINES  ,I6,I8,I12,I10,I11
     1/9H SECTIONS,I5,I8,I12,I10,I11)
   17 FORMAT(/5X,34HSECTIONS SAVED FOR LATER PLOTTING.)
   18 FORMAT(///10X,22HSTORAGE REQUIREMENTS = I6,11H LOCATIONS.)
   20 FORMAT(/15X,18HFRACTIONAL COORDS.,11X,11HMAP COORDS./5H ATOM,8X,1H
     1X,8X,1HY,8X,1HZ,9X,5HPOINT,4X,4HLINE,3X,7HSECTION/)
   22 FORMAT(1X,2A3,1X,3F9.5,3X,3F8.2)
C
	READ(IN,2)TITLE
      READ(IN,4)ITYPE,ISECT,ICENT,KREJ,IPLT,LSPC,IPRNT,IX,IXI,IXF,IY,IYI
     1,IYF,IZ,IZI,IZF,EXXAG,FOOO,TMULT,FMULT
	OPEN(UNIT=S3,TYPE='SCRATCH',FORM='UNFORMATTED')
      READ(S1)NSYM,VOL,(((RMAT(I,J,K),I=1,3),J=1,3),K=1,NSYM)
     1,((TRANS(J,I),J=1,3),I=1,NSYM)
      IF(ITYPE.GE.2) ITYPE=1-ITYPE
      ITEMP=ITYPE+1
      IF(ITYPE.LT.0)ITEMP=20-ITYPE
      WRITE(OUT,14)TITLE,LABEL(ITEMP),LABEL(ISECT+3),FOOO,EXXAG,TMULT
      IF(IPLT.EQ.0)GO TO 100
	WRITE(OUT,17)
100	IF(KREJ)101,102,101
  101 WRITE(OUT,13)
      GO TO 103
  102 WRITE(OUT,12)
  103 WRITE(OUT,15)(LABEL(I),I=9,20),FMULT
      FOOO=FOOO*EXXAG/VOL
      EXXAG=2.0*EXXAG/(VOL*TMULT)
C        TRANSFORM AND CONVERT SYMMETRY OPERATORS
                             DO 113 I=1,NSYM
                             DO 113 J=1,3
      TRANS(J,I)=TWOPI*TRANS(J,I)
                             DO 112 K=1,3
  112 FX(K)=RMAT(J,K,I)
      CALL SECT(ISECT,FX)
                             DO 113 K=1,3
  113 IRMAT(J,K,I)=FX(K)
      DO 120 I=1,6
  120 IEXTR(I)=0
C        READ REFLECTION DATA
	NFBR=504
130	NFBR=NFBR+10
	IF(NFBR.LT.500)GO TO 140
C READ NEW BLOCK
	READ(S1)(FBUF(I),I=1,504)
	NFBR=0
140	IH(1)=IFBUF(NFBR+1)
	IH(2)=IFBUF(NFBR+2)
	IH(3)=IFBUF(NFBR+3)
	FOBS=FBUF(NFBR+4)
	FCALC=FBUF(NFBR+5)
	ACALC=FBUF(NFBR+6)
	BCALC=FBUF(NFBR+7)
      IF(IH(1)-99)160,600,160
C        NOT END, CHECK IF REJECTED
  160 IF(IFBUF(NFBR+8))170,180,170
  170 IF(KREJ)130,180,130
C        REFLECTION NOT REJECTED
  180 IF(ITYPE)190,220,210
  190 IF(ITYPE+2)200,200,205
C	 COEFFICIENT FOR 2*FO-FC MAP
  200 FOBS=2.0*FOBS-FCALC
      GO TO 220
C        COEFFICIENT FOR BETA SYNTHESIS - FO**2/FC
  205 TEMP=FCALC
      IF(ABS(TEMP).LT.1.0) TEMP=FOBS
      FOBS=FOBS**2/TEMP
      GO TO 220
C        DIFFERENCE FOURIER SYNTHESIS
  210 FOBS=FOBS-FCALC
C        CALCULATE PHASE ANGLE
  220 IF(ICENT)230,225,230
  225 BCALC=0.0
  230 PHASE=ATAN2(BCALC,ACALC)
C             DETERMINE REFLECTION TYPE
      IF(IH(1))350,280,350
  280 IF(IH(2))300,290,300
C             00L REFLECTION
  290 IC=2
      GO TO 530
  300 IF(IH(3))320,310,320
C             0K0 REFLECTION
  310 IC=1
      GO TO 530
  320 IF(IH(2)-IH(3))340,330,340
C             0KK REFLECTION
  330 IC=8
      GO TO 530
C             0KL REFLECTION
  340 IC=5
      GO TO 530
  350 IF(IH(2))410,360,410
  360 IF(IH(3))380,370,380
C             H00 REFLECTION
  370 IC=1
      GO TO 530
  380 IF(IH(1)-IH(3))400,390,400
C             H0H REFLECTION
  390 IC=7
      GO TO 530
C             H0L REFLECTION
  400 IC=4
      GO TO 530
  410 IF(IH(3))450,420,450
  420 IF(IH(1)-IH(2))440,430,440
C             HH0 REFLECTION
  430 IC=6
      GO TO 530
C             HK0 REFLECTION
  440 IC=3
      GO TO 530
  450 IF(IH(1)-IH(2))490,460,490
  460 IF(IH(1)-IH(3))480,470,480
C             HHH REFLECTION
  470 IC=9
      GO TO 530
C             HHL REFLECTION
  480 IC=10
      GO TO 530
  490 IF(IH(2)-IH(3))510,500,510
C             HKK OR HKH REFLECTION
  500 IC=11
      GO TO 530
  510 IF(IH(1)-IH(3))520,500,520
C             HKL REFLECTION
  520 IC=12
  530 FMAG=EXXAG*FMULT(IC)*FOBS
	WRITE(S3)IH,FMAG,PHASE
C        EXPAND INDEX LIST
      DO 590 I=1,NSYM
      DO 540 J=1,3
      JH(J)=0.0
      DO 540 K=1,3
  540 JH(J)=JH(J)+IRMAT(K,J,I)*IH(K)
      DO 590 J=1,3
      IF(JH(J)-IEXTR(J))570,570,560
C        SET NEW MAXIMUM
  560 IEXTR(J)=JH(J)
      GO TO 590
  570 IF(JH(J)-IEXTR(J+3))580,590,590
C        SET NEW MINIMUM
  580 IEXTR(J+3)=JH(J)
  590 CONTINUE
      GO TO 130
C        END OF REFLECTION LIST
  600 ENDFILE S3
      WRITE(OUT,16)IXI,IXF,IX,IEXTR(4),IBX,IYI,IYF,IY,IEXTR(5),IBY,
     $IZI,IZF,IZ,IEXTR(6),IBZ
      XX=TWOPI/FLOAT(IX)
      YY=TWOPI/FLOAT(IY)
      ZZ=TWOPI/FLOAT(IZ)
      IXA=IXF-IXI+1
      IYA=IYF-IYI+1
      IZA=IZF-IZI+1
      DO 620 I=1,3
      IEXTR(I+3)=MAX0(IEXTR(I),-IEXTR(I+3))
  620 IEXTR(I)=2*IEXTR(I+3)+1
      IBX=IMAX+1
      IPN2=IBX*IBY
      IPN3=IPN2+IBX*IBY
      IPN4=IPN3+IMAX*IXA
      IPN5=IPN4+IMAX*IXA
      I=MAX0(KMAX,IBZ)
      IPN6=IPN5+I
      IPN7=IPN6+I
      IPN8=IPN7+IBX
      IPN9=IPN8+IBX
      IPN10=IPN9
      IF(IPLT.EQ.0)GO TO 621
C        SECTIONS TO BE SAVED FOR PLOTTING
      IPN10=IPN9+IXA*IYA
  621 IPN10=MAX0(IPN10,1020)
      WRITE(OUT,18)IPN10
      IF(IPN10.GT.MAXST)GO TO 895
C  READ ATOM PARAMETERS FROM THE FILE
 6211 READ(S1,END=624,ERR=624)NATOM,((FBUF(17*I+J-17),J=1,17)
	1,I=1,NATOM)
	WRITE(OUT,20)
	J=-17
	DO 623 I=1,NATOM
	J=J+17
      DO 622 K=1,3
  622 FX(K)=FBUF(J+K+8)
C        TRANSFORM ATOM COORDINATES TO MAP SYSTEM
      CALL SECT(ISECT,FX)
      FX(1)=FX(1)*IX
      FX(2)=FX(2)*IY
      FX(3)=FX(3)*IZ
C        PRINT BOTH ORIGINAL AND TRANSFORMED COORDINATES
  623 WRITE(OUT,22)(FBUF(J+K),K=1,2),(FBUF(J+K),K=9,11)
     1,(FX(K),K=1,3)
C        FILL IN SIN AND COS 2 PI H X
  624 CLOSE(unit=S1)
      K=IPN3
      L=IPN4
      DO 630 I=1,IXA
      HY=XX*FLOAT(I+IXI-1)
      DO 630 J=1,IMAX
      HX=HY*FLOAT(J)
      K=K+1
      L=L+1
      COSHX(K)=COS(HX)
  630 SINHX(L)=SIN(HX)
C        START SECTION LOOP
      DO 890 IZC=1,IZA
	REWIND S3
      J=IZC+IZI-1
C        PRECALCULATE COS AND SIN 2 PI L Z
      HY=ZZ*FLOAT(J)
      DO 650 I=1,IBZ
      HX=HY*FLOAT(I-LMAX-1)
      COSLZ(IPN5+I)=COS(HX)
  650 SINLZ(IPN6+I)=SIN(HX)
C        INITIALIZE SECTION ARRAYS
      DO 660 I=1,IPN2
	DHK(I)=0.0
  660 EHK(I+IPN2)=0.0
      DHK(KMAX+1)=FOOO
C        LOOP THROUGH REFLECTION LIST
C
C		Form Section Sums
C
C			D(hk) = sum(l) (A(hkl)*cos(2pi lz)+B(hkl)*sin(2pi lz))
C			E(hk) = Sum(l) (B(hkl)*cos(2pi lz)-A(hkl)*sin(2pi lz))
C
C
  700 READ(S3,END=770)IH,FMAG,PHASE
C        LOOP THROUGH SYMMETRY OPERATIONS
      DO 760 I=1,NSYM
      HX=PHASE
      DO 720 J=1,3
      HX=HX-IH(J)*TRANS(J,I)
      JH(J)=0
      DO 720 K=1,3
  720 JH(J)=JH(J)+IRMAT(K,J,I)*IH(K)
      IF(JH(1))730,750,750
C        TRANSFORMED H.LT.0, CHANGE REFLECTION TO FRIEDEL EQUIV.
  730 DO 740 J=1,3
  740 JH(J)=-JH(J)
      HX=-HX
  750 J=JH(1)*IBY+JH(2)+KMAX+1
      K=JH(3)+LMAX+1
      FSS=FMAG*SIN(HX)
      FCC=FMAG*COS(HX)
      DHK(J)=DHK(J)+FCC*COSLZ(K+IPN5)+FSS*SINLZ(K+IPN6)
  760 EHK(IPN2+J)=EHK(IPN2+J)+FSS*COSLZ(K+IPN5)-FCC*SINLZ(K+IPN6)
      GO TO 700
C        ALL REFLECTIONS SUMMED, REARRANGE SECTION SUMS
C
C	U(hk) = D(hk) + D(h-k)
C	V(hk) =-D(hk) + D(h-k)
C	W(hk) = E(hk) + E(h-k)
C	X(hk) = E(hk) - E(h-k)
C
  770 DO 780 I=1,IBX
      K=IBY*(I-1)+KMAX+1
      L=K
      DO 780 J=1,KMAX
      K=K+1
      L=L-1
      HX=DHK(K)
      HY=DHK(L)
      UHK(K)=HX+HY
      VHK(L)=HY-HX
      HX=EHK(IPN2+K)
      HY=EHK(IPN2+L)
      WHK(IPN2+K)=HX+HY
  780 XHK(IPN2+L)=HX-HY
C        START LOOP THROUGH STRIPS ALONG X
      DO 880 IXZ=1,IXA,30
      IXX=IXZ+IXI-1
      NXA=MIN0(30,IXF-IXX+1)
      NX=0
      DO 782 I=1,NXA,2
      NX=NX+1
  782 FX(NX)=I+IXX-1
      J=IZC+IZI-1
      IF(IPRNT.EQ.0)WRITE(OUT,8)TITLE,J,IZ,(FX(I),I=1,NX)
C        START LOOP FOR LINES
      DO 880 IAJ=1,IYA
      IPSV=(IAJ-1)*IXA+IXX-IXI+IPN9
      FY=IAJ+IYI-1
      HY=YY*FY
      DO 785 I=1,KMAX
      HX=HY*FLOAT(I)
      COSKY(IPN5+I)=COS(HX)
  785 SINKY(IPN6+I)=SIN(HX)
C        LOOP OVER VALUES OF H
      DO 790 I=1,IBX
      K=IBY*(I-1)+KMAX+1
      L=K
      M=IPN7+I
      N=IPN8+I
      YH(M)=UHK(K)
      ZH(N)=WHK(K+IPN2)
C        SUM OVER K
C
C	Y(h) = Sum(k) (cos(2pi ky)*U(hk) + sin(2pi ky)*X(hk))
C	Z(h) = Sum(k) (sin(2pi ky)*V(hk) + cos(2pi ky)*W(hk))
C
      DO 790 J=1,KMAX
      K=K+1
      L=L-1
      YH(M)=YH(M)+UHK(K)*COSKY(J+IPN5)+XHK(L+IPN2)*SINKY(J+IPN6)
  790 ZH(N)=ZH(N)+WHK(IPN2+K)*COSKY(J+IPN5)+VHK(L)*SINKY(J+IPN6
     1)
C        START POINTS LOOP
      L=0
      K=(IXX-IXI)*IMAX
      DO 870 I=1,NXA
      HX=YH(IPN7+1)
C        SUM OVER H
      DO 800 J=1,IMAX
      K=K+1
  800 HX=HX+YH(IPN7+J+1)*COSHX(IPN3+K)+ZH(IPN8+J+1)*SINHX(IPN4+K)
C        MAP POINT CALCULATED, CONVERT FOR PRINTING AND ADD TO PLOT OUTP
      IF(IPLT.EQ.0)GO TO 805
      IPSV=IPSV+1
      FBUF(IPSV)=HX
  805 IF(IPRNT.NE.0)GO TO 870
      J=HX+SIGN(.5,HX)
      IF(J) 810,860,830
  810 IF(J+1000)820,820,850
  820 J=J+1000
      GO TO 810
  830 IF(J-10000)850,840,840
  840 J=J-10000
      GO TO 830
  850 L=L+1
      IP(L)=J
      FRMT(I+3)=AZRO
      GO TO 870
  860 FRMT(I+3)=ZRO
  870 CONTINUE
      IF(IPRNT.NE.0)GO TO 880
      IF(LSPC.EQ.0)WRITE(OUT,6)
      IF(L)872,872,876
  872 WRITE(OUT,6)FY
      GO TO 880
  876 FRMT(NXA+4)=ENDD
C        PRINT LINE
      WRITE(OUT,FRMT)FY,(IP(I),I=1,L)
  880 CONTINUE
C        SECTION NOW COMPLETE - DUMP PLOT OUTPUT IF NEEDED
      IF(IPLT.EQ.0)GO TO 890
      J=IZC+IZI-1
      WRITE(S2)IXA,IYA,J,IZ,IX,IXI,IXF,IY,IYI,IYF,IZI,IZF
      J=IPN9+1
      WRITE(S2)(FBUF(I),I=J,IPN10)
  890 CONTINUE
C        END OF MAP
      GO TO 900
  895 WRITE(OUT,9)
  900 STOP
      END
      SUBROUTINE SECT(ISECT,IH)
C
      REAL IH(3),ITEMP
C
      IF(ISECT)25,25,5
    5 GO TO (10,30,50,60,70),ISECT
C             XZY SECTIONS
   10 ITEMP=IH(2)
   20 IH(2)=IH(3)
      IH(3)=ITEMP
   25 RETURN
C             YXZ SECTIONS
   30 ITEMP=IH(2)
   40 IH(2)=IH(1)
      IH(1)=ITEMP
      RETURN
C             YZX SECTIONS
   50 ITEMP=IH(1)
      IH(1)=IH(2)
      GO TO 20
C             ZYX SECTIONS
   60 ITEMP=IH(1)
      IH(1)=IH(3)
      IH(3)=ITEMP
      RETURN
C             ZXY SECTIONS
   70 ITEMP=IH(3)
      IH(3)=IH(2)
      GO TO 40
      END
