program  bulles2incc

  integer, parameter                :: szvert  = 100000
  integer, parameter                :: szedges = 100000

  real, dimension(:), allocatable     :: px, py, pz
  integer, dimension(:), allocatable  :: fa, fb
  integer                             :: nbvert, nbedges

  integer                             :: foo
  double precision                    :: cx, cy, cz

  print *, "// sz max vertices  =  ", szvert
  print *, "// sz max edges     =  ", szedges

  allocate (px(szvert))
  allocate (py(szvert))
  allocate (pz(szvert))

  allocate (fa(szedges))
  fa = -666
  allocate (fb(szedges))
  fb = -666

!-
!    first step : load the vertices file
!-
  call load_vertices ("bulles.vertices")
  print *, "// loaded ", nbvert, " vertices"

!-
!    do some fancy computation :)
!-
  cx = 0.0  ;  cy = 0.0  ;  cz = 0.0
  do foo=1, nbvert
    cx = cx + px(foo)
    cy = cy + py(foo)
    cz = cz + pz(foo)
  end do
  cx = cx / nbvert
  cy = cy / nbvert
  cz = cz / nbvert

  print *, "// center ", cx, cy, cz

!-
!    second step : read the edges file and generate
!                  a (may be ?) usable pov include file.
!-
  nbedges = 0
  call load_edges ("bulles.edges")

!-
!    write a POVray include file
!-
  call spit_the_edges_values ("edges.values")

  stop ': BECAUSE WELL DONE BRO'

contains

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

subroutine load_vertices (fname)
  character(len=*), intent(in)  :: fname
  real                          :: rx, ry, rz

  print *, "// loading vertices from: ", fname

  open (unit=42, file=fname, status='old', action='read')
  nbvert = 1
  do
    read (42, *, end=999) rx, ry, rz
    px(nbvert) = rx
    py(nbvert) = ry
    pz(nbvert) = rz
    ! print *, nbvert, rx, ry, rz
    nbvert = nbvert + 1
  enddo

999 continue

  close (unit=42)

end subroutine

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

subroutine load_edges (fname)
  character(len=*), intent(in)  :: fname
  integer                       :: pa, pb

  print *, "// loading edges from: ", fname

  open (unit=42, file=fname, status='old', action='read')

  nbedges = 1
  do
    read (42, *, end=999) pa, pb
    ! print *, "// edge  ", nbedges, pa, pb

!----  same sanity check
    if ( (pa .lt. 0) .or. (pa .gt. nbvert) ) then
      STOP "OVERFLOW PA"
    endif
    if ( (pb .lt. 0) .or. (pb .gt. nbvert) ) then
      STOP "OVERFLOW PB"
    endif

    fa(nbedges) = pa
    fb(nbedges) = pb
    nbedges = nbedges + 1
  enddo

999 continue

  nbedges = nbedges - 1

  close (unit=42)

end subroutine

! ----------------------------------------------------------
subroutine spit_the_edges_values (fname)
  character(len=*), intent(in)  :: fname
  integer                       :: idx

  character(*), parameter  :: fm1 = "('<'(F8.4)', '(F8.4)', '(F8.4)'>')"
  character(*), parameter  :: fm2 = "'(cylinder{"//fm1//","//fm1//"}')"

  print *, "// spit edges to: ", fname
  ! print *, "// fmt 2:", fm1
  ! print *, "// fmt 2:", fm2

  open (unit=42, file=fname, status='replace', action='write')

  do idx = 1, nbedges
    write (42, '(6F12.7)') px(fa(idx)), py(fa(idx)), pz(fa(idx)),    &
                         px(fb(idx)), py(fb(idx)), pz(fb(idx))
  end do

  close (unit=42)

end subroutine

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

end program

