subroutine get_cubic_para(xi,yi,ai,bi,ci,di,ni)
! get parameters for interpolation
integer :: ni,i
double precision :: xi(ni),yi(ni)
double precision :: ai(ni),bi(ni),ci(ni),di(ni)
double precision :: matrix1(ni,ni),ri1(ni)
double precision :: matrix2(ni,ni),ri2(ni)
double precision :: mi(ni),hi(ni)
! build equations
matrix1=0.d0
matrix2=0.d0
do i=1,ni-1
hi(i)=xi(i+1)-xi(i)
enddo
matrix1(1,1)=1.0
matrix1(ni,ni)=1.0
ri1(1)=0.d0
ri1(ni)=0.d0
do i=2,ni-1
matrix1(i-1,i)=hi(i-1)
matrix1(i,i)=2*(hi(i-1)+hi(i))
matrix1(i+1,i)=hi(i)
ri1(i)=(yi(i+1)-yi(i))/hi(i)-(yi(i)-yi(i-1))/hi(i-1)
enddo
ri1=ri1*6
! solve equations
do i=1,ni
matrix2(i,i)=1.d0
enddo
matrix2(2,1)=matrix1(2,1)/matrix1(1,1)
ri2(1)=ri1(1)/matrix1(1,1)
do i=2,ni-1
matrix2(i+1,i)=matrix1(i+1,i)/(matrix1(i,i)-matrix1(i-1,i)*matrix2(i,&
i-1))
enddo
do i=2,ni
ri2(i)=ri1(i)-matrix1(i-1,i)*ri2(i-1)
ri2(i)=ri2(i)/(matrix1(i,i)-matrix1(i-1,i)*matrix2(i,i-1))
enddo
mi(ni)=ri2(ni)
do i=ni-1,1,-1
mi(i)=ri2(i)-matrix2(i+1,i)*mi(i+1)
enddo
! get parameters for interpolation
ai=yi
ci=mi(1:ni)
do i=1,ni-1
bi(i)=(yi(i+1)-yi(i))/hi(i)-hi(i)*mi(i)/2.0-hi(i)*(mi(i+1)-mi(i))/6.0
di(i)=(mi(i+1)-mi(i))/(6.0*hi(i))
enddo
end subroutine get_cubic_para