first dry run

GravityField/essai.f90

@@ -1,62 +1,54 @@
+!
+!       project "gravity field"
+!
 program essai
+  use realfield
   implicit none
 
-  integer                 :: ix, iy
-  real                    :: fx, fy
-  real                    :: foo, bar, maxi, mini
+  type(massbody)          :: planet
+  integer                 :: foo
+  planet%posx = 1337.314
+  planet%posy = 1664.666
+  planet%mass = 1e4
 
-  maxi = 0.0
-  mini = 9e9
+  do foo=1, 10
+    call build_a_field(800, 600, planet)
+    planet%posy = planet%posy + 51.45
+  enddo
 
-  do ix=1, 2000
+  STOP 'YOLO'
+
+contains        !------------------------------------------
+!
+subroutine build_a_field(szx, szy, moon)
+  integer, intent(in)                :: szx, szy
+  type(massbody), intent(in)         :: moon
+
+  integer                            :: ix, iy
+  real                               :: fx, fy
+  real                               :: grav, maxi, mini
+  integer                            :: errcode
+  real, dimension(:,:), allocatable  :: field
+
+  allocate(field(szx, szy), stat=errcode)
+   
+  do ix=1, szx
     fx = real(ix)
-    do iy=1, 2000
+    do iy=1, szy
       fy = real(iy)
-
-      foo = rdist(fx, fy, 222.22, 765.432)
-      bar = gravity(foo, 1337.0)
-
-      maxi = max(maxi, bar)
-      mini = min(mini, bar)
-
+      grav = compute_gravity(ix, iy, moon) / 2180800.0
+      field(ix,iy) = grav
     enddo
   enddo
 
-  print *, "dist : ", mini, maxi
+  maxi = maxval(field)
+  mini = minval(field)
+  print *, "field: ", mini, maxi, maxi-mini
 
-contains        !------------------------------------------
+  deallocate(field)
 
-function gravity(distance, masse)
-  real, intent(in)      :: distance, masse
-  real                  :: gravity
-  real                  :: computed
+end subroutine
 
-  if (distance .LT.  0.010) then
-    computed = 0.0
-  else
-    computed = masse / (distance ** 2)
-  endif
-
-  gravity = computed
-
-end function
-
-                !------------------------------------------
-
-function rdist(ax, ay, bx, by)
-  real, intent(in)       :: ax, ay, bx, by
-  real                      :: rdist
-  real                      :: rx, ry
-
-  rx = real(ax-bx)
-  ry = real(ay-by)
-
-  rdist = sqrt( (rx*rx) + (ry*ry) )
-
-end function
-
-
-                !------------------------------------------
                 !------------------------------------------
 
 end program