Fortraneries/GravityField/realfield.f90

70 lines
1.9 KiB
Fortran
Raw Normal View History

2022-11-28 13:47:44 +01:00
!
! project "gravity field"
!
!-----------------------------------------------------------------------
module realfield
implicit none
!-----------------------------------------------------------------------
! definition of structures
!
2022-11-28 13:47:44 +01:00
type massbody
real :: posx, posy
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-11-28 13:47:44 +01:00
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-11-28 13:47:44 +01:00
dist = sqrt( (rx*rx) + (ry*ry) )
2022-11-30 22:46:44 +01:00
if (dist .LT. 0.11) then
! 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
compute_gravity = body%mass / (dist ** 2)
2022-11-28 13:47:44 +01:00
endif
end function
!-----------------------------------------------------------------------
2022-11-30 03:24:24 +01:00
!-
! 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
!-----------------------------------------------------------------------
2022-11-28 13:47:44 +01:00
!-----------------------------------------------------------------------
end module