program  bulles2incc

  integer, parameter                  :: szvert  = 2000000
  integer, parameter                  :: szedges = 2000000

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

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

  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
    if (nbvert .GE. szvert) then
        STOP "vertex table is full, sorry."
    endif
  enddo
999 continue

  close (unit=42)

end subroutine

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

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

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

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

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

!----  some sanity check
    if ( (pa .lt. 0) .or. (pa .gt. nbvert) )  STOP "OVERFLOW PA"
    if ( (pb .lt. 0) .or. (pb .gt. nbvert) )  STOP "OVERFLOW PB"
    if (pa .eq. pb)                           STOP 'PA == PB'

    fa(nbedges) = pa
    fb(nbedges) = pb

!---- hunting for Degenerate cylinder
    if ( (px(pa) .eq. px(pb))    .and.       &
         (py(pa) .eq. py(pb))    .and.       &
         (pz(pa) .eq. pz(pb)) )   then
        !!! STOP 'DEGENERATED CYLINDER'
        print *, "ignore edge: ", nbedges
        cycle

    endif

    nbedges = nbedges + 1
    if (nbedges .GE. szedges) then
        STOP "edges table is full, sorry"
    endif
  enddo

999 continue

  nbedges = nbedges - 1

  close (unit=42)

end subroutine

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

  print *, "// spit edges to: ", fname

  open (unit=42, file=fname, status='replace', action='write')
  do idx = 1, nbedges
    write (42, '(6F15.9, " ", i7)') px(fa(idx)), py(fa(idx)), pz(fa(idx)),    &
                         px(fb(idx)), py(fb(idx)), pz(fb(idx)), idx
  end do
  close (unit=42)

end subroutine

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

end program

