Fortraneries/GravityField/realfield.f90

290 lines
8.6 KiB
Fortran
Raw Normal View History

2022-12-03 02:25:37 +01:00
!-----------------------------------------------------------------------
!-
! some functions for the project "gravity field"
!-
2022-11-28 13:47:44 +01:00
!-----------------------------------------------------------------------
module realfield
2022-12-03 02:25:37 +01:00
use spitpgm ! XXX
2022-11-28 13:47:44 +01:00
implicit none
!-----------------------------------------------------------------------
2022-12-05 13:10:40 +01:00
!-
! definition of structures
2022-12-03 02:25:37 +01:00
!-
2022-11-28 13:47:44 +01:00
type massbody
2022-12-05 13:10:40 +01:00
real :: posx = 0, posy = 0
real :: heading = 0.29
2022-12-03 20:42:29 +01:00
real :: speed = 1.017
2022-11-29 20:31:45 +01:00
real :: mass = 1.0
integer :: serial = 666
2022-11-28 13:47:44 +01:00
end type
!-----------------------------------------------------------------------
contains
2022-11-29 20:31:45 +01:00
!-----------------------------------------------------------------------
2022-12-05 13:10:40 +01:00
subroutine compute_barycentre_bodies(astres, bcx, bcy)
2022-12-03 02:25:37 +01:00
type(massbody), intent(in) :: astres(:)
2022-12-05 13:10:40 +01:00
real, intent(out) :: bcx, bcy
2022-12-03 02:25:37 +01:00
integer :: foo
2022-12-05 13:10:40 +01:00
real :: cx, cy
2022-12-03 02:25:37 +01:00
2023-02-11 17:00:58 +01:00
! May be we have to use DOUBLE PRECSION here ?
2022-12-03 02:25:37 +01:00
cx = 0.0
cy = 0.0
do foo=1, ubound(astres, 1)
cx = cx + astres(foo)%posx
cy = cy + astres(foo)%posy
enddo
2022-12-05 13:10:40 +01:00
bcx = cx / real(ubound(astres, 1))
bcy = cy / real(ubound(astres, 1))
end subroutine
!-----------------------------------------------------------------------
2023-02-11 17:00:58 +01:00
subroutine print_barycentre_bodies(astres, title)
2022-12-05 13:10:40 +01:00
type(massbody), intent(in) :: astres(:)
2023-02-11 17:00:58 +01:00
character(len=*), intent(in) :: title
2022-12-05 13:10:40 +01:00
2023-02-11 17:00:58 +01:00
real :: cx, cy
2022-12-05 13:10:40 +01:00
call compute_barycentre_bodies(astres, cx, cy)
2023-02-11 17:00:58 +01:00
print *, "barycentre {", title, "} ", cx, cy
2022-12-03 02:25:37 +01:00
end subroutine
!-----------------------------------------------------------------------
!-
2022-12-05 13:10:40 +01:00
! 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
2022-12-03 02:25:37 +01:00
!-
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
2022-12-05 13:10:40 +01:00
!-
! the first planet is the home of Johnny Root
!-
2022-12-03 20:42:29 +01:00
planets(1)%posx = sx / 2
planets(1)%posy = sy / 2
2022-12-12 00:57:00 +01:00
planets(1)%mass = 31e8
2022-12-03 02:25:37 +01:00
planets(1)%serial = 1337
2022-12-05 13:10:40 +01:00
planets(1)%speed = 6.666
2022-12-03 02:25:37 +01:00
else
2022-12-05 13:10:40 +01:00
!-
! others are planets for peones
!-
2022-12-03 02:25:37 +01:00
planets(foo)%posx = rand() * real(sx-1)
planets(foo)%posy = rand() * real(sy-1)
2023-02-11 17:00:58 +01:00
planets(foo)%mass = 7.12e6 + coef*foo
planets(foo)%heading = 2 * 3.141592654 * rand()
2022-12-12 00:57:00 +01:00
if (rand() .LT. 0.15) planets(foo)%speed = 3.14159
2022-12-03 20:42:29 +01:00
planets(foo)%serial = foo*2 + 120
2022-12-03 02:25:37 +01:00
endif
2022-12-12 00:57:00 +01:00
2022-12-03 02:25:37 +01:00
write (*, fmt) foo, planets(foo)
enddo
end subroutine
!-----------------------------------------------------------------------
2022-12-05 13:10:40 +01:00
!-
! the basis of the kluge
!-
function compute_gravity(fx, fy, body)
real, intent(in) :: fx, fy
2022-11-28 13:47:44 +01:00
type(massbody), intent(in) :: body
real :: compute_gravity
real :: rx, ry, dist
2022-11-28 13:47:44 +01:00
rx = fx - body%posx
ry = fy - body%posy
2022-12-13 23:03:54 +01:00
! ??? dist = sqrt( (rx*rx) + (ry*ry) )
dist = (rx*rx) + (ry*ry)
2022-12-05 13:10:40 +01:00
if (dist .LT. 0.08) then
2022-11-30 22:46:44 +01:00
! write (0, *) "dist too small ", dist
2022-11-28 13:47:44 +01:00
compute_gravity = 0e0
2022-11-30 22:46:44 +01:00
else
2022-12-13 23:03:54 +01:00
! ??? compute_gravity = body%mass / (dist ** 2)
compute_gravity = body%mass / dist
2022-11-28 13:47:44 +01:00
endif
end function
2022-12-12 00:57:00 +01:00
!-----------------------------------------------------------------------
!-
! 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
2022-11-28 13:47:44 +01:00
2022-12-12 00:57:00 +01:00
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
2022-11-28 13:47:44 +01:00
!-----------------------------------------------------------------------
2022-11-30 03:24:24 +01:00
!-
! Compute the gravity field in a pre-allocated array relative
2022-12-12 00:57:00 +01:00
! to the massbody 'moon'. Nobody know where the magic numbers
2022-11-30 03:24:24 +01:00
! 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
2022-12-03 02:25:37 +01:00
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))
2022-12-12 00:57:00 +01:00
tmpf = tmpf * 0.018
2022-12-03 02:25:37 +01:00
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()
2022-12-12 00:57:00 +01:00
write(0, '(" dummy", I4, F9.6)') t3, dummy
2022-12-03 02:25:37 +01:00
enddo
2022-12-12 00:57:00 +01:00
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)
2022-11-30 03:24:24 +01:00
end subroutine
!-----------------------------------------------------------------------
2022-11-28 13:47:44 +01:00
!-----------------------------------------------------------------------
end module