! R. J. Purser NOAA/NCEP/EMC 2006
!=============================================================================
module phil
!=============================================================================
!$$$ module documentation block
! . . . .
! module: phil
! prgmmr: purser
!
! abstract: Module procedures pertaining to Hilbert curve transformations and
! representations.
!
! program history log:
! 2009-09-22 lueken - added module doc block
!
! subroutines included:
! sub bsort_s
! sub mergeab_s
! sub hil_to_r_s
! sub r_to_kil_s
! sub xy_to_hil_s
! sub hil_to_xy_s
! sub xc_to_hil_s
! sub hil_to_xc_s
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
use kinds, only: r_kind,i_kind
use constants, only: zero,one,two,four
implicit none
private
public bsort,mergeab,r_to_hil,hil_to_r,xy_to_hil,hil_to_xy,xc_to_hil,hil_to_xc
interface bsort; module procedure bsort_s; end interface
interface mergeab; module procedure mergeab_s; end interface
interface hil_to_r; module procedure hil_to_r_s; end interface
interface r_to_hil; module procedure r_to_hil_s; end interface
interface xy_to_hil; module procedure xy_to_hil_s; end interface
interface hil_to_xy; module procedure hil_to_xy_s; end interface
interface xc_to_hil; module procedure xc_to_hil_s; end interface
interface hil_to_xc; module procedure hil_to_xc_s; end interface
contains
!=============================================================================
recursive subroutine bsort_s(n1,n2,v,next,first)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: bsort_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! n1,n2
! v
!
! output argument list:
! next
! first
!
! note:
! Recursively sort strings by recursive binary "divide and conquer" until
! the divided strings are of length L or less, where L is empirically chosen
! here in the parameter statement (efficiency is found not to be strongly
! sensitive to the L). The values, v, of items in the string are real.
! An integer representation of all the linked-list pointers is used.
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
implicit none
integer(i_kind), intent(IN ) :: n1,n2
real(r_kind), dimension(n1:n2),intent(IN ) :: v
integer(i_kind), dimension(n1:n2),intent( OUT) :: next
integer(i_kind), intent( OUT) :: first
!-----------------------------------------------------------------------------
integer(i_kind),parameter :: L=6
integer(i_kind) :: n,na1,na2,nb1,nb2,i,j,k, &
maxk,maxl,left,firsta,firstb
real(r_kind) :: maxv
!=============================================================================
n=n2+1-n1
if(n<=L)then
! Sort the small number of items by an order (n*n) algorithm:
do i=n1,n2-1
next(i)=i+1
enddo
next(n2)=0
first=n1
do i=n-1,1,-1
k=first
left=0
maxv=v(k)
maxl=left
maxk=k
do j=1,i
left=k
k=next(k)
if(v(k)>=maxv)then
maxv=v(k)
maxl=left
maxk=k
endif
enddo
if(k/=maxk)then
if(maxl==0)then
first=next(maxk)
else
next(maxl)=next(maxk)
endif
next(maxk)=next(k)
next(k)=maxk
endif
enddo
write(11,*) k,next(k)
else
na1=n1; na2=(n1+n2)/2; nb1=na2+1; nb2=n2
call bsort(na1,na2,v(na1:na2),next(na1:na2),firsta)
call bsort(nb1,nb2,v(nb1:nb2),next(nb1:nb2),firstb)
call mergeab(na1,nb2,firsta,firstb, &
v(na1:nb2),next(na1:nb2), first)
write(11,*) firsta,firstb,first
endif
end subroutine bsort_s
!=============================================================================
subroutine mergeab_s(na1,nb2,firsta,firstb, v,next,first)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: mergeab_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! na1,nb2
! firsta,firstb
! v
! next
!
! output argument list:
! next
! first
!
! note:
! Merge a pair (a and b) of individually pre-sorted strings of real values,
! connected as respective linked-lists, into a unified string with ALL the
! items returned in ascending order of values v.
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
implicit none
integer(i_kind), intent(IN ) :: na1,nb2, &
firsta,firstb
real(r_kind),dimension(na1:nb2), intent(IN ) :: v
integer(i_kind), dimension(na1:nb2),intent(INOUT) :: next
integer(i_kind), intent( OUT) :: first
!-----------------------------------------------------------------------------
integer(i_kind),parameter :: hugeint= 1000000
integer(i_kind) :: idum,ia,ib,ic
!=============================================================================
ia=firsta
ib=firstb
if(ia==0)then
first=firstb
return
endif
if(ib==0)then
first=firsta
return
endif
if(v(ia)<v(ib))then
first=ia
ic=ia
ia=next(ia)
else
first=ib
ic=ib
ib=next(ib)
endif
do idum=1,hugeint
if(ia==0)then
next(ic)=ib
return
endif
if(ib==0)then
next(ic)=ia
return
endif
if(v(ia)<v(ib))then
next(ic)=ia
ic=ia
ia=next(ia)
else
next(ic)=ib
ic=ib
ib=next(ib)
endif
enddo
stop 'in mergeab; hugeint too small'
end subroutine mergeab_s
!=============================================================================
subroutine hil_to_r_s(lgen,ngen,hil4,r)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: hil_to_r_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! lgen,ngen
! hil4
!
! output argument list:
! r
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind), intent(IN ) :: lgen,ngen
integer(i_kind),dimension(lgen:ngen),intent(IN ) :: hil4
real(r_kind), intent( OUT) :: r
!-----------------------------------------------------------------------------
real(r_kind),parameter :: f=one/four
real(r_kind) :: p
integer(i_kind) :: i
!=============================================================================
r=zero; if(lgen==0)r=hil4(lgen)
p=f
do i=1,ngen
r=r+p*hil4(i)
p=p*f
enddo
end subroutine hil_to_r_s
!=============================================================================
subroutine r_to_hil_s(lgen,ngen,r,hil4)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: r_to_hil_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! lgen,ngen
! r
!
! output argument list:
! hil4
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind), intent(IN ) :: lgen,ngen
real(r_kind), intent(IN ) :: r
integer(i_kind),dimension(lgen:ngen),intent( OUT) :: hil4
!-----------------------------------------------------------------------------
real(r_kind) :: t
integer(i_kind) :: i
!=============================================================================
t=r
if(lgen==0)then
hil4(0)=t
t=4*(t-hil4(0))
endif
do i=1,ngen
hil4(i)=t
t=4*(t-hil4(i))
enddo
end subroutine r_to_hil_s
!=============================================================================
subroutine xy_to_hil_s(ngen,x,y,hil4)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: xy_to_hil_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram documentation block
!
! input argument list:
! ngen
! x,y
!
! output argument list:
! hil4
!
! note:
! 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.
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
implicit none
integer(i_kind), intent(IN ) :: ngen
real(r_kind), intent(IN ) :: x,y
integer(i_kind),dimension(ngen),intent( OUT) :: hil4
!-----------------------------------------------------------------------------
real(r_kind) :: xr,yr
integer(i_kind),dimension(ngen) :: dig4
integer(i_kind),dimension(0:3, 0:7) :: hil4table, presortable
integer(i_kind) :: igen,presor
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/
!=============================================================================
if(x< zero)stop 'In xy_to_hil; x< 0'
if(x> one) stop 'In xy_to_hil; x> 1'
if(y< zero)stop 'In xy_to_hil; y< 0'
if(y> one) stop 'In xy_to_hil; y> 1'
! Begin by converting the coordinates (x,y) to a single base-4 number in [0,1)
xr=x
yr=y
do igen=1,ngen
xr=xr*2
yr=yr*2
dig4(igen)=0
if(xr>=one)then
dig4(igen)=dig4(igen)+1
xr=xr-one
endif
if(yr>=one)then
dig4(igen)=dig4(igen)+2
yr=yr-one
endif
enddo
! Initialized the present (coarsest scale) orientation code to 0:
presor=0
! 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,ngen
hil4(igen)=hil4table(dig4(igen),presor)
presor =presortable(dig4(igen),presor)
enddo
end subroutine xy_to_hil_s
!=============================================================================
subroutine hil_to_xy_s(ngen,hil4,x,y)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: hil_to_xy_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! ngen
! hil4
!
! output argument list:
! x,y
!
! note:
! Convert a base-4 representation of a Hilbert curve parameter to an (x,y)-
! location in the unit square
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
implicit none
integer(i_kind), intent(IN ) :: ngen
integer(i_kind),dimension(ngen),intent(IN ) :: hil4
real(r_kind), intent( OUT) :: x,y
!-----------------------------------------------------------------------------
real(r_kind) :: frac
integer(i_kind),dimension(0:3, 0:7) :: xtable,ytable,ptable
integer(i_kind) :: igen,presor,h
data xtable/0,0,1,1, 1,0,0,1, 1,1,0,0, 0,1,1,0, &
0,1,1,0, 0,0,1,1, 1,0,0,1, 1,1,0,0/
data ytable/0,1,1,0, 0,0,1,1, 1,0,0,1, 1,1,0,0, &
0,0,1,1, 1,0,0,1, 1,1,0,0, 0,1,1,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/
!=============================================================================
presor=0
x=zero
y=zero
frac=one
do igen=1,ngen
frac=frac/two
h=hil4(igen)
x=x+xtable(h,presor)*frac
y=y+ytable(h,presor)*frac
presor=ptable(h,presor)
enddo
end subroutine hil_to_xy_s
!=============================================================================
subroutine xc_to_hil_s(ngen,xc,hil4)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: xc_to_hil_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! ngen
! xc
!
! output argument list:
! hil4
!
!note:
! hil4(0) will contain the octant index, 0,..,7, according to the rules:
! xc(1)<0: hil4(0) = 0,1, 4,5; xc(1)>0: hil4(0) = 2,3, 6,7
! xc(2)<0: hil4(0) = 1,2, 4,7; xc(2)>0: hil4(0) = 0,3, 5,6
! xc(3)<0: hil4(0) = 4,5,6,7; xc(3)>0: hil4(0) = 0,1,2,3.
!
! Transform from earth-centered cartesians, xc, to an augmented base-4
! representation of the parameter of a surface-filling Hilbert curve. The
! "zeroth" digit of the hilbert curve number is the index [0,7] of the surface
! octant. The remaining ngen digits define the parameter of the portion of the
! Hilbert curve wedged into the denoted octant according to a mapping between
! the unit square and that spherical triangular octant.
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
implicit none
integer(i_kind), intent(IN ) :: ngen
real(r_kind),dimension(3), intent(IN ) :: xc
integer(i_kind), dimension(0:ngen),intent( OUT) :: hil4
!-----------------------------------------------------------------------------
real(r_kind) :: ax,ay,s,t,x,y,z
!=============================================================================
x=xc(1)
y=xc(2)
z=xc(3)
! Choose and perform the appropriate polar stereographic projection:
if(z>=zero)then
x=x/(one+z)
y=y/(one+z)
hil4(0)=0 ! Northern hemisphere
else
x=-x/(one-z)
y=-y/(one-z)
hil4(0)=4 ! Southern hemisphere
endif
ax=abs(x)
ay=abs(y)
! Get tangent, t, of a longitude rotated and reflected into [0,pi/4]:
t=zero
if(ay>ax)then
if(ay>zero)t=ax/ay
else
if(ax>zero)t=ay/ax
endif
s=sqrt(one+t*t) ! <- the corresponding secant
! Expand the unit disk to fill the unit square:
x=x*s
y=y*s
! Subdivide into quadrants of this square = octants of the original sphere:
if(y>zero)then
if(x>zero)then
hil4(0)=hil4(0)+3
else
x=x+one
endif
else
if(x>zero)then
hil4(0)=hil4(0)+2
t=x
x=one+y
y=t
else
hil4(0)=hil4(0)+1
t=x
x=-y
y=-t
endif
endif
call xy_to_hil(ngen,x,y,hil4(1:ngen))
end subroutine xc_to_hil_s
!=============================================================================
subroutine hil_to_xc_s(ngen,hil4,xc)
!=============================================================================
!$$$ subprogram documentation block
! . . . .
! subprogram: hil_to_xc_s
! prgmmr: purser
!
! abstract:
!
! program history log:
! 2009-09-22 lueken - added subprogram doc block
!
! input argument list:
! ngen
! hil4
!
! output argument list:
! xc
!
! note:
! Convert an augmented base-4 representation of the parameter of a
! space-filling Hilbert curve to its corresponding earth-centered cartesian
! respresentation on the unit sphere. The "zeroth" digit of the Hilbert curve
! parameter string, hil4, is the base-8 index of the surface octant, the
! remaining ngen digits are the base-4 parameter representations of the
! portion of the Hilbert curve wedged into this octant via a transformation
! mapping the unit square to the spherical triangle.
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
!=============================================================================
implicit none
integer(i_kind), intent(IN ) :: ngen
integer(i_kind), dimension(0:ngen),intent(IN ) :: hil4
real(r_kind),dimension(3), intent( OUT) :: xc
!-----------------------------------------------------------------------------
real(r_kind) :: ax,ay,s,t,x,y,z,rr
integer(i_kind) :: quad
!=============================================================================
call hil_to_xy(ngen,hil4(1:ngen),x,y)
quad=mod(hil4(0),4)
select case(quad)
case(0)
x=x-one
case(1)
t=x
x=-y
y=-t
case(2)
t=x
x=y
y=t-one
case(3)
end select
ax=abs(x)
ay=abs(y)
! Get tangent, t, of a longitude rotated and reflected into [0,pi/4]:
t=zero
if(ay>ax)then
if(ay>zero)t=ax/ay
else
if(ax>zero)t=ay/ax
endif
s=sqrt(one+t*t)! <- the corresponding secant
! Shrink the unit square to fit the unit disk:
x=x/s
y=y/s
rr=x*x+y*y
z=(one-rr)/(one+rr)
xc(1)=x*(one+z)
xc(2)=y*(one+z)
xc(3)=z
if(hil4(0)>=4)xc=-xc
end subroutine hil_to_xc_s
!=============================================================================
end module phil