explore_spline_akima.f90

Go to the documentation of this file.

program hello
  implicit none
  INTEGER, PARAMETER :: rprec = selected_real_kind(12,100)
  integer, parameter :: dp = rprec
 
  integer, parameter :: l = 91
 
  real(rprec), dimension(1:l) :: am_aux_s = (/ 0.0000000e+00_dp, &
     1.8181818e-02_dp,   3.6363636e-02_dp,   5.4545455e-02_dp,   7.2727273e-02_dp,   9.0909091e-02_dp, &
     1.0909091e-01_dp,   1.2727273e-01_dp,   1.4545455e-01_dp,   1.6363636e-01_dp,   1.8181818e-01_dp, &
     2.0000000e-01_dp,   2.1818182e-01_dp,   2.3636364e-01_dp,   2.5454545e-01_dp,   2.7272727e-01_dp, &
     2.9090909e-01_dp,   3.0909091e-01_dp,   3.2727273e-01_dp,   3.4545455e-01_dp,   3.6363636e-01_dp, &
     3.8181818e-01_dp,   4.0000000e-01_dp,   4.1818182e-01_dp,   4.3636364e-01_dp,   4.5454545e-01_dp, &
     4.7272727e-01_dp,   4.9090909e-01_dp,   5.0909091e-01_dp,   5.2727273e-01_dp,   5.4545455e-01_dp, &
     5.6363636e-01_dp,   5.8181818e-01_dp,   6.0000000e-01_dp,   6.1818182e-01_dp,   6.3636364e-01_dp, &
     6.5454545e-01_dp,   6.7272727e-01_dp,   6.9090909e-01_dp,   7.0909091e-01_dp,   7.2727273e-01_dp, &
     7.4545455e-01_dp,   7.6363636e-01_dp,   7.8181818e-01_dp,   8.0000000e-01_dp,   8.0444444e-01_dp, &
     8.0888889e-01_dp,   8.1333333e-01_dp,   8.1777778e-01_dp,   8.2222222e-01_dp,   8.2666667e-01_dp, &
     8.3111111e-01_dp,   8.3555556e-01_dp,   8.4000000e-01_dp,   8.4444444e-01_dp,   8.4888889e-01_dp, &
     8.5333333e-01_dp,   8.5777778e-01_dp,   8.6222222e-01_dp,   8.6666667e-01_dp,   8.7111111e-01_dp, &
     8.7555556e-01_dp,   8.8000000e-01_dp,   8.8444444e-01_dp,   8.8888889e-01_dp,   8.9333333e-01_dp, &
     8.9777778e-01_dp,   9.0222222e-01_dp,   9.0666667e-01_dp,   9.1111111e-01_dp,   9.1555556e-01_dp, &
     9.2000000e-01_dp,   9.2444444e-01_dp,   9.2888889e-01_dp,   9.3333333e-01_dp,   9.3777778e-01_dp, &
     9.4222222e-01_dp,   9.4666667e-01_dp,   9.5111111e-01_dp,   9.5555556e-01_dp,   9.6000000e-01_dp, &
     9.6444444e-01_dp,   9.6888889e-01_dp,   9.7333333e-01_dp,   9.7777778e-01_dp,   9.8222222e-01_dp, &
     9.8666667e-01_dp,   9.9111111e-01_dp,   9.9555556e-01_dp,   1.0000000e+00_dp, 0.0_dp /)
 
  real(rprec), dimension(1:l) :: am_aux_f = (/ 1.1343073e+05_dp, &
     1.0511264e+05_dp,   9.7719901e+04_dp,   9.1142143e+04_dp,   8.5238372e+04_dp,   7.9893902e+04_dp, &
     7.5021684e+04_dp,   7.0555497e+04_dp,   6.6444168e+04_dp,   6.2647239e+04_dp,   5.9132114e+04_dp, &
     5.5871944e+04_dp,   5.2844250e+04_dp,   5.0029862e+04_dp,   4.7412214e+04_dp,   4.4976739e+04_dp, &
     4.2710495e+04_dp,   4.0601792e+04_dp,   3.8639948e+04_dp,   3.6815020e+04_dp,   3.5117615e+04_dp, &
     3.3538679e+04_dp,   3.2069304e+04_dp,   3.0700518e+04_dp,   2.9423071e+04_dp,   2.8227198e+04_dp, &
     2.7102369e+04_dp,   2.6037002e+04_dp,   2.5018180e+04_dp,   2.4031951e+04_dp,   2.3067082e+04_dp, &
     2.2117542e+04_dp,   2.1181454e+04_dp,   2.0260068e+04_dp,   1.9356991e+04_dp,   1.8477613e+04_dp, &
     1.7628651e+04_dp,   1.6817841e+04_dp,   1.6053656e+04_dp,   1.5345090e+04_dp,   1.4699130e+04_dp, &
     1.4113009e+04_dp,   1.3572603e+04_dp,   1.3055378e+04_dp,   1.2532538e+04_dp,   1.2400160e+04_dp, &
     1.2264938e+04_dp,   1.2126344e+04_dp,   1.1983859e+04_dp,   1.1836908e+04_dp,   1.1684914e+04_dp, &
     1.1527343e+04_dp,   1.1363692e+04_dp,   1.1193491e+04_dp,   1.1016295e+04_dp,   1.0831470e+04_dp, &
     1.0638256e+04_dp,   1.0435957e+04_dp,   1.0224325e+04_dp,   1.0003778e+04_dp,   9.7747634e+03_dp, &
     9.5363620e+03_dp,   9.2816431e+03_dp,   9.0022662e+03_dp,   8.6903284e+03_dp,   8.3409886e+03_dp, &
     7.9527438e+03_dp,   7.5237757e+03_dp,   7.0538828e+03_dp,   6.5554498e+03_dp,   6.0461036e+03_dp, &
     5.5435531e+03_dp,   5.0630772e+03_dp,   4.6140337e+03_dp,   4.2046228e+03_dp,   3.8422577e+03_dp, &
     3.5274335e+03_dp,   3.2543658e+03_dp,   3.0169701e+03_dp,   2.8089885e+03_dp,   2.6245185e+03_dp, &
     2.4577789e+03_dp,   2.3032608e+03_dp,   2.1572214e+03_dp,   2.0204623e+03_dp,   1.8947411e+03_dp, &
     1.7810514e+03_dp,   1.6799853e+03_dp,   1.5917883e+03_dp,   1.5163600e+03_dp, 0.0_dp /)
 
  integer :: i, j, iflag
  integer, parameter :: ns = 25
  real(rprec) :: hs, si
  real(rprec) :: y
 
  j = minloc(am_aux_s(2:), dim=1)
  print *, "length : ", j
 
  hs = 1.0_rprec/(ns-1)
 
  DO i = 2, ns
    si = hs * (i-1.5_rprec)
 
    call spline_akima(si, y, am_aux_s, am_aux_f, j, iflag)
 
    if (iflag .eq. 0) then
      print *, "si = ", si, " => y = ", y
    else
      print *, "error: iflag = ", iflag
    end if
 
  end do
 
end program hello
 
SUBROUTINE spline_akima(x, y, xx, yy, npts, iflag)
      implicit none
 
      INTEGER, PARAMETER :: rprec = selected_real_kind(12,100)
      INTEGER, PARAMETER :: dp = rprec
 
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer,     intent(in)                  :: npts
      real(rprec), intent(in)                  :: x
      real(rprec), intent(out)                 :: y
      real(rprec), dimension(npts), intent(in) :: xx
      real(rprec), dimension(npts), intent(in) :: yy
      integer,     intent(inout)               :: iflag
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      real(rprec), dimension(-1:size(xx)+2) :: a
      real(rprec), dimension(-1:size(xx)+2) :: b
      real(rprec), dimension(-1:size(xx)+2) :: c
      real(rprec), dimension(-1:size(xx)+2) :: d
 
      real(rprec), dimension(-1:size(xx)+2) :: xloc
      real(rprec), dimension(-1:size(xx)+2) :: yloc
 
      real(rprec), dimension(-1:size(xx)+2) :: m
      real(rprec), dimension(-1:size(xx)+2) :: t
      real(rprec), dimension(-1:size(xx)+2) :: dm
      real(rprec), dimension(-1:size(xx)+2) :: p
      real(rprec), dimension(-1:size(xx)+2) :: q
 
      integer :: i
      integer :: ix
      integer :: iv
      integer :: ivm1
 
      real(rprec) :: cl
      real(rprec) :: bl
      real(rprec) :: cr
      real(rprec) :: br
 
      real(rprec) :: dx
 
!-----------------------------------------------
 
      !WRITE(6,*)  x, y
      !WRITE(6,*)npts
      !WRITE(6,*) xx
      !WRITE(6,*) yy
      !stop 'Test!'
 
      iflag = 0  !initialization
 
      ix = size(xx)
!     print *, size(xx)
!     print *, size(yy)
      if(npts > ix) &
         stop'spline_akima: more active points requested than available'
      if(ix /= size(yy)) stop 'size mismatch of xx and yy!'
 
      iv=npts
 
      ivm1  = iv-1
! initialize local variables
      a = 0._dp ; b = 0._dp ; c = 0._dp ; d = 0._dp
      xloc = 0._dp ; yloc = 0._dp
      m = 0._dp ; t = 0._dp ; dm = 0._dp
      p = 0._dp ; q = 0._dp
 
      xloc(1:iv)=xx
      xloc(-1)= 2*xloc(1)-xloc(3)
      xloc( 0)= xloc(1)+xloc(2)-xloc(3)
      xloc(iv+2)= 2*xloc(iv)-xloc(iv-2)
      xloc(iv+1)= xloc(iv)+xloc(iv-1)-xloc(iv-2)
 
      !print *, xloc(-1), xloc(0), xloc(iv+1), xloc(iv+2)
 
      yloc(1:iv)=yy
 
! calculate linear derivatives as far as existent
      m(-1:iv+1) = (yloc(0:iv+2)-yloc(-1:iv+1))/ &
                   (xloc(0:iv+2)-xloc(-1:iv+1))
! values for i=0, -1 and iv, iv+1 by quadratic extrapolation:
      cl = (m(2)-m(1))/(xloc(3)-xloc(1))
      bl = m(1) - cl*(xloc(2)-xloc(1))
      cr = (m(iv-2)-m(iv-1))/(xloc(iv)-xloc(iv-2))
      br = m(iv-2) - cr*(xloc(iv-1)-xloc(iv-2))
 
      !print *, cl, bl, cr, br
 
      yloc( 0)=yloc(1)+bl*(xloc( 0)-xloc(1))+ &
                       cl*(xloc( 0)-xloc(1))**2
      yloc(-1)=yloc(1)+bl*(xloc(-1)-xloc(1))+ &
                       cl*(xloc(-1)-xloc(1))**2
      yloc(iv+1)=yloc(iv)+br*(xloc(iv+1)-xloc(iv))+ &
                          cl*(xloc(iv+1)-xloc(iv))**2
      yloc(iv+2)=yloc(iv)+br*(xloc(iv+2)-xloc(iv))+ &
                          cl*(xloc(iv+2)-xloc(iv))**2
 
      !print *, yloc(-1), yloc(0), yloc(iv+1), yloc(iv+2)
 
! rest of linear derivatives
      m(-1) = (yloc(0)-yloc(-1))/(xloc(0)-xloc(-1))
      m( 0) = (yloc(1)-yloc( 0))/(xloc(1)-xloc( 0))
      m(iv  ) = (yloc(iv+1)-yloc(iv  ))/(xloc(iv+1)-xloc(iv  ))
      m(iv+1) = (yloc(iv+2)-yloc(iv+1))/(xloc(iv+2)-xloc(iv+1))
 
      !print *, m(-1), m(0), m(iv), m(iv+1)
 
! calculate weights for derivatives
      dm(-1:iv)= abs(m(0:iv+1)-m(-1:iv))
      where (dm /= 0._dp) !exclude division by zero
        p(1:iv) = dm(1:iv)/(dm(1:iv)+dm(-1:iv-2))
      end where
      where (dm /= 0._dp) !exclude division by zero
        q(1:iv) = dm(-1:iv-2)/(dm(1:iv)+dm(-1:iv-2))
      end where
 
      print *, "    i           p(i)          q(i)"
      do i=1, iv
        print *, i, p(i), q(i)
      end do
 
      stop
 
      t(1:iv) = p(1:iv)*m(0:iv-1)+q(1:iv)*m(1:iv)
      where ( p(1:iv)+q(1:iv) < tiny(1._dp)) ! in case of two zeros give equal weight
        t(1:iv) = 0.5_dp*m(0:iv-1)+0.5_dp*m(1:iv)
      end where
! fix coefficients
      a = yloc
      b = t
      c(1:iv-1) = (3*m(1:iv-1)-t(2:iv)-2*t(1:iv-1))/ &
                     (xloc(2:iv)-xloc(1:iv-1))
      d(1:iv-1) = (t(2:iv)+t(1:iv-1)-2*m(1:iv-1))/ &
                     (xloc(2:iv)-xloc(1:iv-1))**2
 
! calculation
      if(x<xloc(1) .or. x>xloc(iv)) then
        y=0.0
        iflag=-1
        return
      endif
 
      if(x == xloc(iv)) then
        y = yy(iv)
        iflag = 0
        return
      endif
 
      do i=1,iv-1
        if(x >= xloc(i) .and. x < xloc(i+1))then
          dx=x-xloc(i)
          y=a(i)+dx*(b(i)+dx*(c(i)+d(i)*dx))
          iflag = 0
          return
        endif
      enddo
 
      END SUBROUTINE spline_akima