program daftest
!=========================================================================================
! Performs DAF for 0, 1, 2, 3rd derivatives and compares results
! Compares to bicubic spline
!
use numeric_kinds_module
use bspline
implicit none

integer, parameter :: ngiv=101
integer, parameter :: nask=401

character(len=10) :: date, time, zone
character(len=14) :: tdiff
integer :: time1(8), time2(8), time3(8), time4(8), time5(8)

logical :: eqsp=.true.

integer :: i, j, k, l
integer :: kx=3  ! spline order
integer :: mj
integer :: nodecount(0:2)
integer :: tnodecount(0:2)
real(dp),allocatable :: f(:), v0(:), v1(:), v2(:), v3(:)
real(dp), allocatable :: dif0(:), dif1(:), dif2(:), dif3(:)
real(dp), allocatable :: rel0(:), rel1(:), rel2(:), rel3(:)
real(dp), allocatable :: t0(:), t1(:), t2(:), t3(:)
real(dp), allocatable :: d0(:), d1(:), d2(:), d3(:)
real(dp), allocatable :: delta1(:,:), delta2(:,:), delta0(:,:), delta3(:,:)
real(dp), allocatable :: xgiv(:), xask(:), xtmp(:)
real(dp) :: xmax, xmin, dxa, dxg, pi
real(dp) :: x, ln10, f08
real(dp), allocatable :: sigma(:) ! list of sigma values

real(dp), allocatable :: spl_coef(:)
real(dp), allocatable :: spl_knot(:)
real(dp), allocatable :: spl_val0(:)
real(dp), allocatable :: spl_val1(:)
real(dp), allocatable :: spl_val2(:)
real(dp), allocatable :: spl_val3(:)

real(dp) :: a, b  ! factors in analytical function


integer :: nsig=6
integer :: nm=100

allocate(sigma(nsig))

allocate(xgiv(ngiv), xask(nask))
allocate(f(ngiv))

allocate(v0(nask))
allocate(v1(nask))
allocate(v2(nask))
allocate(v3(nask))

allocate(t0(nask))
allocate(t1(nask))
allocate(t2(nask))
allocate(t3(nask))

allocate(d0(nask))
allocate(d1(nask))
allocate(d2(nask))
allocate(d3(nask))

allocate(dif0(nask))
allocate(dif1(nask))
allocate(dif2(nask))
allocate(dif3(nask))

allocate(rel0(nask))
allocate(rel1(nask))
allocate(rel2(nask))
allocate(rel3(nask))

allocate(delta0(nask,ngiv))
allocate(delta1(nask,ngiv))
allocate(delta2(nask,ngiv))
allocate(delta3(nask,ngiv))

allocate(spl_coef(ngiv))
allocate(spl_knot(ngiv+3))
allocate(spl_val0(nask))
allocate(spl_val1(nask))
allocate(spl_val2(nask))
allocate(spl_val3(nask))


pi=4.d0*atan(1.d0)

! sigmas for testing
sigma(1)=1.41d-1
sigma(2)=1.77d-1
sigma(3)=2.36d-1
sigma(4)=3.54d-1
sigma(5)=3.26676d-1
sigma(6)=0.d0

xmax=10.d0
xmin=-10.d0

dxg=(xmax-xmin)/(ngiv-1)
dxa=(xmax-xmin)/(nask-1)

! linear spacing, square matrices. Make sure the point values match
! well and that the derivatives make sense. 
do i=1,ngiv
   xgiv(i)=xmin+(i-1)*dxg
enddo
do i=1,nask
   xask(i)=xmin+(i-1)*dxa
enddo

a=7.0d0
b=32.d-2

do i=1,ngiv
  x=xgiv(i)
  f(i)=exp(-x*x/b)*cos(a*x)
enddo

! Set up spline
call dbsnak(ngiv, xgiv, kx, spl_knot)
call dbsint(ngiv, xgiv, f, kx, spl_knot, spl_coef)

! Get spline values
do i=1,nask
    x=xask(i)
    spl_val0(i)=dbsval(x,kx,spl_knot, ngiv, spl_coef)
    k=1
    spl_val1(i)=dbsder(k,x,kx,spl_knot, ngiv, spl_coef)
    k=2
    spl_val2(i)=dbsder(k,x,kx,spl_knot, ngiv, spl_coef)
    k=3
    spl_val3(i)=dbsder(k,x,kx,spl_knot, ngiv, spl_coef)
enddo

! get the right answers
do i=1,nask
   x=xask(i)
   t0(i)=exp(-x*x/b)*cos(a*x)
   t1(i)=(-a*sin(a*x)-(2.d0/b)*x*cos(a*x))*exp(-x*x/b)
   t2(i)=exp(-x*x/b)*(-a*a*cos(a*x)-(2.d0/b)*cos(a*x)+4.d0*(x*x/(b*b))*cos(a*x)+(4.d0*a*x/b)*sin(a*x))
   t3(i)=((12.d0*x*Cos(a*x))/(b**2) + (6*a**2*x*Cos(a*x))/(b)  &
     - (8.d0*x*x*x*Cos(a*x))/(B**3) + (a**3*Sin(a*x)) + (6*a*Sin(a*x))/(B) & 
     - (12*a*x**2*Sin(A*x))/(B**2))/exp(x*x/b)
enddo

ln10=log(10.d0)
mj=nm
i=6
print*,"getting daf0"
call date_and_time(date,time,zone,time1)
call daf_deltas(xask,nask,xgiv,ngiv,sigma(i), mj, eqsp, 0, delta0)
print*,"getting daf1"
call date_and_time(date,time,zone,time2)
call daf_deltas(xask,nask,xgiv,ngiv,sigma(i), mj, eqsp, 1, delta1)
print*,"getting daf2"
call date_and_time(date,time,zone,time3)
call daf_deltas(xask,nask,xgiv,ngiv,sigma(i), mj, eqsp, 2, delta2)
call date_and_time(date,time,zone,time4)
call daf_deltas(xask,nask,xgiv,ngiv,sigma(i), mj, eqsp, 3, delta3)
call date_and_time(date,time,zone,time5)
v0=matmul(delta0,f)
v1=matmul(delta1,f)
v2=matmul(delta2,f)
v3=matmul(delta3,f)

do i=1,nask
   d0(i)=-log(abs(v0(i)-t0(i)))/ln10
   d1(i)=-log(abs(v1(i)-t1(i)))/ln10
   d2(i)=-log(abs(v2(i)-t2(i)))/ln10
   d3(i)=-log(abs(v3(i)-t3(i)))/ln10
enddo

open(unit=10,file="val0.csv")
open(unit=11,file="val1.csv")
open(unit=12,file="val2.csv")
open(unit=13,file="val3.csv")

open(unit=20,file="daf_dif0")
open(unit=21,file="daf_dif1")
open(unit=22,file="daf_dif2")
open(unit=23,file="daf_dif3")

open(unit=30,file="daf_sig0")
open(unit=31,file="daf_sig1")
open(unit=32,file="daf_sig2")
open(unit=33,file="daf_sig3")

open(unit=40,file="splined0")
open(unit=41,file="splined1")
open(unit=42,file="splined2")
open(unit=43,file="splined3")

open(unit=50,file="daf_rel0")
open(unit=51,file="daf_rel1")
open(unit=52,file="daf_rel2")
open(unit=53,file="daf_rel3")

open(unit=60,file='spl_dif0')
open(unit=61,file='spl_dif1')
open(unit=62,file='spl_dif2')
open(unit=63,file='spl_dif3')

open(unit=70,file='spl_rel0')
open(unit=71,file='spl_rel1')
open(unit=72,file='spl_rel2')
open(unit=73,file='spl_rel3')

open(unit=80,file='comp_dif0')
open(unit=81,file='comp_dif1')
open(unit=82,file='comp_dif2')
open(unit=83,file='comp_dif3')

do i=1,nask
   ! DAF absolute errors
   dif0(i)=v0(i)-t0(i)
   dif1(i)=v1(i)-t1(i)
   dif2(i)=v2(i)-t2(i)
   dif3(i)=v3(i)-t3(i)
   ! DAF relative errors
   rel0(i)=abs((v0(i)-t0(i))/v0(i))
   rel1(i)=abs((v1(i)-t1(i))/v1(i))
   rel2(i)=abs((v2(i)-t2(i))/v2(i))
   rel3(i)=abs((v3(i)-t3(i))/v3(i))
enddo

do i=1,nask
   write(10,*) xask(i), v0(i), t0(i)
   write(11,*) xask(i), v1(i), t1(i)
   write(12,*) xask(i), v2(i), t2(i)
   write(13,*) xask(i), v3(i), t3(i)
   write(20,*) xask(i), dif0(i)
   write(21,*) xask(i), dif1(i)
   write(22,*) xask(i), dif2(i)
   write(23,*) xask(i), dif3(i)
   write(30,*) xask(i), d0(i)
   write(31,*) xask(i), d1(i)
   write(32,*) xask(i), d2(i)
   write(33,*) xask(i), d3(i)
   write(50,*) xask(i), rel0(i)
   write(51,*) xask(i), rel1(i)
   write(52,*) xask(i), rel2(i)
   write(53,*) xask(i), rel3(i)
enddo

! now get spline errors
do i=1,nask
    d0(i)=v0(i)-spl_val0(i)
    d1(i)=v1(i)-spl_val1(i)
    d2(i)=v2(i)-spl_val2(i)
    d3(i)=v3(i)-spl_val3(i)
    rel0(i)=abs(v0(i)-spl_val0(i)/v0(i))
    rel1(i)=abs(v1(i)-spl_val1(i)/v1(i))
    rel2(i)=abs(v2(i)-spl_val2(i)/v2(i))
    rel3(i)=abs(v3(i)-spl_val3(i)/v3(i))
enddo
    

do i=1,nask
    write(40,*) xask(i), v0(i), spl_val0(i)
    write(41,*) xask(i), v1(i), spl_val1(i)
    write(42,*) xask(i), v2(i), spl_val2(i)
    write(43,*) xask(i), v3(i), spl_val3(i)
    write(60,*) xask(i), d0(i)
    write(61,*) xask(i), d1(i)
    write(62,*) xask(i), d2(i)
    write(63,*) xask(i), d3(i)
    write(70,*) xask(i), rel0(i)
    write(71,*) xask(i), rel1(i)
    write(72,*) xask(i), rel2(i)
    write(73,*) xask(i), rel3(i)
    write(80,*) xask(i), dif0(i), d0(i)
    write(81,*) xask(i), dif1(i), d1(i)
    write(82,*) xask(i), dif2(i), d2(i)
    write(83,*) xask(i), dif3(i), d3(i)
enddo
    

! count nodes
nodecount=0
do i=2,nask
   if (dif0(i)*dif0(i-1) < 0.d0) nodecount(0)=nodecount(0)+1
   if (dif1(i)*dif1(i-1) < 0.d0) nodecount(1)=nodecount(1)+1
   if (dif2(i)*dif2(i-1) < 0.d0) nodecount(2)=nodecount(2)+1
   if (t0(i)*t0(i-1) < 0.d0) tnodecount(0)=tnodecount(0)+1
   if (t1(i)*t1(i-1) < 0.d0) tnodecount(1)=tnodecount(1)+1
   if (t2(i)*t2(i-1) < 0.d0) tnodecount(2)=tnodecount(2)+1
enddo

print*,"nodes 0:",nodecount(0), tnodecount(0)
print*,"nodes 1:",nodecount(1), tnodecount(1)
print*,"nodes 2:",nodecount(2), tnodecount(2)
   
call elapsed_time(time1,time2,tdiff)
print*,"DAF calc time 1 calls ", tdiff
call elapsed_time(time1,time3,tdiff)
print*,"DAF calc time 2 calls ", tdiff
call elapsed_time(time1,time4,tdiff)
print*,"DAF calc time 3 calls ", tdiff
call elapsed_time(time1,time5,tdiff)
print*,"DAF calc time 4 calls ", tdiff

end program