Cylindrical interpolation coefficients and indicies for a 3D grid
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=Float64), | intent(in) | :: | rmin | Minimum R |
||
real(kind=Float64), | intent(in) | :: | dr | R spacing |
||
integer, | intent(in) | :: | nr | Number of R points |
||
real(kind=Float64), | intent(in) | :: | zmin | Minimum Z |
||
real(kind=Float64), | intent(in) | :: | dz | Z spacing |
||
integer, | intent(in) | :: | nz | Number of Z points |
||
real(kind=Float64), | intent(in) | :: | phimin | Minimum phi |
||
real(kind=Float64), | intent(in) | :: | dphi | Phi spacing |
||
integer, | intent(in) | :: | nphi | Number of phi points |
||
real(kind=Float64), | intent(in) | :: | rout | R value to interpolate |
||
real(kind=Float64), | intent(in) | :: | zout | Z value to interpolate |
||
real(kind=Float64), | intent(in) | :: | phiout | Phi value to interpolate |
||
type(InterpolCoeffs3D), | intent(out) | :: | c | Interpolation Coefficients |
||
integer, | intent(out), | optional | :: | err | Error code |
subroutine cyl_interpol3D_coeff(rmin,dr,nr,zmin,dz,nz,phimin,dphi,nphi,rout,zout,phiout,c,err)
!+ Cylindrical interpolation coefficients and indicies for a 3D grid
real(Float64), intent(in) :: rmin
!+ Minimum R
real(Float64), intent(in) :: dr
!+ R spacing
integer, intent(in) :: nr
!+ Number of R points
real(Float64), intent(in) :: zmin
!+ Minimum Z
real(Float64), intent(in) :: dz
!+ Z spacing
integer, intent(in) :: nz
!+ Number of Z points
real(Float64), intent(in) :: phimin
!+ Minimum phi
real(Float64), intent(in) :: dphi
!+ Phi spacing
integer, intent(in) :: nphi
!+ Number of phi points
real(Float64), intent(in) :: rout
!+ R value to interpolate
real(Float64), intent(in) :: zout
!+ Z value to interpolate
real(Float64), intent(in) :: phiout
!+ Phi value to interpolate
type(InterpolCoeffs3D), intent(out) :: c
!+ Interpolation Coefficients
integer, intent(out), optional :: err
!+ Error code
type(InterpolCoeffs2D) :: b
real(Float64) :: r1, r2, phi1, phi2, z1, z2, rp, phip, zp, dV
real(Float64) :: phi
integer :: i, j, k, err_status
err_status = 1
rp = max(rout,rmin)
zp = max(zout,zmin)
phip = max(phiout,phimin)
i = floor((rp-rmin)/dr)+1
j = floor((zp-zmin)/dz)+1
k = floor((phip-phimin)/dphi)+1
if (nphi .eq. 1) then
if (((i.gt.0).and.(i.le.(nr-1))).and.((j.gt.0).and.(j.le.(nz-1)))) then
call interpol2D_coeff(rmin, dr, nr, zmin, dz, nz, rout, zout, b, err_status)
c%b111 = b%b11
c%b121 = b%b12
c%b221 = b%b22
c%b211 = b%b21
c%b212 = 0
c%b222 = 0
c%b122 = 0
c%b112 = 0
c%i = b%i
c%j = b%j
c%k = 1
err_status = 0
endif
else
if ((((i.gt.0).and.(i.le.(nr-1))).and.((j.gt.0).and.(j.le.(nz-1)))).and.((k.gt.0).and.(k.le.(nphi-1)))) then
r1 = rmin + (i-1)*dr
r2 = r1 + dr
z1 = zmin + (j-1)*dz
z2 = z1 + dz
phi1 = phimin + (k-1)*dphi
phi2 = phi1 + dphi
dV = ((r2**2 - r1**2) * (phi2 - phi1) * (z2 - z1))
!! Both volume elements have a factor of 1/2 that cancels out
c%b111 = ((r2**2 - rp**2) * (phi2 - phip) * (z2 - zp)) / dV
c%b121 = ((r2**2 - rp**2) * (phi2 - phip) * (zp - z1)) / dV
c%b221 = ((rp**2 - r1**2) * (phi2 - phip) * (zp - z1)) / dV
c%b211 = ((rp**2 - r1**2) * (phi2 - phip) * (z2 - zp)) / dV
c%b212 = ((rp**2 - r1**2) * (phip - phi1) * (z2 - zp)) / dV
c%b222 = ((rp**2 - r1**2) * (phip - phi1) * (zp - z1)) / dV
c%b122 = ((r2**2 - rp**2) * (phip - phi1) * (zp - z1)) / dV
c%b112 = ((r2**2 - rp**2) * (phip - phi1) * (z2 - zp)) / dV
c%i = i
c%j = j
c%k = k
err_status = 0
endif
endif
if(present(err)) err = err_status
end subroutine cyl_interpol3D_coeff