spline_akima.f

Go to the documentation of this file.

 
      SUBROUTINE spline_akima(x, y, xx, yy, npts, iflag)
      USE stel_kinds
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      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(in)                  :: npts
      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)
 
      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))
      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
! 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))
! 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
 
      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