interp_cubic_spline Subroutine

public subroutine interp_cubic_spline(x_table, y_table, n_table, x_itp, y_itp, n_itp)

Arguments

Type IntentOptional Attributes Name
double precision :: x_table(n_table)
double precision :: y_table(n_table)
integer :: n_table
double precision :: x_itp(n_itp)
double precision :: y_itp(n_itp)
integer :: n_itp

Calls

proc~~interp_cubic_spline~~CallsGraph proc~interp_cubic_spline interp_cubic_spline proc~get_cubic_para get_cubic_para proc~interp_cubic_spline->proc~get_cubic_para proc~mpistop mpistop proc~interp_cubic_spline->proc~mpistop mpi_abort mpi_abort proc~mpistop->mpi_abort

Source Code

  subroutine interp_cubic_spline(x_table,y_table,n_table,x_itp,y_itp,n_itp)
    ! interpolation function fi=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3
    ! first order derivative and second order derivative is continous 
    ! 1D method

    integer :: n_table,n_itp
    double precision :: x_table(n_table),y_table(n_table)
    double precision :: x_itp(n_itp),y_itp(n_itp)

    double precision :: at(n_table),bt(n_table),ct(n_table),dt(n_table)
    double precision :: dxi
    integer :: i,j

    if (x_table(1)>x_itp(1) .or. x_table(n_table)<x_itp(n_itp)) then
      call mpistop("out of interpolation table")
    endif

    call get_cubic_para(x_table,y_table,at,bt,ct,dt,n_table)

    do i=1,n_itp
      LOOPT: do j=1,n_table-1
        if (x_itp(i)>=x_table(j) .and. x_itp(i)<x_table(j+1)) then
          dxi=x_itp(i)-x_table(j)
          y_itp(i)=at(j)+bt(j)*dxi+ct(j)*dxi**2+dt(j)*dxi**3
          exit LOOPT
        endif
      enddo LOOPT
    enddo

    if (x_itp(n_itp)==x_table(n_table)) y_itp(n_itp)=y_table(n_table)

  end subroutine interp_cubic_spline