!-----------------------------------------------------------------------
!-
!      some functions for the project "gravity field" 
!-
!-----------------------------------------------------------------------
module realfield
  use spitpgm               ! XXX

  implicit none

!-----------------------------------------------------------------------
!      definition of structures
!-
type massbody
  real                  :: posx, posy
  real                  :: heading       = 33.21
  real                  :: speed         = 1.007
  real                  :: mass          = 1.0
  integer               :: serial        = 666
end type

!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
subroutine barycentre_bodies(astres)
  type(massbody), intent(in)     :: astres(:)

  real                           :: cx, cy
  integer                        :: foo

  !-
  !   May be we have to use DOUBLE RPECSION here ?
  !-
  cx = 0.0
  cy = 0.0
  do foo=1, ubound(astres, 1)
    cx = cx + astres(foo)%posx
    cy = cy + astres(foo)%posy
  enddo
  cx = cx / real(ubound(astres, 1))
  cy = cy / real(ubound(astres, 1))
  print *, cx, cy

end subroutine
!-----------------------------------------------------------------------
!-
!     make a few solid body to play with...
!-
subroutine create_some_planets(planets, coef, sx, sy)
  type(massbody), intent(inout)     :: planets(:)
  real, intent(in)                  :: coef
  integer, intent(in)               :: sx, sy

  integer                           :: foo
  character(100)                    :: fmt

  fmt = "(I4, '  | ', 2(F10.2, '  '), ' | ', 2F9.3, '  ', e12.3, I7)"

  do foo=1, ubound(planets, 1)
    if (foo .EQ. 1) then
      planets(1)%posx   = 10
      planets(1)%posy   = 10
      planets(1)%mass   = 7e8
      planets(1)%serial = 1337
    else
      planets(foo)%posx     = rand() * real(sx-1)
      planets(foo)%posy     = rand() * real(sy-1)
      planets(foo)%mass     = 7e6 + coef*foo
      planets(foo)%serial   = foo
    endif
    write (*, fmt) foo, planets(foo)
  enddo

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

function compute_gravity(fx, fy, body)
  real, intent(in)              :: fx, fy
  type(massbody), intent(in)    :: body
  real                          :: compute_gravity
  real                          :: rx, ry, dist

  rx = fx - body%posx
  ry = fy - body%posy
  dist = sqrt( (rx*rx) + (ry*ry) )
  if (dist .LT. 0.11) then
    ! write (0, *) "dist too small ", dist
    compute_gravity = 0e0
  else
    compute_gravity = body%mass / (dist ** 2)
  endif

end function

!-----------------------------------------------------------------------
!-
!   Compute the gravity field in a pre-allocated array relative
!   to the massbody 'moon'. Nobody know where the magic number
!   come from, sorry.
!-
subroutine compute_a_field(field, moon)
  real, dimension(:,:), intent(out)  :: field
  type(massbody), intent(in)         :: moon

  integer                            :: ix, iy
  real                               :: fx, fy
  real                               :: grav

  ! print *, "pic size  ", ubound(field, 1), "W", ubound(field, 2), "H"
  ! print *, "mass body ", moon

  do ix=1, ubound(field, 1)
    fx = real(ix)
    do iy=1, ubound(field, 2)
      fy = real(iy)
      grav = compute_gravity(fx, fy, moon)
      field(ix,iy) = grav
    enddo
  enddo

end subroutine
!-----------------------------------------------------------------------
!-
!   compute a field with only one body; and write a pic file
!-
subroutine build_and_write_a_field(szx, szy, moons, fname)
  integer, intent(in)                :: szx, szy
  type(massbody), intent(in)         :: moons(:)
  character(len=*), intent(in)       :: fname

  real                                  :: maxi, mini
  integer                               :: errcode, foo
  real,    dimension(:,:), allocatable  :: field, tmpf
  integer, dimension(:,:), allocatable  :: greymap

  allocate(field(szx, szy), stat=errcode)
  allocate(tmpf(szx, szy), stat=errcode)

  field = 0.0
  do foo=1, ubound(moons, 1)
    call compute_a_field(tmpf, moons(foo))
    tmpf = tmpf * 0.019
    field = field + tmpf
  enddo

  maxi = maxval(field)
  mini = minval(field)
  ! print *, "field: ", mini, maxi, maxi-mini

  allocate(greymap(szx, szy), stat=errcode)
  greymap = 65533
  ! convert from real value to 16 bits int values
  where (field < 65530.0)
    greymap = int(field)
  end where

  call spit_as_pgm_16(greymap, trim(fname))

  ! make valgrind happy
  deallocate(field)
  deallocate(greymap)

end subroutine

!-----------------------------------------------------------------------
!-
!    Yes, I know, this is a disturbing kluge, but I like it :}
!    May be, it's time to read the doc of modern Fortran
!-
subroutine init_random()

  integer, dimension(3) :: tarray
  integer               :: t3
  real                  :: dummy
  call itime(tarray)
  t3 = 8971*tarray(1) + 443*tarray(2) + tarray(3)
  write(0, '(A,3I3,A,I6)') "sranding: ", tarray, " --> ", t3
  call srand(t3)
  ! after initializing the random generator engine,
  ! you MUST use it for initializing the initializer
  do t3=1, 4
    dummy = rand()
    write(0, *) 'dummy ', t3, dummy
  enddo

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

end module