C		PROGRAM BAKOFF
C  
C  TO SUBTRACT BACKGROUND FROM EXAFS SPECTRA.
C  USING VT100 OR CIT 161 COLOUR TERMINAL
C
C  J.C. PHILLIPS SUNY/BNL 2/18/85.
C
	DIMENSION XOLD(512),YOLD(512),Z(128),F(512),W(128),U(256),
     1YNEW(512),XNEW(512),XP(128),YP(128),H(128),A(6),B(6),C(6),G(6),
     2XD(512),YD(512),L(6),WD(512),RD(512),XN(10),FN(10),GN(10),DN(10),
     3THETA(10),W2(126),XCROSS(10),YCROSS(10)
	INTEGER NPCR(10),CITRET
	BYTE FILNAM(10),SCROSS(10)
	COMMON U
	COMMON /HAIRS/ XCROSS,YCROSS,NCROSS,NPCR,SCROSS
	COMMON/VC01B/LP
	COMMON /NORCOM/ XNEW,YNEW
	COMMON /COLOR/ IBACK,IGRAPH
	COMMON /TERTYP/ CITRET
        DATA ZERMON,DSPACE,STPRAD /4.125E9,1.920,1.146E6/
	DATA IBACK,IGRAPH /8,3/
	DATA W/128*1./
	DATA WD/512*1./
	LP=0
C  
C  CONVERT MOTOR AND I0/I FILES
C  TO ENERGY AND LOGE(I0/I) FILES
C
	WRITE(5,200)
	READ(5,201) CITRET
	IF(CITRET.NE.0) CITRET=1
	IF(CITRET.EQ.0) CITRET=2
	WRITE(5,100)
	READ(5,101)FILNAM
	CALL RFILE(XOLD,FILNAM,NDATA)
	WRITE(5,102)
	READ(5,101)FILNAM
	CALL RFILE(YOLD,FILNAM,NDATA)
	WRITE(5,103)
	READ(5,104)OFFSET
	WRITE(5,141)
	READ(5,106) IXTYPE,IYTYPE,ICROSS
C       CALL ASSIGN(1,'MONOC.CON')
C       DEFINE FILE 1 (1,8,U,IFL4)
C       READ(1'1) ZERMON,DSPACE,STPRAD,GAP
C       CALL CLOSE(1)
	ZERMON=ZERMON+OFFSET
	DO 1 I=1,NDATA
	IF(IXTYPE.EQ.0) XOLD(I)=12398./(DSPACE*SIN((ZERMON-XOLD(I))/STPRAD))
	IF(IYTYPE.EQ.1) YOLD(I)=ALOG(YOLD(I))
1	CONTINUE
31	XPER=0
	IGRAPH=3
	CALL XYPLOT(XOLD,YOLD,1,NDATA,2,XMIN,XPER,YMIN,YPER,ICROSS)
C
C  READ IN POINTS FOR PRE-EDGE FIT,
C       ORDER OF POLYNOMIAL,MODE,WEIGHT.
C
	NEXIT=1
	WRITE(5,109)
	READ(5,106) NEXT
	IF(NEXT.EQ.2) GO TO 30
	IF(NCROSS.EQ.0) GO TO 72
	DO 73 I=1,NCROSS
	IF(SCROSS(I).EQ.'L') NLOW=NPCR(I)
	IF(SCROSS(I).EQ.'H') NHI=NPCR(I)
73	CONTINUE
	WRITE(5,142) NLOW,NHI
72	WRITE(5,105)
	READ(5,106) NLOW2,NHI2,NORD,MODE,WEIGHT
	IF(NLOW2.GT.0) NLOW=NLOW2
	IF(NHI2.GT.0) NHI=NHI2
	IF(NLOW.LE.0) NLOW=1
	IF(NHI.LE.0) GO TO 30
	NPNTS=NHI-NLOW+1
	IF(MODE.EQ.1) NPNTS=NPNTS+5
	IF(MODE.EQ.2) NPNTS=NPNTS+1
	X0=XOLD(1)
C	Y0=YOLD(1)
	DO 11 I=1,NDATA
	XNEW(I)=XOLD(I)-X0
	YNEW(I)=YOLD(I)-YMIN
11	CONTINUE
	IF(NPNTS.LE.128) GO TO 12
	WRITE(5,108)
	GO TO 31
C
C MAXIMUM OF 5TH. ORDER POLYNOMIAL ALLOWED
C  DEFAULT TO 1ST. ORDER.
C
12	IF(NORD.GT.0.AND.NORD.LT.6) GO TO 4
	NORD=1
	WRITE(5,107)
4	DO 2 I=NLOW,NHI
	II=I-NLOW+1
	XP(II)=XNEW(I)
	YP(II)=YNEW(I)
	IF(MODE.EQ.3) YP(II)=ALOG(YNEW(I))
2	CONTINUE
C
C  VC01A,PE07A ARE HARWELL LIBRARY ROUTINES FOR
C  FITTING A POLYNOMIAL AND EVALUATING THE POLYNOMIAL.
C
C  MODE=1, STRAIGHT LINE PRELIMINARY FIT IS
C  EXTRAPOLATED TO END OF DATA.
C
	IF(MODE.NE.1) GO TO 13
	NPNTS=NPNTS-5
	CALL VC01A(XP,YP,W,Z,NPNTS,A,B,C,G,H,L,1)
	DO 14 I=1,5
	J=IFIX(FLOAT((NDATA-NHI)*I)/5.)+NHI
	NPNTS=NPNTS+1
	XP(NPNTS)=XNEW(J)
	YP(NPNTS)=PE07A(1,A,B,C,XNEW(J))
	W(NPNTS)=WEIGHT
14	CONTINUE
	GO TO 15
C
C  MODE=2, FIT FORCED THROUGH 0. AT INFINITY.
C
13	IF(MODE.NE.2) GO TO 15
	XP(NPNTS)=XNEW(NLOW)+100.*(XNEW(NDATA)-XNEW(NLOW))
	YP(NPNTS)=0.
	W(NPNTS)=WEIGHT
C
C  FINAL FIT
C
15	CALL VC01A(XP,YP,W,Z,NPNTS,A,B,C,G,H,L,NORD)
	DO 3 I=1,NDATA
	F(I)=PE07A(NORD,A,B,C,XNEW(I))		!+Y0
	IF(MODE.EQ.3) F(I)=EXP(F(I))
	YNEW(I)=YNEW(I)-F(I)
	F(I)=F(I)+YMIN
3	CONTINUE
C
C  PLOT RESULTS OF PRE-EDGE FIT
C
	IGRAPH=2
	CALL XYPLOT(XOLD,F,1,NDATA,1,XMIN,XPER,YMIN,YPER,0)
	WRITE(5,109)
	READ(5,106) NEXT
	NEXT=NEXT+1
	NEXIT=2
	GO TO (34,33,30) NEXT
33	DO 35 I=1,128
35	W(I)=1.
	GO TO 31
34	XPER=0
	IGRAPH=3
	CALL XYPLOT(XOLD,YNEW,1,NDATA,2,XMIN,XPER,YMIN,YPER,0)
	WRITE(5,109)
	READ(5,106) NEXT
	NEXT=NEXT+1
	GO TO (36,33,30) NEXT
C
C  END OF PRE-EDGE FITTING. 
C
C  NOW PLOT EDGE REGION TO CHOOSE POINTS
C  FOR MAIN BACKGOUND SUBTRACTION.
C
36	DO 37 I=1,NDATA
	YOLD(I)=YNEW(I)
37	CONTINUE
38	WRITE(5,110)
	READ(5,106) NLOW2,NHI2
	XPER=0.
	IGRAPH=3
	CALL XYPLOT(XOLD,YOLD,NLOW2,NHI2,2,XMIN,XPER,YMIN,YPER,ICROSS)
	IF(NCROSS.EQ.0) GO TO 74
	NSMTH=0
	NJUMP=0
	DO 75 I=1,NCROSS
	IF(SCROSS(I).EQ.'B') NBELOW=NPCR(I)
	IF(SCROSS(I).EQ.'Z') NZERO=NPCR(I)
	IF(SCROSS(I).EQ.'L') NLOW=NPCR(I)
	IF(SCROSS(I).EQ.'A') NABOVE=NPCR(I)
	IF(SCROSS(I).EQ.'S') NSMTH=IABS(NPCR(I)-NSMTH)
	IF(SCROSS(I).EQ.'J') NJUMP=IABS(NPCR(I)-NJUMP)
	IF(SCROSS(I).EQ.'E') NEND=NPCR(I)
75	CONTINUE
	WRITE(5,143) NBELOW,NZERO,NLOW,NABOVE,NSMTH,NJUMP,NEND
74	WRITE(5,109)
	READ(5,106) NEXT
	NEXIT=3
	IF(NEXT.NE.0) GO TO 38
C
C
C  READ IN RANGE OF POINTS TO BE FITTED,
C       EDGE JUMP POINTS, EZERO POINT.
C
	WRITE(5,111)
	READ(5,112) NBEL2,NZERO2,NLOW2,NABOV2
	IF(NBEL2.GT.0) NBELOW=NBEL2
	IF(NZERO2.GT.0) NZERO=NZERO2
	IF(NLOW2.GT.0) NLOW=NLOW2
	IF(NABOV2.GT.0) NABOVE=NABOV2
	WRITE(5,113)
	READ(5,112) NSMTH2,NJUMP2,NEND2,NOKNOT
	IF(NSMTH2.GT.0) NSMTH=NSMTH2
	IF(NJUMP2.GT.0) NJUMP=NJUMP2
	IF(NEND2.GT.0) NEND=NEND2
	IF(NLOW.LE.0)  NLOW=1
	IF(NJUMP.EQ.0) NEND=NLOW
	IF(NJUMP.EQ.0) NJUMP=1
	EZERO=XOLD(NZERO)
	IF(NOKNOT.LT.1) NOKNOT=2
	WRITE(5,119)
	READ(5,120) THWT,SMOOTH
	IF(THWT.LT.0.0001) THWT=5.
	SMOOTH=SMOOTH/10.0
	IF(SMOOTH.LT.0.0001) SMOOTH=0.3
	NHI2=NDATA
C
C   SMOOTHING OF ABOVE EDGE REGION
C
	IF(NSMTH.LT.1) GO TO 20
	WRITE(5,114)
	NJ2=NJUMP
	K=NLOW-NJUMP
	NHI2=NLOW-1
	DO 42 I=NLOW,NDATA
	IF(K.GE.NEND) GO TO 55
	DO 56 J=1,100
	ITEST=NLOW+NJUMP*(J-1)-I
	IF(ITEST) 56,57,42
56	CONTINUE
55	NJ2=1
57	NHI2=NHI2+1
	K=K+NJ2
	YNEW(NHI2)=YOLD(I)
	XNEW(NHI2)=XOLD(I)
42	CONTINUE
	NR2=NSMTH/2
	DO 40 I=1,NR2
	J=NLOW-I
	K=NLOW+I
	YNEW(J)=2.*YNEW(NLOW)-YNEW(K)
40	CONTINUE
	NBOT=NLOW+NR2
	NTOP=NHI2-NR2
	NCYCLE=0
43	NSTOP=0
	NCYCLE=NCYCLE+1
	IF(NCYCLE.GT.10) GO TO 41
	DO 44 I=NLOW,NTOP
	F(I)=0.
	DO 45 J=1,NSMTH
	K=I+J-NR2-1
	F(I)=F(I)+YNEW(K)
45	CONTINUE
	F(I)=F(I)/NSMTH
	IF(I.LE.NBOT) GO TO 44
	TEST=ABS(100.*(F(I)-YNEW(I))/YNEW(I))
	IF(TEST.GT.SMOOTH) NSTOP=1
44	CONTINUE
	IF(NSTOP.EQ.0) GO TO 41
C
C   BOTTOM END FIX
C
	DO 8 I=1,NR2
	J=NLOW-I
	K=NLOW+I
	YNEW(J)=2.*F(NLOW)-F(K)
8	CONTINUE
	IF(NCYCLE.GT.1) GO TO 46
C
C  TOP END OF REGION
C
	NTOP2=NTOP+1
	DO 7 I=NTOP2,NHI2
	F(I)=YNEW(I)
7	CONTINUE
46	DO 9 I=NLOW,NHI2
	YNEW(I)=F(I)
9	CONTINUE
	GO TO 43
41	XPER=0.
	IGRAPH=3
	CALL XYPLOT(XOLD,YOLD,NLOW,NDATA,CITRET,XMIN,XPER,YMIN,YPER,0)
	IGRAPH=2
	CALL XYPLOT(XNEW,F,NLOW,NHI2,1,XMIN,XPER,YMIN,YPER,0)
	WRITE(5,115) NCYCLE
	WRITE(5,109)
	READ(5,106) NEXT
20	NPNTS=0
	DO 52 I=NLOW,NHI2
	NPNTS=NPNTS+1
	IF(NSMTH.LT.1) XD(NPNTS)=XOLD(I)
	IF(NSMTH.GT.0) XD(NPNTS)=XNEW(I)
	IF(NSMTH.LT.1) YD(NPNTS)=YOLD(I)
	IF(NSMTH.GT.0) YD(NPNTS)=F(I)
52	CONTINUE
	WRITE(5,116)
	NSPL=NPNTS
	NOKNOT=NOKNOT+1
	NKJUMP=NSPL/NOKNOT
	DO 64 I=2,NOKNOT
	J=(I-1)*NKJUMP
	XN(I)=XD(J)
	THETA(I)=THWT
64	CONTINUE
	XN(1)=XD(1)
	THETA(1)=THWT
	NOKNOT=NOKNOT+1
	XN(NOKNOT)=XD(NSPL)
	THETA(NOKNOT)=THWT
C
C  VB06A,TG01B ARE HARWELL LIBRARY ROUTINES FOR
C  FITTING A CUBIC SPLINE AND EVALUATING THE SPLINE
C
	CALL VB06A(NSPL,NOKNOT,XD,YD,WD,RD,XN,FN,GN,DN,THETA,0,W)
	DO 54 I=1,NLOW
	F(I)=0.
54	CONTINUE
	EJUMP=YOLD(NABOVE)-YOLD(NBELOW)
	NPNTS=0
	NHI2=NDATA-NSMTH/2
	DO 65 I=NLOW,NDATA
	NPNTS=NPNTS+1
	J=I-NLOW+1
	K=J-2
	XNEW(J)=SQRT(XOLD(I)-EZERO)/1.95193
	F(I)=TG01B(K,NOKNOT,XN,FN,GN,XOLD(I))
	YNEW(J)=10000.*(YOLD(I)-F(I))/EJUMP
	XNEW(J)=10.*XNEW(J)
65	CONTINUE
C
C  PLOT RESULTS
C
	XPER=0.
	IGRAPH=3
	CALL XYPLOT(XOLD,YOLD,NLOW,NDATA,CITRET,XMIN,XPER,YMIN,YPER,0)
	IGRAPH=2
	CALL XYPLOT(XOLD,F,NLOW,NDATA,1,XMIN,XPER,YMIN,YPER,0)
	WRITE(5,109)
	READ(5,106) NEXT
	IF(NEXT.NE.0) GO TO 38
	IPASS=0
78	XPER=0.
	IPASS=IPASS+1
	IGRAPH=3
	IF(IPASS.EQ.1) CALL XYPLOT(XNEW,YNEW,1,NPNTS,1,XMIN,XPER,YMIN,YPER,0)
	IF(IPASS.EQ.2) CALL XYPLOT(XNEW,F,1,NPNTS,1,XMIN,XPER,YMIN,YPER,0)
	A(1)=XNEW(1)
	A(2)=XNEW(NPNTS)
	B(1)=0.
	B(2)=0.
	IGRAPH=5
	CALL XYPLOT(A,B,1,2,1,XMIN,XPER,YMIN,YPER,0)	! PLOT ZERO LINE
	WRITE(5,109)
	READ(5,106) NEXT
	IF(NEXT.NE.0) GO TO 38
C
C  OUTPUT RESULTS
C
C	WRITE(5,130)
C	CALL NORMAL(NPNTS)
	IF(IPASS.EQ.2) GO TO 76
	DO 77 I=1,NPNTS
	F(I)=YNEW(I)*XNEW(I)**2/(XNEW(1)**2)
77	CONTINUE
	GO TO 78
76	FILNAM(8)='E'
	FILNAM(9)='K'
	FILNAM(10)='V'
	CALL WFILE(XNEW,FILNAM,NPNTS)
	WRITE(5,117) FILNAM
	FILNAM(8)='E'
	FILNAM(9)='X'
	FILNAM(10)='A'
	CALL WFILE(YNEW,FILNAM,NPNTS)
	WRITE(5,118) FILNAM
	GO TO 80
30	FILNAM(8)='N'
	FILNAM(9)='R'
	FILNAM(10)='G'
	CALL WFILE(XOLD,FILNAM,NDATA)
	WRITE(5,117) FILNAM
	IF(NEXIT.GT.1) GO TO 70
	FILNAM(8)='L'
	FILNAM(9)='S'
	FILNAM(10)='P'
	GO TO 71
70	FILNAM(8)='P'
	FILNAM(9)='B'
	FILNAM(10)='S'
71	IF(NEXIT.EQ.2) CALL SPSCAL(YNEW,1,NDATA)
	IF(NEXIT.NE.2) CALL SPSCAL(YOLD,1,NDATA)
	IF(NEXIT.EQ.2) CALL WFILE(YNEW,FILNAM,NDATA)
	IF(NEXIT.NE.2) CALL WFILE(YOLD,FILNAM,NDATA)
	WRITE(5,118) FILNAM
80	CALL VTPLOT(0,0,6)
	CALL EXIT
100	FORMAT('$ XAXIS FILENAME    ? ')
101	FORMAT(10A1)
102	FORMAT('$ SPECTRUM FILENAME ? ')
103	FORMAT('$ MONOC. OFFSET     ? ')
104	FORMAT(F10.0)
105	FORMAT('$ NLOW,NHI,NORD,MODE,WEIGHT ? ')
106	FORMAT(4I5,F12.5)
107	FORMAT(/,'  ORDER OF POLYNOMIAL DEFAULTED TO 1 ')
108	FORMAT(' TOO MANY POINTS, TRY AGAIN ')
109	FORMAT(/,'$ NEXT [CR=CONTINUE, 1=GO BACK, 2=EXIT] ? ')
110	FORMAT('$ PLOT BETWEEN POINTS ? ')
111	FORMAT('$ NBELOW,NZERO,NLOW,NABOVE    ? ')
112	FORMAT(10I5)
113	FORMAT('$ NSMTH,NJUMP,NEND,NOKNOT     ? ')
114	FORMAT(//,'  SMOOTHING, HAVE PATIENCE !!! ')
115	FORMAT('  SMOOTHING FINISHED AFTER',I5,' CYCLES')
116	FORMAT(//,'  SPLINING, HAVE PATIENCE !!! ')
117	FORMAT(/,' X-AXIS OUTPUT FILE =  ',10A1)
118	FORMAT(/,' Y-AXIS OUTPUT FILE =  ',10A1)
119	FORMAT('$ KNOT WEIGHT, SMOOTHING PAR. [DEF. 5.0, 3.0] ? ')
120	FORMAT(2F12.5)
C130	FORMAT(' NOW PERFORM NORMALIZATION FOR ATOMIC ',
C     1'ABSORPTION FALLOFF',/)
141	FORMAT('$ IXTYPE,IYTYPE,ICROSS  ? ')
142	FORMAT(/' NLOW, NHI =',2I5)
143	FORMAT(/' NBELOW  NZERO   NLOW NABOVE  NSMTH  NJUMP   NEND',/,
     1' ',7(2X,I5))
200	FORMAT('$ TYPE <CR> FOR VT100, OR 1 FOR CIT-161 ')
201	FORMAT(I5)
	END
