3th version only has C++ codes.
I rewrite the Interpolation class for Fortran, and add a Akima method.
needed:
use nrutil ! From Numerical Recipes 2th version
use stdlib_optval ! stdlib
!插值基类 默认线性插值
module Base_InterpMod
implicit none
private
!integer(8),public::num1=0,num2=0,num11=0,num22=0
type,public::Base_Interp !插值基类 线性插值
integer::n__ !一维数组元素个数
integer::mm__ !参与插值节点数
integer::jsav__ !节点序号中间保存值 为下一次搜索作准备
integer::dj__ !节点序号间隔
logical::cor__ !0False使用locate 1True使用hunt
real,allocatable::xx__(:),yy__(:) !一维数组
contains
procedure::create !构造方法
procedure::interp !插值函数 对外接口
procedure,non_overridable::locate,hunt !两种搜索方法
procedure::rawinterp !具体的插值算法实现 不同插值方法重写该实例方法
final::destroy !析构方法
end type
contains
subroutine destroy(this) !析构方法
type(Base_Interp)::this
if(allocated(this%xx__))deallocate(this%xx__)
if(allocated(this%yy__))deallocate(this%yy__)
write(*,*)"Base_Interp类已完成析构!"
end subroutine
subroutine create(this,x,y,m,yp1,ypn) !接口统一的构造方法
use nrutil,only:assert_eq
use stdlib_optval,only:optval !20251027增加
class(Base_Interp)::this
real,intent(in)::x(:),y(:) !一维数组
integer,optional,intent(in)::m !一般情况下为参与插值节点数 重心有理插值时为分子多项式P(x)的最高次数
real,optional,intent(in)::yp1,ypn !只在样条插值时使用 首尾端点导数值
this%n__=assert_eq(size(x),size(y),'create')
this%mm__=optval(m,2) !不传入时默认为2
select type(this)
type is(Base_Interp)
this%mm__=2 !若this是父类线性插值 强制参与插值节点数为2
end select
this%jsav__=0 !C++为0
this%cor__=.false. !第一次用locate搜索 False使用locate True使用hunt
!这里自动释放旧 分配新 无须手动deallocate与allocate 前提：x与this%xx__同维度
this%xx__=x !复制数据 x是实参的别名
this%yy__=y
this%dj__=min(1,int(real(this%n__)**0.25)) !C++:MIN(1,(int)pow((Doub)n,0.25))
end subroutine
function interp(this,x) !插值函数
class(Base_Interp)::this
real,intent(in)::x !插值点
real::interp !插值计算值
integer::j !插值位置
!如果之前成功搜索且距离不远，使用hunt，否则使用locate
if(this%cor__)then
j=this%hunt(x) !True使用hunt
else
j=this%locate(x) !False使用locate
end if
!此时j为参与插值第1个节点的序号
interp=this%rawinterp(j,x) !调用具体的插值算法 j用于数组切片
!print *,num1,num2,num11,num22
end function
function locate(this,x) !二分法检索
class(Base_Interp)::this
real,intent(in)::x !插值点
integer::locate
!局部变量
integer::jl,jm,ju !下、中、上节点序号
logical::ascnd
!检查是否在范围内
if((this%n__<2) .or. (this%mm__<2) .or. (this%mm__>this%n__))then
!print *,this%mm__,this%n__
error stop "locate size error"
end if
ascnd = (this%xx__(this%n__) >= this%xx__(1)) !判断数组是否升序
jl=0 !下界
ju=this%n__+1 !上界
do
if (ju-jl <= 1) exit !n<=0或差距为1
jm=shiftr(ju+jl,1) !位右移1位 第3版修改 第2版为jm=(ju+jl)/2 !下标中点
!jm取值范围为[0,n] 0=(0+1)/2;n=(n+n+1)/2
if (ascnd .EQV. (x >= this%xx__(jm))) then !解位于中点或中点偏上
jl=jm !中点作为新的下界
else !解位于中点偏下 x<xx(jm)
ju=jm !中点作为新的上界
end if
!num1=num1+1
end do
!num11=num11+1
this%cor__=abs(jl-this%jsav__)<=this%dj__ !下一次计算采用的搜索方法 False使用locate True使用hunt
this%jsav__=jl !保存该值 下一次搜索作准备
!C++第三版P115:MAX(0,MIN(n-mm,jl-((mm-2)>>1))) C++第二版P89:MIN(MAX(jl-((mm-1)>>1),0),n-mm))
!确保索引在有效范围内
locate=max(1,min(this%n__+1-this%mm__,jl-(this%mm__-1)/2)) !C++第三版P115
end function
function hunt(this,x) !猎寻搜索
class(Base_Interp)::this
real,intent(in)::x !插值点
integer::hunt
integer::jl,jm,ju,inc !下、中、上节点序号
logical::ascnd
jl=this%jsav__ !利用上次计算的节点
inc=1
!print *,this%n__,this%mm__
!检查是否在范围内
if((this%n__<2) .or. (this%mm__<2) .or. (this%mm__>this%n__))error stop "hunt size error"
ascnd = (this%xx__(this%n__) >= this%xx__(1)) !判断数组是否升序
if (jl <= 0 .or. jl > this%n__) then !预测值不准 无参考价值 直接转到二分法处理
jl=0
ju=this%n__+1
else
inc=1 !搜索步长
!决定搜索方向
if (x >= this%xx__(jl) .eqv. ascnd) then !向右搜索
do
ju=jl+inc
!num2=num2+1
if (ju > this%n__) then !搜索结束 已至表尾
ju=this%n__+1
exit
else
if (x < this%xx__(ju) .eqv. ascnd) exit
jl=ju
inc=inc+inc !步长加倍
end if
end do
else !向左搜索
ju=jl
do
jl=ju-inc
!num2=num2+1
if (jl < 1) then !搜索结束 已至表尾
jl=0
exit
else
if (x >= this%xx__(jl) .eqv. ascnd) exit
ju=jl
inc=inc+inc !步长加倍
end if
end do
end if
end if
!在找到的范围内进行二分查找精确定位
do
if (ju-jl <= 1) then
exit
else
jm=shiftr(ju+jl,1)
if (x >= this%xx__(jm) .eqv. ascnd) then
jl=jm
else
ju=jm
end if
end if
!num2=num2+1
end do
!num22=num22+1
this%cor__=abs(jl-this%jsav__)<=this%dj__ !下一次计算采用的搜索方法 False使用locate True使用hunt
this%jsav__=jl !保存该值 下一次搜索作准备
!C++第三版P115:MAX(0,MIN(n-mm,jl-((mm-2)>>1))) C++第二版P89:MIN(MAX(jl-((mm-1)>>1),0),n-mm))
!确保索引在有效范围内
hunt=max(1,min(this%n__+1-this%mm__,jl-(this%mm__-1)/2)) !C++第三版P115
end function
function rawinterp(this,j,x) !不同插值方法重写该实例方法
class(Base_Interp)::this
integer,intent(in)::j !参与插值第1节点 j用于数组切片
real,intent(in)::x !插值点
real::rawinterp
!线性插值
if(this%xx__(j)==this%xx__(j+1))then !避免除0
rawinterp=this%yy__(j)
else
rawinterp=this%yy__(j)+((x-this%xx__(j))/(this%xx__(j+1)-this%xx__(j)))*(this%yy__(j+1)-this%yy__(j))
end if
end function
end module Base_InterpMod
!Akima插值
module Akima_InterpMod
use Base_InterpMod
implicit none
private
type,public,extends(Base_Interp)::Akima_Interp
integer::k__ !子区间序号
real::s__(4)=0.0 !三次多项式系数
contains
!无须重写构造方法
procedure::rawinterp !插值实现重写
end type
contains
function rawinterp(this,j,x) !源码摘自徐士良常用算法集
use nrutil,only:iminloc
class(Akima_Interp)::this
integer,intent(in)::j !参与插值第1节点 j用于数组切片
real,intent(in)::x !插值点
real::rawinterp
real::u(5),p,q !u为各点斜率
!这里this%k__为子区间序号 与基类中的this%mm__含义不同
rawinterp=0.0; this%s__=0.0
if (this%n__<1)then
this%k__ = 0
return
end if
if (this%n__==1) then
this%k__ = 0
this%s__(1)=this%yy__(1)
rawinterp=this%yy__(1)
return
end if
if (this%n__==2)then
this%k__ = 1
this%s__(1)=this%yy__(1)
this%s__(2)=(this%yy__(2)-this%yy__(1))/(this%xx__(2)-this%xx__(1))
rawinterp=(this%yy__(1)*(x-this%xx__(2))-this%yy__(2)*(x-this%xx__(1)))/(this%xx__(1)-this%xx__(2))
return
end if
if (x<=this%xx__(2)) then
this%k__=1
else if (x>=this%xx__(this%n__)) then
this%k__=this%n__-1
else
!方法1 interp方法不需重写
this%k__=this%jsav__ !不采用形参值j 因为j受原始this%mm__影响 这里默认全部参与插值
end if
if(this%k__>=this%n__)this%k__=this%n__-1 !增加
u(3)=(this%yy__(this%k__+1)-this%yy__(this%k__))/(this%xx__(this%k__+1)-this%xx__(this%k__))
if (this%n__==3)then
if (this%k__==1)then
u(4)=(this%yy__(3)-this%yy__(2))/(this%xx__(3)-this%xx__(2))
u(5)=2.0*u(4)-u(3)
u(2)=2.0*u(3)-u(4)
u(1)=2.0*u(2)-u(3)
else if(this%k__==2)then
u(2)=(this%yy__(2)-this%yy__(1))/(this%xx__(2)-this%xx__(1))
u(1)=2.0*u(2)-u(3)
u(4)=2.0*u(3)-u(2)
u(5)=2.0*u(4)-u(3)
end if
else
if (this%k__<=2)then
u(4)=(this%yy__(this%k__+2)-this%yy__(this%k__+1))/(this%xx__(this%k__+2)-this%xx__(this%k__+1))
if (this%k__==2)then
u(2)=(this%yy__(2)-this%yy__(1))/(this%xx__(2)-this%xx__(1))
u(1)=2.0*u(2)-u(3)
if (this%n__==4) then
u(5)=2.0*u(4)-u(3)
else
u(5)=(this%yy__(5)-this%yy__(4))/(this%xx__(5)-this%xx__(4))
endif
else
u(2)=2.0*u(3)-u(4)
u(1)=2.0*u(2)-u(3)
u(5)=(this%yy__(4)-this%yy__(3))/(this%xx__(4)-this%xx__(3))
end if
else if (this%k__>=(this%n__-2)) then
u(2)=(this%yy__(this%k__)-this%yy__(this%k__-1))/(this%xx__(this%k__)-this%xx__(this%k__-1))
if (this%k__==(this%n__-2))then
u(4)=(this%yy__(this%n__)-this%yy__(this%n__-1))/(this%xx__(this%n__)-this%xx__(this%n__-1))
u(5)=2.0*u(4)-u(3)
if (this%n__==4) then
u(1)=2.0*u(2)-u(3)
else
u(1)=(this%yy__(this%k__-1)-this%yy__(this%k__-2))/(this%xx__(this%k__-1)-this%xx__(this%k__-2))
end if
else
u(4)=2.0*u(3)-u(2);
u(5)=2.0*u(4)-u(3);
u(1)=(this%yy__(this%k__-1)-this%yy__(this%k__-2))/(this%xx__(this%k__-1)-this%xx__(this%k__-2));
end if
else !5个点
u(2)=(this%yy__(this%k__)-this%yy__(this%k__-1))/(this%xx__(this%k__)-this%xx__(this%k__-1))
u(1)=(this%yy__(this%k__-1)-this%yy__(this%k__-2))/(this%xx__(this%k__-1)-this%xx__(this%k__-2))
u(4)=(this%yy__(this%k__+2)-this%yy__(this%k__+1))/(this%xx__(this%k__+2)-this%xx__(this%k__+1))
u(5)=(this%yy__(this%k__+3)-this%yy__(this%k__+2))/(this%xx__(this%k__+3)-this%xx__(this%k__+2))
end if
end if
this%s__(1)=abs(u(4)-u(3))
this%s__(2)=abs(u(1)-u(2))
if ((this%s__(1)+1.0==1.0).and.(this%s__(2)+1.0==1.0))then
p=(u(2)+u(3))/2.0
else
p=(this%s__(1)*u(2)+this%s__(2)*u(3))/(this%s__(1)+this%s__(2))
end if
this%s__(1)=abs(u(4)-u(5))
this%s__(2)=abs(u(3)-u(2))
if ((this%s__(1)+1.0==1.0).and.(this%s__(2)+1.0==1.0))then
q=(u(3)+u(4))/2.0
else
q=(this%s__(1)*u(3)+this%s__(2)*u(4))/(this%s__(1)+this%s__(2))
end if
this%s__(1)=this%yy__(this%k__)
this%s__(2)=p
this%s__(4)=this%xx__(this%k__+1)-this%xx__(this%k__)
this%s__(3)=(3.0*u(3)-2.0*p-q)/this%s__(4)
this%s__(4)=(q+p-2.0*u(3))/(this%s__(4)*this%s__(4))
p=x-this%xx__(this%k__)
rawinterp=this%s__(1)+p*(this%s__(2)+p*(this%s__(3)+this%s__(4)*p)) !改动
end function
end module Akima_InterpMod
!多项式插值
module Polint_InterpMod
use Base_InterpMod
implicit none
private
type,public,extends(Base_Interp)::Polint_Interp
real::dy__=0.0 !误差估计
contains
!无须重写构造方法
procedure::rawinterp !插值实现重写
end type
contains
function rawinterp(this,j,x)
use nrutil,only:iminloc
class(Polint_Interp)::this
integer,intent(in)::j !参与插值第1节点 j用于数组切片
real,intent(in)::x !插值点
real::rawinterp
integer::m,n,ns
real,allocatable,dimension(:)::xa,ya,c,d,den,ho
allocate(xa(this%mm__),ya(this%mm__),c(this%mm__),d(this%mm__),den(this%mm__),ho(this%mm__))
xa(1:this%mm__)=this%xx__(j:j+this%mm__-1) !只有mm__个元素参与插值
ya(1:this%mm__)=this%yy__(j:j+this%mm__-1)
n=this%mm__
c=ya !初始化表c与表d
d=ya
ho=xa-x
ns=iminloc(abs(x-xa)) !表中最近入口的指针
rawinterp=ya(ns) !rawinterp的初始逼近
ns=ns-1
do m=1,n-1 !对表中的每一列
den(1:n-m)=ho(1:n-m)-ho(1+m:n)
if (any(den(1:n-m) == 0.0)) error stop 'Polint: calculation failure'
!该错误仅当两个输入的xx__相同时才发生
den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
d(1:n-m)=ho(1+m:n)*den(1:n-m) !c和d在此处被更新
c(1:n-m)=ho(1:n-m)*den(1:n-m)
if (2*ns < n-m) then
this%dy__=c(ns+1)
else
this%dy__=d(ns)
ns=ns-1
end if
!表中的一列处理完后 判断要把修正值c还是d加到rawinterp的累积值上
rawinterp=rawinterp+this%dy__
end do
deallocate(xa,ya,c,d,den,ho)
end function
end module Polint_InterpMod
!有理函数插值
module Ratint_InterpMod
use Base_InterpMod
implicit none
private
type,public,extends(Base_Interp)::Ratint_Interp
real::dy__=0.0 !误差估计
contains
!无须重写构造方法
procedure::rawinterp !插值实现重写
end type
contains
function rawinterp(this,j,x)
use nrutil,only:iminloc
class(Ratint_Interp)::this
integer,intent(in)::j !参与插值第1节点 j用于数组切片
real,intent(in)::x !插值点
real::rawinterp
integer::m,n,ns
real,allocatable,dimension(:)::xa,ya,c,d,dd,h,t
real,parameter::TINY=1.0e-25 !一个小的数
allocate(xa(this%mm__),ya(this%mm__),c(this%mm__),d(this%mm__),dd(this%mm__),h(this%mm__),t(this%mm__))
xa(1:this%mm__)=this%xx__(j:j+this%mm__-1) !只有mm__个元素参与插值
ya(1:this%mm__)=this%yy__(j:j+this%mm__-1)
n=this%mm__
h=xa-x
ns=iminloc(abs(h))
rawinterp=ya(ns)
if (x == xa(ns)) then
this%dy__=0.0
return
end if
c=ya
d=ya+TINY !防止少见的0除0情况
ns=ns-1
do m=1,n-1
t(1:n-m)=(xa(1:n-m)-x)*d(1:n-m)/h(1+m:n)
dd(1:n-m)=t(1:n-m)-c(2:n-m+1)
if (any(dd(1:n-m) == 0.0)) error stop 'failure in ratint'
!这种错误情况表明该内插函数在欲求的x处存在极点
dd(1:n-m)=(c(2:n-m+1)-d(1:n-m))/dd(1:n-m)
d(1:n-m)=c(2:n-m+1)*dd(1:n-m)
c(1:n-m)=t(1:n-m)*dd(1:n-m)
if (2*ns < n-m) then
this%dy__=c(ns+1)
else
this%dy__=d(ns)
ns=ns-1
end if
rawinterp=rawinterp+this%dy__
end do
deallocate(xa,ya,c,d,dd,h,t)
end function
end module Ratint_InterpMod
!三次样条插值
module Spline_InterpMod
use Base_InterpMod
implicit none
private
type,public,extends(Base_Interp)::Spline_Interp
real,allocatable::y2__(:) !二阶导数数组
contains
procedure::create !重写构造方法
procedure::rawinterp !插值实现重写
procedure::sety2 !计算二阶导数
final::destroy !析构方法
end type
contains
subroutine destroy(this) !析构方法
type(Spline_Interp)::this
if(allocated(this%y2__))deallocate(this%y2__)
write(*,*)"Spline_Interp类已完成析构!"
end subroutine
subroutine create(this,x,y,m,yp1,ypn) !重写构造方法
use stdlib_optval,only:optval !20251027增加
class(Spline_Interp)::this
real,intent(in)::x(:),y(:) !一维数组
integer,optional,intent(in)::m !一般情况下为参与插值节点数 重心有理插值时为分子多项式P(x)的最高次数
real,optional,intent(in)::yp1,ypn !样条插值时首尾两端点导数值
real,parameter::C=1.0e99
call this%Base_Interp%create(x,y,m) !调用父类方法 m实际没有使用
if(allocated(this%y2__))deallocate(this%y2__) !释放上一次计算的内存 20251115增加
allocate(this%y2__(size(x)))
call this%sety2(x,y,optval(yp1,C),optval(ypn,C)) !计算二阶导数
end subroutine
!计算二阶导数
subroutine sety2(this,x,y,yp1,ypn)
use nr,only:tridag
class(Spline_Interp)::this
real,intent(in)::x(:),y(:) !一维数组
real,intent(in)::yp1,ypn !首尾两端点导数值
integer::n
real,dimension(size(x))::a,b,c,r
n=size(this%y2__)
c(1:n-1)=x(2:n)-x(1:n-1) !sig
r(1:n-1)=6.0*((y(2:n)-y(1:n-1))/c(1:n-1)) !u[i]
r(2:n-1)=r(2:n-1)-r(1:n-2)
a(2:n-1)=c(1:n-2)
b(2:n-1)=2.0*(c(2:n-1)+a(2:n-1))
b(1)=1.0
b(n)=1.0
if (yp1 > 0.99e99) then !左边界条件设成自然的
r(1)=0.0 !u[0]
c(1)=0.0 !y2[0]
else !否则有特定的一阶导数
r(1)=(3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) !u[0]
c(1)=0.5 !y2[0]
end if
if (ypn > 0.99e99) then !右边界条件设成自然的
r(n)=0.0 !un
a(n)=0.0 !qn
else !否则有特定的一阶导数
r(n)=(-3.0/(x(n)-x(n-1)))*((y(n)-y(n-1))/(x(n)-x(n-1))-ypn) !un
a(n)=0.5 !qn
end if
!三对角算法 分解循环+回代循环
call tridag(a(2:n),b(1:n),c(1:n-1),r(1:n),this%y2__(1:n))
end subroutine
function rawinterp(this,j,x)
class(Spline_Interp)::this
integer,intent(in)::j !参与插值第1节点 j用于数组切片 下界节点
real,intent(in)::x
real::rawinterp
integer::khi
real::a,b,h
khi=j+1 !上界节点
h=this%xx__(khi)-this%xx__(j)
if (h == 0.0) error stop 'Bad xa input in splint'
a=(this%xx__(khi)-x)/h
b=(x-this%xx__(j))/h
rawinterp=a*this%yy__(j)+b*this%yy__(khi)+((a*a*a-a)*this%y2__(j)+(b*b*b-b)*this%y2__(khi))*(h*h)/6.0
end function
end module Spline_InterpMod
!重心有理插值
module BaryRat_InterpMod
use Base_InterpMod
implicit none
private
type,public,extends(Base_Interp)::BaryRat_Interp
real,allocatable::w__(:) !权重
integer::d__=1 !d分子多项式P(x)的最高次数 0<=d<=n-1 n-1-d为分母多项式Q(x)的最高次数 推荐d=2-4
contains
procedure::create !重写构造方法
procedure::rawinterp !插值实现重写
procedure::interp
final::destroy
end type
contains
subroutine destroy(this) !析构方法
type(BaryRat_Interp)::this
if(allocated(this%w__))deallocate(this%w__)
write(*,*)"BaryRat_Interp类已完成析构!"
end subroutine
subroutine create(this,x,y,m,yp1,ypn) !重写构造方法
class(BaryRat_Interp)::this
real,intent(in)::x(:),y(:) !一维数组
integer,optional,intent(in)::m !一般情况下为参与插值节点数 重心有理插值时为分子多项式P(x)的最高次数
real,optional,intent(in)::yp1,ypn !样条插值时首尾两端点导数值 没有使用
integer::k,imin,imax,i,j,jmax
real::temp,sum,term
call this%Base_Interp%create(x,y,size(x)) !调用父类方法
if(allocated(this%w__))deallocate(this%w__) !释放上一次计算的内存 20251115增加
allocate(this%w__(this%n__))
if(present(m))this%d__=m !这里m为分子多项式P(x)的最高次数
if(this%n__<=this%d__) error stop "d too large for number of points in BaryRat_interp"
do k=0,this%n__-1 !为了与C++一致
imin=max(k-this%d__,0)
if(k>=this%n__-this%d__)then
imax=this%n__-this%d__-1
else
imax=k
end if
if(mod(imin,2)==1)then !奇数
temp=-1.0
else
temp=1.0 !偶数
end if
sum=0.0
do i=imin,imax
jmax=min(i+this%d__,this%n__-1)
term=1.0
do j=i,jmax
if(j==k)cycle
term=term*(this%xx__(k+1)-this%xx__(j+1)) !Fortran从1开始
end do
term=temp/term
temp=-temp
sum=sum+term
end do
this%w__(k+1)=sum !Fortran从1开始
end do
end subroutine
function rawinterp(this,j,x)
class(BaryRat_Interp)::this
integer,intent(in)::j !没有使用 参与插值第1节点 j用于数组切片
real,intent(in)::x
real::rawinterp
real::num,den,h,temp
integer::i
num=0.0;den=0.0
do i=1,this%n__
h=x-this%xx__(i)
if(h==0.0)then
rawinterp=this%yy__(i)
return
else
temp=this%w__(i)/h
num=num+temp*this%yy__(i)
den=den+temp
end if
end do
rawinterp=num/den
end function
function interp(this,x) !插值函数
class(BaryRat_Interp)::this
real,intent(in)::x !插值点
real::interp !插值计算值
!不调用locate与hunt函数
interp=this%rawinterp(1,x)
end function
end module BaryRat_InterpMod
module Interpolation
use Base_InterpMod !线性插值
use Polint_InterpMod !多项式插值
use Ratint_InterpMod !有理插值
use Spline_InterpMod !三次样条插值
use BaryRat_InterpMod !重心有理插值
use Akima_InterpMod !Akima插值
implicit none
end module
program main
use Interpolation !插值
implicit none
type(BaryRat_Interp)::a !重心有理插值
real,parameter::x(*)=[1.0,2.0,3.0,4.0,5.0,6.0],y(*)=[0.0,3.0,8.0,15.0,24.0,35.0]
call a%create(x,y,3) !y=x*x-1 y'=2x
write(*,*)a%interp(0.9),a%cor__,a%jsav__ !T 0
!d=0 -0.51284645785584748
!d=1 -0.30388530279487141
!d=2 -0.18999999999999997
!d=3 -0.18999999999999934
!d=4 -0.19000000000000086
!d=5 -0.18999999999999967
write(*,*)a%interp(1.5),a%cor__,a%jsav__ !T 1
!d=0 2.0674740484429064
!d=1 1.4454094292803972
!d=2 1.2500000000000004
!d=3 1.2499999999999996
!d=4 1.2500000000000018
!d=5 1.2500000000000004
write(*,*)a%interp(3.1),a%cor__,a%jsav__ !F 3
write(*,*)a%interp(3.5),a%cor__,a%jsav__ !T 3
write(*,*)a%interp(6.5),a%cor__,a%jsav__ !F 5
!d=0 39.233902249806050
!d=1 40.386340977068791
!d=2 41.250000000000000
!d=3 41.250000000000028
!d=4 41.250000000000000
!d=5 41.250000000000000
end program