      PROGRAM MINDIR
C  ----------------------------------------------------------------
C                    USER PLEASE NOTE
C  THIS PROGRAMME CONTAINS A RANDOM NUMBER GENERATOR WHICH USES
C  INTEGERS.   SOME COMPILERS HAVE DEFAULT INTEGERS WITH ONLY 
C  16 BITS WHICH IS NOT SUFFICIENT.   HOWEVER, IN THAT CASE THERE
C  WILL ALWAYS BE A COMPILER OPTION WHICH GIVES LONG INTEGERS
C  AND THIS SHOULD BE USED.
C  ----------------------------------------------------------------  
C  THIS PROGRAMME GIVES A SIMULATION OF A DIRECT-METHODS SOLUTION 
C  FOR EITHER A CENTROSYMMETRIC OR A NON-CENTROSYMMETRIC STRUCTURE.
C  IT FINDS A NUMBER OF SETS OF PHASES WHICH SATISFY THE TANGENT 
C  FORMULA, THEN APPLIES SIMPLE FIGURE-OF-MERIT TESTS AND THEN
C  COMPUTES THE E-MAP WITH THE HIGHEST COMBINED FIGURE OF MERIT. 
C  IF THIS DOES NOT GIVE A SOLUTION THEN IT SUBSEQUENTLY COMPUTES
C  AND PRINT E-MAPS FOR ANY PHASE SET SELECTED BY THE USER.      
C  THE STEPS IN THE PROCESS ARE:
C	(1) GENERATE RANDOM PHASES FOR ALL REFLECTIONS AND GIVE 
C           INITIAL WEIGHTS 0.25.
C	(2) CALCULATE A DENSITY MAP, WITH WEIGHTS AND SQUARE.
C	(3) FOURIER TRANSFORM THE SQUARED MAP AND FIND NEW PHASES
C           AND WEIGHTS. ACCEPT NEW PHASE AND WEIGHT ONLY IF NEW
C           WEIGHT IS GREATER THAN PREVIOUS WEIGHT.
C	(4) COMPARE WITH PREVIOUS PHASES.  IF MEAN CHANGE IS LESS
C	    THAN 2 DEGREES OR IF THIS IS SIXTEENTH REFINEMENT CYCLE
C	    THEN GOTO (5).   OTHERWISE RETURN TO 2.
C	(5) RETURN TO (1) UNLESS NUMBER OF TRIALS EQUALS NTRIAL.
C       (6) DO ONE LAST CYCLE ACCEPTING NEW PHASES REGARDLESS OF
C           WEIGHTS.      
C	(6) OUTPUT FIGURES OF MERIT INCLUDING COMBINED FOM.
C       (7) THE MAP FOR THE BEST FOM IS OUTPUT
C	(8) THEREAFTER USER REQUESTS E-MAP FOR DESIGNATED PHASE SETS
C  THE INITIAL INPUT IS AN OUTPUT FILE FROM FTOUE WITH E'S.
      REAL RO(0:60,0:60),E(1000),PHI(1000),FP(1000),PHIP(1000),
     +CFOM(500),A(500),Z(500),WAIT(1000),EP(1000),WATE(1000)
      INTEGER H(1000),K(1000),HMAX,LIST(0:60)
      CHARACTER ANS*3,FNAME*10,DIR*1,BT*6
      CHARACTER PGG*(*),PG*4
      PARAMETER (PGG='p1  p2  pm  pg  cm  pmm pmg pgg cmm p4  p4m p4g p3
     +  p3m1p31mp6  p6m ') 
      COMMON RO,E,PHI,H,K,NX,NY,HMAX,KMAX,NOREF,WAIT
      DATA NPAGE,NTRIAL/15,50/
      OPEN(UNIT=9,FILE='LPT1')
      PI=4.0*ATAN(1.0)
      CV=PI/180.0
      VC=1/CV
      MARK=0
   80 WRITE(6,'('' IF YOU ARE RUNNING THE PROGRAMME FROM THE'')')
      WRITE(6,'('' BEGINNING THEN INPUT "B"'')')
      WRITE(6,'('' IF YOU WANT TO PRINT A MAP FROM A PREVIOUS'')')
      WRITE(6,'('' RUN THEN INPUT "P" '')')
      READ(5,50)ANS
      IF(ANS.EQ.'B'.OR.ANS.EQ.'b')GOTO 79
      IF(ANS.EQ.'P'.OR.ANS.EQ.'p')THEN
      MARK=1
      GOTO 78
      ENDIF
      GOTO 80
   78 WRITE(6,'('' READ IN THE NUMBER OF THE PHASE SET FROM THE'')')
      WRITE(6,'('' PREVIOUS RUN YOU WISH TO PRINT '')')
      READ(5,*)IBEST
      GOTO 76
C ENTER SPACE GROUP AND CELL DATA
   79 WRITE(6,'('' INPUT X-RAY WAVELENGTH IN ANGSTROMS '')')
      READ(5,*)XLAM
      WRITE(6,'('' INPUT CELL DIMENSIONS A,B [IN ANGSTROMS '')')
      WRITE(6,'('' AND GAMMA IN DEGREES '')')
      READ(5,*)AC,BC,GAMMA
C  READ IN PLANE GROUP SYMBOL AND TRANSLATE INTO IDENTIFYING INTEGER
  107 WRITE(6,'('' INPUT PLANE GROUP - ONE OF THE FOLLOWING '')')
      WRITE(6,'('' p1, p2, pm, pg, pgg, cm, pmm, pmg, cmm, p4,'')')
      WRITE(6,'('' p4m, p4g, p3, p3m1, p31m, p6, p6m '')')
      WRITE(6,'('' Note: use lower case symbols as shown '')')
      READ(5,51)PG
   51 FORMAT(A4)
C  FORM IDENTIFYING INTEGER 1 TO 17 FOR PLANE GROUP
      NSG=INDEX(PGG,PG)+3
      NPG=NSG/4      
      IF(NPG*4.NE.NSG)GOTO 107
      WRITE(6,'('' READ IN NUMBER OF ATOMS IN THE UNIT CELL'')')
      READ(5,*)NATOM
      FAC=2.0/SQRT(FLOAT(NATOM))
      FNAME='UE.DAT'
    3 WRITE(6,'(''IT IS ASSUMED THAT THE INPUT FILE FROM FTOUE'')')
      WRITE(6,'(''IS NAMED UE.DAT. DO YOU WANT TO CHANGE IT? Y/N'')')        
      READ(5,50)ANS
   50 FORMAT(A3)
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 1
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 2
      GOTO 3
    2 WRITE(6,'(''INPUT THE REQUIRED NAME.'')')
      READ(5,100)FNAME
  100 FORMAT(A10)
    1 OPEN(UNIT=10,FILE=FNAME)
   76 FNAME='DIRECT.DAT'
  103 WRITE(6,'(''IT IS ASSUMED THAT THE PROGRAMME OUTPUT FILE IS'')')
      WRITE(6,'(''NAMED DIRECT.DAT. DO YOU WANT TO CHANGE IT? Y/N'')')        
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 101
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 102
      GOTO 103
  102 WRITE(6,'(''INPUT THE REQUIRED NAME.'')')
      READ(5,100)FNAME
  101 OPEN(UNIT=11,FILE=FNAME)
      REWIND(11)
      IF(MARK.EQ.1)GOTO 177 
      READ(10,*)IUE,NOREF,(H(I),K(I),E(I),I=1,NOREF)  
      IF(IUE.EQ.1)THEN
      WRITE(6,'(''U VALUES HAVE BEEN INPUT INSTEAD OF E'')')
      STOP
      ENDIF
C  END OF DATA INPUT.      
C  STORE ALL VALUES OF E IN EP AND SELECT SUBSET OF LARGE E'S
C  ALSO LOCATE POSITION OF E(0)
      EMAX=0
      DO 129 I=1,NOREF
      EP(I)=E(I)
      IF(E(I).GT.EMAX)THEN
      EMAX=E(I)
      M0=I
      ENDIF
      IF(E(I).LT.1.0)E(I)=0
  129 CONTINUE
      E(M0)=0.5*E(M0)
C  CALCULATE THE INTERVALS FOR E-MAP
      HMAX=0
      KMAX=0
      DO 5 I=1,NOREF
      IF(ABS(H(I)).GT.HMAX)HMAX=ABS(H(I))
      IF(ABS(K(I)).GT.KMAX)KMAX=ABS(K(I))
    5 CONTINUE
      IF(HMAX.GT.15)HMAX=15
      IF(KMAX.GT.15)KMAX=15
      NX=4*HMAX
      NY=4*KMAX
C  ALLOCATE RANDOM PHASES. A SEED IS REQUIRED FOR THE RANDOM
C  NUMBER GENERATOR.
      WRITE(6,'('' INPUT A SEED FOR THE RANDOM NUMBER GENERATOR'')')
      WRITE(6,'('' AS AN INTEGER BETWEEN 0 AND 10,000'')')
      READ(5,*)IDUM
      ISEED=IDUM
      IDUM=-IDUM
C  NTRIAL SETS OF PHASES WILL BE GENERATED      
      DO 10 ITRIAL=1,NTRIAL
      IEND=0
C  GENERATE INITIAL SET OF RANDOM PHASES
      CALL RANFAZ(NPG,H,K,NOREF,PHI,IDUM,M0) 
C  INITIALISE WEIGHTS TO 0.25
      DO 555 I=1,NOREF
      WAIT(I)=0.25
  555 CONTINUE
      ICYCLE=0
   70 CALL EMAP
C  THE E-MAP HAS BEEN CALCULATED.  IT IS NOW SQUARED.
      DO 12 I=0,NX
      DO 12 J=0,NY
      RO(I,J)=RO(I,J)**2
   12 CONTINUE 
      CALL FTMAP(FP,PHIP)        
C  THE STRUCTURE AMPLITUDES OF THE SQUARED MAP ARE IN ARRAY FP
C  AND THE PHASES IN ARRAY PHIP.  CALCULATE WEIGHTS FOR THE
C  NEW PHASE ESTIMATES.   FIRST A FACTOR MUST BE FOUND WHICH
C  CONVERTS COEFFICIENTS FP TO THOSE OF SQUARED DENSITY.
C  NOW CALCULATE NEXT PHASES UNLESS IT IS FINAL CYCLE
      SUM=0
      DO 536 I=1,NOREF
      SUM=SUM+(WAIT(I)*E(I))**2
  536 CONTINUE
      SCALE=SUM/FP(M0)
      DO 537 I=1,NOREF
      FP(I)=SCALE*FP(I)
  537 CONTINUE
      RF=4.0
C  ADJUST THE CONSTANT RF UNTIL AT LEAST ONE PHASE IS CHANGED.
C  HOWEVER, RF HAS A MINIMUM VALUE OF 2.5 
  571 INEW=0
      DO 538 I=1,NOREF    
      WATE(I)=FAC*E(I)*FP(I)/RF
      IF(WATE(I).GT.1.0)WATE(I)=1.0
      IF(WATE(I).GT.WAIT(I))INEW=INEW+1
  538 CONTINUE
      IF(INEW.EQ.0.AND.RF.GT.2.51)THEN
      RF=RF-0.5
      GOTO 571
      ENDIF 
C  PHASES CHANGED IF WATE > WAIT. FOR FINAL CYCLE ALL WEIGHTS AND
C  PHASES ARE CHANGED.
      DO 667 I=1,NOREF
      IF(IEND.EQ.1)THEN
      WAIT(I)=WATE(I)
      GOTO 667
      ENDIF
      IF(WATE(I).GT.WAIT(I).OR.WATE(I).GT.0.999999)THEN
      WAIT(I)=WATE(I)
      ELSE
      PHIP(I)=PHI(I)
      ENDIF
  667 CONTINUE
      IF(IEND.EQ.1)GOTO 115      
C  CHECK THAT THERE HAVE BEEN AT LEAST FIVE CYCLES
      IF(ICYCLE.LE.5)THEN
      DO 159 I=1,NOREF
      PHI(I)=PHIP(I)
  159 CONTINUE
      GOTO 370
      ENDIF      
C  CHECK FOR CONDITION THAT EITHER THERE HAVE BEEN 16 CYCLES
C  OF REFINEMENT OR THAT MEAN PHASE SHIFT IS LESS THAN 2 DEGREE
      SUM=0.0
      DO 13 I=1,NOREF
      DEL=ABS(PHI(I)-PHIP(I))
      IF(DEL.GT.180.0)DEL=360.0-DEL
      SUM=SUM+DEL
      PHI(I)=PHIP(I)
   13 CONTINUE
      IF(SUM/NOREF.LE.2.0)GOTO 60
  370 ICYCLE=ICYCLE+1
      IF(ICYCLE.EQ.16)GOTO 60       
      GOTO 70
C  NOW DO A FINAL CYCLE WHERE ALL THE WEIGHTS AND PHASES ARE CHANGED.
   60 IEND=1
      GOTO 70   
C  CALCULATE THE FIGURES OF MERIT FOR THIS TRIAL.   
  115 A(ITRIAL)=0.0
      Z(ITRIAL)=0
      DO 637 I=1,NOREF
      IF(E(I).GT.1.3)A(ITRIAL)=A(ITRIAL)+FP(I)*WAIT(I)
      IF(EP(I).LT.0.5)Z(ITRIAL)=Z(ITRIAL)+FP(I)
  637 CONTINUE
C  WRITE THE SOLUTION ON AN OUTPUT TAPE. 
      WRITE(11,*)NOREF,(H(I),K(I),E(I),WAIT(I),PHI(I),I=1,NOREF)
   10 CONTINUE
C  ALL SETS OF PHASES HAVE NOW BEEN GENERATED TOGETHER WITH FOMS.
C  NORMALIZE THE FOMS TO MAXIMUM OF 1000
      AMAX=0
      ZMAX=0
      AMIN=1.0E8
      ZMIN=1.0E8
      DO 737 I=1,NTRIAL
      IF(A(I).GT.AMAX)AMAX=A(I)
      IF(Z(I).GT.ZMAX)ZMAX=Z(I)
      IF(A(I).LT.AMIN)AMIN=A(I)
      IF(Z(I).LT.ZMIN)ZMIN=Z(I)
  737 CONTINUE
      DO 837 I=1,NTRIAL
      A(I)=1000.0*A(I)/AMAX
      Z(I)=1000.0*Z(I)/ZMAX
  837 CONTINUE            
      AMIN=1000.0*AMIN/AMAX
      ZMIN=1000.0*ZMIN/ZMAX
C  CALCULATE COMBINED FIGURE OF MERIT
      DO 937 I=1,NTRIAL
      CFOM(I)=(A(I)-AMIN)/(1000.0-AMIN)+(1000.0-Z(I))/(1000.0-ZMIN)
  937 CONTINUE
      WRITE(9,231)ISEED
  231 FORMAT(17H THE SEED USED IS,I7)
      WRITE(9,200)
  200 FORMAT(2(28H  SOLN    A       Z     CFOM,4X))
      WRITE(9,300)(I,A(I),Z(I),CFOM(I),I=1,NTRIAL)      
  300 FORMAT(2(I6,2F8.1,F6.3,4X))
      WRITE(9,'('' '')')
      WRITE(9,'('' '')')
      CBEST=0
      DO 65 I=1,NTRIAL
      IF(CFOM(I).GT.CBEST)THEN
      CBEST=CFOM(I)
      IBEST=I
      ENDIF
   65 CONTINUE
  177 REWIND(11)
      DO 66 I=1,IBEST
      READ(11,*)NOREF,(H(II),K(II),E(II),WAIT(II),PHI(II),II=1,NOREF)
   66 CONTINUE
      DO 715 I=1,NOREF
      IF(ABS(H(I)).GT.HMAX)HMAX=ABS(H(I))
      IF(ABS(K(I)).GT.KMAX)KMAX=ABS(K(I))
  715 CONTINUE
      IF(HMAX.GT.15)HMAX=15
      IF(KMAX.GT.15)KMAX=15
      NX=4*INT(0.75*HMAX+0.5)
      NY=4*INT(0.75*KMAX+0.5)
      CALL EMAP
C  SCALE THE E-MAP READY FOR OUTPUT
      ROMAX=0
      DO 67 I=0,NX
      DO 67 J=0,NY
      IF(ABS(RO(I,J)).GT.ROMAX)ROMAX=ABS(RO(I,J))
   67 CONTINUE
      DO 68 I=0,NX
      DO 68 J=0,NY
      RO(I,J)=RO(I,J)*100.0/ROMAX
   68 CONTINUE
      DO 166 I=0,60
      LIST(I)=I
  166 CONTINUE
C  OUTPUT BEST OR CHOSEN MAP
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
      DO 444 IL=1,6
  444 WRITE(9,'('' '')')
      IF(MARK.EQ.1)THEN
      WRITE(9,651)IBEST
      GOTO 652
      ENDIF
  651 FORMAT(16H SOLUTION NUMBER,I4)
      WRITE(9,650)IBEST,A(IBEST),Z(IBEST),CFOM(IBEST)
  650 FORMAT(16H SOLUTION NUMBER,I4,3X,14H A, Z, CFOM = ,2F8.1,F9.3)    
  652 WRITE(9,500)DIR,BT
  500 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)
  600 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,600)IX,(RO(IX,IY),IY=N1,N2)
      WRITE(9,'(''  '')')
   26 CONTINUE
   45 STOP
      END






      SUBROUTINE EMAP
C  EMAP CALCULATES A TWO-DIMENSIONAL E-MAP FROM THE OUTPUT OF 
C  DIRECT WITH PHASE CHANGES CONTROLLED BY WEIGHTS. VALUES OF 
C  HMAX AND HMIN MUST BE <= 15 AND THE MINIMUM GRID SPACING IS 1/60.
      DIMENSION A(-15:15,0:15),B(-15:15,0:15),AC(0:15,0:30),
     +AS(0:15,0:30),BC(0:15,0:30),BS(0:15,0:30),CSX(0:60),
     +CSY(0:60),RO(0:60,0:60),H(1000),K(1000),WAIT(1000),
     +E(1000),PHI(1000)
      INTEGER H,HMAX,HH
      COMMON RO,E,PHI,H,K,NX,NY,HMAX,KMAX,NOREF,WAIT
      TPI=6.2831853
C  CONSTANT FOR CONVERTING DEGREES TO RADIANS
      DTR=0.0174533
C  CLEAR A AND B ARRAYS
      DO 1 HH=-15,15
      DO 1 KK=0,15
      A(HH,KK)=0
      B(HH,KK)=0
    1 CONTINUE
      DO 77 I=1,NOREF
      A(H(I),K(I))=E(I)*WAIT(I)*COS(DTR*PHI(I))
      B(H(I),K(I))=E(I)*WAIT(I)*SIN(DTR*PHI(I))
   77 CONTINUE   
      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,15
      DO 12 J=0,30
      AC(I,J)=0
      AS(I,J)=0
      BC(I,J)=0
      BS(I,J)=0
   12 CONTINUE
      DO 13 KK=0,KMAX
      DO 13 IX=0,NX/2
      DO 13 HH=-HMAX,HMAX
      IF(HH.EQ.0.AND.KK.EQ.0)GOTO 13
      COZ=CSX(MOD(HH*IX+NX16,NX))
      ZIN=CSX(MOD(HH*IX+NX16-NX4,NX))
      AC(KK,IX)=AC(KK,IX)+A(HH,KK)*COZ
      AS(KK,IX)=AS(KK,IX)+A(HH,KK)*ZIN            
      BC(KK,IX)=BC(KK,IX)+B(HH,KK)*COZ
      BS(KK,IX)=BS(KK,IX)+B(HH,KK)*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 E(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 E-MAP.
      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 KK=0,KMAX
      COZ=CSY(MOD(KK*IY+NY16,NY))
      ZIN=CSY(MOD(KK*IY+NY16-NY4,NY))
      RO(IX,IY)=RO(IX,IY)+W*(AC(KK,IX)+BS(KK,IX))*COZ
     ++W*(BC(KK,IX)-AS(KK,IX))*ZIN
      RO(NX-IX,IY)=RO(NX-IX,IY)+W*(AC(KK,IX)-BS(KK,IX))*COZ
     ++W*(BC(KK,IX)+AS(KK,IX))*ZIN
   16 CONTINUE
C  THE DENSITY TABLE IS NOW COMPLETE. 
      RETURN
      END
   
    
      SUBROUTINE FTMAP(FP,PHIP)
C  THIS CALCULATES THE FOURIER TRANSFORM OF THE SQUARED E-MAP - G(h).
C  NO ATTEMPT IS MADE TO SCALE THE AMPLITUDES
      REAL RO(0:60,0:60),E(1000),PHI(1000),PHIP(1000),FP(1000),
     +CSX(0:60),CSY(0:60),A(1:60,0:15),B(1:60,0:15),AF(-15:15,0:15),
     +BF(-15:15,0:15),WAIT(1000)
      INTEGER H(1000),HH,HMAX,K(1000)
      COMMON RO,E,PHI,H,K,NX,NY,HMAX,KMAX,NOREF,WAIT
      NX16=16*NX
      NX4=NX/4
      NY16=16*NY
      NY4=NY/4
      TPI=8.0*ATAN(1.0)
      CONV=360.0/TPI
C  SET UP COSINE TABLES FOR USE IN SUMMATIONS
      DO 1 I=0,NX
    1 CSX(I)=COS(TPI*I/NX)      
      DO 2 I=0,NY
    2 CSY(I)=COS(TPI*I/NY)
C  SUMMATION OVER THE Y DIRECTION FOR PARTICULAR K AND X VALUES
C  FIRST CLEAR ARRAYS
      DO 11 IX=1,NX
      DO 11 KK=0,KMAX
      A(IX,KK)=0
      B(IX,KK)=0
   11 CONTINUE
      DO 12 IX=1,NX
      DO 12 KK=0,KMAX
      DO 12 IY=1,NY
      COZ=CSY(MOD(KK*IY+NY16,NY))
      ZIN=CSY(MOD(KK*IY+NY16-NY4,NY))
      A(IX,KK)=A(IX,KK)+RO(IX,IY)*COZ
      B(IX,KK)=B(IX,KK)+RO(IX,IY)*ZIN
   12 CONTINUE
C  NOW FOR EACH H AND K DO THE SUMMATIONS OVER X
C  FIRST CLEAR ARRAYS
      DO 14 HH=-HMAX,HMAX
      DO 14 KK=0,KMAX
      AF(HH,KK)=0
      BF(HH,KK)=0
   14 CONTINUE
      DO 15 HH=-HMAX,HMAX
      DO 15 KK=0,KMAX
      DO 16 IX=1,NX
      COZ=CSX(MOD(HH*IX+NX16,NX))
      ZIN=CSX(MOD(HH*IX+NX16-NX4,NX))
      AF(HH,KK)=AF(HH,KK)+A(IX,KK)*COZ-B(IX,KK)*ZIN
      BF(HH,KK)=BF(HH,KK)+A(IX,KK)*ZIN+B(IX,KK)*COZ
   16 CONTINUE
      FF=SQRT(AF(HH,KK)**2+BF(HH,KK)**2)
      PHIF=ATAN2(BF(HH,KK),AF(HH,KK))
      AF(HH,KK)=FF
      BF(HH,KK)=PHIF*CONV
   15 CONTINUE
C  TRANSFER MAGNITUDES AND PHASES TO ARRAYS FP AND PHIP
      DO 20 I=1,NOREF             
      FP(I)=AF(H(I),K(I))
      PHIP(I)=BF(H(I),K(I))
   20 CONTINUE
      RETURN
      END







 
      SUBROUTINE RANFAZ(NPG,H,K,NOREF,PHI,IDUM,M0)
      REAL PHI(1000)
      INTEGER H(1000),K(1000),KEY(1000)
C  RANDOM PHASES ARE ALLOCATED TO ALL REFLECTIONS.  SYMMETRY-RELATED
C  REFLECTIONS ARE ALLOCATED PHASES WITH THE PROPER RELATIONSHIPS.
C  THE PROGRAMME CAN BE APPLIED TO ALL 17 TWO-DIMENSIONAL SPACE
C  GROUPS.
      PI=4.0*ATAN(1.0)
      GOTO (1,2,3,4,6,6,7,8,6,10,11,11,13,14,14,13,14)NPG 
    1 DO 101 I=1,NOREF
      XX=RAN2(IDUM)
      PHI(I)=360.00*XX
  101 CONTINUE    
      GOTO 500
    2 DO 102 I=1,NOREF
      XX=RAN2(IDUM)
      PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180
  102 CONTINUE
      GOTO 500
    3 DO 103 I=1,NOREF
      KEY(I)=0
  103 CONTINUE
      DO 203 I=1,NOREF
      IF(KEY(I).EQ.1)GOTO 203
      XX=RAN2(IDUM)
      IF(K(I).EQ.0)GOTO 303
      PHI(I)=360.0*XX
      GOTO 403
  303 PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180
  403 DO 503 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 503
      IF(H(J).EQ.-H(I).AND.K(J).EQ.K(I))GOTO 603
      GOTO 503
  603 KEY(J)=1
      PHI(J)=PHI(I)
  503 CONTINUE
  203 CONTINUE    
      GOTO 500
    4 DO 104 I=1,NOREF
      KEY(I)=0
  104 CONTINUE            
      DO 204 I=1,NOREF
      IF(KEY(I).EQ.1)GOTO 204
      XX=RAN2(IDUM)
      IF(K(I).EQ.0)GOTO 304
      PHI(I)=360.0*XX
      GOTO 404
  304 PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180
  404 DO 504 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 504
      IF(H(J).EQ.-H(I).AND.K(J).EQ.K(I))GOTO 604
      GOTO 504
  604 KEY(J)=1
      PHI(J)=PHI(I)+(1.0-1.0*(-1.0)**K(I))*90.0
  504 CONTINUE
  204 CONTINUE    
      GOTO 500
    6 DO 106 I=1,NOREF
      KEY(I)=0
  106 CONTINUE
      DO 206 I=1,NOREF
      IF(KEY(I).EQ.1)GOTO 206
      XX=RAN2(IDUM)
      PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180.0
      DO 306 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 306
      IF(H(J).EQ.-H(I).AND.K(J).EQ.K(I))GOTO 406
      GOTO 306
  406 KEY(J)=1
      PHI(J)=PHI(I)      
  306 CONTINUE
  206 CONTINUE
      GOTO 500
    7 DO 107 I=1,NOREF
      KEY(I)=0
  107 CONTINUE
      DO 207 I=1,NOREF
      IF(KEY(I).EQ.1)GOTO 207
      XX=RAN2(IDUM)
      PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180.0
      DO 307 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 307
      IF(H(J).EQ.-H(I).AND.K(J).EQ.K(I))GOTO 407
      GOTO 307
  407 KEY(J)=1
      PHI(J)=PHI(I)+(1.0-1.0*(-1.0)**H(I))*90.0      
  307 CONTINUE
  207 CONTINUE
      GOTO 500
    8 DO 108 I=1,NOREF
      KEY(I)=0
  108 CONTINUE
      DO 208 I=1,208
      IF(KEY(I).EQ.1)GOTO 208
      XX=RAN2(IDUM)
      PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180.0
      DO 308 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 308
      IF(H(J).EQ.-H(I).AND.K(J).EQ.K(I))GOTO 408
      GOTO 308
  408 KEY(J)=1
      PHI(J)=PHI(I)+(1.0-1.0*(-1.0)**(H(I)+K(I)))*90.0      
  308 CONTINUE
  208 CONTINUE
      GOTO 500 
   10 DO 110 I=1,NOREF
      KEY(I)=0
  110 CONTINUE
      DO 210 I=1,210
      IF(KEY(I).EQ.1)GOTO 210
      XX=RAN2(IDUM)
      PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180.0
      DO 310 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 310
      IF(H(J).EQ.-K(I).AND.K(J).EQ.H(I))GOTO 410
      IF(H(J).EQ.K(I).AND.K(J).EQ.-H(I))GOTO 410
      GOTO 310
  410 KEY(J)=1
      PHI(J)=PHI(I)  
  310 CONTINUE
  210 CONTINUE
      GOTO 500
   11 DO 111 I=1,NOREF
      KEY(I)=0
  111 CONTINUE
      DO 211 I=1,NOREF
      IF(KEY(I).EQ.1)GOTO 211
      XX=RAN2(IDUM)
      PHI(I)=0
      IF(XX.GT.0.5)PHI(I)=180.0
      DO 311 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 311
      IF(H(J).EQ.K(I).AND.K(J).EQ.H(I))GOTO 511
      IF(H(J).EQ.-K(I).AND.K(J).EQ.-H(I))GOTO 511
      IF(H(J).EQ.-H(I).AND.K(J).EQ.K(I))GOTO 411
      IF(H(J).EQ.-K(I).AND.K(J).EQ.H(I))GOTO 411
      IF(H(J).EQ.K(I).AND.K(J).EQ.-H(I))GOTO 411
      GOTO 311
  511 KEY(J)=1
      PHI(J)=PHI(I)
      GOTO 311
  411 KEY(J)=1
      PHI(J)=PHI(I)+(1.0-1.0*(-1.0)**(H(I)+K(I)))*90.0 
  311 CONTINUE
  211 CONTINUE
      GOTO 500
   13 DO 113 I=1,NOREF
      KEY(I)=0
  113 CONTINUE
      DO 213 I=1,NOREF
      IF(KEY(I).EQ.1)GOTO 213
      XX=RAN2(IDUM)
      IF(NPG.EQ.13)PHI(I)=XX*360.0
      IF(NPG.EQ.16)PHI(I)=0
      IF(NPG.EQ.16.AND.XX.GT.0.5)PHI(I)=180.0
      DO 313 J=I+1,NOREF
      IF(KEY(J).EQ.1)GOTO 313
      IF(H(J).EQ.K(I).AND.K(J).EQ.-H(I)-K(I))GOTO 413
      IF(H(J).EQ.-K(I).AND.K(J).EQ.H(I)+K(I))GOTO 513
      IF(H(J).EQ.-H(I)-K(I).AND.K(J).EQ.H(I))GOTO 413
      IF(H(J).EQ.H(I)+K(I).AND.K(J).EQ.-H(I))GOTO 513
      GOTO 313
  413 PHI(J)=PHI(I)
      GOTO 313
  513 PHI(J)=360.0-PHI(I)
  313 CONTINUE
  213 CONTINUE
      GOTO 500 
   14 DO 114 I=1,NOREF
      KEY(I)=0
  114 CONTINUE
      DO 214 I=1,NOREF      
      IF(KEY(I).EQ.1)GOTO 214
      XX=RAN2(IDUM)
      PHI(I)=XX*360.0
      IF(NPG.EQ.19)PHI(I)=0
      IF(NPG.EQ.19.AND.XX.GT.0.5)PHI(I)=180.0
      DO 314 J=I+1,NOREF
      IF(H(J).EQ.-K(I).AND.K(J).EQ.-H(I))GOTO 414
      IF(H(J).EQ.K(I).AND.K(J).EQ.H(I))GOTO 514
      IF(H(J).EQ.K(I).AND.K(J).EQ.-H(I)-K(I))GOTO 614
      IF(H(J).EQ.-K(I).AND.K(J).EQ.H(I)+K(I))GOTO 714
      IF(H(J).EQ.H(I)+K(I).AND.K(J).EQ.-K(I))GOTO 414
      IF(H(J).EQ.-H(I)-K(I).AND.K(J).EQ.K(I))GOTO 514
      IF(H(J).EQ.H(I).AND.K(J).EQ.-H(I)-K(I))GOTO 514
      IF(H(J).EQ.-H(I).AND.K(J).EQ.H(I)+K(I))GOTO 414
      IF(H(J).EQ.H(I)+K(I).AND.K(J).EQ.-H(I))GOTO 714
      IF(H(J).EQ.-H(I)-K(I).AND.K(J).EQ.H(I))GOTO 614
      GOTO 314
  414 PHI(J)=PHI(I)
      IF(NPG.EQ.15)PHI(J)=360-PHI(I)
      GOTO 314
  514 PHI(J)=360.0-PHI(I)
      IF(NPG.EQ.15)PHI(J)=PHI(I)
      GOTO 314
  614 PHI(J)=PHI(I)
      GOTO 314
  714 PHI(J)=PHI(I)
  314 CONTINUE
  214 CONTINUE
  500 PHI(M0)=0
      RETURN
      END

C  THE FOLLOWING FUNCTION ROUTINE WHICH GENERATES RANDOM NUMBERS IN
C  THE RANGE 0.0 TO 1.0 COMES FROM "NUMERICAL RECIPES" BY W.H.PRESS, 
C  B.P.FLANNERY, S.A.TEUKOLSKY AND W.T.VETTERLING, CAMBRIDGE 
C  UNIVERSITY PRESS, 1986.     
      FUNCTION RAN2(IDUM)
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1.0/M)
      DIMENSION IR(97)
      DATA IFF/0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
      IFF=1
      IDUM=MOD(IC-IDUM,M)
      DO 1 J=1,97
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
    1 CONTINUE
      IDUM=MOD(IA*IDUM+IC,M)
      ENDIF
      IY=IDUM
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1)PAUSE
      IY=IR(J)
      RAN2=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
      

