lag_interp.F90

Go to the documentation of this file.

!> Fortran module to prepare for Lagrange polynomial interpolation.
!> based on GSI: lagmod.f90
 
module lag_interp_mod
 
use kinds, only: kind_real
use gnssro_mod_constants,only: zero, one
 
implicit none
 
! set default to private
private
! set subroutines to public
public :: lag_interp_const
public :: lag_interp_const_tl
public :: lag_interp_weights
public :: lag_interp_weights_tl
public :: lag_interp_smthweights
public :: lag_interp_smthweights_tl
 
contains
 
subroutine lag_interp_const(q,x,n)
! Precompute the N constant denominator factors of the N-point 
! Lagrange polynomial interpolation formula.
!
! input argument list:
!   X -    The N abscissae.
!   N -    The number of points involved.
!
! output argument list:
!   Q -    The N denominator constants.
IMPLICIT NONE
 
INTEGER          ,INTENT(in   ) :: n
REAL(kind_real),INTENT(  out) :: q(n)
REAL(kind_real),INTENT(in   ) :: x(n)
!-----------------------------------------------------------------------------
INTEGER                      :: i,j
!=============================================================================
DO i=1,n
   q(i)=one
   DO j=1,n
      IF(j /= i)q(i)=q(i)/(x(i)-x(j))
   ENDDO
ENDDO
end subroutine lag_interp_const
 
 
!============================================================================
subroutine lag_interp_const_tl(q,q_TL,x,x_TL,n)
!============================================================================
IMPLICIT NONE
 
INTEGER          ,INTENT(in   ) :: n !number of points involved
REAL(kind_real),DIMENSION(n),INTENT(  out) :: q,q_tl
REAL(kind_real),DIMENSION(n),INTENT(in   ) :: x,x_tl
!-----------------------------------------------------------------------------
INTEGER                      :: i,j
REAL(kind_real)                         :: rat
!=============================================================================
DO i=1,n
   q(i)=one
   q_tl(i)=zero
   DO j=1,n
      IF(j /= i) THEN
         rat=one/(x(i)-x(j))
         q_tl(i)=(q_tl(i)-q(i)*(x_tl(i)-x_tl(j))*rat)*rat
         q(i)=q(i)*rat
      ENDIF
   ENDDO
ENDDO
end subroutine lag_interp_const_tl
 
!============================================================================
subroutine lag_interp_weights(x,xt,q,w,dw,n)
! Construct the Lagrange weights and their derivatives when 
! target abscissa is known and denominators Q have already 
! been precomputed
!
! input argument list:
!   X   - Grid abscissae
!   XT  - Target abscissa
!   Q   - Q factors (denominators of the Lagrange weight formula)
!   N   - Number of grid points involved in the interpolation
!
! output argument list:
!   W   - Lagrange weights
!   DW  - Derivatives, dW/dX, of Lagrange weights W
!============================================================================
IMPLICIT NONE
 
INTEGER          ,INTENT(IN   ) :: n
REAL(kind_real)             ,INTENT(IN   ) :: xt
REAL(kind_real),INTENT(IN   ) :: x(n),q(n)
REAL(kind_real),INTENT(  OUT) :: w(n),dw(n)
!-----------------------------------------------------------------------------
REAL(kind_real)            :: d(n),pa(n),pb(n),dpa(n),dpb(n)
INTEGER                      :: j
!============================================================================
pa(1)=one
dpa(1)=zero
do j=1,n-1
   d(j)=xt-x(j)
   pa(j+1)=pa(j)*d(j)
   dpa(j+1)=dpa(j)*d(j)+pa(j)
enddo
d(n)=xt-x(n)
 
pb(n)=one
dpb(n)=zero
do j=n,2,-1
   pb(j-1)=pb(j)*d(j)
   dpb(j-1)=dpb(j)*d(j)+pb(j)
enddo
do j=1,n
   w(j)= pa(j)*pb(j)*q(j)
   dw(j)=(pa(j)*dpb(j)+dpa(j)*pb(j))*q(j)
enddo
end subroutine lag_interp_weights
 
 
!============================================================================
subroutine lag_interp_weights_tl(x,x_TL,xt,q,q_TL,w,w_TL,dw,dw_TL,n)
!============================================================================
IMPLICIT NONE
 
INTEGER          ,INTENT(IN   ) :: n
REAL(kind_real)             ,INTENT(IN   ) :: xt
REAL(kind_real),DIMENSION(n),INTENT(IN   ) :: x,q,x_tl,q_tl
REAL(kind_real),DIMENSION(n),INTENT(  OUT) :: w,dw,w_tl,dw_tl
!-----------------------------------------------------------------------------
REAL(kind_real),DIMENSION(n)            :: d,pa,pb,dpa,dpb
REAL(kind_real),DIMENSION(n)            :: d_tl,pa_tl,pb_tl,dpa_tl,dpb_tl
INTEGER                      :: j
!============================================================================
pa(1)=one
dpa(1)=zero
pa_tl(1)=zero
dpa_tl(1)=zero
 
do j=1,n-1
   d(j)=xt-x(j)
   d_tl(j)=-x_tl(j)
   pa(j+1)=pa(j)*d(j)
   pa_tl(j+1)=pa_tl(j)*d(j)+pa(j)*d_tl(j)
   dpa(j+1)=dpa(j)*d(j)+pa(j)
   dpa_tl(j+1)=dpa_tl(j)*d(j)+dpa(j)*d_tl(j)+pa_tl(j)
enddo
d(n)=xt-x(n)
d_tl(n)=-x_tl(n)
 
pb(n)=one
dpb(n)=zero
pb_tl(n)=zero
dpb_tl(n)=zero
do j=n,2,-1
   pb(j-1)=pb(j)*d(j)
   pb_tl(j-1)=pb_tl(j)*d(j)+pb(j)*d_tl(j)
   dpb(j-1)=dpb(j)*d(j)+pb(j)
   dpb_tl(j-1)=dpb_tl(j)*d(j)+dpb(j)*d_tl(j)+pb_tl(j)
enddo
do j=1,n
   w(j)= pa(j)*pb(j)*q(j)
   dw(j)=(pa(j)*dpb(j)+dpa(j)*pb(j))*q(j)
   w_tl(j)=(pa_tl(j)*pb(j)+pa(j)*pb_tl(j))*q(j)+pa(j)*pb(j)*q_tl(j)
   dw_tl(j)=(pa_tl(j)*dpb(j)+pa(j)*dpb_tl(j)+dpa_tl(j)*pb(j)+dpa(j)*pb_tl(j))*q(j)+ &
            (pa(j)*dpb(j)+dpa(j)*pb(j))*q_tl(j)
enddo
end subroutine lag_interp_weights_tl
 
!============================================================================
subroutine lag_interp_smthweights(x,xt,aq,bq,w,dw,n)
! Construct weights and their derivatives for interpolation 
! to a given target based on a linearly weighted mixture of 
! the pair of Lagrange interpolators which omit the respective 
! end points of the source nodes provided. The number of source 
! points provided must be even and the nominal target interval 
! is the unique central one. The linear weighting pair is 
! determined by the relative location of the target within 
! this central interval, or else the extreme values, 0 and 1, 
! when target lies outside this interval.  The objective is to 
! provide an interpolator whose derivative is continuous.
!============================================================================
IMPLICIT NONE
 
INTEGER            ,INTENT(IN   ) :: n
REAL(kind_real)               ,INTENT(IN   ) :: xt
REAL(kind_real),INTENT(IN   ) :: x(n)
REAL(kind_real),INTENT(IN   ) :: aq(n-1),bq(n-1)
REAL(kind_real),INTENT(  OUT) :: w(n),dw(n)
!-----------------------------------------------------------------------------
REAL(kind_real)               :: aw(n),bw(n),daw(n),dbw(n)
REAL(kind_real)               :: xa,xb,dwb,wb
INTEGER                       :: na
!============================================================================
CALL lag_interp_weights(x(1:n-1),xt,aq,aw(1:n-1),daw(1:n-1),n-1)
CALL lag_interp_weights(x(2:n  ),xt,bq,bw(2:n  ),dbw(2:n  ),n-1)
aw(n)=zero
daw(n)=zero
bw(1)=zero
dbw(1)=zero
na=n/2
IF(na*2 /= n)stop 'In lag_interp_smthWeights; n must be even'
xa =x(na     )
xb =x(na+1)
dwb=one/(xb-xa)
wb =(xt-xa)*dwb
IF(wb>one )THEN
   wb =one
   dwb=zero
ELSEIF(wb<zero)THEN
   wb =zero
   dwb=zero
ENDIF
 
bw =bw -aw
dbw=dbw-daw
 
w =aw +wb*bw
dw=daw+wb*dbw+dwb*bw
end subroutine lag_interp_smthweights
 
!============================================================================
subroutine lag_interp_smthweights_tl(x,x_TL,xt,aq,aq_TL,bq,bq_TL,dw,dw_TL,n)
!============================================================================
IMPLICIT NONE
 
INTEGER            ,INTENT(IN   ) :: n
REAL(kind_real)               ,INTENT(IN   ) :: xt
REAL(kind_real),DIMENSION(n)  ,INTENT(IN   ) :: x,x_tl
REAL(kind_real),DIMENSION(n-1),INTENT(IN   ) :: aq,bq,aq_tl,bq_tl
REAL(kind_real),DIMENSION(n)  ,INTENT(  OUT) :: dw_tl,dw
!-----------------------------------------------------------------------------
REAL(kind_real),DIMENSION(n)               :: aw,bw,daw,dbw
REAL(kind_real),DIMENSION(n)               :: aw_tl,bw_tl,daw_tl,dbw_tl
REAL(kind_real)                            :: xa,xb,dwb,wb
REAL(kind_real)                            :: xa_tl,xb_tl,dwb_tl,wb_tl
INTEGER                         :: na
!============================================================================
CALL lag_interp_weights_tl(x(1:n-1),x_tl(1:n-1),xt,aq,aq_tl,aw(1:n-1),aw_tl(1:n-1),&
             daw(1:n-1),daw_tl(1:n-1),n-1)
CALL lag_interp_weights_tl(x(2:n     ),x_tl(2:n     ),xt,bq,bq_tl,bw(2:n     ),bw_tl(2:n     ),&
             dbw(2:n     ),dbw_tl(2:n     ),n-1)
aw(n) =zero
daw(n)=zero
bw(1) =zero
dbw(1)=zero
!
aw_tl(n) =zero
daw_tl(n)=zero
bw_tl(1) =zero
dbw_tl(1)=zero
na=n/2
IF(na*2 /= n)stop 'In lag_interp_smthWeights; n must be even'
xa   =x(na)
xa_tl=x_tl(na)
xb   =x(na+1)
xb_tl=x_tl(na+1)
dwb   = one          /(xb-xa)
dwb_tl=-(xb_tl-xa_tl)/(xb-xa)**2
wb   =             (xt-xa)*dwb
wb_tl=(-xa_tl)*dwb+(xt-xa)*dwb_tl
IF(wb > one)THEN
   wb    =one
   dwb   =zero
   wb_tl =zero
   dwb_tl=zero
ELSEIF(wb < zero)THEN
   wb    =zero
   dwb   =zero
   wb_tl =zero
   dwb_tl=zero
ENDIF
 
bw    =bw    -aw
bw_tl =bw_tl -aw_tl
dbw   =dbw   -daw
dbw_tl=dbw_tl-daw_tl
 
!w=aw+wb*bw
dw   =daw   + wb   *dbw           + dwb   *bw
dw_tl=daw_tl+(wb_tl*dbw+wb*dbw_tl)+(dwb_tl*bw+dwb*bw_tl)
end subroutine lag_interp_smthweights_tl
 
end module lag_interp_mod