Fortraneries/GravityField/realfield.f90

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

  implicit none

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

!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
subroutine compute_barycentre_bodies(astres, bcx, bcy)
  type(massbody), intent(in)     :: astres(:)
  real, intent(out)              :: bcx, bcy
  integer                        :: foo
  real                           :: cx, cy

  !   May be we have to use DOUBLE PRECSION here ?
  cx = 0.0
  cy = 0.0
  do foo=1, ubound(astres, 1)
    cx = cx + astres(foo)%posx
    cy = cy + astres(foo)%posy
  enddo
  bcx = cx / real(ubound(astres, 1))
  bcy = cy / real(ubound(astres, 1))
end subroutine
!-----------------------------------------------------------------------
subroutine print_barycentre_bodies(astres, title)
  type(massbody), intent(in)     :: astres(:)
  character(len=*), intent(in)   :: title

  real                           :: cx, cy
  call compute_barycentre_bodies(astres, cx, cy)
  print *, "barycentre {", title, "}  ", cx, cy

end subroutine
!-----------------------------------------------------------------------
!-
!              make a few solid body to play with...
!-     
!       planets    :   an array of type(massbody) to be filled
!       coef       :   for setting the mass of the body
!       sx, sy     :   borders of the universe
!-
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
      !-
      !   the first planet is the home of Johnny Root
      !-
      planets(1)%posx   = sx / 2
      planets(1)%posy   = sy / 2
      planets(1)%mass   = 31e8
      planets(1)%serial = 1337
      planets(1)%speed  = 6.666
    else
      !-
      !   others are planets for peones
      !-
      planets(foo)%posx     = rand() * real(sx-1)
      planets(foo)%posy     = rand() * real(sy-1)
      planets(foo)%mass     = 7.12e6 + coef*foo*3
      planets(foo)%heading  = 2 * 3.141592654 * rand()
      if (rand() .LT. 0.15) planets(foo)%speed = 3.14159
      planets(foo)%serial   = foo*2 + 120
    endif

    write (*, fmt) foo, planets(foo)
  enddo
end subroutine
!-----------------------------------------------------------------------
!-
!           the basis of the kluge
!-
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) )
  dist = (rx*rx) + (ry*ry)
  if (dist .LT. 0.08) then
    ! write (0, *) "dist too small ", dist
    compute_gravity = 0e0
  else
    ! ??? compute_gravity = body%mass / (dist ** 2)
    compute_gravity = body%mass / dist
  endif

end function
!-----------------------------------------------------------------------
!-
!   Export a massbody area to a text file. no error check, wtf ?
!-
subroutine save_bodies_to_txt_file (astres, fname)
  type(massbody), intent(in)        :: astres(:)
  character(len=*), intent(in)      :: fname

  character(50)                     :: fmt
  integer                           :: io, idx

  write(0, "('saving planets to ', A20)") fname
  fmt = "( 2(F9.3, ' ')   2(F9.3, '  '), F14.3,  I8)"
  open(newunit=io, file=fname)
  do idx = 1, ubound(astres, 1)
    write(io, fmt) astres(idx)
  enddo
  close(io)
end subroutine
!-----------------------------------------------------------------------
!-
!   Compute the gravity field in a pre-allocated array relative
!   to the massbody 'moon'. Nobody know where the magic numbers
!   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.018
    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", I4, F9.6)') t3, dummy
  enddo

end subroutine
!-----------------------------------------------------------------------
!-
!    dump a field of reals numbers to disk - preliminary version
!-
subroutine dump_a_field_to_file(field, fname)
  real, dimension(:,:), intent(in)  :: field
  character(len=*), intent(in)      :: fname

  integer                         :: header(8)
  integer                         :: io

  print *, "D) field size  ", ubound(field, 1), "W", ubound(field, 2), "H"
  print *, "D) filename    ", fname

  header     = 0
  header(1)  = 574908040             ! magic number
  header(2)  = 1                     ! this is a dump of real field
  header(3)  = ubound(field, 1)
  header(4)  = ubound(field, 2)
  header(5)  = 666

  open(newunit=io, file=fname, form='unformatted')
  
  write(io) header
  write(io) field

  close(io)

end subroutine
!-----------------------------------------------------------------------
!-
!    load a real field from file  - preliminary version
!-
subroutine load_a_field_from_file(field, fname)
  real, dimension(:,:), intent(in)  :: field
  character(len=*), intent(in)      :: fname

  integer                         :: header(8)
  integer                         :: io, foo

  print *, "L) field size  ", ubound(field, 1), "W", ubound(field, 2), "H"

  !-
  ! how to check if the field array was valid ?
  !-

  open(newunit=io, file=fname, form='unformatted', status='old',     &
                                                   action='read')
  read(io) header
  do foo=1, 8
    print *, foo, header(foo)
  enddo

  STOP ' --- FUCKED UP BEYOND ALL REPAIR ---'

  close(io)

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

end module