C		PGM PREFIT
C
C  TO COMPARE EXPERIMENTAL AND THEORETICAL EXAFS SPECTRA. 
C  USING CIT-161 COLOR GRAPHICS TERMINAL OR VT100
C
C  THEORY IS TEO AND LEE PARAMETRISED PHASES AND AMPLITUDES.
C  TEO ET AL. JACS 99 3854 (1977) AND LEE ET AL. JACS 99 3856 (1977)
C
C  J.C. PHILLIPS SUNY/BNL VERSION 2/19/85.
C
	DIMENSION W(512),X(512),Y(512),Z(512),XOLD(512),YOLD(512)
	DIMENSION XZ(2),YZ(2)
	BYTE CHAR,FNAM1(10)
	INTEGER CITRET
	COMMON /COLOR/ IBACK,IGRAPH
	COMMON /TERTYP/ CITRET
	DATA IBACK,IGRAPH /8,3/
	DATA YZ /0.,0./
	DATA PI /3.1416/
	WRITE(5,150)
	READ(5,151) CITRET
	IF(CITRET.GT.0) CITRET=1	! CIT-161
	IF(CITRET.EQ.0) CITRET=2	! VT100
	WRITE(5,100)
	READ(5,101) FNAM1
	CALL RFILE(XOLD,FNAM1,NPNTS)
	WRITE(5,102)
	READ(5,101) FNAM1
	CALL RFILE(YOLD,FNAM1,NPNTS)
	WRITE(5,103)
	READ(5,101) FNAM1
	CALL RFILE(Z,FNAM1,NPAR)
	DO 3 I=1,NPAR
	Z(I)=Z(I)/1000.
3	CONTINUE
4	WRITE(5,104) (Z(I), I=1,10)
	WRITE(5,120)
	WRITE(5,121) (Z(I), I=11,NPAR)
	WRITE(5,105)
	READ(5,106) EZERO,POWER,EKVMIN,EKVMAX,EFPATH,NSHELL
	IF(ABS(EZERO).GT.0.000001) Z(1)=EZERO
	IF(POWER.GT.0.00001) Z(2)=POWER
	IF(EKVMIN.GT.0.00001) Z(3)=EKVMIN
	IF(EKVMAX.GT.0.00001) Z(4)=EKVMAX
	IF(EFPATH.GT.0.00001) Z(5)=EFPATH
	IF(NSHELL.GT.0) NPAR=10*(NSHELL+1)
	WRITE(5,107)
	READ(5,106) SHAKE,A0,A1,A2,A3
	IF(SHAKE.GT.0.00001) Z(6)=SHAKE
	IF(ABS(A0).GT.0.000001) Z(7)=A0
	IF(ABS(A1).GT.0.000001) Z(8)=A1
	IF(ABS(A2).GT.0.000001) Z(9)=A2
	IF(ABS(A3).GT.0.000001) Z(10)=A3
	NSHELL=NPAR/10-1
	DO 5 I=1,NSHELL
	K=10*I+1
	WRITE(5,108) I
	READ(5,106) OCCUP,RADIUS,DEBWAL
	IF(OCCUP.GT.0.1) Z(K)=OCCUP
	IF(RADIUS.GT.0.00001) Z(K+1)=RADIUS
	IF(DEBWAL.GT.0.00001) Z(K+2)=DEBWAL
	WRITE(5,109)
	READ(5,106) A,B,C
	WRITE(5,110)
	READ(5,106) B0,B1,B2,B3
	IF(ABS(B0).GT.0.000001) Z(K+6)=B0
	IF(ABS(B1).GT.0.000001) Z(K+7)=B1
	IF(ABS(B2).GT.0.000001) Z(K+8)=B2
	IF(ABS(B3).GT.0.000001) Z(K+9)=B3
	IF(ABS(A).GT.0.000001) Z(K+3)=A
	IF(ABS(B).GT.0.000001) Z(K+4)=B
	IF(ABS(C).GT.0.000001) Z(K+5)=C
5	CONTINUE
	WRITE(5,122)
	NLOW=1
	NHI=NPNTS
	DO 1 I=1,NPNTS
	E=0.038102*XOLD(I)**2-Z(1)
	IF(E.LT.0.00001) X(I)=0.001
	IF(E.GT.0.) X(I)=SQRT(E/3.8102)
	Y(I)=YOLD(I)*(X(I)**Z(2))
	W(I)=0.
1	CONTINUE
	DO 2 J=1,NSHELL
	K=J*10+1
	IF(Z(K).LT.0.1) GO TO 2
	DO 2 I=1,NPNTS
	Q=X(I)
	IF(Q.LT.Z(3)) NLOW=I
	IF(Q.LE.Z(4)) NHI=I
	AMP=Z(K+3)/(1.+Z(K+4)*Z(K+4)*(Q-Z(K+5))**2)
	PHA=Z(K+6)+Z(7)-PI+(Z(K+7)+Z(8))*Q
	PHA=PHA+(Z(K+8)+Z(9))*Q**2+(Z(K+9)+Z(10))/(Q**3)
	S=(Q**(Z(2)-1.))*Z(K)*AMP*SIN(2.*Q*Z(K+1)+PHA)/(Z(K+1)**2)
	S=S*Z(6)*EXP(-2.0*(Z(K+2)*Q*Q+Z(K+1)/Z(5)))
	W(I)=W(I)+10000.*S
2	CONTINUE
	XPER=0.
	IGRAPH=3
	IPEN=1
	CALL XYPLOT(X,W,NLOW,NHI,IPEN,XMIN,XPER,YMIN,YPER,ICROSS)
	IGRAPH=2
	IF(CITRET.EQ.2) IPEN=2
	CALL XYPLOT(X,Y,NLOW,NHI,IPEN,XMIN,XPER,YMIN,YPER,ICROSS)
	IGRAPH=4
	XZ(1)=X(NLOW)
	XZ(2)=X(NHI)
	CALL XYPLOT(XZ,YZ,1,2,1,XMIN,XPER,YMIN,YPER,ICROSS)
	WRITE(5,111)
	READ(5,101) CHAR
	IF(CHAR.EQ.'E') GO TO 7
	CALL VTPLOT(0,799,6)
	GO TO 4
7	WRITE(5,113)
	READ(5,101) CHAR
	IF(CHAR.NE.'Y') GO TO 12
	DO 8 I=1,NPAR
	Z(I)=1000.*Z(I)
8	CONTINUE
	CALL WFILE(Z,FNAM1,NPAR)
12	WRITE(5,114)
	READ(5,101) FNAM1
	IF(FNAM1(1).EQ.' ') GO TO 9
	NPNTS=NHI-NLOW+1
	DO 13 I=1,NPNTS
	W(I)=W(I+NLOW-1)
	X(I)=X(I+NLOW-1)
	Y(I)=Y(I+NLOW-1)
13	CONTINUE
	CALL WFILE(W,FNAM1,NPNTS)
	FNAM1(8)='E'
	FNAM1(9)='E'
	FNAM1(10)='X'
	CALL WFILE(Y,FNAM1,NPNTS)
	FNAM1(9)='K'
	FNAM1(10)='V'
	DO 29 I=1,NPNTS
29	X(I)=10.*X(I)
	CALL WFILE(X,FNAM1,NPNTS)
9	CALL VTPLOT(0,799,6)
	CALL EXIT
100	FORMAT('$ X AXIS FILE ? ')
101	FORMAT(10A1)
102	FORMAT('$ Y AXIS FILE ? ')
103	FORMAT('$ SHELL FILE  ? ')
104	FORMAT(/,' ',6X,'EZERO',7X,'POWER',8X,'KMIN',8X,'KMAX',6X,'LAMBDA',/,
     15F12.5,//,' ',6X,'SHAKE',10X,'A0',10X,'A1',10X,'A2',10X,
     2'A3',/,5F12.5)
120	FORMAT(/,' OCC  RADIUS  DEBWAL     AMPLITUDE PARAMETERS     ',
     1'        PHASE PARAMETERS',/)
121	FORMAT(F4.0,2F8.4,4X,3F7.4,4X,F7.3,F9.4,F9.5,F7.2,/)
122	FORMAT(/,' CALCULATING, BE PATIENT !!! ')
105	FORMAT(/,'$ EZERO, POWER, KMIN, KMAX, LAMBDA, NSHELL ? ')
106	FORMAT(5F12.5,I5)
107	FORMAT('$ SHAKE, CENTRAL ATOM PHASE PARAMETERS [A0-4] ? ')
108	FORMAT(/,' SHELL ',I2,//,'$ OCCUPANCY, RADIUS, DEBYE-WALLER ? ')
109	FORMAT('$ AMPLITUDE PARAMETERS [A,B,C] ? ')
110	FORMAT('$ PHASE PARAMETERS [B0-4] ? ')
111	FORMAT('$ CR TO CONTINUE, E TO EXIT ')
113	FORMAT('$ DO YOU WANT TO SAVE THE NEW PARAMETERS ? ')
114	FORMAT('$ TO SAVE CALCULATED SPECTRUM INPUT A FILENAME - ')
150	FORMAT('$ TYPE <CR> FOR VT100, 1 FOR CIT-161 ')
151	FORMAT(I5)
	END
