C		PROGRAM SPLINE
C  
C  TO SPLINE A SPECTRUM
C
C  J.C. PHILLIPS SUNY/BNL 9/30/85
C
C  BASED ON PGM XSPLIN BY
C  J.C. PHILLIPS, EMBL OUTSTATION, HAMBURG
C		16/2/81
C
	DIMENSION X(512),Y(512),WD(512),RD(512),
     1XN(60),FN(60),GN(60),DN(60),THETA(60),W(548)
	INTEGER*2 CITRET
	BYTE FILNAM(10)
	COMMON /COLOR/ IBACK,IGRAPH
	COMMON /TERTYP/ CITRET
	DATA WD/512*1./
	DATA THETA/60*0.0/
	DATA IBACK /8/
C
C  READ IN INPUT FILENAMES, RANGE OF POINTS TO BE FITTED,
C
C	WRITE(5,107)
C	READ(5,108) CITRET
	IF(CITRET.NE.0) CITRET=1
	IF(CITRET.EQ.0) CITRET=2
	WRITE(5,100)
	READ(5,101)FILNAM
	CALL RFILE(X,FILNAM,NCHAN)
	WRITE(5,102)
	READ(5,101) FILNAM
	CALL RFILE(Y,FILNAM,NCHAN)
C	WRITE(5,103)
C	READ(5,101) FILNAM
3	WRITE(5,105)
	READ(5,108) INTER
	IF(INTER.EQ.0) GO TO 7
	IF(INTER.LT.6) INTER=6
	NOKNOT=0
	NSPL=NCHAN
	DO 5 I=1,512
	WD(I)=1.
5	CONTINUE
	DO 6 I=1,60
	THETA(I)=0.
6	CONTINUE
	DO 2 I=1,NCHAN,INTER
	NOKNOT=NOKNOT+1
	XN(NOKNOT)=X(I)
2	CONTINUE
	IF(X(NCHAN).EQ.XN(NOKNOT)) GO TO 4
	NOKNOT=NOKNOT+1
	XN(NOKNOT)=X(NCHAN)
C
C  VB06A,TG01B ARE HARWELL LIBRARY ROUTINES FOR
C  FITTING A CUBIC SPLINE AND EVALUATING THE SPLINE
C
4	CALL VB06A(NSPL,NOKNOT,X,Y,WD,RD,XN,FN,GN,DN,THETA,0,W)
	IF=-1
	W(1)=TG01B(IF,NOKNOT,XN,FN,GN,X(1))
	IF=1
	DO 1 I=2,NCHAN
	W(I)=TG01B(IF,NOKNOT,XN,FN,GN,X(I))
1	CONTINUE
	FILNAM(8)='S'
	FILNAM(9)='P'
	FILNAM(10)='L'
	CALL WFILE(W,FILNAM,NCHAN)
	IGRAPH=3
	XPER=0.
C	CALL XYPLOT(X,Y,1,NCHAN,CITRET,XMIN,XPER,YMIN,YPER,0)
	IGRAPH=2
C	CALL XYPLOT(X,W,1,NCHAN,1,XMIN,XPER,YMIN,YPER,0)
	GO TO 3
7	DO 8 I=1,NCHAN
	X(I)=ABS(Y(I)-W(I))
8	CONTINUE
	NHI=NCHAN-3
	DO 9 I=1,3
	W(I)=X(I)
	W(NHI+I)=X(NHI+I)
9	CONTINUE
	DO 10 I=4,NHI
	W(I)=0.
	J=I-3
	K=I+3
	DO 11 L=J,K
	W(I)=W(I)+X(L)*X(L)
11	CONTINUE
	W(I)=SQRT(W(I)/6.)
10	CONTINUE
	FILNAM(9)='I'
	FILNAM(10)='G'
	CALL WFILE(W,FILNAM,NCHAN)
C	CALL VTPLOT(0,0,6)
	CALL EXIT
100	FORMAT('$ X-AXIS FILENAME   ? ')
101	FORMAT(10A1)
102	FORMAT('$ Y-AXIS FILENAME   ? ')
103	FORMAT('$  SIGMA FILENAME   ? ')
104	FORMAT('$ OUTPUT FILENAME   ? ')
105	FORMAT('$ NO. OF POINTS BETWEEN KNOTS  ? ')
107	FORMAT('$ TYPE 1 FOR CIT-161 TERMINAL, CR FOR VT100 ')
108	FORMAT(I5)
	END
