OctreeBitPlus.f90

!-------------------------------------------

! NEW OCTREEBIT LIBRARY (WRITTEN 26-27/02/2004) BY JEAN BRAUN

! This version addresses problems that arose in the computation
! of the icon array for very large octrees. This new version is much
! more efficient and uses the same amount of memory (sometimes less).
! Many of the querry routines have been greatly improved by storing
! additional information in the octree structure. Data storage and
! navigation within the octree are based on the methods described by
! Frisken and Perry from Mitsubishi Electric Research Lab.
! The computation of the icon array is based on a three stage method:
! (1) Compute icon and a nodal position array disregarding the nodal connectivity
! (2) Sort the nodal arrays by ordering along x, y and z coordinates
! (3) Remove the spurious nodes and modify the icon array accordingly
! Because a fast sorting method is used (O(n log log n)), this approach
! is rather efficient.


!=======================!
!=====[OCTREE_INIT]=====!
!=======================!

subroutine octree_init (octree,noctree)

! This routine must be called before any other routine when an octree is created
! the basic structure of the octree is 
! octree(1)=maximum level (unit cube is level 0)
! octree(2)=number of leaves
! octree(3)=total length of octree
! for each cube in the octree (at location loc)
! octree(loc)=level
! octree(loc+1)=address of parent
! octree(loc+2 to loc+9)=address of children (if negative the child is a leaf and the
!         value is the leaf number in the sequence of leaves)

implicit none

integer noctree,octree(noctree)

integer levelmax,nleaves,length,loc,k

levelmax=1
nleaves=8
length=13

octree(1)=levelmax
octree(2)=nleaves
octree(3)=length

loc=4

octree(loc)=1
octree(loc+1)=0
  do k=1,8
  octree(loc+1+k)=-k
  enddo

return
end

!====================================!
!=====[FIND_INTEGER_COORDINATES]=====!
!====================================!

subroutine find_integer_coordinates (x,y,z,ix,iy,iz,levelmax)

! returns the integer coordinates of point (x,y,z) in (ix,iy,iz)
! the integer coordinates are determined by the maximum level in the octree
! levelmax

! the integer coordinate is comprise between 1 and 2**levelmax
! and corresponds to the location of the leaf containing the
! point of given coordinates

implicit none

double precision x,y,z
integer ix,iy,iz,levelmax

double precision powermax

powermax=2.d0**levelmax
ix=int(x*powermax)
iy=int(y*powermax)
iz=int(z*powermax)

return
end

!=================================!
!=====[FIND_REAL_COORDINATES]=====!
!=================================!

subroutine find_real_coordinates (ix,iy,iz,x,y,z,levelmax)

! returns the real coordinates of a point of integer coordinates (ix,iy,iz)
! the real coordinates are determined by the maximum level in the octree

implicit none

double precision x,y,z
integer ix,iy,iz,levelmax

double precision powermax

powermax=2**levelmax
x=dfloat(ix)/powermax
y=dfloat(iy)/powermax
z=dfloat(iz)/powermax

return
end

!========================================!
!=====[OCTREE_CREATE_FROM_PARTICLE]=====!
!========================================!

subroutine octree_create_from_particle (octree,noctree,x,y,z,np,level)

! updates the octree by creating a leaf at point (x,y,z)
! of level level
! if the leaf (or a cube of smaller level) exists, the routine has no effect on the octree
! note that x,y,z must belong to [0,1[

implicit none

integer noctree,octree(noctree),np
double precision x,y,z
integer level
double precision xp,yp,zp

integer levelin,ix,iy,iz,levelmax,loc,ip

levelmax=octree(1)

xp=x
yp=y
zp=z
levelin=0
loc=4
  if (xp.eq.1.d0) xp=1.d0-1.d-20
  if (yp.eq.1.d0) yp=1.d0-1.d-20
  if (zp.eq.1.d0) zp=1.d0-1.d-20
  if (xp.eq.0.d0) xp=1.d-20
  if (yp.eq.0.d0) yp=1.d-20
  if (zp.eq.0.d0) zp=1.d-20
  if (xp*(xp-1.d0).lt.0.d0 .and. yp*(yp-1.d0).lt.0.d0 .and. zp*(zp-1.d0).lt.0.d0) &
    call update (octree,noctree,xp,yp,zp,level,levelin,loc)

return
end

!========================================!
!=====[OCTREE_CREATE_FROM_PARTICLES]=====!
!========================================!

subroutine octree_create_from_particles (octree,noctree,x,y,z,np,level)

! updates the octree by creating a leaf at points (x(1:np),y(1:np),z(1:np))
! of level level(1:np)
! if the leaf (or a cube of smaller level) exists, the routine has no effect on the octree
! note that x,y,z must belong to [0,1[

implicit none

integer noctree,octree(noctree),np
double precision x(np),y(np),z(np)
integer level(np)
double precision xp,yp,zp

integer levelin,ix,iy,iz,levelmax,loc,ip

levelmax=octree(1)

do ip=1,np
xp=x(ip)
yp=y(ip)
zp=z(ip)
levelin=0
loc=4
  if (xp.eq.1.d0) xp=1.d0-1.d-20
  if (yp.eq.1.d0) yp=1.d0-1.d-20
  if (zp.eq.1.d0) zp=1.d0-1.d-20
  if (xp.eq.0.d0) xp=1.d-20
  if (yp.eq.0.d0) yp=1.d-20
  if (zp.eq.0.d0) zp=1.d-20
  if (xp*(xp-1.d0).lt.0.d0 .and. yp*(yp-1.d0).lt.0.d0 .and. zp*(zp-1.d0).lt.0.d0) &
    call update (octree,noctree,xp,yp,zp,level(ip),levelin,loc)
enddo

return
end

!==================!
!=====[UPDATE]=====!
!==================!

recursive subroutine update (octree,noctree,x,y,z,level,levelin,loc)

! DO NOT USE
! internal routine called by OctreeUpdate
! this routine divides a cube of the octree in 8
! if the requested level is greater than0
! and if the requested level is stricly greater than levelin
! loc is the location in the octree of the current cube (to be divided)

implicit none

integer noctree,octree(noctree),ix,iy,iz,levelin,level,loc
double precision x,y,z
integer ibits_jean
external ibits_jean

integer levelmax,ibitx,ibity,ibitz,ichild,iaddress,newaddress

if (level.eq.0) return

if (levelin+1.eq.level) return

levelmax=octree(1)
call find_integer_coordinates (x,y,z,ix,iy,iz,levelmax)
ibitx=ibits_jean(ix,levelmax-levelin-1)
ibity=ibits_jean(iy,levelmax-levelin-1)
ibitz=ibits_jean(iz,levelmax-levelin-1)

ichild=loc+2+ibitz*4+ibity*2+ibitx

if (ichild.gt.octree(3)) then
print*,ichild,' requested ',octree(3),' available'
stop 'octree needs to grow...'
endif

iaddress=octree(ichild)

  if (iaddress.lt.0) then
!adds a level
  newaddress=octree(3)+1
! updates current level
  octree(ichild)=newaddress
! updates octree length
  octree(3)=octree(3)+10
! updates number of leaves
  octree(2)=octree(2)+7
! updates maximum level
  if (levelin+2.gt.octree(1)) octree(1)=levelin+2
! creates new division
  octree(newaddress)=levelin+2
  octree(newaddress+1)=ichild
  octree(newaddress+2)=iaddress
  octree(newaddress+3)=-(octree(2)-6)
  octree(newaddress+4)=-(octree(2)-5)
  octree(newaddress+5)=-(octree(2)-4)
  octree(newaddress+6)=-(octree(2)-3)
  octree(newaddress+7)=-(octree(2)-2)
  octree(newaddress+8)=-(octree(2)-1)
  octree(newaddress+9)=-(octree(2))
  call update (octree,noctree,x,y,z,level,levelin+1,newaddress)
  else
  call update (octree,noctree,x,y,z,level,levelin+1,iaddress)
  endif

return
end


!=================!
!=====[IBITS]=====!
!=================!

integer function ibits_jean(i,ipos)

!internal function, do not modify

integer j
j=i/2**ipos
ibits_jean=j-(j/2)*2
return
end

!============================!
!=====[OCTREE_FIND_LEAF]=====!
!============================!

subroutine octree_find_leaf (octree,noctree,x,y,z,leaf,level,loc,x0,y0,z0,dxyz)

! given an octree of size noctree
! this routine returns the leaf number in which a point (x,y,z) resides
! the level of the leaf (0 is unit cube)
! the location in the octree of the part describing the parent of the leaf (loc)
! the centroid of the leaf (x0,y0,z0) and its size (dxyz)

implicit none

integer noctree,octree(noctree)
double precision x,y,z,x0,y0,z0,dxyz
integer leaf,level,loc

integer ix,iy,iz,levelin,locin,levelmax

levelmax=octree(1)
call find_integer_coordinates (x,y,z,ix,iy,iz,levelmax)

levelin=0
locin=4
leaf=0
call find_leaf (octree,noctree,ix,iy,iz,levelin,locin,leaf,level,loc)

call find_integer_coordinates (x,y,z,ix,iy,iz,level)
call find_real_coordinates (ix,iy,iz,x0,y0,z0,level)
dxyz=1.d0/2.d0**level

return
end

!=====================!
!=====[FIND_LEAF]=====!
!=====================!

recursive subroutine find_leaf (octree,noctree,ix,iy,iz,levelin,locin,leaf,level,loc)

! DO NOT USE
! internal routine used by octree_find_leaf

implicit none

integer noctree,octree(noctree),ix,iy,iz,levelin,leaf,level,loc,locin
integer levelmax,ibitx,ibity,ibitz,ichild,iaddress
integer ibits_jean
external ibits_jean

levelmax=octree(1)
!ibitx=ibits(ix,levelmax-levelin-1,1)
!ibity=ibits(iy,levelmax-levelin-1,1)
!ibitz=ibits(iz,levelmax-levelin-1,1)
ibitx=ibits_jean(ix,levelmax-levelin-1)
ibity=ibits_jean(iy,levelmax-levelin-1)
ibitz=ibits_jean(iz,levelmax-levelin-1)

ichild=locin+2+ibitz*4+ibity*2+ibitx

if (ichild.gt.octree(3)) stop 'octree needs to grow...'
iaddress=octree(ichild)

  if (iaddress.lt.0) then
  leaf=-iaddress
  level=levelin+1
  loc=locin
  return
  else
  call find_leaf (octree,noctree,ix,iy,iz,levelin+1,iaddress,leaf,level,loc)
  endif

return
end

!===============================!
!=====[OCTREETHROUGHLEAVES]=====!
!===============================!

subroutine OctreeThroughLeaves (octree,noctree)

! DO NOT USE
! this is a general routine that simply goes through the leaves
! this routine is used as a template to build other routines.
! it should not be used

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz

loc=4
ix=0
iy=0
iz=0
call throughleaves (octree,noctree,loc,ix,iy,iz)

return
end

!=========================!
!=====[THROUGHLEAVES]=====!
!=========================!

recursive subroutine throughleaves (octree,noctree,loc,ix,iy,iz)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
    else
    call throughleaves (octree,noctree,locn,ixn,iyn,izn)
    endif
  enddo
  enddo
  enddo

return
end

!===========================!
!=====[OCTREE_SMOOTHEN]=====!
!===========================!

subroutine octree_smoothen (octree,noctree)

! as it names indicates this routine smoothens the octree
! ie it ensures that no two adjacent leaves are more than one level apart

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz,nleaves,nleaves0
integer ioctree_number_of_elements

nleaves=ioctree_number_of_elements (octree,noctree)

do while (nleaves.ne.nleaves0)
nleaves0=nleaves
loc=4
ix=0
iy=0
iz=0
call smooth (octree,noctree,loc,ix,iy,iz)
nleaves=ioctree_number_of_elements (octree,noctree)
enddo

return
end

!==================!
!=====[SMOOTH]=====!
!==================!

recursive subroutine smooth (octree,noctree,loc,ix,iy,iz)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer k,idx,idy,idz,ip
integer ipmax,iddx,iddy,iddz,ixp,iyp,izp,levelin,locp
double precision xp,yp,zp

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! the binary coordinate of their bottom corner (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
! the address of their parent (loc)

! first check that `right' neighbours are not down by more than 1 level
    ipmax=2**levelmax
    ip=2**level
    izp=izn+ipower
    iyp=iyn+ipower
    ixp=ixn+ipower
      if (izp.lt.ipmax) then
      zp=dfloat(izp)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (iyp.lt.ipmax) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyp)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (ixp.lt.ipmax) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixp)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif

! second check if the 'left' neighbours are not down by more than one level
    izp=izn-1
    iyp=iyn-1
    ixp=ixn-1
      if (izp.ge.0) then
      zp=dfloat(izp)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (iyp.ge.0) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyp)/ipmax
      xp=dfloat(ixn)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif
      if (ixp.ge.0) then
      zp=dfloat(izn)/ipmax
      yp=dfloat(iyn)/ipmax
      xp=dfloat(ixp)/ipmax
      levelin=0
      locp=4
      call update (octree,noctree,xp,yp,zp,level-1,levelin,locp)
      endif

    else
    call smooth(octree,noctree,locn,ixn,iyn,izn)
    endif
  enddo
  enddo
  enddo

return
end

!=================================!
!=====[OCTREE_SUPER_SMOOTHEN]=====!
!=================================!

subroutine octree_super_smoothen (octree,noctree)

! as it names indicates this routine further smoothens the octree
! ie it ensures that for any given leaf and its six closest neighbouring
! leaves, ther is no difference in level larger than 1

! Note that contrary to octree_smoothen, octree_super_smoothen can
! be used iteratively to increase the smoothness of the octree

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz,nleaves,nleaves0
integer ioctree_number_of_elements

nleaves=ioctree_number_of_elements (octree,noctree)

do while (nleaves.ne.nleaves0)
nleaves0=nleaves
loc=4
ix=0
iy=0
iz=0
call super_smooth (octree,noctree,loc,ix,iy,iz,nleaves)
nleaves=ioctree_number_of_elements (octree,noctree)
enddo

call octree_smoothen (octree,noctree)

return
end

!========================!
!=====[SUPER_SMOOTH]=====!
!========================!

recursive subroutine super_smooth (octree,noctree,loc,ix,iy,iz,nleaves0)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer k,idx,idy,idz,ip
integer ipmax,iddx,iddy,iddz,ixp,iyp,izp,levelin,locp
double precision xp,yp,zp
integer levelout(6),levelset,leaf,locin,nleaves0

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! the binary coordinate of their bottom corner (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
! the address of their parent (loc)

! first find `right' neighbours levels
    ipmax=2**levelmax
    ip=2**level
    izp=izn+ipower
    iyp=iyn+ipower
    ixp=ixn+ipower
    levelout=level
      if (izp.lt.ipmax) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyn,izp,levelin,locin,leaf,levelout(1),locp)
      if (leaf.gt.nleaves0) levelout(1)=level
      endif
      if (iyp.lt.ipmax) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyp,izn,levelin,locin,leaf,levelout(2),locp)
      if (leaf.gt.nleaves0) levelout(2)=level
      endif
      if (ixp.lt.ipmax) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixp,iyn,izn,levelin,locin,leaf,levelout(3),locp)
      if (leaf.gt.nleaves0) levelout(3)=level
      endif

! second find 'left' neighbours levels
    izp=izn-1
    iyp=iyn-1
    ixp=ixn-1
      if (izp.ge.0) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyn,izp,levelin,locin,leaf,levelout(4),locp)
      if (leaf.gt.nleaves0) levelout(4)=level
      endif
      if (iyp.ge.0) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixn,iyp,izn,levelin,locin,leaf,levelout(5),locp)
      if (leaf.gt.nleaves0) levelout(5)=level
      endif
      if (ixp.ge.0) then
      levelin=0
      locin=4
      leaf=0
      call find_leaf (octree,noctree,ixp,iyn,izn,levelin,locin,leaf,levelout(6),locp)
      if (leaf.gt.nleaves0) levelout(6)=level
      endif

! if the difference in level between the leaf and any of its neighbours is greater
! than 1, all leaves are set to the maximum level minus one at least

      if (maxval(levelout)-minval(levelout).gt.1) then
      levelset=maxval(levelout)-1
! now check that `right' neighbours are at least at this level
      izp=izn+ipower
      iyp=iyn+ipower
      ixp=ixn+ipower
        if (izp.lt.ipmax) then
        zp=dfloat(izp)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (iyp.lt.ipmax) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyp)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (ixp.lt.ipmax) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixp)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
  ! now check that 'left' neighbours are at least at this level
      izp=izn-1
      iyp=iyn-1
      ixp=ixn-1
        if (izp.ge.0) then
        zp=dfloat(izp)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (iyp.ge.0) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyp)/ipmax
        xp=dfloat(ixn)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
        if (ixp.ge.0) then
        zp=dfloat(izn)/ipmax
        yp=dfloat(iyn)/ipmax
        xp=dfloat(ixp)/ipmax
        levelin=0
        locp=4
        call update (octree,noctree,xp,yp,zp,levelset,levelin,locp)
        endif
  
      endif

    else
    call super_smooth(octree,noctree,locn,ixn,iyn,izn,nleaves0)
    endif
  enddo
  enddo
  enddo

return
end

!===============================!
!=====[OCTREE_SHOW_AS_MESH]=====!
!===============================!

subroutine octree_show_as_mesh (octree,noctree)

! routine to show the octree as a mesh in vrml
! produces a VRML file called mesh.wrl
! this is mostly used for debugging purposes

implicit none

integer noctree,octree(noctree)

integer loc,ix,iy,iz

open (7,file='mesh.wrl',status='unknown')
write (7,'(a)') '#VRML V2.0 utf8'
write (7,'(a)') 'Transform { children ['
write (7,'(a)') 'NavigationInfo { '
write (7,'(a)') 'type ["EXAMINE"]'
write (7,'(a)') 'headlight FALSE}'
write (7,'(a)') 'Background{groundColor 1 1 1 skyColor 1 1 1}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  .8 1 .5}'
write (7,'(a)') 'DirectionalLight {ambientIntensity  0.2'
write (7,'(a)') '                  color      1 1 1'
write (7,'(a)') '                  direction  -.8 -1 -.5}'
write (7,'(a)') 'Transform { children Viewpoint {'
write (7,'(a)') '      description "Starting"'
write (7,'(a)') '      fieldOfView 1'
write (7,'(a,3f12.8,a)') '      position    ',-.41885125637054443, -.8311104774475098, 1.5406757593154907
write (7,'(a,4f12.8,a)') '      orientation ', .7352051138877869, -.10698327422142029, -.669348955154419, 1.3260198831558228,'}}'
write (7,'(a,a)') 'DEF Node0 Shape{geometry Sphere{radius 0.0075', &
                ' }appearance Appearance{material Material{diffuseColor 1 0 0}}}'
write (7,'(a,3f10.6,a)') 'Transform{translation',0.,0.,0., &
                         ' children [USE Node0]}'
loc=4
ix=0
iy=0
iz=0
call show (octree,noctree,loc,ix,iy,iz)
write (7,'(a)') ']}'
close (7)

return
end

!================!
!=====[SHOW]=====!
!================!

recursive subroutine show (octree,noctree,loc,ix,iy,iz)

! DO NOT USE
! internal routine
! to show the octree as a mesh

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k
double precision x1,y1,z1,x2,y2,z2,dxyz

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
    call find_real_coordinates (ixn,iyn,izn,x1,y1,z1,levelmax)
    dxyz=1.d0/(2.d0**level)
    x2=x1+dxyz
    y2=y1+dxyz
    z2=z1+dxyz
    write (7,'(a,24f10.3,a,a,a)') &
         'Shape { geometry IndexedLineSet { coord Coordinate { point [', &
          x1,y1,z1,x2,y1,z1,x2,y2,z1,x1,y2,z1,x1,y1,z2,x2,y1,z2,x2,y2,z2,x1,y2,z2, &
          ']} coordIndex [0 1 2 3 0 -1 4 5 6 7 4 -1 0 4 -1 1 5 -1 2 6 -1 3 7 -1', &
          ']}appearance Appearance { material Material { emissiveColor 0 0 0}}}'
    else
    call show (octree,noctree,locn,ixn,iyn,izn)
    endif
  enddo
  enddo
  enddo

return
end

!=====================================!
!=====[OCTREE_FIND_ELEMENT_LEVEL]=====!
!=====================================!

subroutine octree_find_element_level (octree,noctree,levs,nleaves)

! routine to return the level of each leaf
! the result is returned in the array levs of dimension nleaves

! the function ioctree_number_of_elements should be called first to
! find nleaves and dimension levs accordingly

implicit none

integer noctree,octree(noctree),nleaves,levs(nleaves)

integer loc,ix,iy,iz

loc=4
ix=0
iy=0
iz=0
call levels (octree,noctree,loc,ix,iy,iz,levs,nleaves)

return
end

!==================!
!=====[LEVELS]=====!
!==================!

recursive subroutine levels (octree,noctree,loc,ix,iy,iz,levs,nleaves)

! DO NOT USE
! internal routine
! on entry we have the address of the current cube
! and the binary coordinates of its bottom corner

implicit none

integer noctree,octree(noctree),loc,ix,iy,iz
integer nleaves,levs(nleaves)

integer level,levelmax,ipower,ixn,iyn,izn,locn
integer idx,idy,idz,k

level=octree(loc)
levelmax=octree(1)

  do idz=0,1
  do idy=0,1
  do idx=0,1
  k=idx+idy*2+idz*4
  locn=octree(loc+2+k)
  ipower=2**(levelmax-level)
  ixn=ix+idx*ipower
  iyn=iy+idy*ipower
  izn=iz+idz*ipower
    if (locn.lt.0) then
! here i am going through the leaves one by one and i know
! their coordinates (ixn,iyn,izn)
! their address (loc+2+k)
! their leaf number (-locn)
! their level (level)
    if (-locn.gt.nleaves) stop 'array level needs to grow'
    levs(-locn)=level
    else
    call levels (octree,noctree,locn,ixn,iyn,izn,levs,nleaves)
    endif
  enddo
  enddo
  enddo

return
end

!======================================!
!=====[IOCTREE_NUMBER_OF_ELEMENTS]=====!
!======================================!

integer function ioctree_number_of_elements (octree,noctree)

! function that returns (in ioctree_number_of_elements) the number of leaves in the octree

implicit none

integer noctree,octree(noctree)

ioctree_number_of_elements=octree(2)

return
end

!=================================!
!=====[IOCTREE_MAXIMUM_LEVEL]=====!
!=================================!

integer function ioctree_maximum_level (octree,noctree)

! function that returns (in ioctree_maximum_level) the maximum level in the octree

implicit none

integer noctree,octree(noctree)

ioctree_maximum_level=octree(1)

return
end

!========================!
!=====[IOCTREE_SIZE]=====!
!========================!

integer function ioctree_size (octree,noctree)

! function that returns (in ioctree_size) the size of the octree

implicit none