      PROGRAM ISOCOEFF
C  THIS PROGRAMME CALCULATES THE FOURIER COEFFICIENTS AND PHASES FOR
C  A DOUBLE FOURIER SUMMATION GIVEN THAT THE POSITIONS OF I-R ATOMS 
C  AND THE STRUCTURE AMPLITUDES FOR THE TWO ISOMORPHOUS STRUCTURES 
C  ARE KNOWN.   THE SIMULATION REQUIRES:
C     (i)   STRUCFAC OUTPUT FILE FOR STRUCTURE A
C     (ii)  STRUCFAC OUTPUT FILE FOR STRUCTURE B
C     (iii) STRUCFAC OUTPUT FILE FOR I-R ATOMS A
C     (iv)  STRUCFAC OUTPUT FILE FOR I-R ATOMS B
C  THE OUTPUT FILE WILL BE SUITABLE AS AN INPUT FILE FOR FOUR2, THE
C  PRINT-OUT OF WHICH WILL BE THE DOUBLE FOURIER SUMMATION.
      REAL FA(1000),FB(1000),FIRA(1000),FIRB(1000),PHI(1000)
      INTEGER H(1000),K(1000)
      CHARACTER ANS*3,FNAME*10
C  SET NAME OF INPUT FILE FOR STRUCTURE A
      FNAME='SF1.DAT'
    3 WRITE(6,'('' IT IS ASSUMED THAT THE STRUCTURE FACTOR '')')
      WRITE(6,'('' INFORMATION FOR STRUCTURE A IS IN FILE SF1.DAT.'')')
      WRITE(6,'('' DO YOU WANT TO CHANGE THIS?   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,'(''READ IN NAME OF FILE [<= 10 CHARACTERS] '')')
      READ(5,100)FNAME
  100 FORMAT(A10)
    1 OPEN(UNIT=10,FILE=FNAME)    
C  SET NAME OF INPUT FILE FOR STRUCTURE B 
      FNAME='SF2.DAT'
   13 WRITE(6,'('' IT IS ASSUMED THAT THE STRUCTURE FACTOR '')')
      WRITE(6,'('' INFORMATION FOR STRUCTURE B IS IN FILE SF2.DAT.'')')
      WRITE(6,'('' DO YOU WANT TO CHANGE THIS?   Y/N  '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 11
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 12
      GOTO 13
   12 WRITE(6,'(''READ IN NAME OF FILE [<= 10 CHARACTERS] '')')
      READ(5,100)FNAME
   11 OPEN(UNIT=11,FILE=FNAME)    
C  SET NAME OF INPUT FILE FOR I-R ATOMS A
      FNAME='SFIRA.DAT'
   23 WRITE(6,'('' IT IS ASSUMED THAT THE STRUCTURE FACTOR '')')
      WRITE(6,'('' INFORMATION FOR I-R ATOMS A IS IN FILE '')')
      WRITE(6,'('' SFIRA.DAT. DO YOU WANT TO CHANGE THIS?   Y/N  '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 21
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 22
      GOTO 23
   22 WRITE(6,'(''READ IN NAME OF FILE [<= 10 CHARACTERS] '')')
      READ(5,100)FNAME
   21 OPEN(UNIT=12,FILE=FNAME)    
C  SET NAMES OF INPUT FILE FOR I-R ATOMS B
      FNAME='SFIRB.DAT'
   33 WRITE(6,'('' IT IS ASSUMED THAT THE STRUCTURE FACTOR '')')
      WRITE(6,'('' INFORMATION FOR I-R ATOMS B IS IN FILE '')')
      WRITE(6,'('' SFIRB.DAT. DO YOU WANT TO CHANGE THIS?   Y/N  '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 31
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 32
      GOTO 33
   32 WRITE(6,'(''READ IN NAME OF FILE [<= 10 CHARACTERS] '')')
      READ(5,100)FNAME
   31 OPEN(UNIT=13,FILE=FNAME)    
C  SET NAME OF OUTPUT FILE
      FNAME='SDF.DAT'
   43 WRITE(6,'('' IT IS ASSUMED THAT THE OUTPUT FILE FOR THE  '')')
      WRITE(6,'('' DOUBLE FOURIER SUMMATION WILL BE IN SDF.DAT. '')')
      WRITE(6,'('' DO YOU WANT TO CHANGE THIS?   Y/N  '')')
      READ(5,50)ANS
      IF(ANS.EQ.'N'.OR.ANS.EQ.'n')GOTO 41
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y')GOTO 42
      GOTO 43
   42 WRITE(6,'(''READ IN NAME OF FILE [<= 10 CHARACTERS] '')')
      READ(5,100)FNAME
   41 OPEN(UNIT=14,FILE=FNAME)    
C  READ IN ALL DATA FROM INPUT FILES
      READ(10,*)NOREF,(H(I),K(I),FA(I),PHI(I),I=1,NOREF)
      READ(11,*)NOREF,(H(I),K(I),FB(I),PHI(I),I=1,NOREF)
      READ(12,*)NOREF,(H(I),K(I),FIRA(I),PHI(I),I=1,NOREF)
      READ(13,*)NOREF,(H(I),K(I),FIRB(I),PHI(I),I=1,NOREF)
C  CALCULATE VALUES OF K*COS(DELTA PHI)
      DO 10 I=1,NOREF
C  APPLY EQUATION (8.22)
      CPSI=((FIRB(I)-FIRA(I))**2+FB(I)**2-FA(I)**2)/2.0/FB(I)
     +/(FIRB(I)-FIRA(I))
C  CHECK FOR COS(PSI) WITHIN PROPER LIMITS.  IF NOT BRING TO LIMIT
      IF(ABS(CPSI).GT.1.0)CPSI=SIGN(1.0,CPSI)    
C  APPLY EQUATION (8.23)      
      XK=SQRT(FB(I)**2+FIRB(I)**2-2.0*FB(I)*FIRB(I)*CPSI)
C  APPLY EQUATION (8.24) TO FIND COS(DELTA PHI)
      CDPHI=(FB(I)**2-XK**2-FIRB(I)**2)/2.0/XK/FIRB(I)
      FA(I)=XK*CDPHI
   10 CONTINUE
C  WRITE OUTPUT FILE
      WRITE(14,*)NOREF,(H(I),K(I),FA(I),PHI(I),I=1,NOREF)
      STOP
      END

      


