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

  implicit none

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

  write(0, *) "--- start of voxelize"

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

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

  coefs(1) =  2.23     ;     coefs(2) =  0.42
  coefs(3) = -0.64     ;     coefs(4) = -2.42

  KB(1) =  1.51        ;     KB(2) = -1.89
  KB(3) =  1.69        ;     KB(4) =  0.79
  call compute_pickover(points, KB)

  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, *) "--- maximum of the cube= ", 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

  dval = DBLE(MAXVAL(cube))
  write(0, *) "--- end of voxelize"

!-----------------------------------------------------
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
!------------------------------------------------------------
subroutine spit_cube_as_union(fname, datas, notused)
  character(*)                                 :: fname
  type(integer), dimension(:,:,:), intent(in)  :: datas
  integer                                      :: notused

  integer             :: fd, errcode, foo

  open (newunit=fd, file='WS/k-pick.txt',            &
        status='unknown', position='append',         &
        action='write', iostat=errcode)
  if (0 .NE. errcode) then
    STOP ' : SPITUNION : FAIL OPEN OUTPUT FILE'
  endif


  close(fd)

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