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

  implicit none

  integer, parameter                          :: DIMC = 320
  integer, dimension(:,:,:), allocatable      :: cube
  type(t_point3d), dimension(:), allocatable  :: points
  integer                                     :: errcode, foo, argc
  integer                                     :: ix, iy, iz
  integer                                     :: nbr_points, maxcube
  double precision, dimension(4)              :: KA, KB, KM
  double precision                            :: dmaxcube, delta
  character(200)                              :: filename, string

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

  argc = IARGC()
  if (2 .NE. argc) then
    STOP ": VOXELIZE NEED PARAMETERS !"
  endif

  call getarg(1, filename)
  call getarg(2, string)       ;   read (string, *) delta
  write(0, "( ' --- delta = ', F11.6)") delta

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

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

  KA(1) =  -1.3402     ;     KA(2) =   1.5245
  KA(3) =   1.0966     ;     KA(4) =  -2.3423
  KB(1) =  -1.2100     ;     KB(2) =   1.3685
  KB(3) =   1.3237     ;     KB(4) =  -2.3992
  call interp4dp(KA, KB, KM, delta)
  write(0, "(' --- coefs = ', 4F11.6)") KM
  call compute_pickover(points, KM)
  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

  maxcube = MAXVAL(cube)
  dmaxcube = DBLE(maxcube)
  write(0, *) "--- maxval(cube) = ", maxcube

  call spit_cube_as_union(filename,   cube,       &
                          maxcube/2000, dble(9000.00))

  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(DIMC/2))

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

end subroutine
!------------------------------------------------------------
!           USELESS USE OF LOOPS !

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 print_cube(cube, points, scaling)
  type(integer), dimension(:,:,:), intent(in)  :: cube
  type(t_point3d), dimension(:), intent(in)    :: points
  double precision, intent(in)                 :: scaling

  integer                     :: foo

  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)) / scaling
  enddo

end subroutine

!------------------------------------------------------------
!           generation Povray STL source file !

subroutine spit_cube_as_union(fname, voxels, limit, scaling)
  character(*), intent(in)                     :: fname
  type(integer), dimension(:,:,:), intent(in)  :: voxels
  integer, intent(in)                          :: limit
  double precision, intent(in)                 :: scaling

  integer             :: fd, errcode
  integer             :: ix, iy, iz, maxv
  integer             :: nbrvox = 0
  double precision    :: bx, by, bz, valeur
  character(200)      :: chaine

  !   molly-guard
  maxv = limit
  if (maxv .LT. 2)   maxv = 2

  open (newunit=fd, file=trim(fname),           &
        status='replace',                       &
        action='write', iostat=errcode)
  if (0 .NE. errcode) then
    STOP ' : SPIT UNION, FAIL OPEN OUTPUT FILE'
  endif

  write(fd, *) "// generated file, don't touch it bastard !"
  write(fd, *) "// version 2.09"
  write(fd, *) "#declare DIMC = ", DIMC, ";"
  write(fd, *) "#declare Voxels = object {"
  write(fd, *) "union {"

  bx = 0.0   ;   by = 0.0   ;   bz = 0.0

  do ix=lbound(voxels,1), ubound(voxels,1)
    do iy=lbound(voxels,2), ubound(voxels,2)
      do iz=lbound(voxels,3), ubound(voxels,3)

        if (cube(ix,iy,iz) .LT. maxv) then
          ! print *, "foo = ", foo, cube(ix,iy,iz)
          cycle                                   ! REDO FROM START
        endif

        nbrvox  = nbrvox + 1
        bx = bx + dble(ix)
        by = by + dble(iy)
        bz = bz + dble(iz)
        valeur = DBLE(cube(ix,iy,iz)) / scaling
        ! XXX
        if (valeur .GT. 1.5) then
          valeur = 1.5
        endif

        write(chaine, "( 'translate <', I4, ',', I4, ',', I4, '> ' )")   &
                            ix, iy, iz
        write(unit=fd,                                                   &
              fmt="( 'object { VOXEL scale ', F11.6, 1X, A, ' }'  )",    &
              iostat=errcode)                                            &
              valeur, trim(chaine)
      enddo
    enddo
  enddo

  write(fd, *) "} }"
  write(fd, *) "// limit = ", limit
  write(fd, "( '#declare NbrVox = ', I9,    ';' )")   nbrvox
  write(fd, "( '#declare Bary_X = ', F11.6, ';' )")   bx / dble(nbrvox)
  write(fd, "( '#declare Bary_Y = ', F11.6, ';' )")   by / dble(nbrvox)
  write(fd, "( '#declare Bary_Z = ', F11.6, ';' )")   bz / dble(nbrvox)

  close(fd)

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