      PROGRAM FTOUE
C  THIS PROGRAMME ACCEPTS THE OUTPUT FROM STRUCFAC AND CONVERTS THE
C  F's INTO EITHER U's OR E's.   THE F's ARE FIRST SORTED INTO GROUPS
C  OF ROUGHLY 40 IN ORDER OF INCREASING VALUE OF SIN(THETA).   FOR 
C  EACH GROUP A CONVERSION FACTOR IS FOUND WHICH WILL, FOR U's, MAKE
C  <|U**2|>= SIGMA(f**2)/(SIGMA(f))**2 OR, FOR E's MAKE <|E**2|> = 1.
C  THESE FACTORS ARE THEORETICALLY WHAT IS REQUIRED TO CONVERT F**2 TO 
C  U**2 OR E**2 AT <SIN(THETA)> FOR EACH GROUP.   A BEST FIT OF THE
C  KIND   " FACTOR = C(1)+C(2)*s+C(3)*s**2 " [s = SIN(THETA)]  IS FOUND
C  TO FIT THE CONVERSION FACTORS AND <SIN(THETA)> FOR EACH GROUP.   
C  THE QUADRATIC FORMULA IS THEN USED WITH THE S FOR EACH REFLECTION
C  TO CONVERT THE VALUE OF F TO U OR E.   THESE ARE THEN PLACED IN A
C  FILE UE.DAT FOR SUBSEQUENT USE.
      REAL F(1000),THETLIST(100),F2LIST(100),AVSUMF2(30),AVTHET(30),
     +EPS(30),C(3),UE(1000)
      INTEGER ISYMB(13),NATOM(13),HT(1000),KT(1000),NUMLIST(100),
     +IEND(0:30)
      CHARACTER Q*4,FNAME*10,ASYMB*3,CU*1,ATOM*(*)
      PARAMETER (ATOM='H  LI C  N  O  NA AL P  S  CL FE BR HG END')
      OPEN(UNIT=9,FILE='LPT1')
      FNAME='UE.DAT'
      WRITE(6,'('' THE OUTPUT FILE WHICH CONTAINS H,K,U (or E) IS'')')
      WRITE(6,'('' NAMED UE.DAT.   DO YOU WANT TO CHANGE IT? (Y/N)'')')
    3 READ(5,50)Q
      IF(Q.EQ.'Y'.OR.Q.EQ.'y')GOTO 1
      IF(Q.EQ.'N'.OR.Q.EQ.'n')GOTO 2
      GOTO 3
    1 WRITE(6,'('' READ IN THE FILENAME <= 10 CHARACTERS] '')')
      READ(5,75)FNAME
   50 FORMAT(A4)
   75 FORMAT(A10)  
    2 OPEN(UNIT=10,FILE=FNAME)
C  INPUT THE UNIT CELL DATA
      WRITE(6,'('' INPUT CELL DIMENSIONS A, B [IN ANGSTROMS] AND'')')
      WRITE(6,'('' GAMMA IN DEGREES. '')')
      READ(5,*)A,B,GAMMA
C  NOW CALCULATE FACTORS WHICH WILL BE USED LATER TO FIND VALUES
C  OF SIN[THETA]/LAMDA
      CO1=0.25/A/A
      CO2=0.25/B/B
      CO3=0.5/A/B*COS(GAMMA*0.0174533)
      WRITE(6,'('' INPUT THE RADIATION WAVELENGTH IN ANGSTROMS'')')
      READ(5,*)XLAM 
C  DECIDE WHETHER TO OUTPUT U's OR E's.
      WRITE(6,'('' DO YOU WANT U OR E?   (U/E) '')')
      READ(5,50)Q
      IF(Q.EQ.'U'.OR.Q.EQ.'u')THEN
      IUE=1
      GOTO 4
      ENDIF
      IF(Q.EQ.'E'.OR.Q.EQ.'e')THEN
      IUE=2
      GOTO 4
      ENDIF
C  INPUT THE ATOMIC CONTENTS OF THE WHOLE UNIT CELL
    4 WRITE(6,'('' THE CONTENTS OF THE WHOLE UNIT CELL [NOT JUST'')')
      WRITE(6,'('' ONE ASYMMETRIC UNIT] IS INPUT IN THE FORM X1,N1'')')
      WRITE(6,'('' X2,N2,ETC WHERE X1 IS THE TYPE OF ATOM, WHICH '')')
      WRITE(6,'('' MUST BE ONE OF H,LI,C,N,O,NA,AL,P,S,CL,FE,BR,HG,'')')
      WRITE(6,'('' AND N1 IS THE TOTAL NUMBER OF THAT TYPE. '')')
      WRITE(6,'(''  '')')
      ITOT=0
    6 WRITE(6,'('' INPUT ATOMIC TYPE'')')
      WRITE(6,'('' IF ALL ATOMS HAVE BEEN INPUT THEN INPUT END'')')
      READ(5,100)ASYMB
  100 FORMAT(A3)
      NX=INDEX(ATOM,ASYMB)+2
      NASYMB=NX/3
      IF(NASYMB*3.NE.NX)GOTO 6               
      IF(NASYMB.EQ.14)GOTO 7
      ITOT=ITOT+1
      ISYMB(ITOT)=NASYMB
      WRITE(6,'('' INPUT NUMBER OF ATOMS OF THIS TYPE IN WHOLE CELL'')')
      READ(5,*)NATOM(ITOT)
      GOTO 6
C THE DATA FROM STRUCFAC CAN NOW BE READ IN
    7 FNAME='SF1.DAT'
   10 WRITE(6,'('' THE OUTPUT FILE FROM STRUCFAC IS ASSUMED TO BE'')')
      WRITE(6,'('' NAMED SF1.DAT. DO YOU WANT TO CHANGE THIS? (Y/N)'')')
      READ(5,50)Q
      IF(Q.EQ.'Y'.OR.Q.EQ.'y')GOTO 8
      IF(Q.EQ.'N'.OR.Q.EQ.'n')GOTO 9
      GOTO 10
    8 WRITE(6,'('' INPUT THE FILE NAME YOU WANT [<= 10 CHARACTERS]'')')
      READ(5,75)FNAME
    9 OPEN(UNIT=11,FILE=FNAME)
      READ(11,*)NOREF,(HT(I),KT(I),F(I),DUMMY,I=1,NOREF)
C  EACH REFLECTION IS READ IN AND VALUES OF NUMBER OF REFLECTIONS
C  SIGMA(F**2) AND SIGMA(SIN[THETA]) ARE FOUND FOR SIN[THETA] RANGES
C  OF 0.01. REFLECTIONS WITH SIN[THETA] LESS THAN O.O5 ARE EXCLUDED
C  SINCE THEY DISTURB THE STATISTICS.  FIRST ALL ARRAYS ARE CLEARED
      DO 66 I=1,100
      THETLIST(I)=0
      F2LIST(I)=0
      NUMLIST(I)=0
   66 CONTINUE
      DO 20 I=1,NOREF
      STHET=XLAM*SQRT(CO1*HT(I)**2+CO2*KT(I)**2-CO3*HT(I)*KT(I))         
      IF(STHET.LT.0.05)GOTO 20
      ITG=INT(100*STHET)+1
      IF(ITG.GT.100)ITG=100
      THETLIST(ITG)=THETLIST(ITG)+STHET
      F2LIST(ITG)=F2LIST(ITG)+F(I)**2
      NUMLIST(ITG)=NUMLIST(ITG)+1
   20 CONTINUE
C  NOW FIND GROUPS OF 40 OR SO FOR FINDING AVERAGES
C  THESE GROUPS WILL SLIGHTLY OVERLAP
      DO 21 I=1,30
      IEND(I)=0
   21 CONTINUE  
      IEND(0)=1
      IX=1
      IN=1
   31 SUM=0
      DO 30 I=IN,100
      SUM=SUM+NUMLIST(I)
      IF(SUM.GE.40)THEN
      IEND(IX)=I
      IX=IX+1
      IN=I
      IF(IN.EQ.100)GOTO 30
      GOTO 31
      ENDIF
   30 CONTINUE
C  IF THE FINAL GROUP HAS MORE THAN 20 MEMBERS THEN IT IS ACCEPTED
C  OTHERWISE IT IS ADDED TO THE PREVIOUS GROUP
      IF(SUM.LT.20)IX=IX-1
      IEND(IX)=100 
C  THE TOTAL NUMBER OF GROUPS IS IX.   THE BEGINNINGS AND ENDS OF THE 
C  GROUPS ARE IEND(0) TO IEND(1), IEND(1) TO IEND(2) ...........
C  ...., IEND(IX-1) TO IEND(IX).  NOW THERE ARE CALCULATED THE AVERAGE
C  F**2 IN EACH GROUP AND <SIN[THETA]>.
      DO 40 I=1,IX
      SUMF2=0
      SUMTHET=0
      NSUM=0
      DO 41 J=IEND(I-1),IEND(I)
      SUMF2=SUMF2+F2LIST(J)
      SUMTHET=SUMTHET+THETLIST(J)
      NSUM=NSUM+NUMLIST(J)
   41 CONTINUE
      AVSUMF2(I)=SUMF2/NSUM
      AVTHET(I)=SUMTHET/NSUM
   40 CONTINUE
C  NOW CALCULATE CONVERSION FACTORS FOR CALCULATION OF U OR E
      DO 60 I=1,IX
      STLAM=AVTHET(I)/XLAM
      IF(IUE.EQ.2)THEN
      FACTOR=SQRT(1.0/AVSUMF2(I))
      GOTO 62
      ENDIF
      SUM1=0
      SUM2=0
      DO 61 J=1,ITOT
      IA=ISYMB(ITOT)        
      CALL SCAT(IA,STLAM,SCF)
      SUM1=SUM1+NATOM(J)*SCF
      SUM2=SUM2+NATOM(J)*SCF*SCF
   61 CONTINUE
      FACTOR=SQRT(SUM2/AVSUMF2(I))/SUM1
   62 EPS(I)=FACTOR
   60 CONTINUE 
C  A BEST PARABOLIC FIT IS NOW TO BE FOUND TO GIVE EPS AS A FUNCTION 
C  OF SIN[THETA]/LAMDA [=s].
      CALL PARAFIT(AVTHET,EPS,IX,C)
C  C(1)+C(2)*s+C(3)*s*s GIVES THE CONVERSION FACTOR FROM F TO U OR E
C  THE PROPER VALUE OF E(0,0) OR U(0,0) IS KNOW AND IS NOW CALCULATED
   94 WRITE(6,'('' DO YOU WANT TO OUTPUT THE VALUES OF a, b AND c '')')
      WRITE(6,'('' GIVEN IN EQUATION 2.48?    Y/N     '')')
      READ(5,50)Q
      IF(Q.EQ.'Y'.OR.Q.EQ.'y')THEN
      IABC=1
      GOTO 93
      ENDIF
      IF(Q.EQ.'N'.OR.Q.EQ.'n')THEN
      IABC=0
      GOTO 93
      ENDIF
      GOTO 94
   93 IF(IUE.EQ.1)THEN
      Z=1.0
      GOTO 89
      ENDIF
      SUM1=0
      SUM2=0
      DO 87 I=1,ITOT
      NN=ISYMB(I)
      CALL SCATFAC(NN,0.0,SCF)
      SUM1=SUM1+NATOM(I)*SCF
      SUM2=SUM2+NATOM(I)*SCF**2
   87 CONTINUE
      Z=SUM1/SQRT(SUM2)
C  FIND E(0)
   89 DO 90 I=1,NOREF
      IF(HT(I).EQ.0.AND.KT(I).EQ.0)THEN
      UE(I)=Z
      GOTO 90
      ENDIF
   90 CONTINUE
      DO 91 I=1,NOREF
      IF(HT(I).EQ.0.AND.KT(I).EQ.0)GOTO 91
      STHET=XLAM*SQRT(CO1*HT(I)**2+CO2*KT(I)**2-CO3*HT(I)*KT(I))
      CONV=C(1)+C(2)*STHET+C(3)*STHET**2
      UE(I)=CONV*F(I)
C  CHECK THAT THE CALCULATED E IS NOT BIGGER THAN E(0) WHICH CAN
C  HAPPEN FOR A LARGE STRUCTURE FACTOR TAKING INTO ACCOUNT THE 
C  INACCURACIES IN SCALING.
      IF(UE(I).GT.0.95*Z)UE(I)=0.95*Z
   91 CONTINUE
C  ALL F's HAVE NOW BEEN CONVERTED TO U OR E.   THEY ARE NOW WRITTEN ON 
C  A FILE FOR FUTURE USE.
      WRITE(10,*)IUE,NOREF,(HT(I),KT(I),UE(I),I=1,NOREF)
C  THE VALUES OF a, b AND c ARE WRITTEN IF REQUESTED
      IF(IABC.EQ.0)GOTO 95
      WRITE(9,500)C
  500 FORMAT(5H a = ,E9.3,5H b = ,E9.3,5H c = ,E9.3)    
   95 CU='U'
      IF(IUE.EQ.2)CU='E'
      WRITE(9,200)CU
  200 FORMAT(11H VALUES OF ,A1)
      WRITE(9,300)CU,CU,CU
  300 FORMAT(3(16H       H    K   ,A1,3X))
      WRITE(9,400)(HT(I),KT(I),UE(I),I=1,NOREF)
  400 FORMAT(3(5X,I3,I5,F7.3))      
      STOP
      END              

      SUBROUTINE SCATFAC(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(13,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)/11.00,9.76,8.34,6.89,5.47,4.29,3.40,
     +2.76,2.31,2.00,1.78/
      DATA (SF(7,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(8,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(9,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(10,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(11,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(12,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(13,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
      
      SUBROUTINE PARAFIT(X,Y,NDATA,C)
C  THIS SUBROUTINE FINDS THE COEFFICIENTS C(1), C(2), C(3) TO GIVE
C  THE BEST FIT OF Y(I)=C(1)+C(2)*X(I)+C(3)*X(I)**2 TO A SET OF 
C  POINTS X(I),Y(I), I=1 TO NDATA.  IT USES THE SUBROUTINE GAUSSJ
C  WHICH IS DERIVED FROM "NUMERICAL RECIPES" BY W.H.Press,
C  B.P.Flannery, S.A.Taulkolsky &amp; W.T.Vetterling, PUBLISHED BY
C  Cambridge University Press, 1987.
      DIMENSION X(NDATA),Y(NDATA),C(3),A(3,3),B(3),FUN(3)
      DO 1 J=1,3
      DO 2 K=1,3
      A(J,K)=0.
    2 CONTINUE
      B(J)=0.
    1 CONTINUE
      DO 3 I=1,NDATA
      FUN(1)=1
      FUN(2)=X(I)
      FUN(3)=X(I)*X(I)
      DO 4 J=1,3
      DO 5 K=1,J
      A(J,K)=A(J,K)+FUN(K)*FUN(J)
    5 CONTINUE
      B(J)=B(J)+Y(I)*FUN(J)
    4 CONTINUE
    3 CONTINUE
      DO 6 J=1,3
      DO 7 K=1,J
      A(K,J)=A(J,K)
    7 CONTINUE
    6 CONTINUE
      CALL GAUSSJ(A,3,3,B,1,1)
      DO 8 I=1,3
      C(I)=B(I)
    8 CONTINUE
      RETURN
      END

      SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
      PARAMETER (NMAX=50)
      DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
      DO 11 J=1,N
      IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
      BIG=0.
      DO 13 J=1,N
      IF(IPIV(J).NE.1)THEN
      DO 12 K=1,N
      IF (IPIV(K).EQ.0) THEN
      IF (ABS(A(J,K)).GE.BIG)THEN
      BIG=ABS(A(J,K))
      IROW=J
      ICOL=K
      ENDIF
      ELSE IF (IPIV(K).GT.1) THEN
      PAUSE 'Singular matrix'
      ENDIF
12    CONTINUE
      ENDIF
13    CONTINUE
      IPIV(ICOL)=IPIV(ICOL)+1
      IF (IROW.NE.ICOL) THEN
      DO 14 L=1,N
      DUM=A(IROW,L)
      A(IROW,L)=A(ICOL,L)
      A(ICOL,L)=DUM
14    CONTINUE
      DO 15 L=1,M
      DUM=B(IROW,L)
      B(IROW,L)=B(ICOL,L)
      B(ICOL,L)=DUM
15    CONTINUE
      ENDIF
      INDXR(I)=IROW
      INDXC(I)=ICOL
      IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.'
      PIVINV=1./A(ICOL,ICOL)
      A(ICOL,ICOL)=1.
      DO 16 L=1,N
      A(ICOL,L)=A(ICOL,L)*PIVINV
16    CONTINUE
      DO 17 L=1,M
      B(ICOL,L)=B(ICOL,L)*PIVINV
17    CONTINUE
      DO 21 LL=1,N
      IF(LL.NE.ICOL)THEN
      DUM=A(LL,ICOL)
      A(LL,ICOL)=0.
      DO 18 L=1,N
      A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18    CONTINUE
      DO 19 L=1,M
      B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19    CONTINUE
      ENDIF
21    CONTINUE
22    CONTINUE
      DO 24 L=N,1,-1
      IF(INDXR(L).NE.INDXC(L))THEN
      DO 23 K=1,N
      DUM=A(K,INDXR(L))
      A(K,INDXR(L))=A(K,INDXC(L))
      A(K,INDXC(L))=DUM
23    CONTINUE
      ENDIF
24    CONTINUE
      RETURN
      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(13,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)/11.00,9.76,8.34,6.89,5.47,4.29,3.40,
     +2.76,2.31,2.00,1.78/
      DATA (SF(7,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(8,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(9,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(10,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(11,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(12,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(13,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
      
