-
Antoine Cyril David Hoffmann authoredAntoine Cyril David Hoffmann authored
lag_interp_mod.F90 6.67 KiB
!! This source is taken from GENE https://genecode.org/ !!
!>lagrange_interpolation contains subroutines to perform
!!a mid-point lagrange interpolation of order 3
MODULE lagrange_interpolation
USE prec_const
IMPLICIT NONE
PUBLIC:: lag3interp, lag3deriv, lag3interp_2d
PUBLIC:: lag3interp_complex
PRIVATE
INTERFACE lag3interp
MODULE PROCEDURE lag3interp_scalar, lag3interp_array
END INTERFACE
INTERFACE lag3deriv
MODULE PROCEDURE lag3deriv_scalar, lag3deriv_array
END INTERFACE
CONTAINS
!> Third order lagrange interpolation
SUBROUTINE lag3interp_scalar(y_in,x_in,n_in,y_out,x_out)
INTEGER, INTENT(IN) :: n_in
REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
REAL(dp), INTENT(IN) :: x_out
REAL(dp), INTENT(OUT) :: y_out
REAL(dp), DIMENSION(1) :: xout_wrap, yout_wrap
xout_wrap = x_out
call lag3interp_array(y_in,x_in,n_in,yout_wrap,xout_wrap,1)
y_out = yout_wrap(1)
END SUBROUTINE lag3interp_scalar
!> Third order lagrange interpolation
subroutine lag3interp_array(y_in,x_in,n_in,y_out,x_out,n_out)
INTEGER, INTENT(IN) :: n_in,n_out
REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
REAL(dp), DIMENSION(n_out), INTENT(OUT) :: y_out
REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
IF (x_in(n_in) > x_in(1)) THEN
jstart=3
jfirst=1
jlast=n_out
jstep=1
ELSE
jstart=n_in-2
jfirst=n_out
jlast=1
jstep=-1
END IF
j1=jstart
DO j=jfirst,jlast,jstep
x=x_out(j)
DO WHILE (x >= x_in(j1) .AND. j1 < n_in-1 .AND. j1 > 2)
j1=j1+jstep
END DO
j2=j1+jstep
j0=j1-jstep
jm=j1-2*jstep
!... extrapolate inside or outside
x2=x_in(j2)
x1=x_in(j1)
x0=x_in(j0)
xm=x_in(jm)
aintm=(x-x0)*(x-x1)*(x-x2)/((xm-x0)*(xm-x1)*(xm-x2))
aint0=(x-xm)*(x-x1)*(x-x2)/((x0-xm)*(x0-x1)*(x0-x2))
aint1=(x-xm)*(x-x0)*(x-x2)/((x1-xm)*(x1-x0)*(x1-x2))
aint2=(x-xm)*(x-x0)*(x-x1)/((x2-xm)*(x2-x0)*(x2-x1))
y_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
+aint1*y_in(j1)+aint2*y_in(j2)
END DO
END SUBROUTINE Lag3interp_array
!> Third order lagrange interpolation for complex arrays
SUBROUTINE lag3interp_complex(y_in,x_in,n_in,y_out,x_out,n_out)
INTEGER, INTENT(IN) :: n_in,n_out
COMPLEX, DIMENSION(n_in), INTENT(IN) :: y_in
REAL(dp), DIMENSION(n_in), INTENT(IN) :: x_in
COMPLEX, DIMENSION(n_out), INTENT(OUT) :: y_out
REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
IF (x_in(n_in) > x_in(1)) THEN
jstart=3
jfirst=1
jlast=n_out
jstep=1
ELSE
jstart=n_in-2
jfirst=n_out
jlast=1
jstep=-1
END IF
j1=jstart
DO j=jfirst,jlast,jstep
x=x_out(j)
DO WHILE (x >= x_in(j1) .AND. j1 < n_in-1 .AND. j1 > 2)
j1=j1+jstep
END DO
j2=j1+jstep
j0=j1-jstep
jm=j1-2*jstep
!... extrapolate inside or outside
x2=x_in(j2)
x1=x_in(j1)
x0=x_in(j0)
xm=x_in(jm)
aintm=(x-x0)*(x-x1)*(x-x2)/((xm-x0)*(xm-x1)*(xm-x2))
aint0=(x-xm)*(x-x1)*(x-x2)/((x0-xm)*(x0-x1)*(x0-x2))
aint1=(x-xm)*(x-x0)*(x-x2)/((x1-xm)*(x1-x0)*(x1-x2))
aint2=(x-xm)*(x-x0)*(x-x1)/((x2-xm)*(x2-x0)*(x2-x1))
y_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
+aint1*y_in(j1)+aint2*y_in(j2)
END DO
END SUBROUTINE lag3interp_complex
!>2D interpolation
!\TODO check whether a "REAL(dp)" 2D interpolation would
!! be more appropriate
SUBROUTINE lag3interp_2d(y_in,x1_in,n1_in,x2_in,n2_in,&
&y_out,x1_out,n1_out,x2_out,n2_out)
INTEGER, INTENT(IN) :: n1_in,n2_in,n1_out,n2_out
REAL(dp), DIMENSION(n1_in,n2_in), INTENT(IN) :: y_in
REAL(dp), DIMENSION(n1_in) :: x1_in
REAL(dp), DIMENSION(n2_in) :: x2_in
REAL(dp), DIMENSION(n1_out), INTENT(IN) :: x1_out
REAL(dp), DIMENSION(n2_out), INTENT(IN) :: x2_out
REAL(dp), DIMENSION(n1_out,n2_out), INTENT(OUT) :: y_out
!local variables
REAL(dp), DIMENSION(n2_in) :: y2_in_tmp
REAL(dp), DIMENSION(n2_out) :: y2_out_tmp
REAL(dp), DIMENSION(n1_in,n2_out) :: y_tmp
INTEGER :: i
DO i=1,n1_in
y2_in_tmp = y_in(i,:)
call lag3interp(y2_in_tmp,x2_in,n2_in,&
y2_out_tmp,x2_out,n2_out)
y_tmp(i,:) = y2_out_tmp
ENDDO
DO i=1,n2_out
call lag3interp(y_tmp(:,i),x1_in,n1_in,&
y_out(:,i),x1_out,n1_out)
END DO
END SUBROUTINE lag3interp_2d
!> Third order lagrange interpolation
SUBROUTINE lag3deriv_scalar(y_in,x_in,n_in,dydx_out,x_out)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n_in
REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
REAL(dp), INTENT(IN) :: x_out
REAL(dp), INTENT(OUT) :: dydx_out
REAL(dp), DIMENSION(1) :: xout_wrap, dydxout_wrap
xout_wrap = x_out
call lag3deriv_array(y_in,x_in,n_in,dydxout_wrap,xout_wrap,1)
dydx_out = dydxout_wrap(1)
END SUBROUTINE lag3deriv_scalar
!>Returns Derivative based on a 3rd order lagrange interpolation
subroutine lag3deriv_array(y_in,x_in,n_in,dydx_out,x_out,n_out)
INTEGER :: n_in,n_out
REAL(dp), DIMENSION(n_in), INTENT(IN) :: y_in,x_in
REAL(dp), DIMENSION(n_out), INTENT(IN) :: x_out
REAL(dp), DIMENSION(n_out), INTENT(OUT) :: dydx_out
REAL(dp) :: x,aintm,aint0,aint1,aint2,xm,x0,x1,x2
INTEGER :: j,jm,j0,j1,j2
INTEGER :: jstart,jfirst,jlast,jstep
IF (x_in(n_in) > x_in(1)) THEN
jstart=3 jfirst=1
jlast=n_out
jstep=1
ELSE
jstart=n_in-2
jfirst=n_out
jlast=1
jstep=-1
END IF
j1=jstart
DO j=jfirst,jlast,jstep
x=x_out(j)
DO WHILE (x >= x_in(j1) .AND. j1 < n_in-1 .AND. j1 > 2)
j1=j1+jstep
END DO
j2=j1+jstep
j0=j1-jstep
jm=j1-2*jstep
!... extrapolate inside or outside
x2=x_in(j2)
x1=x_in(j1)
x0=x_in(j0)
xm=x_in(jm)
aintm=((x-x1)*(x-x2)+(x-x0)*(x-x2)+(x-x0)*(x-x1)) &
/((xm-x0)*(xm-x1)*(xm-x2))
aint0=((x-x1)*(x-x2)+(x-xm)*(x-x2)+(x-xm)*(x-x1)) &
/((x0-xm)*(x0-x1)*(x0-x2))
aint1=((x-x0)*(x-x2)+(x-xm)*(x-x2)+(x-xm)*(x-x0)) &
/((x1-xm)*(x1-x0)*(x1-x2))
aint2=((x-x0)*(x-x1)+(x-xm)*(x-x1)+(x-xm)*(x-x0)) &
/((x2-xm)*(x2-x0)*(x2-x1))
dydx_out(j)=aintm*y_in(jm)+aint0*y_in(j0) &
+aint1*y_in(j1)+aint2*y_in(j2)
END DO
end subroutine Lag3deriv_array
end module lagrange_interpolation