!-----------------------------------------------------
!           VOXELIZE
!           ========
!    this is the main program, see also mkvoxvidz.sh
!    showvoxels.pov and vox2inc.awk
!-----------------------------------------------------
program voxelize
  use fraktals

  integer, parameter                          :: DIM = 100
  integer, dimension(:,:,:), allocatable      :: cube
  type(t_point3d), dimension(:), allocatable  :: points
  integer                                     :: errcode, foo, maxi
  integer                                     :: ix, iy, iz
  double precision, dimension(4)              :: coefs
  double precision                            :: dval

  ! XXX foo = (DIM*DIM*DIM) / (1024)
  ! XXX PRINT *, "memory request for cube (in Kwords) ", foo

  allocate (cube(DIM,DIM,DIM), stat=errcode)
  if (0 .NE. errcode) then
    STOP " : NO ENOUGH MEMORY FOR CUBE"
  endif

  nbr_points = 3500000
  allocate(points(nbr_points), stat=errcode)
  if (0 .NE. errcode) then
    STOP " : NO ENOUGH MEMORY FOR POINTS"
  endif

  coefs(1) =  2.24     ;     coefs(2) =  0.43
  coefs(3) = -0.65     ;     coefs(4) = -2.43
  call compute_pickover(points, coefs)

  call clear_cube(cube)

  ! and now, we loop over all the pre-computed
  ! points of the attractor
  ! 
  do foo=1, nbr_points
    call fcoor2icoor(points(foo)%x, ix)
    call fcoor2icoor(points(foo)%y, iy)
    call fcoor2icoor(points(foo)%z, iz)
    cube(ix,iy,iz) = cube(ix,iy,iz) + 1
  enddo

  dval = DBLE(MAXVAL(cube))

  write(0, *) "--- cube maximum = ", dval

  do foo=1, nbr_points
    call fcoor2icoor(points(foo)%x, ix)
    call fcoor2icoor(points(foo)%y, iy)
    call fcoor2icoor(points(foo)%z, iz)
    print *, ix, iy, iz,                   &
             cube(ix,iy,iz),               &
             DBLE(cube(ix,iy,iz)) / dval
  enddo

!-----------------------------------------------------
contains
!-----------------------------------------------------
  ! or maybe, we can write a function ?
subroutine fcoor2icoor(in, out)
  double precision, intent(in)       :: in
  integer, intent(out)               :: out
  double precision                   :: invalue
  integer                            :: outvalue

  invalue = (in + 2.0) / 2.0
  outvalue = int(invalue * real(DIM/2))

  ! add molly-guard here
  out = outvalue
  if (outvalue .LT. 1)      out = 1
  if (outvalue .GE. DIM)    out = DIM-1

end subroutine
!-----------------------------------------------------

subroutine clear_cube(cube)
  type(integer), dimension(:,:,:), intent(out)  :: cube

  integer        :: i, j, k

  do i=lbound(cube, 1), ubound(cube, 1)
    do j=lbound(cube, 2), ubound(cube, 2)
      do k=lbound(cube, 3), ubound(cube, 3)
        cube(i, j, k) = 0
      enddo
    enddo
  enddo

end subroutine
!-----------------------------------------------------
end program voxelize
!-----------------------------------------------------