      PROGRAM HEAVY
C  THIS PROGRAMME ALLOWS A SIMULATION OF A HEAVY-ATOM-METHOD SOLUTION
C  OF EITHER A CENTROSYMMETRIC OR NON-CENTROSYMMETRIC STRUCTURE.  
C  DATA FILES ARE READ IN [SF1.DAT & SF2.DAT OR OTHER CHOSEN NAMES]
C  CONTAINING THE OUTPUT FROM "STRUCFAC" FOR THE WHOLE STRUCTURE AND
C  FOR THE HEAVY ATOMS RESPECTIVELY.   THE OUTPUT FILE, SUITABLE FOR 
C  INPUT TO "FOUR2", HAS THE PHASE OF THE HEAVY ATOM CONTRIBUTION AND 
C  MAGNITUDE OF THE WHOLE STRUCTURE MODIFIED BY EITHER THE SIM WEIGHT
C  FOR A NON-CENTROSYMMETRIC STRUCTURE OR THAT GIVEN BY WOOLFSON FOR
C  A CENTROSYMMETRIC STRUCTURE.
C  THERE IS ALSO AN OPTION FOR ACCEPTING THE HEAVY ATOM PHASE WITHOUT
C  WEIGHTING THE MAGNITUDE.
      REAL F1(1000),F2(1000),PHI(1000)
      INTEGER H(1000),K(1000),NATOM(14),MTYPE(14)
      CHARACTER FNAME*10,ANS*1,ATYPE*3,ATOM*(*)
      PARAMETER (ATOM='H  LI C  N  O  F  NA AL P  S  CL FE BR HG END')
      FNAME='SF1.DAT'
    3 WRITE(6,'('' IT IS ASSUMED THAT THE FILE CONTAINING THE DATA'')')
      WRITE(6,'('' FOR THE WHOLE STRUCTURE IS CALLED SF1.DAT.  DO '')')
      WRITE(6,'('' YOU WANT TO CHANGE IT?   Y/N '')')
      READ(5,50)ANS
   50 FORMAT(A1)
      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,'('' READ IN THE REQUIRED FILE NAME '')')
      READ(5,100)FNAME
    1 OPEN(UNIT=10,FILE=FNAME)
      FNAME='SF2.DAT'
    6 WRITE(6,'('' IT IS ASSUMED THAT THE FILE CONTAINING THE DATA'')')
      WRITE(6,'('' FOR THE HEAVY ATOM STRUCTURE IS CALLED SF2.DAT.'')')
      WRITE(6,'('' DO YOU WANT TO CHANGE IT?   Y/N '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 4
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 5
      GOTO 6
    5 WRITE(6,'('' READ IN THE REQUIRED FILE NAME '')')
      READ(5,100)FNAME
    4 OPEN(UNIT=11,FILE=FNAME)  
      FNAME='SF.DAT'
    9 WRITE(6,'('' IT IS ASSUMED THAT THE FILE CONTAINING OUTPUT'')')
      WRITE(6,'('' DATA FOR INPUT TO FOUR2 IS CALLED SF.DAT.  DO '')')
      WRITE(6,'('' YOU WANT TO CHANGE IT?   Y/N '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 7
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 8
      GOTO 9
    8 WRITE(6,'('' READ IN THE REQUIRED FILE NAME '')')
      READ(5,100)FNAME
  100 FORMAT(A10)
    7 OPEN(UNIT=12,FILE=FNAME)
   46 WRITE(6,'('' DO YOU WANT TO TAKE THE STRUCTURE AMPLITUDE'')')
      WRITE(6,'('' WITHOUT WEIGHT?  Y/N '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 12
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')THEN
      KEY=1
      GOTO 13
      ENDIF
      GOTO 46 
   12 WRITE(6,'('' INPUT C FOR A CENTROSYMMETRIC STRUCTURE '')')
      WRITE(6,'('' OR N FOR A NON-CENTROSYMMETRIC STRUCTURE '')')
      READ(5,50)ANS
      IF(ANS.EQ.'C'.OR.ANS.EQ.'c')GOTO 10
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 11
      GOTO 12
   10 KEY=2
      GOTO 13
   11 KEY=3
   13 READ(10,*)NOREF,(H(I),K(I),F1(I),PHI(I),I=1,NOREF)
      READ(11,*)NOREF,(H(I),K(I),F2(I),PHI(I),I=1,NOREF)
      WRITE(6,'('' INPUT CELL DIMENSIONS A, B IN ANGSTROMS '')')
      WRITE(6,'('' AND GAMMA IN DEGREES '')')
      READ(5,*)A,B,GAMMA
C  CO1, CO2 AND CO3 ARE USED FOR CALCULATING SIN(THETA)/LAMBDA
      CO1=0.25/A/A
      CO2=0.25/B/B
      CO3=0.5/A/B*COS(GAMMA*0.0174533)
      WRITE(6,'(''THERE ARE NOW INPUT THE TYPES AND NUMBERS '')')
      WRITE(6,'(''OF ATOMS WHOSE COORDINATES ARE NOT KNOWN. '')')
      WRITE(6,'(''NOTE-THE NUMBERS ARE FOR THE WHOLE UNIT CELL!'')')      
      NOTYPE=0
   30 WRITE(6,'('' INPUT ATOMIC TYPE - ONE OF H, LI, C, N, O, F,'')')
      WRITE(6,'('' NA, AL, P, S, CL, FE, BR, HG OR END IF ALL '')')
      WRITE(6,'('' ATOMS HAVE BEEN READ IN '')')
      READ(5,150)ATYPE
  150 FORMAT(A3)
C  GENERATE A NUMBER FROM 1 TO 14 TO IDENTIFY THE TYPE OF ATOM 
C  OR 15 IF ALL ATOMS HAVE BEEN READ IN
      NX=INDEX(ATOM,ATYPE)+2
      NASYMB=NX/3
      IF(NASYMB*3.NE.NX)GOTO 30                
      IF(NASYMB.EQ.15)GOTO 40
      NOTYPE=NOTYPE+1
      MTYPE(NOTYPE)=NASYMB
      WRITE(6,'('' INPUT THE NUMBER OF THIS TYPE OF ATOM IN THE'')')
      WRITE(6,'('' WHOLE UNIT CELL '')')
      READ(5,*)NATOM(NOTYPE)
      GOTO 30
   40 DO 31 I=1,NOREF
C  CALCULATE THE VALUE OF (SIN(THETA)/LAMBDA)
      SL=SQRT(H(I)*H(I)*CO1+K(I)*K(I)*CO2-H(I)*K(I)*CO3)
      SIG=0
      DO 32 J=1,NOTYPE
      CALL SCAT(MTYPE(J),SL,SCF)
      SIG=SIG+NATOM(J)*SCF*SCF
   32 CONTINUE
      X=2*F1(I)*F2(I)/SIG
      GOTO(33,34,35)KEY
C  THIS GIVES UNWEIGHTED STRUCTURE AMPLITUDE UNLESS THE HEAVY ATOM
C  CONTRIBUTION IS ZERO WHEN THE TERM IS REMOVED FROM THE SYNTHESIS.
   33 FACTOR=1
      IF(F2(I).LT.1.0E-6)FACTOR=0
      GOTO 36
C  THIS GIVES THE WEIGHT GIVEN BY WOOLFSON
   34 FACTOR=EXP(X)
      FACTOR=FACTOR/(1.0+FACTOR)
      GOTO 36
C  THERE FOLLOWS A CALCULATION OF THE SIM WEIGHT.  IF X&gt;6 THEN IT 
C  IS EFFECTIVELY 1.  FOR SMALLER VALUES OF X A POLYNOMIAL
C  APPROXIMATION IS USED.
   35 IF(X.GT.6.0)THEN
      FACTOR=1
      GOTO 36
      ENDIF
      FACTOR=0.5658*X-0.1304*X*X+0.0106*X*X*X
   36 F1(I)=F1(I)*FACTOR
   31 CONTINUE
      WRITE(12,*)NOREF,(H(I),K(I),F1(I),PHI(I),I=1,NOREF)
      STOP
      END
   


      SUBROUTINE SCAT(I,SL,SCF)
C  THERE NOW FOLLOW THE SCATTERING FACTORS FROM INTERNATIONAL TABLES
C  FOR ATOMS H, LI, C, N, O, NA, AL, P, S, CL, FE, BR AND HG AT 
C  INTERVALS OF 0.1 FOR SIN(THETA)/LAMBDA. PARABOLIC INTERPOLATION IS
C  USED TO OBTAIN THE SCATTERING FACTOR FOR SIN(THETA)/LAMBDA = SL.
      DIMENSION SF(14,0:10)
      DATA (SF(1,I),I=0,10)/1.000,0.811,0.481,0.251,0.130,0.071,
     +0.040,0.024,0.015,0.010,0.007/
      DATA (SF(2,I),I=0,10)/3.000,2.215,1.741,1.521,1.269,1.032,
     +0.823,0.650,0.513,0.404,0.320/
      DATA (SF(3,I),I=0,10)/6.000,5.126,3.581,2.502,1.950,1.685,
     +1.536,1.426,1.322,1.218,1.114/
      DATA (SF(4,I),I=0,10)/7.000,6.203,4.600,3.241,2.397,1.944,
     +1.698,1.550,1.444,1.350,1.263/
      DATA (SF(5,I),I=0,10)/8.000,7.250,5.634,4.094,3.010,2.338,
     +1.944,1.714,1.566,1.462,1.374/
      DATA (SF(6,I),I=0,10)/9.000,8.293,6.691,5.044,3.760,2.878,
     +2.312,1.958,1.735,1.587,1.481/
      DATA (SF(7,I),I=0,10)/11.00,9.76,8.34,6.89,5.47,4.29,3.40,
     +2.76,2.31,2.00,1.78/
      DATA (SF(8,I),I=0,10)/13.00,11.23,9.16,7.88,6.72,5.69,4.71,
     +3.88,3.21,2.71,2.32/
      DATA (SF(9,I),I=0,10)/15.00,13.17,10.34,8.59,7.54,6.67,5.83,
     +5.02,4.28,3.64,3.11/
      DATA (SF(10,I),I=0,10)/16.00,14.33,11.21,8.99,7.83,7.05,6.31,
     +5.56,4.82,4.15,3.56/
      DATA (SF(11,I),I=0,10)/17.00,15.33,12.00,9.44,8.07,7.29,6.64,
     +5.96,5.27,4.60,4.00/
      DATA (SF(12,I),I=0,10)/26.00,23.68,20.09,16.77,13.84,11.47,
     +9.71,8.47,7.60,6.99,6.51/
      DATA (SF(13,I),I=0,10)/35.00,32.43,27.70,23.82,20.84,18.27,
     +15.91,13.78,11.93,10.41,9.19/
      DATA (SF(14,I),I=0,10)/80.00,75.48,67.14,59.31,52.65,47.04,
     +42.31,38.22,34.64,31.43,28.59/
      J=INT(SL/0.1)
      IF(J.EQ.10)J=9
      IF(J.EQ.0)J=1
      A=(SF(I,J+1)+SF(I,J-1)-2*SF(I,J))/0.02
      B=(SF(I,J+1)-SF(I,J-1))/0.2
      SCF=A*(SL-0.1*J)**2+B*(SL-0.1*J)+SF(I,J)
      RETURN
      END
      

