! **************************************
! * Module phil0 *
! * R. J. Purser NOAA/NCEP/EMC 2009 *
! * jim.purser@noaa.gov *
! **************************************
!
! Module procedures pertaining to Hilbert curve transformations and
! representations.
!
! Direct dependencies
! Modules: kinds, pietc
!
!=============================================================================
module phil0
!=============================================================================
! Several different ways of expressing location are indicated in the names
! of the routines that perform the transformations from one representation
! to another. Single and double precsion versions are also indicated by the
! suffixes _s or _d, but mostly, only the _d versions are supported.
!
! The location representations are indicated as follows:
! r: Real parameter of the Hilbert curve
! hil#: Base-# expansion of parameter r, # =4, 8, 16, or 48 for mixed base 4,8
! dig#: Base-# digits in standard expansion of cartesian location in sq or cube
! xy*: Index (optional) of sq or cube, and Cartesian location within it.
! xs: 3D unit-Cartesian vector of point on the unit sphere
! ea: Equal-area coordinates of square map panel of cubed sphere
! gn: Gnomonic coordinates of square map panel of cubed sphere.
!
! The PUBLIC interfaces relate only to module procedures explicitly required
! by the general user; the PRIVATE ones are internally used as intermediate
! steps.
! PUBLIC:
! r_to_hil#, where # = 4,8,16: Encode the Hilbert parameter in base-#,
! or, when # = 48: mixed base 4 and 8.
! hil#_to_r: Inverse of r_to_hil# within the allowed precision.
! hil#_to_x*: Transform from Hilbert curve code to Cartesian location in
! unit square (#=4, x* = xy), cube (#=8, x*=xyz) or hypercube
! (#=16, x*=xyza).
! xy*_to_hil#: Inverse of hil#_to_xy*
! hil4_to_xs: Transform base-4 Hilbert parameter expansion to Cartesian
! 3-vector of an equal-area cubed sphere. The 0th digit of the
! Hilbert code expansion, hil, denotes the index, from 0 to 23,
! of the principal square of this mapping, these square being
! grouped in triplets around each of the 8 cube-vertices. (This
! convention was designed so that one hemisphere is filled by
! the first half of the Hilbert curve, the other hemisphere by
! the second half.)
! xs_to_hil4: The inverse of hil4_to_xs
! hil48_to_xs: Like hil4_to_xs except the mixed basis, 4, then 8, allows
! for the case of a spherical shell whose thickness, or 3rd
! dimension is invoked once the base becomes 8 instead of 4.
! But, in addition to unit 3-vector output, xs, there is also
! a real scalar output variable, r, that denotes a scaled
! radial, or vertical, coordinate variable whose resolution is
! the same as the horizontal resolution of xs when the latter
! is regarded as being composed of 24 deformed squares around
! the spherical surface.
! xs_to_hil48: Inverse of hil48_to_xs.
!
! PRIVATE:
! hil#_to_dig#: Convert from hilbert base-# expansion to ordinary Cartesian
! base-# expansion within a square (# =4), cube (# = 8) or
! hypercube (# = 16).
! dig#_to_hil#: Exact inverse of hil#_to_dig#
! dig#_to_x*: Base-# expansion (#=4, 8, 16) to real Cartesians within unit
! square (x* = x), cube (x* =xyz), hypercube (x* =xyza).
! x*_to_dig#: Inverse of dig#_to_x*.
! gn_to_ea: Cube face's gnomonic coordinates to equal-area coordinates
! ea_to_gn: Inverse of gn_to_ea.
!=============================================================================
use kinds, only: sp,dp,i_kind
use pietc, only: T,F,u0,u1,u2,o2,u4,o4,pi
implicit none
real(dp),parameter:: u8=8.0_dp,o8=u1/8.0_dp,u16=16.0_dp,o16=u1/16.0_dp,pio6=pi/6.0_dp
private
public r_to_hil4, r_to_hil8, r_to_hil16, r_to_hil48, &
hil4_to_r, hil8_to_r, hil16_to_r, hil48_to_r, &
hil4_to_xy, hil8_to_xyz, hil16_to_xyza, hil48_to_xyz, &
xy_to_hil4, xyz_to_hil8, xyza_to_hil16, xyz_to_hil48, &
hil4_to_xs, hil48_to_xs, &
xs_to_hil4, xs_to_hil48, &
! hil4_to_dig4, hil8_to_dig8, hil16_to_dig16, &
! dig4_to_hil4, dig8_to_hil8, dig16_to_hil16, &
! dig4_to_xy, dig8_to_xyz, dig16_to_xyza, &
! xy_to_dig4, xyz_to_dig8, xyza_to_dig16, &
! gn_to_ea,ea_to_gn, &
hil4_to_rz, hil8_to_rz, hil16_to_rz
!-- Public interfaces:
interface r_to_hil4; module procedure r_to_hil4_s,r_to_hil4_d; end interface
interface r_to_hil8; module procedure r_to_hil8_d; end interface
interface r_to_hil16; module procedure r_to_hil16_d; end interface
interface r_to_hil48; module procedure r_to_hil48_d; end interface
interface hil4_to_r; module procedure hil4_to_r_d; end interface
interface hil8_to_r; module procedure hil8_to_r_d; end interface
interface hil16_to_r; module procedure hil16_to_r_d; end interface
interface hil48_to_r; module procedure hil48_to_r_d; end interface
!..
interface hil4_to_xy; module procedure hil4_to_xy_d; end interface
interface hil8_to_xyz; module procedure hil8_to_xyz_d; end interface
interface hil16_to_xyza;module procedure hil16_to_xyza_d; end interface
interface hil48_to_xyz; module procedure hil48_to_xyz_d; end interface
!..
interface xy_to_hil4
module procedure xy_to_hil4z_s,xy_to_hil4z_d,xy_to_hil4_s,xy_to_hil4_d
end interface
interface xyz_to_hil8
module procedure xyz_to_hil8z_d,xyz_to_hil8_d
end interface
interface xyza_to_hil16
module procedure xyza_to_hil16z_d,xyza_to_hil16_d
end interface
interface xyz_to_hil48
module procedure xyz_to_hil48_d
end interface
interface hil4_to_xs; module procedure hil4_to_xs_d; end interface
interface hil48_to_xs; module procedure hil48_to_xs_d; end interface
interface xs_to_hil4; module procedure xs_to_hil4_d; end interface
interface xs_to_hil48; module procedure xs_to_hil48_d; end interface
!--
!-- private interfaces:
interface hil4_to_dig4; module procedure hil4_to_dig4; end interface
interface hil8_to_dig8; module procedure hil8_to_dig8; end interface
interface hil16_to_dig16; module procedure hil16_to_dig16; end interface
!..
interface dig4_to_hil4; module procedure dig4_to_hil4; end interface
interface dig8_to_hil8; module procedure dig8_to_hil8; end interface
interface dig16_to_hil16; module procedure dig16_to_hil16; end interface
!..
interface dig4_to_xy; module procedure dig4_to_xy_d; end interface
interface dig8_to_xyz; module procedure dig8_to_xyz_d; end interface
interface dig16_to_xyza; module procedure dig16_to_xyza_d; end interface
!..
interface xy_to_dig4
module procedure xy_to_dig4_s,xy_to_dig4_d
end interface
interface xyz_to_dig8
module procedure xyz_to_dig8_s,xyz_to_dig8_d
end interface
interface xyza_to_dig16
module procedure xyza_to_dig16_d
end interface
!..
interface gn_to_ea ; module procedure gn_to_ea_s,gn_to_ea_d; end interface
interface ea_to_gn ; module procedure ea_to_gn_s,ea_to_gn_d; end interface
!
!== Earlier versions of hil*_to_r, still public, but deprecated
! (scheduled for future deletion):
interface hil4_to_rz; module procedure hil4_to_rz_s,hil4_to_rz_d;end interface
interface hil8_to_rz; module procedure hil8_to_rz_d; end interface
interface hil16_to_rz; module procedure hil16_to_rz_d; end interface
contains
!=============================================================================
subroutine r_to_hil4_s(lgen,ngen,r,hil4)! [r_to_hil4]
!=============================================================================
! Take a real number r and peel off the base-4 digits from place lgen to ngen
! putting them into hil4 and leaving r as the remainder in [0,1). Note that
! by doing things this way, we can concatenate the similar operations
! and even change bases from one link to the next.
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
real(sp), intent(inout):: r
integer(i_kind),dimension(lgen:ngen),intent( out):: hil4
!-----------------------------------------------------------------------------
integer(i_kind):: i,j
!=============================================================================
do i=lgen,ngen
if(i>0)r=4.0_sp*r
j=r; r=r-j; hil4(i)=j
enddo
end subroutine r_to_hil4_s
!=============================================================================
subroutine r_to_hil4_d(lgen,ngen,r,hil4)! [r_to_hil4]
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
real(dp), intent(inout):: r
integer(i_kind),dimension(lgen:ngen),intent( out):: hil4
!-----------------------------------------------------------------------------
integer(i_kind):: i,j
!=============================================================================
do i=lgen,ngen
if(i>0)r=u4*r
j=r; r=r-j; hil4(i)=j
enddo
end subroutine r_to_hil4_d
!=============================================================================
subroutine r_to_hil8_d(lgen,ngen,r,hil8)! [r_to_hil8]
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
real(dp), intent(inout):: r
integer(i_kind),dimension(lgen:ngen),intent( out):: hil8
!-----------------------------------------------------------------------------
integer(i_kind):: i,j
!=============================================================================
do i=lgen,ngen
if(i>0)r=u8*r
j=r; r=r-j; hil8(i)=j
enddo
end subroutine r_to_hil8_d
!=============================================================================
subroutine r_to_hil16_d(lgen,ngen,r,hil16)! [r_to_hil16]
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
real(dp), intent(inout):: r
integer(i_kind),dimension(lgen:ngen),intent( out):: hil16
!-----------------------------------------------------------------------------
integer(i_kind):: i,j
!=============================================================================
do i=lgen,ngen
if(i>0)r=u16*r
j=r; r=r-j; hil16(i)=j
enddo
end subroutine r_to_hil16_d
!=============================================================================
subroutine r_to_hil48_d(lgen,ngen4,ngen48,r,hil)! [r_to_hil48]
!=============================================================================
implicit none
integer(i_kind), intent(in ):: lgen,ngen4,ngen48
real(dp), intent(inout):: r
integer(i_kind),dimension(lgen:ngen48),intent( out):: hil
!-----------------------------------------------------------------------------
call r_to_hil4(lgen, ngen4, r,hil(lgen:ngen4))
call r_to_hil8(ngen4+1,ngen48,r,hil(ngen4+1:ngen48))
end subroutine r_to_hil48_d
!=============================================================================
subroutine hil4_to_r_d(lgen,ngen,hil,r)! [hil4_to_r]
!=============================================================================
! Be sure to define r on input !
!=============================================================================
implicit none
integer(i_kind), intent(in ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(in ):: hil
real(dp), intent(inout):: r
!-----------------------------------------------------------------------------
integer(i_kind):: i
!=============================================================================
do i=ngen,lgen,-1
r=r+hil(i)
if(i==0)return
r=r*o4
enddo
end subroutine hil4_to_r_d
!=============================================================================
subroutine hil8_to_r_d(lgen,ngen,hil,r)! [hil8_to_r]
!=============================================================================
! Be sure to define r on input !
!=============================================================================
implicit none
integer(i_kind), intent(in ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(in ):: hil
real(dp), intent(inout):: r
!-----------------------------------------------------------------------------
integer(i_kind):: i
!=============================================================================
do i=ngen,lgen,-1
r=r+hil(i)
if(i==0)return
r=r*o8
enddo
end subroutine hil8_to_r_d
!=============================================================================
subroutine hil16_to_r_d(lgen,ngen,hil,r)! [hil16_to_r]
!=============================================================================
! Be sure to define r on input !
!=============================================================================
implicit none
integer(i_kind), intent(in ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(in ):: hil
real(dp), intent(inout):: r
!-----------------------------------------------------------------------------
integer(i_kind) :: i
!=============================================================================
do i=ngen,lgen,-1
r=r+hil(i)
if(i==0)return
r=r*o16
enddo
end subroutine hil16_to_r_d
!=============================================================================
subroutine hil48_to_r_d(lgen,ngen4,ngen48,hil,r)! [hil48_to_r]
!=============================================================================
! Be sure to define r on input !
!=============================================================================
implicit none
integer(i_kind), intent(in ):: lgen,ngen4,ngen48
integer(i_kind),dimension(lgen:ngen48),intent(in ):: hil
real(dp), intent(inout):: r
!=============================================================================
call hil8_to_r(ngen4+1,ngen48,hil(ngen4+1:ngen48),r)
call hil4_to_r(lgen, ngen4, hil(lgen:ngen4), r)
end subroutine hil48_to_r_d
!=============================================================================
subroutine hil4_to_xy_d(hil4,x,y)! [hil4_to_xy]
!=============================================================================
implicit none
integer(i_kind),dimension(:),intent(in ):: hil4
real(dp), intent(out):: x,y
!-----------------------------------------------------------------------------
real(dp) :: frac
integer(i_kind),dimension(size(hil4)):: dig4
integer(i_kind) :: presor,ngen
!============================================================================
dig4=hil4
presor=0; call hil4_to_dig4(presor,dig4)
call dig4_to_xy(dig4,x,y)
ngen=size(hil4)
frac=o2**(ngen+1); x=x+frac; y=y+frac ! <- to ensure final result is unbiased
end subroutine hil4_to_xy_d
!=============================================================================
subroutine hil8_to_xyz_d(hil8,x,y,z)! [hil8_to_xyz]
!=============================================================================
implicit none
integer(i_kind),dimension(:),intent(in ):: hil8
real(dp), intent(out):: x,y,z
!-----------------------------------------------------------------------------
real(dp) :: frac
integer(i_kind),dimension(size(hil8)):: dig8
integer(i_kind) :: presor,ngen
!============================================================================
dig8=hil8
presor=0; call hil8_to_dig8(presor,dig8)
call dig8_to_xyz(dig8,x,y,z)
ngen=size(hil8)
frac=o2**(ngen+1); x=x+frac; y=y+frac; z=z+frac ! <- final result is unbiased
end subroutine hil8_to_xyz_d
!=============================================================================
subroutine hil16_to_xyza_d(hil16,x,y,z,a)! [hil16_to_xyza]
!=============================================================================
implicit none
integer(i_kind),dimension(:),intent(in ):: hil16
real(dp), intent(out):: x,y,z,a
!-----------------------------------------------------------------------------
real(dp) :: frac
integer(i_kind),dimension(size(hil16)):: dig16
integer(i_kind) :: presor,ngen
!============================================================================
dig16=hil16
presor=0; call hil16_to_dig16(presor,dig16)
call dig16_to_xyza(dig16,x,y,z,a)
ngen=size(hil16)
frac=o2**(ngen+1); x=x+frac; y=y+frac; z=z+frac; a=a+frac ! result is unbiased
end subroutine hil16_to_xyza_d
!=============================================================================
subroutine hil48_to_xyz_d(ngen4,hil,x,y,z)! [hil48_to_xyz]
!=============================================================================
! Take a mixed (4 and 8) radix hilbert curve parameter, hil48, whose first
! ngen4 digits are the base-4 part, while the rest are base-8, and output
! the implied unbiased x,y,z coordinates, where x and y are in the unit square
! while z lies within [0,1/2**ngen4]. The resolution of the curve is
! 1/2**(ngen48) where ngen48=size(hil48), and to ensure the results are
! unbiased, the raw output from the final stage, dig8_to_xyz, has all the
! coordinates incremented by half the final resolution.
!=============================================================================
implicit none
integer(i_kind), intent(in ):: ngen4
integer(i_kind),dimension(:),intent(in ):: hil
real(dp), intent(out):: x,y,z
!-----------------------------------------------------------------------------
real(dp) :: frac8,frac,x4,y4
integer(i_kind),dimension(ngen4) :: dig4
integer(i_kind),dimension(size(hil)-ngen4):: dig8
integer(i_kind),dimension(0:7) :: p4to8
integer(i_kind) :: presor,ngen48
data p4to8/0,4,9,7,1,6,10,3/ ! Convert h4 orientation code to h8 code
!=============================================================================
ngen48=size(hil)
presor=0
if(ngen4>0)then ! Treat the radix-4 part (if there is one)
dig4=hil(1:ngen4)
call hil4_to_dig4(presor,dig4)
call dig4_to_xy(dig4,x4,y4)
else
x4=0.0_dp; y4=0.0_dp
endif
if(ngen48>ngen4)then ! Treat the radix-8 part (if there is one)
dig8=hil(ngen4+1:ngen48)
presor=p4to8(presor)
call hil8_to_dig8(presor,dig8)
call dig8_to_xyz(dig8,x,y,z)
else
x=0.0_dp; y=0.0_dp; z=0.0_dp
endif
frac8=o2**ngen4
frac =o2**(ngen48+1)
x=x4+frac8*x+frac; y=y4+frac8*y+frac; z= frac8*z+frac
end subroutine hil48_to_xyz_d
!=============================================================================
subroutine xy_to_hil4_s(x,y,hil4)! [xy_to_hil4]
!=============================================================================
! Convert an (x,y)-representation of a point in the proper interior of the
! unit square to an ngen-digit base-4 representation of the parameter of
! a space-filling Hilbert curve.
!=============================================================================
implicit none
real(sp), intent(IN ):: x,y
integer(i_kind),dimension(:),intent(OUT):: hil4
!-----------------------------------------------------------------------------
real(sp):: xr,yr
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y
call xy_to_dig4(xr,yr,hil4)
presor=0; call dig4_to_hil4(presor,hil4)
end subroutine xy_to_hil4_s
!=============================================================================
subroutine xy_to_hil4_d(x,y,hil4)! [xy_to_hil4]
!=============================================================================
! Convert an (x,y)-representation of a point in the proper interior of the
! unit square to an ngen-digit base-4 representation of the parameter of
! a space-filling Hilbert curve.
!=============================================================================
implicit none
real(dp), intent(IN ):: x,y
integer(i_kind),dimension(:),intent(OUT):: hil4
!-----------------------------------------------------------------------------
real(dp):: xr,yr
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y
call xy_to_dig4(xr,yr,hil4)
presor=0; call dig4_to_hil4(presor,hil4)
end subroutine xy_to_hil4_d
!=============================================================================
subroutine xyz_to_hil8_d(x,y,z,hil8)! [xyz_to_hil8]
!=============================================================================
implicit none
real(dp), intent(IN ):: x,y,z
integer(i_kind),dimension(:),intent(OUT):: hil8
!-----------------------------------------------------------------------------
real(dp):: xr,yr,zr
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y; zr=z
call xyz_to_dig8(xr,yr,zr,hil8)
presor=0; call dig8_to_hil8(presor,hil8)
end subroutine xyz_to_hil8_d
!=============================================================================
subroutine xyza_to_hil16_d(x,y,z,a,hil16)! [xyza_to_hil16]
!=============================================================================
implicit none
real(dp), intent(IN ):: x,y,z,a
integer(i_kind),dimension(:),intent(OUT):: hil16
!-----------------------------------------------------------------------------
real(dp):: xr,yr,zr,ar
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y; zr=z; ar=a
call xyza_to_dig16(xr,yr,zr,ar,hil16)
presor=0; call dig16_to_hil16(presor,hil16)
end subroutine xyza_to_hil16_d
!=============================================================================
subroutine xyz_to_hil48_d(ngen4,x,y,z,hil)! [xyz_to_hil48]
!=============================================================================
implicit none
integer(i_kind), intent(in ):: ngen4
real(dp), intent(in ):: x,y,z
integer(i_kind),dimension(:),intent(out):: hil
!-----------------------------------------------------------------------------
real(dp) :: xr,yr,zr
integer(i_kind),dimension(0:7):: p4to8
integer(i_kind) :: presor,ngen48
data p4to8/0,4,9,7,1,6,10,3/ ! Convert h4 orientation code to h8 code
!=============================================================================
ngen48=size(hil)
xr=x; yr=y; zr=z*u2**ngen4
presor=0
if(ngen4 >0 )then ! Treat radix-4 part (if any):
call xy_to_dig4 (xr,yr, hil(1:ngen4) )
call dig4_to_hil4(presor, hil(1:ngen4) )
endif
presor=p4to8(presor)
if(ngen48>ngen4)then ! Treat radix-8 part (if any):
call xyz_to_dig8(xr,yr,zr,hil(ngen4+1:ngen48))
call dig8_to_hil8(presor, hil(ngen4+1:ngen48))
endif
end subroutine xyz_to_hil48_d
!=============================================================================
subroutine hil4_to_xs_d(ngen,hil,xs)! [hil4_to_xs]
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: ngen
integer(i_kind),dimension(0:ngen),intent(IN ):: hil
real(dp),dimension(3), intent(OUT):: xs
!-----------------------------------------------------------------------------
integer(i_kind) :: m,m6
real(dp):: x,y,q,xx,yy
!=============================================================================
call hil4_to_xy(hil(1:ngen),x,y)
m=hil(0)
m6=mod(m,6)
select case(m6)
case(0); q=x; x=y; y=q
case(1); q=x; x=1.0_dp-y; y=q
case(2); x=1.0_dp-x; y=1.0_dp-y
case(3); x=x; y=1.0_dp-y
case(4); q=x; x=1.0_dp-y; y=1.0_dp-q
case(5); q=x; x=y; y=1.0_dp-q
end select
call ea_to_gn(x,y,xx,yy)
x=xx; y=yy
select case(m)
case(1,2,9,10,13,14,21,22); xs(1)= x
case(3,4,7, 8,15,16,19,20); xs(1)=-x
case(0,11,12,23); xs(1)= 1.0_dp
case(5, 6,17,18); xs(1)=-1.0_dp
end select
select case(m)
case(6,11,12,17); xs(2)= x
case(0,5,18,23); xs(2)=-x
case(7,10,13,16); xs(2)= y
case(1,4,19,22); xs(2)=-y
case(8,9,14,15); xs(2)= 1.0_dp
case(2,3,20,21); xs(2)=-1.0_dp
end select
select case(m)
case(0,2,3,5,6,8,9,11); xs(3)= y
case(12,14,15,17,18,20,21,23); xs(3)=-y
case(1,4,7,10); xs(3)= 1.0_dp
case(13,16,19,22); xs(3)=-1.0_dp
end select
q=sqrt( dot_product(xs,xs) ); xs=xs/q
end subroutine hil4_to_xs_d
!=============================================================================
subroutine hil48_to_xs_d(ngen4,ngen48,hil,xs,r)! [hil48_to_xs]
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: ngen4,ngen48
integer(i_kind),dimension(0:ngen48),intent(IN ):: hil
real(dp),dimension(3), intent(OUT):: xs
real(dp), intent(out):: r
!-----------------------------------------------------------------------------
integer(i_kind) :: m,m6
real(dp):: x,y,q,xx,yy
!=============================================================================
call hil48_to_xyz(ngen4,hil(1:ngen48),x,y,r)
m=hil(0); m6=mod(m,6)
select case(m6)
case(0); q=x; x=y; y=q
case(1); q=x; x=1.0_dp-y; y=q
case(2); x=1.0_dp-x; y=1.0_dp-y
case(3); x=x; y=1.0_dp-y
case(4); q=x; x=1.0_dp-y; y=1.0_dp-q
case(5); q=x; x=y; y=1.0_dp-q
end select
call ea_to_gn(x,y,xx,yy)
x=xx; y=yy
select case(m)
case(1,2,9,10,13,14,21,22); xs(1)= x
case(3,4,7, 8,15,16,19,20); xs(1)=-x
case(0,11,12,23); xs(1)= 1.0_dp
case(5, 6,17,18); xs(1)=-1.0_dp
end select
select case(m)
case(6,11,12,17); xs(2)= x
case(0,5,18,23); xs(2)=-x
case(7,10,13,16); xs(2)= y
case(1,4,19,22); xs(2)=-y
case(8,9,14,15); xs(2)= 1.0_dp
case(2,3,20,21); xs(2)=-1.0_dp
end select
select case(m)
case(0,2,3,5,6,8,9,11); xs(3)= y
case(12,14,15,17,18,20,21,23); xs(3)=-y
case(1,4,7,10); xs(3)= 1.0_dp
case(13,16,19,22); xs(3)=-1.0_dp
end select
q=sqrt( dot_product(xs,xs) ); xs=xs/q
end subroutine hil48_to_xs_d
!=============================================================================
subroutine xs_to_hil4_d(ngen,xs,hil)! [xs_to_hil4]
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: ngen
real(dp),dimension(3), intent(IN ):: xs
integer(i_kind), dimension(0:ngen),intent(OUT):: hil
!-----------------------------------------------------------------------------
real(dp) :: x,y,z,ax,ay,az,q,xx,yy
integer(i_kind) :: k,L,m,m6
integer(i_kind),dimension(0:7) :: Ltab,ktab
data Ltab/1,0,2,0,1,1,2,1/
data ktab/6,7,5,4,1,0,2,3/
!=============================================================================
x=xs(1); ax=abs(x)
y=xs(2); ay=abs(y)
z=xs(3); az=abs(z)
k=0; if( x>0 )k=k+1; if( y>0 )k=k+2; if( z>0 )k=k+4
L=0; if(ax>ay)L=L+1; if(ay>az)L=L+2; if(az>ax)L=L+4
L=Ltab(L)
k=ktab(k)
if(mod(k,2)==0)then; m=k*3+L
else; m=k*3+2-L
endif
hil(0)=m
select case(m)
case(0,5,6,11,12,17,18,23); x=ay/ax; y=az/ax
case(2,3,8,9,14,15,20,21); x=ax/ay; y=az/ay
case(1,4,7,10,13,16,19,22); x=ax/az; y=ay/az
end select
m6=mod(m,6)
call gn_to_ea(x,y,xx,yy); x=xx; y=yy
select case(m6)
case(0); q=x; x=y; y=q
case(1); q=x; x=y; y=1.0_dp-q
case(2); x=1.0_dp-x; y=1.0_dp-y
case(3); x=x; y=1.0_dp-y
case(4); q=x; x=1.0_dp-y; y=1.0_dp-q
case(5); q=x; x=1.0_dp-y; y=q
end select
call xy_to_hil4(x,y,hil(1:ngen))
end subroutine xs_to_hil4_d
!=============================================================================
subroutine xs_to_hil48_d(ngen4,ngen48,xs,r,hil)! [xs_to_hil48]
!=============================================================================
implicit none
integer(i_kind), intent(in ):: ngen4,ngen48
real(dp),dimension(3), intent(in ):: xs
real(dp), intent(in ):: r
integer(i_kind), dimension(0:ngen48),intent(out):: hil
!-----------------------------------------------------------------------------
real(dp) :: x,y,z,ax,ay,az,q,xx,yy
integer(i_kind) :: k,L,m,m6
integer(i_kind),dimension(0:7) :: Ltab,ktab
data Ltab/1,0,2,0,1,1,2,1/
data ktab/6,7,5,4,1,0,2,3/
!=============================================================================
x=xs(1); ax=abs(x)
y=xs(2); ay=abs(y)
z=xs(3); az=abs(z)
k=0; if( x>0 )k=k+1; if( y>0 )k=k+2; if( z>0 )k=k+4
L=0; if(ax>ay)L=L+1; if(ay>az)L=L+2; if(az>ax)L=L+4
L=Ltab(L)
k=ktab(k)
if(mod(k,2)==0)then; m=k*3+L
else; m=k*3+2-L
endif
hil(0)=m
select case(m)
case(0,5,6,11,12,17,18,23); x=ay/ax; y=az/ax
case(2,3,8,9,14,15,20,21); x=ax/ay; y=az/ay
case(1,4,7,10,13,16,19,22); x=ax/az; y=ay/az
end select
! Insert equal-area transformation of (x,y) here
m6=mod(m,6)
call gn_to_ea(x,y,xx,yy) ; x=xx; y=yy
select case(m6)
case(0); q=x; x=y; y=q
case(1); q=x; x=y; y=1.0_dp-q
case(2); x=1.0_dp-x; y=1.0_dp-y
case(3); x=x; y=1.0_dp-y
case(4); q=x; x=1.0_dp-y; y=1.0_dp-q
case(5); q=x; x=1.0_dp-y; y=q
end select
call xyz_to_hil48(ngen4,x,y,r,hil(1:ngen48))
end subroutine xs_to_hil48_d
!-----------------
! Private routines:
!=============================================================================
subroutine hil4_to_dig4(presor,hil4)! [hil4_to_dig4]
!=============================================================================
implicit none
integer(i_kind), intent(inout):: presor
integer(i_kind),dimension(:),intent(inout):: hil4
!-----------------------------------------------------------------------------
integer(i_kind),dimension(0:3, 0:7):: dtable,ptable
integer(i_kind) :: igen,h
data dtable/0,2,3,1, 1,0,2,3, 3,1,0,2, 2,3,1,0, &
0,1,3,2, 2,0,1,3, 3,2,0,1, 1,3,2,0/
data ptable/4,0,0,6, 7,1,1,5, 6,2,2,4, 5,3,3,7, &
0,4,4,2, 3,5,5,1, 2,6,6,0, 1,7,7,3/
!=============================================================================
do igen=1,size(hil4)
h=hil4(igen)
hil4(igen)=dtable(h,presor)
presor =ptable(h,presor)
enddo
end subroutine hil4_to_dig4
!=============================================================================
subroutine hil8_to_dig8(presor,hil8)! [hil8_to_dig8]
!=============================================================================
implicit none
integer(i_kind), intent(inout):: presor
integer(i_kind),dimension(:),intent(inout):: hil8
!-----------------------------------------------------------------------------
integer(i_kind),dimension(0:7,0:23):: dtable,ptable
integer(i_kind) :: igen,h
data dtable/ &
0,2,6,4,5,7,3,1, 0,4,5,1,3,7,6,2, 0,1,3,2,6,7,5,4, &
1,5,7,3,2,6,4,0, 1,0,4,5,7,6,2,3, 1,3,2,0,4,6,7,5, &
2,6,4,0,1,5,7,3, 2,3,7,6,4,5,1,0, 2,0,1,3,7,5,4,6, &
3,1,5,7,6,4,0,2, 3,7,6,2,0,4,5,1, 3,2,0,1,5,4,6,7, &
4,0,2,6,7,3,1,5, 4,5,1,0,2,3,7,6, 4,6,7,5,1,3,2,0, &
5,7,3,1,0,2,6,4, 5,1,0,4,6,2,3,7, 5,4,6,7,3,2,0,1, &
6,4,0,2,3,1,5,7, 6,2,3,7,5,1,0,4, 6,7,5,4,0,1,3,2, &
7,3,1,5,4,0,2,6, 7,6,2,3,1,0,4,5, 7,5,4,6,2,0,1,3/
data ptable/ &
1, 2, 2,18,18,17,17,10, 2, 0, 0,16,16, 9, 9,20, 0, 1, 1,11,11,19,19,15, &
5, 4, 4,21,21, 7, 7,14, 3, 5, 5,13,13,23,23, 6, 4, 9, 9, 8, 8,12,12,22, &
8, 7, 7,12,12, 4, 4,23, 6, 8, 8,22,22,14,14, 3, 7, 6, 6, 5, 5,21,21,13, &
10,11,11,15,15,20,20, 1, 11, 9, 9,19,19, 0, 0,17, 9,10,10, 2, 2,16,16,18, &
14,13,13, 6, 6,22,22, 5, 12,14,14, 4, 4, 8, 8,21, 13,12,12,23,23,33, 3, 7, &
16,17,17, 9, 9, 2, 2,19, 17,15,15, 1, 1,18,18,11, 15,16,16,20,20,10,10, 0, &
19,20,20, 0, 0,11,11,16, 20,18,18,10,10,15,15, 2, 18,19,19,17,17, 1, 1, 9, &
23,22,22, 3, 3,13,13, 8, 21,23,23, 7, 7, 5, 5,12, 22,21,21,14,14, 6, 6, 4/
!=============================================================================
do igen=1,size(hil8)
h=hil8(igen)
hil8(igen)=dtable(h,presor)
presor =ptable(h,presor)
enddo
end subroutine hil8_to_dig8
!=============================================================================
subroutine hil16_to_dig16(presor,hil16)! [hil16_to_dig16]
!=============================================================================
implicit none
integer(i_kind), intent(inout):: presor
integer(i_kind),dimension(:),intent(inout):: hil16
!-----------------------------------------------------------------------------
integer(i_kind),dimension(0:15,0:63):: dtable,ptable
integer(i_kind) :: igen,h
data dtable/ &
0, 8,12, 4, 6,14,10, 2, 3,11,15, 7, 5,13, 9, 1, &
0, 8, 9, 1, 5,13,12, 4, 6,14,15, 7, 3,11,10, 2, &
0, 8,10, 2, 3,11, 9, 1, 5,13,15, 7, 6,14,12, 4, &
0, 1, 3, 2, 6, 7, 5, 4,12,13,15,14,10,11, 9, 8, &
1, 9,11, 3, 7,15,13, 5, 4,12,14, 6, 2,10, 8, 0, &
1, 9,13, 5, 4,12, 8, 0, 2,10,14, 6, 7,15,11, 3, &
1, 9, 8, 0, 2,10,11, 3, 7,15,14, 6, 4,12,13, 5, &
1, 3, 2, 0, 4, 6, 7, 5,13,15,14,12, 8,10,11, 9, &
2,10, 8, 0, 4,12,14, 6, 7,15,13, 5, 1, 9,11, 3, &
2,10,14, 6, 7,15,11, 3, 1, 9,13, 5, 4,12, 8, 0, &
2,10,11, 3, 1, 9, 8, 0, 4,12,13, 5, 7,15,14, 6, &
2, 0, 1, 3, 7, 5, 4, 6,14,12,13,15,11, 9, 8,10, &
3,11,15, 7, 5,13, 9, 1, 0, 8,12, 4, 6,14,10, 2, &
3,11,10, 2, 6,14,15, 7, 5,13,12, 4, 0, 8, 9, 1, &
3,11, 9, 1, 0, 8,10, 2, 6,14,12, 4, 5,13,15, 7, &
3, 2, 0, 1, 5, 4, 6, 7,15,14,12,13, 9, 8,10,11, &
4,12,14, 6, 2,10, 8, 0, 1, 9,11, 3, 7,15,13, 5, &
4,12, 8, 0, 1, 9,13, 5, 7,15,11, 3, 2,10,14, 6, &
4,12,13, 5, 7,15,14, 6, 2,10,11, 3, 1, 9, 8, 0, &
4, 6, 7, 5, 1, 3, 2, 0, 8,10,11, 9,13,15,14,12, &
5,13, 9, 1, 3,11,15, 7, 6,14,10, 2, 0, 8,12, 4, &
5,13,12, 4, 0, 8, 9, 1, 3,11,10, 2, 6,14,15, 7, &
5,13,15, 7, 6,14,12, 4, 0, 8,10, 2, 3,11, 9, 1, &
5, 4, 6, 7, 3, 2, 0, 1, 9, 8,10,11,15,14,12,13, &
6,14,10, 2, 0, 8,12, 4, 5,13, 9, 1, 3,11,15, 7, &
6,14,15, 7, 3,11,10, 2, 0, 8, 9, 1, 5,13,12, 4, &
6,14,12, 4, 5,13,15, 7, 3,11, 9, 1, 0, 8,10, 2, &
6, 7, 5, 4, 0, 1, 3, 2,10,11, 9, 8,12,13,15,14, &
7,15,13, 5, 1, 9,11, 3, 2,10, 8, 0, 4,12,14, 6, &
7,15,11, 3, 2,10,14, 6, 4,12, 8, 0, 1, 9,13, 5, &
7,15,14, 6, 4,12,13, 5, 1, 9, 8, 0, 2,10,11, 3, &
7, 5, 4, 6, 2, 0, 1, 3,11, 9, 8,10,14,12,13,15, &
8, 0, 2,10,14, 6, 4,12,13, 5, 7,15,11, 3, 1, 9, &
8, 0, 4,12,13, 5, 1, 9,11, 3, 7,15,14, 6, 2,10, &
8, 0, 1, 9,11, 3, 2,10,14, 6, 7,15,13, 5, 4,12, &
8,10,11, 9,13,15,14,12, 4, 6, 7, 5, 1, 3, 2, 0, &
9, 1, 5,13,15, 7, 3,11,10, 2, 6,14,12, 4, 0, 8, &
9, 1, 0, 8,12, 4, 5,13,15, 7, 6,14,10, 2, 3,11, &
9, 1, 3,11,10, 2, 0, 8,12, 4, 6,14,15, 7, 5,13, &
9, 8,10,11,15,14,12,13, 5, 4, 6, 7, 3, 2, 0, 1, &
10, 2, 6,14,12, 4, 0, 8, 9, 1, 5,13,15, 7, 3,11, &
10, 2, 3,11,15, 7, 6,14,12, 4, 5,13, 9, 1, 0, 8, &
10, 2, 0, 8, 9, 1, 3,11,15, 7, 5,13,12, 4, 6,14, &
10,11, 9, 8,12,13,15,14, 6, 7, 5, 4, 0, 1, 3, 2, &
11, 3, 1, 9,13, 5, 7,15,14, 6, 4,12, 8, 0, 2,10, &
11, 3, 7,15,14, 6, 2,10, 8, 0, 4,12,13, 5, 1, 9, &
11, 3, 2,10, 8, 0, 1, 9,13, 5, 4,12,14, 6, 7,15, &
11, 9, 8,10,14,12,13,15, 7, 5, 4, 6, 2, 0, 1, 3, &
12, 4, 0, 8,10, 2, 6,14,15, 7, 3,11, 9, 1, 5,13, &
12, 4, 5,13, 9, 1, 0, 8,10, 2, 3,11,15, 7, 6,14, &
12, 4, 6,14,15, 7, 5,13, 9, 1, 3,11,10, 2, 0, 8, &
12,13,15,14,10,11, 9, 8, 0, 1, 3, 2, 6, 7, 5, 4, &
13, 5, 7,15,11, 3, 1, 9, 8, 0, 2,10,14, 6, 4,12, &
13, 5, 1, 9, 8, 0, 4,12,14, 6, 2,10,11, 3, 7,15, &
13, 5, 4,12,14, 6, 7,15,11, 3, 2,10, 8, 0, 1, 9, &
13,15,14,12, 8,10,11, 9, 1, 3, 2, 0, 4, 6, 7, 5, &
14, 6, 4,12, 8, 0, 2,10,11, 3, 1, 9,13, 5, 7,15, &
14, 6, 2,10,11, 3, 7,15,13, 5, 1, 9, 8, 0, 4,12, &
14, 6, 7,15,13, 5, 4,12, 8, 0, 1, 9,11, 3, 2,10, &
14,12,13,15,11, 9, 8,10, 2, 0, 1, 3, 7, 5, 4, 6, &
15, 7, 3,11, 9, 1, 5,13,12, 4, 0, 8,10, 2, 6,14, &
15, 7, 6,14,10, 2, 3,11, 9, 1, 0, 8,12, 4, 5,13, &
15, 7, 5,13,12, 4, 6,14,10, 2, 0, 8, 9, 1, 3,11, &
15,14,12,13, 9, 8,10,11, 3, 2, 0, 1, 5, 4, 6, 7/
data ptable/ &
3, 2, 2,49,49,26,26,40,40,14,14,61,61,22,22,39,&
3, 0, 0,38,38,20,20,49,49,24,24,62,62,12,12,43,&
3, 1, 1,40,40,13,13,38,38,21,21,60,60,25,25,51,&
0, 1, 1,14,14,25,25,23,23,49,49,62,62,41,41,36,&
7, 5, 5,46,46,29,29,52,52,17,17,58,58, 9, 9,35,&
7, 6, 6,52,52,18,18,33,33,10,10,56,56,30,30,47,&
7, 4, 4,33,33, 8, 8,46,46,28,28,57,57,16,16,55,&
5, 4, 4,10,10,16,16,31,31,52,52,58,58,32,32,45,&
11, 9, 9,34,34,17,17,56,56,29,29,54,54, 5, 5,47,&
11,10,10,56,56,30,30,45,45, 6, 6,52,52,18,18,35,&
11, 8, 8,45,45, 4, 4,34,34,16,16,53,53,28,28,59,&
9, 8, 8, 6, 6,28,28,19,19,56,56,54,54,44,44,33,&
15,14,14,61,61,22,22,36,36, 2, 2,49,49,26,26,43,&
15,12,12,42,42,24,24,61,61,20,20,50,50, 0, 0,39,&
15,13,13,36,36, 1, 1,42,42,25,25,48,48,21,21,63,&
12,13,13, 2, 2,21,21,27,27,61,61,50,50,37,37,40,&
19,17,17,58,58, 9, 9,32,32, 5, 5,46,46,29,29,55,&
19,18,18,32,32, 6, 6,53,53,30,30,44,44,10,10,59,&
19,16,16,53,53,28,28,58,58, 8, 8,45,45, 4, 4,35,&
17,16,16,30,30, 4, 4,11,11,32,32,46,46,52,52,57,&
23,22,22,37,37,14,14,60,60,26,26,41,41, 2, 2,51,&
23,20,20,50,50, 0, 0,37,37,12,12,42,42,24,24,63,&
23,21,21,60,60,25,25,50,50, 1, 1,40,40,13,13,39,&
20,21,21,26,26,13,13, 3, 3,37,37,42,42,61,61,48,&
27,26,26,41,41, 2, 2,48,48,22,22,37,37,14,14,63,&
27,24,24,62,62,12,12,41,41, 0, 0,38,38,20,20,51,&
27,25,25,48,48,21,21,62,62,13,13,36,36, 1, 1,43,&
24,25,25,22,22, 1, 1,15,15,41,41,38,38,49,49,60,&
31,29,29,54,54, 5, 5,44,44, 9, 9,34,34,17,17,59,&
31,30,30,44,44,10,10,57,57,18,18,32,32, 6, 6,55,&
31,28,28,57,57,16,16,54,54, 4, 4,33,33, 8, 8,47,&
29,28,28,18,18, 8, 8, 7, 7,44,44,34,34,56,56,53,&
35,33,33,10,10,57,57,16,16,53,53,30,30,45,45, 7,&
35,34,34,16,16,54,54, 5, 5,46,46,28,28,58,58,11,&
35,32,32, 5, 5,44,44,10,10,56,56,29,29,52,52,19,&
33,32,32,46,46,52,52,59,59,16,16,30,30, 4, 4, 9,&
39,38,38,21,21,62,62,12,12,42,42,25,25,50,50, 3,&
39,36,36, 2, 2,48,48,21,21,60,60,26,26,40,40,15,&
39,37,37,12,12,41,41, 2, 2,49,49,24,24,61,61,23,&
36,37,37,42,42,61,61,51,51,21,21,26,26,13,13, 0,&
43,42,42,25,25,50,50, 0, 0,38,38,21,21,62,62,15,&
43,40,40,14,14,60,60,25,25,48,48,22,22,36,36, 3,&
43,41,41, 0, 0,37,37,14,14,61,61,20,20,49,49,27,&
40,41,41,38,38,49,49,63,63,25,25,22,22, 1, 1,12,&
47,45,45, 6, 6,53,53,28,28,57,57,18,18,33,33,11,&
47,46,46,28,28,58,58, 9, 9,34,34,16,16,54,54, 7,&
47,44,44, 9, 9,32,32, 6, 6,52,52,17,17,56,56,31,&
45,44,44,34,34,56,56,55,55,28,28,18,18, 8, 8, 5,&
51,50,50, 1, 1,42,42,24,24,62,62,13,13,38,38,23,&
51,48,48,22,22,36,36, 1, 1,40,40,14,14,60,60,27,&
51,49,49,24,24,61,61,22,22,37,37,12,12,41,41, 3,&
48,49,49,62,62,41,41,39,39, 1, 1,14,14,25,25,20,&
55,53,53,30,30,45,45, 4, 4,33,33,10,10,57,57,19,&
55,54,54, 4, 4,34,34,17,17,58,58, 8, 8,46,46,31,&
55,52,52,17,17,56,56,30,30,44,44, 9, 9,32,32, 7,&
53,52,52,58,58,32,32,47,47, 4, 4,10,10,16,16,29,&
59,57,57,18,18,33,33, 8, 8,45,45, 6, 6,53,53,31,&
59,58,58, 8, 8,46,46,29,29,54,54, 4, 4,34,34,19,&
59,56,56,29,29,52,52,18,18,32,32, 5, 5,44,44,11,&
57,56,56,54,54,44,44,35,35, 8, 8, 6, 6,28,28,17,&
63,62,62,13,13,38,38,20,20,50,50, 1, 1,42,42,27,&
63,60,60,26,26,40,40,13,13,36,36, 2, 2,48,48,23,&
63,61,61,20,20,49,49,26,26,41,41, 0, 0,37,37,15,&
60,61,61,50,50,37,37,43,43,13,13, 2, 2,21,21,24/
!=============================================================================
do igen=1,size(hil16)
h=hil16(igen)
hil16(igen)=dtable(h,presor)
presor =ptable(h,presor)
enddo
end subroutine hil16_to_dig16
!=============================================================================
subroutine dig4_to_hil4(presor,hil4)! [dig4_to_hil4]
!=============================================================================
! On input, hil4 contains the ngen base-4 digits defining the location of the
! target point in the unit square, and on output, it will contain the first
! ngen base-4 digits of the corresponding parameter, in [0,1), of the
! Hilbert curve. On input, presor contains the orientation code of the whole
! curve; on output, it contains the orientation code of the final segment.
! The orientation code is identified with the oriented edge segments of the
! unit square in the following way:
! ORIENTATION CODE EDGE SEGMENT
! 0 (0,0)-->(1,0)
! 1 (1,0)-->(1,1)
! 2 (1,1)-->(0,1)
! 3 (0,1)-->(1,0)
! 4 (0,0)-->(0,1)
! 5 (0,1)-->(1,1)
! 6 (1,1)-->(1,0)
! 7 (1,0)-->(0,0)
!=============================================================================
implicit none
integer(i_kind), intent(inout):: presor
integer(i_kind),dimension(:),intent(inout):: hil4
!-----------------------------------------------------------------------------
integer(i_kind),dimension(0:3, 0:7) :: hil4table, presortable
integer(i_kind) :: igen,j
data hil4table /0,3,1,2, 1,0,2,3, 2,1,3,0, 3,2,0,1, &
0,1,3,2, 1,2,0,3, 2,3,1,0, 3,0,2,1/
data presortable/4,6,0,0, 1,7,1,5, 2,2,4,6, 7,3,5,3, &
0,4,2,4, 5,5,3,1, 6,0,6,2, 3,1,7,7/
!=============================================================================
! At successive refinements, update the present orientation and refine the
! segment on the hilbert curve according to the quadrant of the present
! square occupied by the next generation of refinement's square:
do igen=1,size(hil4)
j=hil4(igen)
hil4(igen)=hil4table (j,presor)
presor =presortable(j,presor)
enddo
end subroutine dig4_to_hil4
!=============================================================================
subroutine dig8_to_hil8(presor,hil8)! [dig8_to_hil8]
!=============================================================================
! Like dig4_to_hil4, but for a cube. Orientation codes now run from 0 thru
! 23 and are grouped in triplets according to the vertex the associated edge
! begins with; within each triplet, the orientations are in the x, y, z order.
! The path of 8 stations through the 2*2*2 cube associated with overall
! orientation p is given by the array, hil8table(:,p); by finding the "0"
! and the "7" in this path, one can locate the path start and endpoint, and
! hence interpret the orientation that p refers to.
!=============================================================================
implicit none
integer(i_kind), intent(inout):: presor
integer(i_kind),dimension(:),intent(inout):: hil8
!-----------------------------------------------------------------------------
integer(i_kind),dimension(0:7, 0:23) :: hil8table, presortable
integer(i_kind) :: igen,j
data hil8table / &
0,7,1,6,3,4,2,5, 0,3,7,4,1,2,6,5, 0,1,3,2,7,6,4,5, &
7,0,4,3,6,1,5,2, 1,0,6,7,2,3,5,4, 3,0,2,1,4,7,5,6, &
3,4,0,7,2,5,1,6, 7,6,0,1,4,5,3,2, 1,2,0,3,6,5,7,4, &
6,1,7,0,5,2,4,3, 4,7,3,0,5,6,2,1, 2,3,1,0,5,4,6,7, &
1,6,2,5,0,7,3,4, 3,2,4,5,0,1,7,6, 7,4,6,5,0,3,1,2, &
4,3,5,2,7,0,6,1, 2,1,5,6,3,0,4,7, 6,7,5,4,1,0,2,3, &
2,5,3,4,1,6,0,7, 6,5,1,2,7,4,0,3, 4,5,7,6,3,2,0,1, &
5,2,6,1,4,3,7,0, 5,4,2,3,6,7,1,0, 5,6,4,7,2,1,3,0/
data presortable/ &
1,10, 2,17,18,18, 2,17, 2,16,20,16, 0, 0, 9, 9, 0, 1,11, 1,15,19,11,19, &
14, 5,21,21, 7, 4, 7, 4, 5, 3,23, 6, 5,13,23,13, 8, 4, 9, 9, 8,22,12,12, &
12,12, 8,23, 7, 4, 7, 4, 3,14, 6, 8,22,14,22, 8, 6, 6, 7, 5,21,21,13, 5, &
20,11, 1,10,20,11,15,15, 19,17,19,11, 0, 0, 9, 9, 10, 2,10, 9,16, 2,16,18, &
13,22,13,22,14, 5, 6, 6, 4,14, 4, 8,12,14,21, 8, 7,23, 3, 3,13,23,12,12, &
9, 9, 2,17,19,16, 2,17, 15,15,18,18, 1,17, 1,11, 10, 0,10,20,16,15,16,20, &
20,11, 0, 0,20,11,19,16, 15,15,18,18, 2,10,20,10, 17, 1, 9, 1,17,19,18,19, &
13,22,13,22, 3, 3, 8,23, 5, 7,23, 7, 5,12,23,21, 6, 6,14, 4,21,21,14,22/
!=============================================================================
! At successive refinements, update the present orientation and refine the
! segment on the hilbert curve according to the quadrant of the present
! square occupied by the next generation of refinement's cube:
do igen=1,size(hil8)
j=hil8(igen)
hil8(igen)=hil8table (j,presor)
presor =presortable(j,presor)
enddo
end subroutine dig8_to_hil8
!=============================================================================
subroutine dig16_to_hil16(presor,hil16)! [dig16_to_hil16]
!=============================================================================
! Like dig4_to_hil4, but for a 4D hypercube.
!=============================================================================
implicit none
integer(i_kind), intent(inout):: presor
integer(i_kind),dimension(:),intent(inout):: hil16
!-----------------------------------------------------------------------------
integer(i_kind),dimension(0:15, 0:63) :: hil16table,presortable
integer(i_kind) :: igen,j
data hil16table / &
0,15,7,8,3,12,4,11,1,14,6,9,2,13,5,10, 0,3,15,12,7,4,8,11,1,2,14,13,6,5,9,10,&
0,7,3,4,15,8,12,11,1,6,2,5,14,9,13,10, 0,1,3,2,7,6,4,5,15,14,12,13,8,9,11,10,&
15,0,12,3,8,7,11,4,14,1,13,2,9,6,10,5, 7,0,8,15,4,3,11,12,6,1,9,14,5,2,10,13,&
3,0,4,7,12,15,11,8,2,1,5,6,13,14,10,9, 3,0,2,1,4,7,5,6,12,15,13,14,11,8,10,9,&
3,12,0,15,4,11,7,8,2,13,1,14,5,10,6,9, 15,8,0,7,12,11,3,4,14,9,1,6,13,10,2,5,&
7,4,0,3,8,11,15,12,6,5,1,2,9,10,14,13, 1,2,0,3,6,5,7,4,14,13,15,12,9,10,8,11,&
8,7,15,0,11,4,12,3,9,6,14,1,10,5,13,2, 12,15,3,0,11,8,4,7,13,14,2,1,10,9,5,6,&
4,3,7,0,11,12,8,15,5,2,6,1,10,13,9,14, 2,3,1,0,5,4,6,7,13,12,14,15,10,11,9,8,&
7,8,4,11,0,15,3,12,6,9,5,10,1,14,2,13, 3,4,12,11,0,7,15,8,2,5,13,10,1,6,14,9,&
15,12,8,11,0,3,7,4,14,13,9,10,1,2,6,5, 7,4,6,5,0,3,1,2,8,11,9,10,15,12,14,13,&
12,3,11,4,15,0,8,7,13,2,10,5,14,1,9,6, 4,7,11,8,3,0,12,15,5,6,10,9,2,1,13,14,&
8,15,11,12,7,0,4,3,9,14,10,13,6,1,5,2, 6,7,5,4,1,0,2,3,9,8,10,11,14,15,13,12,&
4,11,3,12,7,8,0,15,5,10,2,13,6,9,1,14, 8,11,7,4,15,12,0,3,9,10,6,5,14,13,1,2,&
12,11,15,8,3,4,0,7,13,10,14,9,2,5,1,6, 4,5,7,6,3,2,0,1,11,10,8,9,12,13,15,14,&
11,4,8,7,12,3,15,0,10,5,9,6,13,2,14,1, 11,12,4,3,8,15,7,0,10,13,5,2,9,14,6,1,&
11,8,12,15,4,7,3,0,10,9,13,14,5,6,2,1, 5,6,4,7,2,1,3,0,10,9,11,8,13,14,12,15,&
1,14,2,13,6,9,5,10,0,15,3,12,7,8,4,11, 1,6,14,9,2,5,13,10,0,7,15,8,3,4,12,11,&
1,2,6,5,14,13,9,10,0,3,7,4,15,12,8,11, 15,12,14,13,8,11,9,10,0,3,1,2,7,4,6,5,&
14,1,9,6,13,2,10,5,15,0,8,7,12,3,11,4, 2,1,13,14,5,6,10,9,3,0,12,15,4,7,11,8,&
6,1,5,2,9,14,10,13,7,0,4,3,8,15,11,12, 14,15,13,12,9,8,10,11,1,0,2,3,6,7,5,4,&
6,9,1,14,5,10,2,13,7,8,0,15,4,11,3,12, 14,13,1,2,9,10,6,5,15,12,0,3,8,11,7,4,&
2,5,1,6,13,10,14,9,3,4,0,7,12,11,15,8, 12,13,15,14,11,10,8,9,3,2,0,1,4,5,7,6,&
13,2,14,1,10,5,9,6,12,3,15,0,11,4,8,7, 9,14,6,1,10,13,5,2,8,15,7,0,11,12,4,3,&
5,6,2,1,10,9,13,14,4,7,3,0,11,8,12,15, 13,14,12,15,10,9,11,8,2,1,3,0,5,6,4,7,&
2,13,5,10,1,14,6,9,3,12,4,11,0,15,7,8, 6,5,9,10,1,2,14,13,7,4,8,11,0,3,15,12,&
14,9,13,10,1,6,2,5,15,8,12,11,0,7,3,4, 8,9,11,10,15,14,12,13,7,6,4,5,0,1,3,2,&
9,6,10,5,14,1,13,2,8,7,11,4,15,0,12,3, 5,2,10,13,6,1,9,14,4,3,11,12,7,0,8,15,&
13,14,10,9,2,1,5,6,12,15,11,8,3,0,4,7, 11,8,10,9,12,15,13,14,4,7,5,6,3,0,2,1,&
5,10,6,9,2,13,1,14,4,11,7,8,3,12,0,15, 13,10,2,5,14,9,1,6,12,11,3,4,15,8,0,7,&
9,10,14,13,6,5,1,2,8,11,15,12,7,4,0,3, 9,10,8,11,14,13,15,12,6,5,7,4,1,2,0,3,&
10,5,13,2,9,6,14,1,11,4,12,3,8,7,15,0, 10,9,5,6,13,14,2,1,11,8,4,7,12,15,3,0,&
10,13,9,14,5,2,6,1,11,12,8,15,4,3,7,0, 10,11,9,8,13,12,14,15,5,4,6,7,2,3,1,0/
data presortable/ &
3,39,40,40,49,61,49,61, 2,22,26,14, 2,22,26,14,&
3,38,43,62,49,38,49,62, 0, 0,12,12,20,20,24,24,&
3,38,40,40,51,38,60,60, 1,13, 1,13,25,21,25,21,&
0, 1,14, 1,23,25,14,25,36,41,62,41,23,49,62,49,&
35, 7,58,46,52,52,58,46, 9, 5, 9, 5,17,29,17,29,&
33, 7,33,47,52,52,56,56,18, 6,10,30,18, 6,10,30,&
33, 7,33,46,57,55,57,46, 4, 4, 8, 8,16,16,28,28,&
10, 5, 4, 4,10,31,16,16,58,45,32,32,58,31,52,52,&
34,54,11,47,34,54,56,56, 9, 5, 9, 5,17,29,17,29,&
35,45,11,45,52,52,56,56,18, 6,10,30,18, 6,10,30,&
34,45,11,45,34,53,59,53, 4, 4, 8, 8,16,16,28,28,&
8, 8, 9, 6,28,28,19, 6,44,44,33,54,56,56,19,54,&
36,36,43,15,49,61,49,61, 2,22,26,14, 2,22,26,14,&
50,39,42,15,50,61,42,61, 0, 0,12,12,20,20,24,24,&
36,36,42,15,48,48,42,63, 1,13, 1,13,25,21,25,21,&
13, 2,13,12,21, 2,21,27,37,50,37,40,61,50,61,27,&
32,32,58,46,19,55,58,46, 9, 5, 9, 5,17,29,17,29,&
32,32,44,44,19,53,59,53,18, 6,10,30,18, 6,10,30,&
35,45,58,45,19,53,58,53, 4, 4, 8, 8,16,16,28,28,&
11,30, 4, 4,17,30,16,16,11,46,32,32,57,46,52,52,&
41,37,41,37,51,23,60,60, 2,22,26,14, 2,22,26,14,&
50,37,42,37,50,23,42,63, 0, 0,12,12,20,20,24,24,&
50,39,40,40,50,23,60,60, 1,13, 1,13,25,21,25,21,&
13, 3,13,26,21,20,21,26,37, 3,37,42,61,48,61,42,&
41,37,41,37,48,48,27,63, 2,22,26,14, 2,22,26,14,&
41,38,41,62,51,38,27,62, 0, 0,12,12,20,20,24,24,&
36,36,43,62,48,48,27,62, 1,13, 1,13,25,21,25,21,&
22, 1,15, 1,22,25,24,25,38,41,15,41,38,49,60,49,&
34,54,44,44,34,54,59,31, 9, 5, 9, 5,17,29,17,29,&
32,32,44,44,57,55,57,31,18, 6,10,30,18, 6,10,30,&
33,54,33,47,57,54,57,31, 4, 4, 8, 8,16,16,28,28,&
8, 8,18, 7,28,28,18,29,44,44,34, 7,56,56,34,53,&
33,45,33,45,57,53,57,53,35, 7,10,30,16,16,10,30,&
34,54,58,46,34,54,58,46,35, 5,11, 5,16,16,28,28,&
32,32,44,44,52,52,56,56,35, 5,10, 5,19,29,10,29,&
9,30, 4, 4,59,30,16,16,33,46,32,32,59,46,52,52,&
50,38,42,62,50,38,42,62, 3,39,12,12,25,21,25,21,&
36,36,40,40,48,48,60,60, 2,39,26,15, 2,21,26,21,&
41,37,41,37,49,61,49,61, 2,39,12,12, 2,23,24,24,&
13, 0,13,26,21,51,21,26,37,36,37,42,61,51,61,42,&
50,38,42,62,50,38,42,62, 0, 0,43,15,25,21,25,21,&
36,36,40,40,48,48,60,60, 3,22,43,14,25,22,25,14,&
41,37,41,37,49,61,49,61, 0, 0,43,14,20,20,27,14,&
22, 1,12, 1,22,25,63,25,38,41,40,41,38,49,63,49,&
33,45,33,45,57,53,57,53,18, 6,11,47,18, 6,28,28,&
34,54,58,46,34,54,58,46, 9, 7, 9,47,16,16,28,28,&
32,32,44,44,52,52,56,56, 9, 6, 9,47,17, 6,17,31,&
8, 8,18, 5,28,28,18,55,44,44,34,45,56,56,34,55,&
50,38,42,62,50,38,42,62, 1,13, 1,13,51,23,24,24,&
36,36,40,40,48,48,60,60, 1,22, 1,14,51,22,27,14,&
41,37,41,37,49,61,49,61, 3,22,12,12,51,22,24,24,&
39, 1,14, 1,20,25,14,25,39,41,62,41,48,49,62,49,&
33,45,33,45,57,53,57,53, 4, 4,10,30,19,55,10,30,&
34,54,58,46,34,54,58,46, 4, 4, 8, 8,17,55,17,31,&
32,32,44,44,52,52,56,56, 9, 7, 9,30,17,55,17,30,&
10,47, 4, 4,10,29,16,16,58,47,32,32,58,53,52,52,&
33,45,33,45,57,53,57,53,18, 6, 8, 8,18, 6,59,31,&
34,54,58,46,34,54,58,46, 4, 4, 8, 8,19,29,59,29,&
32,32,44,44,52,52,56,56,18, 5,11, 5,18,29,59,29,&
8, 8,35, 6,28,28,17, 6,44,44,35,54,56,56,57,54,&
50,38,42,62,50,38,42,62, 1,13, 1,13,20,20,27,63,&
36,36,40,40,48,48,60,60, 2,13,26,13, 2,23,26,63,&
41,37,41,37,49,61,49,61, 0, 0,26,15,20,20,26,63,&
13, 2,13,43,21, 2,21,24,37,50,37,43,61,50,61,60/
!=============================================================================
! At successive refinements, update the present orientation and refine the
! segment on the hilbert curve according to the quadrant of the present
! square occupied by the next generation of refinement's hypercube:
do igen=1,size(hil16)
j=hil16(igen)
hil16(igen)=hil16table (j,presor)
presor =presortable(j,presor)
enddo
end subroutine dig16_to_hil16
!=============================================================================
subroutine dig4_to_xy_d(dig,x,y)! [dig4_to_xy]
!=============================================================================
implicit none
integer(i_kind),dimension(:),intent(in ):: dig
real(dp), intent(out):: x,y
!-----------------------------------------------------------------------------
real(dp) :: s
logical,dimension(0:3):: xofd4,yofd4
integer(i_kind) :: igen,d
data xofd4/F,T,F,T/
data yofd4/F,F,T,T/
!=============================================================================
x=u0; y=u0; s=u1
do igen=size(dig),1,-1
d=dig(igen); if(xofd4(d))x=x+s; if(yofd4(d))y=y+s
s=s*u2
enddo
x=x/s; y=y/s
end subroutine dig4_to_xy_d
!=============================================================================
subroutine dig8_to_xyz_d(dig,x,y,z)! [dig8_to_xyz]
!=============================================================================
implicit none
integer(i_kind),dimension(:),intent(in ):: dig
real(dp), intent(out):: x,y,z
!-----------------------------------------------------------------------------
real(dp) :: s
logical,dimension(0:7):: xofd8,yofd8,zofd8
integer(i_kind) :: igen,d
data xofd8/F,T,F,T,F,T,F,T/
data yofd8/F,F,T,T,F,F,T,T/
data zofd8/F,F,F,F,T,T,T,T/
!=============================================================================
x=u0; y=u0; z=u0; s=u1
do igen=size(dig),1,-1
d=dig(igen); if(xofd8(d))x=x+s; if(yofd8(d))y=y+s; if(zofd8(d))z=z+s
s=s*u2
enddo
x=x/s; y=y/s; z=z/s
end subroutine dig8_to_xyz_d
!=============================================================================
subroutine dig16_to_xyza_d(dig,x,y,z,a)! [dig16_to_xyza]
!=============================================================================
implicit none
integer(i_kind),dimension(:),intent(in ):: dig
real(dp), intent(out):: x,y,z,a
!-----------------------------------------------------------------------------
real(dp) :: s
logical,dimension(0:15):: xofd16,yofd16,zofd16,aofd16
integer(i_kind) :: igen,d
data xofd16/F,T,F,T,F,T,F,T,F,T,F,T,F,T,F,T/
data yofd16/F,F,T,T,F,F,T,T,F,F,T,T,F,F,T,T/
data zofd16/F,F,F,F,T,T,T,T,F,F,F,F,T,T,T,T/
data aofd16/F,F,F,F,F,F,F,F,T,T,T,T,T,T,T,T/
!=============================================================================
x=u0; y=u0; z=u0; a=u0; s=u1
do igen=size(dig),1,-1
d=dig(igen)
if(xofd16(d))x=x+s;if(yofd16(d))y=y+s;if(zofd16(d))z=z+s;if(aofd16(d))a=a+s
s=s*u2
enddo
x=x/s; y=y/s; z=z/s; a=a/s
end subroutine dig16_to_xyza_d
!=============================================================================
subroutine xy_to_dig4_s(x,y,dig4)! [xy_to_dig4]
!=============================================================================
! Convert an (x,y)-representation of a point in the square, [0,1]*[0,1]
! to an ngen-digit base-4 number, dig4, where ngen is the size of array dig4.
!=============================================================================
implicit none
real(sp), intent(inout):: x,y
integer(i_kind),dimension(:),intent( out):: dig4
!-----------------------------------------------------------------------------
integer(i_kind):: igen
!=============================================================================
if(x< 0.0_sp)stop 'In xy_to_dig4; x< 0.0_sp'
if(x> 1.0_sp)stop 'In xy_to_dig4; x> 1.0_sp'
if(y< 0.0_sp)stop 'In xy_to_dig4; y< 0.0_sp'
if(y> 1.0_sp)stop 'In xy_to_dig4; y> 1.0_sp'
dig4=0
do igen=1,size(dig4)
x=x*2.0_sp; y=y*2.0_sp
if(x>=1.0_sp)then; dig4(igen)=dig4(igen)+1; x=x-1.0_sp; endif
if(y>=1.0_sp)then; dig4(igen)=dig4(igen)+2; y=y-1.0_sp; endif
enddo
end subroutine xy_to_dig4_s
!=============================================================================
subroutine xy_to_dig4_d(x,y,dig4)! [xy_to_dig4]
!=============================================================================
implicit none
real(dp), intent(inout):: x,y
integer(i_kind),dimension(:),intent( out):: dig4
!-----------------------------------------------------------------------------
integer(i_kind):: igen
!=============================================================================
if(x< u0)stop 'In xy_to_dig4; x< 0.0_dp'
if(x> u1)stop 'In xy_to_dig4; x> 1.0_dp'
if(y< u0)stop 'In xy_to_dig4; y< 0.0_dp'
if(y> u1)stop 'In xy_to_dig4; y> 1.0_dp'
dig4=0
do igen=1,size(dig4)
x=x*u2; y=y*u2
if(x>=u1)then; dig4(igen)=dig4(igen)+1; x=x-u1; endif
if(y>=u1)then; dig4(igen)=dig4(igen)+2; y=y-u1; endif
enddo
end subroutine xy_to_dig4_d
!=============================================================================
subroutine xyz_to_dig8_s(x,y,z,dig8)! [xyz_to_dig8]
!=============================================================================
! Convert an (x,y,z)-representation of a point in the cube, [0,1]*[0,1]*[0,1]
! to an ngen-digit base-8 number, dig8.
!=============================================================================
implicit none
real(sp), intent(inout):: x,y,z
integer(i_kind),dimension(:),intent( out):: dig8
!-----------------------------------------------------------------------------
integer(i_kind):: igen
!=============================================================================
if(x< 0.0_sp)stop 'In xyz_to_dig8; x< 0.0_sp'
if(x> 1.0_sp)stop 'In xyz_to_dig8; x> 1.0_sp'
if(y< 0.0_sp)stop 'In xyz_to_dig8; y< 0.0_sp'
if(y> 1.0_sp)stop 'In xyz_to_dig8; y> 1.0_sp'
if(z< 0.0_sp)stop 'In xyz_to_dig8; z< 0.0_sp'
if(z> 1.0_sp)stop 'In xyz_to_dig8; z> 1.0_sp'
dig8=0
do igen=1,size(dig8)
x=x*2.0_sp; y=y*2.0_sp; z=z*2.0_sp
if(x>=1.0_sp)then; dig8(igen)=dig8(igen)+1; x=x-1.0_sp; endif
if(y>=1.0_sp)then; dig8(igen)=dig8(igen)+2; y=y-1.0_sp; endif
if(z>=1.0_sp)then; dig8(igen)=dig8(igen)+4; z=z-1.0_sp; endif
enddo
end subroutine xyz_to_dig8_s
!=============================================================================
subroutine xyz_to_dig8_d(x,y,z,dig8)! [xyz_to_dig8]
!=============================================================================
implicit none
real(dp), intent(inout):: x,y,z
integer(i_kind),dimension(:),intent( out):: dig8
!-----------------------------------------------------------------------------
integer(i_kind):: igen
!=============================================================================
if(x< u0)stop 'In xyz_to_dig8; x< 0.0_dp'
if(x> u1)stop 'In xyz_to_dig8; x> 1_0_dp'
if(y< u0)stop 'In xyz_to_dig8; y< 0_0_dp'
if(y> u1)stop 'In xyz_to_dig8; y> 1_0_dp'
if(z< u0)stop 'In xyz_to_dig8; z< 0_0_dp'
if(z> u1)stop 'In xyz_to_dig8; z> 1_0_dp'
dig8=0
do igen=1,size(dig8)
x=x*u2; y=y*u2; z=z*u2
if(x>=u1)then; dig8(igen)=dig8(igen)+1; x=x-u1; endif
if(y>=u1)then; dig8(igen)=dig8(igen)+2; y=y-u1; endif
if(z>=u1)then; dig8(igen)=dig8(igen)+4; z=z-u1; endif
enddo
end subroutine xyz_to_dig8_d
!=============================================================================
subroutine xyza_to_dig16_d(x,y,z,a,dig16)! [xyza_to_dig16]
!=============================================================================
! Convert an (x,y,z,t)-representation of a point in the hypercube,
! [0,1]*[0,1]*[0,1]*[0,1]
! to an ngen-digit base-16 number, dig16.
!=============================================================================
implicit none
real(dp), intent(inout):: x,y,z,a
integer(i_kind),dimension(:),intent( out):: dig16
!-----------------------------------------------------------------------------
integer(i_kind):: igen
!=============================================================================
if(x< u0)stop 'In xyza_to_dig16; x< 0_0_dp'
if(x> u1)stop 'In xyza_to_dig16; x> 1_0_dp'
if(y< u0)stop 'In xyza_to_dig16; y< 0_0_dp'
if(y> u1)stop 'In xyza_to_dig16; y> 1_0_dp'
if(z< u0)stop 'In xyza_to_dig16; z< 0_0_dp'
if(z> u1)stop 'In xyza_to_dig16; z> 1_0_dp'
if(a< u0)stop 'In xyza_to_dig16; a< 0_0_dp'
if(a> u1)stop 'In xyza_to_dig16; a> 1_0_dp'
dig16=0
do igen=1,size(dig16)
x=x*u2; y=y*u2; z=z*u2; a=a*u2
if(x>=u1)then; dig16(igen)=dig16(igen)+1; x=x-u1; endif
if(y>=u1)then; dig16(igen)=dig16(igen)+2; y=y-u1; endif
if(z>=u1)then; dig16(igen)=dig16(igen)+4; z=z-u1; endif
if(a>=u1)then; dig16(igen)=dig16(igen)+8; a=a-u1; endif
enddo
end subroutine xyza_to_dig16_d
!=============================================================================
subroutine gn_to_ea_s(x1,y1, x2,y2)! [gn_to_ea]
!=============================================================================
! Gnomonic to equal-area cubic
!=============================================================================
implicit none
real(sp),intent(IN ):: x1,y1
real(sp),intent(OUT):: x2,y2
!-----------------------------------------------------------------------------
integer(i_kind) :: iquad
real(sp) :: x,q,xx,xxp,rxxp,p
!=============================================================================
iquad=1
if(y1 > x1)iquad=iquad+1
if(x1 < -y1)iquad=iquad+2
if(x1== 0.0_sp .and. y1==0.0_sp )then; x2=0.0_sp; y2=0.0_sp; return; endif
select case(iquad)
case(1,4); x=abs(x1); q=y1/x
case(2,3); x=abs(y1); q=x1/x
end select
xx=x*x; xxp=xx+1.0_sp; rxxp=sqrt(xxp)
q=q*sqrt((xx+xxp)/(xxp+q*q*xx))
p=sqrt(asin(xx/xxp)/pio6)
select case(iquad)
case(1); x2= p; y2=p*q
case(2); y2= p; x2=p*q
case(3); y2=-p; x2=p*q
case(4); x2=-p; y2=p*q
end select
end subroutine gn_to_ea_s
!=============================================================================
subroutine gn_to_ea_d(x1,y1, x2,y2)! [gn_to_ea]
!=============================================================================
! Gnomonic to equal-area cubic
!============================================================================
implicit none
real(dp),intent(IN ):: x1,y1
real(dp),intent(OUT):: x2,y2
!-----------------------------------------------------------------------------
integer(i_kind) :: iquad
real(dp):: x,q,xx,xxp,rxxp,p
!=============================================================================
iquad=1
if(y1 > x1)iquad=iquad+1
if(x1 < -y1)iquad=iquad+2
if(x1==0.0_dp .and. y1==0.0_dp)then; x2=0.0_dp; y2=0.0_dp; return; endif
select case(iquad)
case(1,4); x=abs(x1); q=y1/x
case(2,3); x=abs(y1); q=x1/x
end select
xx=x*x; xxp=xx+1; rxxp=sqrt(xxp)
q=q*sqrt((xx+xxp)/(xxp+q*q*xx))
p=sqrt(asin(xx/xxp)/pio6)
select case(iquad)
case(1); x2= p; y2=p*q
case(2); y2= p; x2=p*q
case(3); y2=-p; x2=p*q
case(4); x2=-p; y2=p*q
end select
end subroutine gn_to_ea_d
!=============================================================================
subroutine ea_to_gn_s(x2,y2, x1,y1)! [ea_to_gn]
!=============================================================================
! Equal-area cubic to gnomonic
!============================================================================
implicit none
real(sp),intent(IN ):: x2,y2
real(sp),intent(OUT):: x1,y1
!-----------------------------------------------------------------------------
integer(i_kind) :: iquad
real(sp):: x,q,xx,xxp,p,pp,s
!=============================================================================
iquad=1
if(y2 > x2)iquad=iquad+1
if(x2 < -y2)iquad=iquad+2
if(x2==0.0_sp .and. y2==0.0_sp)then; x1=0.0_sp; y1=0.0_sp; return; endif
select case(iquad)
case(1); p= x2; q=y2/p
case(2); p= y2; q=x2/p
case(3); p=-y2; q=x2/p
case(4); p=-x2; q=y2/p
end select
pp=p*p
s=sin(pio6*pp)
xx=max(0._sp,s/(1.0_sp-s))
xxp=xx+1.0_sp
x=sqrt(xx)
q=q * sqrt( xxp/(xxp+xx*(1.0_sp-q*q)) )
select case(iquad)
case(1); x1= x; y1=x*q
case(2); y1= x; x1=x*q
case(3); y1=-x; x1=x*q
case(4); x1=-x; y1=x*q
end select
end subroutine ea_to_gn_s
!=============================================================================
subroutine ea_to_gn_d(x2,y2, x1,y1)! [ea_to_gn]
!=============================================================================
! equal-area cubic to gnomonic
!=============================================================================
implicit none
real(dp),intent(IN ):: x2,y2
real(dp),intent(OUT):: x1,y1
!-----------------------------------------------------------------------------
integer(i_kind) :: iquad
real(dp):: x,q,xx,xxp,p,pp,s
!=============================================================================
iquad=1
if(y2 > x2)iquad=iquad+1
if(x2 < -y2)iquad=iquad+2
if(x2==0.0_dp .and. y2==0.0_dp)then; x1=0.0_dp; y1=0.0_dp; return; endif
select case(iquad)
case(1); p= x2; q=y2/p
case(2); p= y2; q=x2/p
case(3); p=-y2; q=x2/p
case(4); p=-x2; q=y2/p
end select
pp=p*p
s=sin(pio6*pp)
xx=max(u0,s/(1.0_dp-s))
xxp=xx+u1
x=sqrt(xx)
q=q * sqrt( xxp/(xxp+xx*(1.0_dp-q*q)) )
select case(iquad)
case(1); x1= x; y1=x*q
case(2); y1= x; x1=x*q
case(3); y1=-x; x1=x*q
case(4); x1=-x; y1=x*q
end select
end subroutine ea_to_gn_d
!=============================================================================
!=============================================================================
! Routines whose usage is deprecated (better versions are now available)
!=============================================================================
subroutine hil4_to_rz_s(lgen,ngen,hil4,r)! [hil4_to_rz]
!=============================================================================
! Deprecated; replace by hil4_to_r with r also predefined
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(IN ):: hil4
real(sp), intent(OUT):: r
!-----------------------------------------------------------------------------
real(sp),parameter:: o4=1._sp/4_sp
real(sp) :: p
integer(i_kind) :: i
!=============================================================================
r=0.0_sp; if(lgen==0)r=hil4(lgen)
p=o4
do i=1,ngen
r=r+p*hil4(i)
p=p*o4
enddo
end subroutine hil4_to_rz_s
!=============================================================================
subroutine hil4_to_rz_d(lgen,ngen,hil4,r)! [hil4_to_rz]
!=============================================================================
! Deprecated; replace by hil4_to_r with r also predefined
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(IN ):: hil4
real(dp), intent(OUT):: r
!-----------------------------------------------------------------------------
real(dp):: p
integer(i_kind) :: i
!=============================================================================
r=0.0_dp; if(lgen==0)r=hil4(lgen)
p=o4
do i=1,ngen
r=r+p*hil4(i)
p=p*o4
enddo
end subroutine hil4_to_rz_d
!=============================================================================
subroutine hil8_to_rz_d(lgen,ngen,hil8,r)! [hil8_to_rz]
!=============================================================================
! Deprecated; replace by hil8_to_r with r also predefined
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(IN ):: hil8
real(dp), intent(OUT):: r
!-----------------------------------------------------------------------------
real(dp):: p
integer(i_kind) :: i
!=============================================================================
r=0.0_dp; if(lgen==0)r=hil8(lgen)
p=o8
do i=1,ngen
r=r+p*hil8(i)
p=p*o8
enddo
end subroutine hil8_to_rz_d
!=============================================================================
subroutine hil16_to_rz_d(lgen,ngen,hil16,r)! [hil16_to_rz]
!=============================================================================
! Deprecated; replace by hil16_to_r with r also predefined
implicit none
integer(i_kind), intent(IN ):: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(IN ):: hil16
real(dp), intent(OUT):: r
!-----------------------------------------------------------------------------
real(dp):: p
integer(i_kind) :: i
!=============================================================================
r=0.0_dp; if(lgen==0)r=hil16(lgen)
p=o16
do i=1,ngen
r=r+p*hil16(i)
p=p*o16
enddo
end subroutine hil16_to_rz_d
!=============================================================================
subroutine xy_to_hil4z_s(ngen,x,y,hil4)! [xy_to_hil4]
!=============================================================================
! DEPRECATED (since ngen is a redundant variable)
! Convert an (x,y)-representation of a point in the proper interior of the
! unit square to an ngen-digit base-4 representation of the parameter of
! a space-filling Hilbert curve.
!=============================================================================
implicit none
integer(i_kind), intent(IN ):: ngen
real(sp), intent(IN ):: x,y
integer(i_kind),dimension(ngen),intent(OUT):: hil4
!-----------------------------------------------------------------------------
real(sp):: xr,yr
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y
call xy_to_dig4(xr,yr,hil4)
presor=0
call dig4_to_hil4(presor,hil4)
end subroutine xy_to_hil4z_s
!=============================================================================
subroutine xy_to_hil4z_d(ngen,x,y,hil4)! [xy_to_hil4]
!=============================================================================
! DEPRECATED (since ngen is a redundant variable)
implicit none
integer(i_kind), intent(IN ):: ngen
real(dp), intent(IN ):: x,y
integer(i_kind),dimension(ngen),intent(OUT):: hil4
!-----------------------------------------------------------------------------
real(dp):: xr,yr
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y
call xy_to_dig4(xr,yr,hil4)
presor=0
call dig4_to_hil4(presor,hil4)
end subroutine xy_to_hil4z_d
!=============================================================================
subroutine xyz_to_hil8z_d(ngen,x,y,z,hil8)! [xyz_to_hil8]
!=============================================================================
! DEPRECATED (since ngen is a redundant variable)
implicit none
integer(i_kind), intent(IN ):: ngen
real(dp), intent(IN ):: x,y,z
integer(i_kind),dimension(ngen),intent(OUT):: hil8
!-----------------------------------------------------------------------------
real(dp):: xr,yr,zr
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y; zr=z
call xyz_to_dig8(xr,yr,zr,hil8)
presor=0
call dig8_to_hil8(presor,hil8)
end subroutine xyz_to_hil8z_d
!=============================================================================
subroutine xyza_to_hil16z_d(ngen,x,y,z,a,hil16)! [xyza_to_hil16]
!=============================================================================
! DEPRECATED (since ngen is a redundant variable)
implicit none
integer(i_kind), intent(IN ):: ngen
real(dp), intent(IN ):: x,y,z,a
integer(i_kind),dimension(ngen),intent(OUT):: hil16
!-----------------------------------------------------------------------------
real(dp):: xr,yr,zr,ar
integer(i_kind) :: presor
!=============================================================================
xr=x; yr=y; zr=z; ar=a
call xyza_to_dig16(xr,yr,zr,ar,hil16)
presor=0
call dig16_to_hil16(presor,hil16)
end subroutine xyza_to_hil16z_d
end module phil0