!-------------------------------------------------------------------!
subroutine header(nc,ni,a,b)
implicit none
integer, 	intent(out)	:: nc,ni
real(kind=8),	intent(out)	:: a,b
!
open(1234,file='ph15.dat')
read(1234,*) a,b,nc,ni
close(1234)
end subroutine header
!-------------------------------------------------------------------!
subroutine load(nc,ni,c,dat)
implicit none
integer, 				intent(in)	:: nc,ni
real(kind=8),	dimension(nc,ni,6), 	intent(out) 	:: c
real(kind=8),	dimension(ni,2), 	intent(out)	:: dat
integer :: i,j,p
!
open(1234,file='ph15.dat')
read(1234,*) 
do p=1,ni
	read(1234,*) dat(p,1),dat(p,2)
	do i=1,nc
		read(1234,*) j,c(i,p,:)
	enddo
enddo
close(1234)
end subroutine load
!-------------------------------------------------------------------!
subroutine pos_ph(a,b,nc,ni,cf,dat,t,x)
implicit none
!
real(kind=8), 				intent(in)	:: a,b
integer, 				intent(in) 	:: nc,ni
real(kind=8),	dimension(nc,ni,6), 	intent(in)	:: cf
real(kind=8),	dimension(ni,2), 	intent(in)	:: dat
real(kind=8), 				intent(in)	:: t
real(kind=8),	dimension(6), 		intent(out) 	:: x
!
integer 			:: it,p,j
real(kind=8), 	dimension(nc,6)	:: c
real(kind=8) 			:: d,dd,sv,y,y2
!
if ((t-a)*(t-b).gt.0d0) then 
	write(*,*) 'time is out of range'
	stop
endif
it=int((t-a)/(b-a)*ni)+1
if (it==(ni+1)) it=it-1
if ((dat(it,1).gt.t).or.(dat(it,2).lt.t)) then
	write(*,'(a13,f9.1,a12,f9.1,a1,f9.1,a1)') 'time error : ' &
		,t,' is not in [',dat(it,1),':',dat(it,2),']'
	stop
endif
c(:,:)=cf(:,it,:)
do p=1,6
	sv=0d0
	d=0.d0
	dd=0.d0
	y=(2d0*t-dat(it,1)-dat(it,2))/(dat(it,2)-dat(it,1)) 
	y2=2.d0*y
	do j=nc,2,-1 
		sv=d
		d=y2*d-dd+c(j,p)
		dd=sv
	enddo
	x(p)=y*d-dd+0.5d0*c(1,p)
enddo
end subroutine pos_ph
!-------------------------------------------------------------------!