      PROGRAM FOUR2
C  FOUR2 CALCULATES A TWO-DIMENSIONAL PROJECTED ELECTRON DENSITY MAP
C  FROM STRUCTURE FACTORS INPUT FROM THE OUTPUT FILE OF STRUCFAC.
C  THE DENSITY IS SCALED TO A MAXIMUM VALUE OF 100 AND THE SCALING
C  FACTOR TO CONVERT TO ELECTRONS/ANGSTROM**2 IS GIVEN.   TO TAKE 
C  ACCOUNT OF THE LIMITED WIDTH OF PRINTER PAPER THE OUTPUT IS IN 
C  SECTIONS WHICH HAVE TO BE JOINED SIDEWAYS.   FOR ECONOMY THE 
C  OUTPUT IS ON AN ORTHOGONAL GRID AND THE SCALING WILL NOT BE 
C  ACCURATE ALONG THE AXES.   THE VALUES OF HMAX AND HMIN MUST BE
C  <= 15 AND THE MINIMUM GRID SPACING IS 1/60.
C  THERE IS A PATTERSON FUNCTION OPTION IN WHICH CASE THE SCALING 
C  FACTOR IS ELECTRONS**2/ANGSTROM**2.    
      DIMENSION A(-25:25,0:25),B(-25:25,0:25),AC(0:25,0:50),
     +AS(0:25,0:50),BC(0:25,0:50),BS(0:25,0:50),CSX(0:60),
     +CSY(0:60),RO(0:60,0:60),LIST(0:60),HH(1000),KK(1000),
     +F(1000),PHI(1000)
      INTEGER H,HH
      CHARACTER ANS*1,FNAME*10,DIR*1,BT*6
      DATA NPAGE/15/
      OPEN(UNIT=9,FILE='LPT1')
      TPI=6.2831853
C  CONSTANT FOR CONVERTING DEGREES TO RADIANS
      DTR=0.0174533
C  CLEAR A AND B ARRAYS
      DO 1 H=-25,25
      DO 1 K=0,25
      A(H,K)=0
      B(H,K)=0
    1 CONTINUE
      FNAME='SF1.DAT'
    4 WRITE(6,'('' IT IS ASSUMED THAT THE INPUT DATA FILE IS NAMED'')')
      WRITE(6,'('' SF1.DAT.   DO YOU WISH TO CHANGE THE NAME?  Y/N'')')
      READ(5,50)ANS
   50 FORMAT(A1)
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 2
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 3
      GOTO 4
    3 WRITE(6,'('' READ IN FILENAME [<= 10 CHARACTERS] '')')
      READ(5,55)FNAME
   55 FORMAT(A10)
    2 OPEN(UNIT=10,FILE=FNAME)
      HMAX=0
      KMAX=0
C  INPUT THE CELL DIMENSIONS A &amp; B AND THE CELL ANGLE GAMMA.
      WRITE(6,'('' INPUT CELL DIMENSIONS A,B AND GAMMA (degrees)'')')
      WRITE(6,'('' MAKE SURE THAT THEY ARE CONSISTENT WITH THE '')')
      WRITE(6,'('' X,Y COORDINATES USED IN PROGRAMME STRUCFAC. '')')
      READ(5,*)AA,BB,GAMMA
      AREA=AA*BB*SIN(DTR*GAMMA)
      READ(10,*)NOREF,(HH(I),KK(I),F(I),PHI(I),I=1,NOREF)
   95 WRITE(6,'('' DO YOU WANT THE PATTERSON OPTION?  Y/N '')')
      WRITE(6,'('' REMEMBER - IF THE INPUT TAPE IS FROM PROGRAMME'')')
      WRITE(6,'('' ISOFILE THEN YOU SHOULD INPUT N '')')
      READ(5,50)ANS
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 93      
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 94
      GOTO 95
   93 DO 96 I=1,NOREF
      F(I)=F(I)*F(I)
      PHI(I)=0
   96 CONTINUE
   94 DO 77 I=1,NOREF
      IF(ABS(HH(I)).GT.HMAX)HMAX=ABS(HH(I))
      IF(KK(I).GT.KMAX)KMAX=KK(I)
      A(HH(I),KK(I))=F(I)*COS(DTR*PHI(I))
      B(HH(I),KK(I))=F(I)*SIN(DTR*PHI(I))
   77 CONTINUE   
C  INPUT OF DATA COMPLETE
      WRITE(6,'('' INPUT NX AND NY, THE NUMBER OF DIVISIONS IN THE'')')
      WRITE(6,'('' X & Y DIRECTIONS. THEY SHOULD BE DIVISIBLE BY 4.'')')
      READ(5,*)NX,NY
      NX4=NX/4
      NY4=NY/4
      NX16=16*NX
      NY16=16*NY
C  TABLES ARE NOW SET UP OF COS(2*PI*N/NX) AND COS(2*PI*M/NY) WITH N
C  FROM 0 TO NX AND M FROM 0 TO NY. THESE TABLES ARE USED IN THE
C  SUMMATIONS WHICH FOLLOW AND THEIR USE IS FAR MORE EFFICIENT THAN
C  CALCULATING EACH SINE AND COSINE AS IT IS REQUIRED.
      DO 10 I=0,NX
   10 CSX(I)=COS(TPI*I/NX)
      DO 11 I=0,NY
   11 CSY(I)=COS(TPI*I/NY)
C  NOW DO THE SUMMATIONS AC, AS, BC AND BS AFTER CLEARING ARRAYS
      DO 12 I=0,25
      DO 12 J=0,50
      AC(I,J)=0
      AS(I,J)=0
      BC(I,J)=0
      BS(I,J)=0
   12 CONTINUE
      DO 13 K=0,KMAX
      DO 13 IX=0,NX/2
      DO 13 H=-HMAX,HMAX
      IF(H.EQ.0.AND.K.EQ.0)GOTO 13
      COZ=CSX(MOD(H*IX+NX16,NX))
      ZIN=CSX(MOD(H*IX+NX16-NX4,NX))
      AC(K,IX)=AC(K,IX)+A(H,K)*COZ
      AS(K,IX)=AS(K,IX)+A(H,K)*ZIN            
      BC(K,IX)=BC(K,IX)+B(H,K)*COZ
      BS(K,IX)=BS(K,IX)+B(H,K)*ZIN
   13 CONTINUE
C  THE SUMMATIONS OVER X FOR AC,AS,BC AND BS ARE FOR X = 0 TO 1/2.
C  THE REMAINING REGION FROM X = 1/2 TO 1 WILL BE GIVEN BY SYMMETRY.
C  THE F(0,0) TERM HAS BEEN EXCLUDED.   IT HAS HALF THE WEIGHT OF THE
C  OTHER TERMS AND ITS CONTRIBUTION IS NOW LOADED INTO ARRAY RO WHICH
C  WILL CONTAIN THE SCALED DENSITY.
      DO 15 IX=0,NX
      DO 15 IY=0,NY
   15 RO(IX,IY)=A(0,0)/2.0
      DO 16 IX=0,NX/2
      W=1.0
      IF(IX.EQ.NX/2)W=0.5
      DO 16 IY=0,NY
      DO 16 K=0,KMAX
      COZ=CSY(MOD(K*IY+NY16,NY))
      ZIN=CSY(MOD(K*IY+NY16-NY4,NY))
      RO(IX,IY)=RO(IX,IY)+W*(AC(K,IX)+BS(K,IX))*COZ
     ++(BC(K,IX)-AS(K,IX))*ZIN
      RO(NX-IX,IY)=RO(NX-IX,IY)+W*(AC(K,IX)-BS(K,IX))*COZ
     ++W*(BC(K,IX)+AS(K,IX))*ZIN
   16 CONTINUE
C  THE DENSITY TABLE IS NOW COMPLETE AND IT IS SCALED TO A MAXIMUM 
C  MAGNITUDE OF 100.
      RMAX=0
      DO 20 I=0,NX-1
      DO 20 J=0,NY-1
      IF(ABS(RO(I,J)).GT.RMAX)RMAX=ABS(RO(I,J))
   20 CONTINUE
      FAC=100/RMAX
      DO 21 I=0,NX
      DO 21 J=0,NY
      RO(I,J)=FAC*RO(I,J)
   21 CONTINUE
C  CALCULATE FACTOR WHICH GIVES PROJECTED DENSITY FROM PRINTED VALUES
      CAF=1/FAC/AREA
C  OUTPUT THE FACTOR
      WRITE(9,700)CAF
  700 FORMAT(42H FACTOR CONVERTING TO DENSITY/UNIT AREA IS,F8.4)
      DO 66 I=0,60
      LIST(I)=I
   66 CONTINUE
C  FIND A CONVENIENT ORIENTATION FOR THE OUTPUT WHICH IS THAT WITH 
C  THE LEAST NUMBER OF DIVISIONS ACROSS THE PAGE.
      DIR='X'
      BT='BOTTOM'
      M=NX/NPAGE+1
      L=0
      IF(NY.LT.NX)THEN
      DIR='Y'
      BT='   TOP'
      M=NY/NPAGE+1
      L=1
      ENDIF
      WRITE(9,100)DIR,BT
  100 FORMAT(13H OUTPUT WITH ,A1,23H HORIZONTAL AND ORIGIN ,A6,5H LEFT)
      IF(L.EQ.1)GOTO 30
C  THE MAP IS PRINTED IN SECTIONS WHICH MUST BE JOINED AT SIDE EDGES.
C  THE NUMBER OF GRIDPOINTS ACROSS THE PAGE IS NPAGE WHICH IS SET AT
C  15.   FOR OTHER AVAILABLE PAGE WIDTHS THIS CAN BE ALTERED BY THE
C  DATA STATEMENT NEAR THE BEGINNING OF THE PROGRAMME.
      DO 25 N0=1,M
      WRITE(9,'(''  '')')
      WRITE(9,'(''  '')')
      N1=(N0-1)*NPAGE
      N2=NPAGE*N0
      IF(N0.EQ.M)N2=NX
      WRITE(9,150)(LIST(I),I=N1,N2)
  150 FORMAT(5X,21I4)
      WRITE(9,'(''  '')')
      DO 25 IY=NY,0,-1
      WRITE(9,200)IY,(RO(IX,IY),IX=N1,N2)
  200 FORMAT(I4,2X,21F4.0)
      WRITE(9,'(''  '')')
   25 CONTINUE
      GOTO 45
   30 DO 26 N0=1,M
      WRITE(9,'(''  '')')
      WRITE(9,'(''  '')')
      N1=(N0-1)*NPAGE
      N2=NPAGE*N0
      IF(N0.EQ.M)N2=NY
      WRITE(9,150)(LIST(I),I=N1,N2)
      WRITE(9,'(''  '')')
      DO 26 IX=0,NX
      WRITE(9,200)IX,(RO(IX,IY),IY=N1,N2)
      WRITE(9,'(''  '')')
   26 CONTINUE
   45 STOP
      END
   
    

