273 lines
4.8 KiB
Fortran
273 lines
4.8 KiB
Fortran
module simplespline
|
|
|
|
use const_and_precisions, only : wp_
|
|
implicit none
|
|
|
|
contains
|
|
|
|
function spli(cspli,n,k,dx)
|
|
implicit none
|
|
integer, intent(in) :: n, k
|
|
real(wp_), intent(in) :: cspli(n,4), dx
|
|
real(wp_) :: spli
|
|
spli=cspli(k,1)+dx*(cspli(k,2)+dx*(cspli(k,3)+dx*cspli(k,4)))
|
|
end function spli
|
|
|
|
function splid(cspli,n,k,dx)
|
|
implicit none
|
|
integer, intent(in) :: n, k
|
|
real(wp_), intent(in) :: cspli(n,4), dx
|
|
real(wp_) :: splid
|
|
splid=cspli(k,2)+dx*(2.0_wp_*cspli(k,3)+3.0_wp_*dx*cspli(k,4))
|
|
end function splid
|
|
|
|
subroutine difcs(x,y,n,iopt,c,ier)
|
|
implicit none
|
|
integer, intent(in) :: n, iopt
|
|
real(wp_), intent(in) :: x(n), y(n)
|
|
real(wp_), intent(inout) :: c(n*4)
|
|
integer :: ier
|
|
integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3
|
|
real(wp_) :: xb,xc,ya,yb,h,a,r,dya,dyb,dy2
|
|
jmp =1
|
|
if (n <= 1) return
|
|
!
|
|
! initialization
|
|
!
|
|
xc =x(1)
|
|
yb =y(1)
|
|
h =0.0_wp_
|
|
a =0.0_wp_
|
|
r =0.0_wp_
|
|
dyb =0.0_wp_
|
|
!
|
|
! iol=0 - given derivative at first point
|
|
! ioh=0 - given derivative at last point
|
|
!
|
|
iol =iopt-1
|
|
ioh =iopt-2
|
|
if (ioh == 1) then
|
|
iol =0
|
|
ioh =0
|
|
end if
|
|
dy2 =c(2)
|
|
!
|
|
! form the system of linear equations
|
|
! and eliminate subsequentially
|
|
!
|
|
j =1
|
|
do i=1,n
|
|
j2 =n+i
|
|
j3 =j2+n
|
|
a =h*(2.0_wp_-a)
|
|
dya =dyb+h*r
|
|
if (i>=n) then
|
|
!
|
|
! set derivative dy2 at last point
|
|
!
|
|
dyb =dy2
|
|
h =0.0_wp_
|
|
if (ioh/=0) then
|
|
dyb =dya
|
|
goto 13
|
|
end if
|
|
else
|
|
j =j+jmp
|
|
xb =xc
|
|
xc =x(j)
|
|
h =xc-xb
|
|
!
|
|
! ii=0 - increasing abscissae
|
|
! ii=1 - decreasing abscissae
|
|
!
|
|
ii =0
|
|
if (h==0) return
|
|
if (h<0) ii =1
|
|
ya =yb
|
|
yb =y(j)
|
|
dyb =(yb-ya)/h
|
|
if (i<=1) then
|
|
j1 =ii
|
|
if (iol/=0) goto 13
|
|
dya =c(1)
|
|
end if
|
|
end if
|
|
if (j1-ii /= 0) return
|
|
a =1.0_wp_/(h+h+a)
|
|
13 continue
|
|
r =a*(dyb-dya)
|
|
c(j3)=r
|
|
a =h*a
|
|
c(j2)=a
|
|
c(i) =dyb
|
|
end do
|
|
!
|
|
! back substitution of the system of linear equations
|
|
! and computation of the other coefficients
|
|
!
|
|
a =1.0_wp_
|
|
j1 =j3+n+ii-ii*n
|
|
i =n
|
|
do iol=1,n
|
|
xb =x(j)
|
|
h =xc-xb
|
|
xc =xb
|
|
a =a+h
|
|
yb =r
|
|
r =c(j3)-r*c(j2)
|
|
ya =r+r
|
|
c(j3)=ya+r
|
|
c(j2)=c(i)-h*(ya+yb)
|
|
c(j1)=(yb-r)/a
|
|
c(i) =y(j)
|
|
a =0.0_wp_
|
|
j =j-jmp
|
|
i =i-1
|
|
j2 =j2-1
|
|
j3 =j3-1
|
|
j1 =j3+n+ii
|
|
end do
|
|
ier =0
|
|
end subroutine difcs
|
|
|
|
subroutine difcsn(xx,yy,nmx,n,iopt,cc,ier)
|
|
!
|
|
! same as difcs but with dimension(xx,yy) = nmx > n
|
|
!
|
|
implicit none
|
|
integer, intent(in) :: nmx, n, iopt
|
|
real(wp_), intent(in) :: xx(nmx), yy(nmx)
|
|
real(wp_), intent(inout) :: cc(nmx,4)
|
|
integer :: ier
|
|
integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3
|
|
real(wp_) :: x(n),y(n),c(n*4),xb,xc,ya,yb,h,a,r,dya,dyb,dy2
|
|
!
|
|
do i=1,n
|
|
x(i)=xx(i)
|
|
y(i)=yy(i)
|
|
end do
|
|
ii=0
|
|
do j=1,4
|
|
do i=1,n
|
|
ii=ii+1
|
|
c(ii)=cc(i,j)
|
|
end do
|
|
end do
|
|
!
|
|
jmp =1
|
|
if (n>1) then
|
|
!
|
|
! initialization
|
|
!
|
|
xc =x(1)
|
|
yb =y(1)
|
|
h =0.0_wp_
|
|
a =0.0_wp_
|
|
r =0.0_wp_
|
|
dyb =0.0_wp_
|
|
!
|
|
! iol=0 - given derivative at first point
|
|
! ioh=0 - given derivative at last point
|
|
!
|
|
iol =iopt-1
|
|
ioh =iopt-2
|
|
if (ioh==1) then
|
|
iol =0
|
|
ioh =0
|
|
end if
|
|
dy2 =c(2)
|
|
!
|
|
! form the system of linear equations
|
|
! and eliminate subsequentially
|
|
!
|
|
j =1
|
|
do i=1,n
|
|
j2 =n+i
|
|
j3 =j2+n
|
|
a =h*(2.0_wp_-a)
|
|
dya =dyb+h*r
|
|
if (i>=n) then
|
|
!
|
|
! set derivative dy2 at last point
|
|
!
|
|
dyb =dy2
|
|
h =0.0_wp_
|
|
if (ioh/=0) then
|
|
dyb =dya
|
|
goto 13
|
|
end if
|
|
else
|
|
j =j+jmp
|
|
xb =xc
|
|
xc =x(j)
|
|
h =xc-xb
|
|
!
|
|
! ii=0 - increasing abscissae
|
|
! ii=1 - decreasing abscissae
|
|
!
|
|
ii =0
|
|
if (h==0) goto 16
|
|
if (h<0) ii =1
|
|
ya =yb
|
|
yb =y(j)
|
|
dyb =(yb-ya)/h
|
|
if (i<=1) then
|
|
j1 =ii
|
|
if (iol/=0) goto 13
|
|
dya =c(1)
|
|
end if
|
|
end if
|
|
if (j1/=ii) goto 16
|
|
a =1.0_wp_/(h+h+a)
|
|
13 continue
|
|
r =a*(dyb-dya)
|
|
c(j3)=r
|
|
a =h*a
|
|
c(j2)=a
|
|
c(i) =dyb
|
|
end do
|
|
!
|
|
! back substitution of the system of linear equations
|
|
! and computation of the other coefficients
|
|
!
|
|
a =1.0_wp_
|
|
j1 =j3+n+ii-ii*n
|
|
i =n
|
|
do iol=1,n
|
|
xb =x(j)
|
|
h =xc-xb
|
|
xc =xb
|
|
a =a+h
|
|
yb =r
|
|
r =c(j3)-r*c(j2)
|
|
ya =r+r
|
|
c(j3)=ya+r
|
|
c(j2)=c(i)-h*(ya+yb)
|
|
c(j1)=(yb-r)/a
|
|
c(i) =y(j)
|
|
a =0.0_wp_
|
|
j =j-jmp
|
|
i =i-1
|
|
j2 =j2-1
|
|
j3 =j3-1
|
|
j1 =j3+n+ii
|
|
end do
|
|
ier =0
|
|
end if
|
|
!
|
|
16 continue
|
|
ii=0
|
|
do j=1,4
|
|
do i=1,nmx
|
|
if(i<=n) then
|
|
ii=ii+1
|
|
cc(i,j)=c(ii)
|
|
else
|
|
cc(i,j)=0.0_wp_
|
|
end if
|
|
end do
|
|
end do
|
|
!
|
|
end subroutine difcsn
|
|
|
|
end module simplespline |