c     ES361
c     fortran program to determine the m values for relaxed curve
c
      implicit real *8 (A-H,O-Z)
	dimension x(1001),f(1001),a(1001),b(1001),c(1001),r(1001)
      OPEN(5,FILE='sp.dat')
	open(6,FILE='output')
	read(5,*)N
	write(6,*)'N is the number of data points, N=',N
c	stop
	Call main1(N,x,f,a,b,c,r)
      stop
	end
	subroutine main1(N,x,f,a,b,c,r)
      implicit real *8 (a-h,o-z)
	dimension x(N),f(N),a(N),b(N),c(N),r(N)
	read(5,*) (x(i),f(i), i=1,N)
	write(6,101) (i,x(i),f(i), i=1,N)
101   format(//,4x,'i',15x,'x-values',25x,'f-values',/,
     *3x,'---',7x,'--------------------',10x,'--------------------',3x,/
     *,(I5,9x,(1pe12.6),13x,(1pe12.6)))
c	stop
      do 100 i=1,N-2
	a(i)=x(i+1)-x(i)
	b(i)=2*(x(i+1)-x(i)+x(i+2)-x(i+1))
	c(i)=x(i+2)-x(i+1)
      r(i)=6*((f(i+2)-f(i+1))/(x(i+2)-x(i+1))-
     *(f(i+1)-f(i))/(x(i+1)-x(i)))
100   continue
c      boundary condition for relax curve
c      write(*,*) (r(i),a(i),b(i),c(i), i=1,N)
      a(1)=0.
	c(N-2)=0.
	call tridia(N-2,a,b,c,r)
	write(6,102)
	write(6,*) ' m values'
	write(6,*) '----------'
      write(6,103) (r(i), i=1,N-2)
102   format(//)
103   format(f12.6)
      return
	end
      subroutine tridia(N,a,b,c,f)
      IMPLICIT REAL * 8 (A-H,O-Z)
      DIMENSION f(N),a(N),b(N),c(N)
      Do 10 i=2,N
      a(i)=a(i)/b(i-1)
      b(i)=b(i)-a(i)*c(i-1)
10    f(i)=f(i)-a(i)*f(i-1)
      f(n)=f(n)/b(n)
      do 20 i=1,n-1
      ii=n-i
20    f(ii)=(f(ii)-c(ii)*f(ii+1))/b(ii)
      return
      end
