problema gravitatie

[florin_try]

Int-un post anterior cineva a mentionat ca ar dori sa aiba un program care sa simuleze miscarea planetelor si animatie pe aceasta tema.

M-am gindit sa scriu un programel pe aceasta tema, bineinteles am inceput cu ceva simplu si static. Programul e static in sensul ca primeste date de intrare dintr-un fisier la inceputul executiei si tipareste in fisiere datele de ierire (energii si coordonate) si nu permite iteractia utilizatorului in timpul executiei. Bineinteles, un profesionist se va amuza de chestiile ‘dinamice’ si va rula static iar animatiile le va face din fisierele obtinute la finalul executiei.
Am sa il fac si dinamic in flash sau java (urasc java btw) sau poate chiar html5/java-script mai ales ca e timpul sa ma joc un pic cu ‘canvas’ din html-5 (atit cit e acum standardizata) sau vre-un webkit – vad eu cind o sa am timp.

Pentru pasul 0 , programul l-am scris in fortran90/95 din urmatoarele motive: i) sunt o multime de compilatori free (l-am testat pe vechiul compaq v6 si noul intel ifc) ii) este de 10 ori mai rapid decit java iii) este posibil ca problemele pe care le-am avut si le voi prezenta sa fie usor observate de cineva cu experienta in gravitatie computationala, caci probabil foloseste fortran deasemeni. iv)  e usor de inteles si poate fi usor tradus in limbaj C-like fie manual fie automat (intotdeauna sunt sceptic in conversiile automate).

Datele de intrare sunt in 2 fisiere locazilate in acelasi director ca si  executabilul.
Un fisier este ‘prms.gvr’ si contine:
1)   pasul de integrare (in ore) si numarul de pasi de integrare [pe prima linie]
2)   cit de des se tiparesc datele [pe  a doua linie]
3)   un flag care acum nu conteaza

Iata un exemplu de fisier de intrare prms.grv
0.002   4000000  time_step N_steps (pasul de integrare este 0.002 ore )
10000   Nprint (la fiecare 10000 pasi printeaza date)
1      integrator

Un alt fisier este cel care contine configuratia initiala a planetelor.
Pe prima linie e numarul de plaanete
Linia a 2-a e ignorata
De la linia a treia pina la linia ce contine vxyz la inceputul liniei sunt puse masa si coordonatele carteziene (la urma e si un string ce contine numele plaentei):
masa   x   y    z  numele_planetei 
pentru fienare obiect in parte
dupa care urmeaza vxyz
iar de la vxyz se citesc vitezele initiale ale planetelor (cele 3 componente)
Orice linie care incepe cu ! este ignorata (asta da flexibilitate, voi arata mai departe)

Iata un exemplu de fisier de intrare config.grv cu Soarele si primele planete pina la Jupiter.

6
xyz
330.0d3     0.0  0.0 0.0    Sun
0.0553d0    0.39 0.0  0.0   Mercury
0.815d0     0.72  0.0  0.0  Venus
1.0d0       1.0 0.0 0.0     Earth
0.107       1.52  0.0  0.0   Mars
317.83      5.2   0.0  0.0   Jupiter
vxyz
 0.0d0 0.0d0 0.0d0  Sun
 0.0 -172.332d3 0.0d0  Mercury
 0.0  -126.936d3  0.0d0 Venus
 0.0 -107.300d3 0.0d0  Earth
 0.0  -86.868d3 0.0d0  Mars
 0.0  -49.392d3  0.0d0  Jupiter


unitatile de intrare sunt:
masa = unitati masa pamint
distanta = unitati astronomice
viteza = km/h


Programul contine:
1)   O parte de definitie a constantelor si variabilelor in modulul data_module
2)   O parte de intrate in modulul read_input, subrutina get_indata
3)   O parte de output in modulul out_module, subrutina write_them
4)   O parte de calcul a fortelor si de integrare in modulul compute. Fortele sunt calculate in subrutina forces iar integrarea e facuta in integrate_VV cu un algoritm reversibil in timp si de ordin 2, insa fara pas adaptiv (cel mai simplu mod de a integra). Algoritmul de integrare e in wikipedia.

Programul e general in sensul ca poate lucra cu un numar oarecare N de planete. Suficient de rapid pentru N de ordinul sutelor. Cu totul alte metode trebuie insa folosite pentru galaxii (miliarde de obiecte).


Ca iesire progtamul genereaza 3 fisiere
1) eng.hist care contine
timp (zile)  ,  energie_kinetica (MJ)  ,  Energie_potentiala (MJ) , Energie Totala (MJ), impuls (in SI)

2) xyz.hist care contine pe fiecare linie timpul si coordonatele x y z ale fiecarei plannete
iar =ultima inregistrare e distanta |r1-r2| (sic) . distantele sunt in unitati astronomice si la intrare si la iesire.

3) vxyz.hist e similar cu xyz.hist dar contine viteze in km/h

[florin_try]

Iata si codul:



module precision
 integer,parameter :: rp = kind(0.0d0) !SELECTED_REAL_KIND(14,200)
end module precision

module reduced_units
 use precision
 real(rp), parameter :: unit_mass = 5.9742e24_rp ! Earth mass
 real(rp), parameter :: unit_length = 149597871e3_rp ! distance Earth - Sun at its max.
 real(rp), parameter :: unit_time = 3600.0e0_rp  ! 1 hour

 real(rp), parameter :: unit_force = unit_mass * unit_length / unit_time**2
 real(rp), parameter :: unit_energy = unit_force * unit_length
 real(rp), parameter :: unit_speed = unit_length/unit_time
 real(rp), parameter :: unit_linmom = unit_mass *  unit_speed
end module reduced_units

module data_module
 use reduced_units
  implicit none;
  integer Nprint ! how often print (in timesteps)
  integer N_objects ! number of planets (satelites etc)
  real(rp), parameter :: KAPPA = 6.67428e-11_rp !m^3/Kg/s^2;
  real(rp), allocatable :: xyz(:,:),vxyz(:,:),fxyz(:,:),mass(,axyz(:,:),pxyz(:,:)
  real(rp) Enkin, Enpot, Entot; ! kinetic, potential and total energy
  real(rp) linmom(3) !linear momentum
  character(250), allocatable :: name_object(; ! ignore it for now
  character(250) :: file_name_cfg = 'config.grv'; ! the name of the input file the name of the input file where the config. of planets is
  character(250) :: file_name_prm = 'prms.grv';   ! the name of the input file with whatever parameters
  logical  close_contact ;   ! true if objects colide
  real(rp), parameter :: CONTACT=7000.0e3_rp/unit_length   ! I set this to be  a distance of colision (for now)
  ! normally the distance of collision depends on the radius of each object.
  real(rp) time_step; ! the interval of time over which the equations of motion are asdvanced
  integer N_steps    ! the number of steps to be performed
  integer flag_integration_method; ! just go with VV for now.
  ! Clearly time addaptive integration methods are required for precision.

end module data_module

module allocate_arrays
implicit none ! for "obvious" reasons
contains
  subroutine xyz_alloc    ! memory management
   use data_module, only : N_objects, xyz,vxyz,fxyz,mass,axyz,pxyz
    integer N
    N=N_objects;
    allocate(xyz(N,3),vxyz(N,3),fxyz(N,3),mass(N),axyz(N,3),pxyz(N,3))
  end subroutine xyz_alloc
  subroutine xyz_DEalloc  ! free the memory
   use data_module, only : N_objects, xyz,vxyz,fxyz,mass,axyz,pxyz
     deallocate(xyz,vxyz,fxyz,mass,axyz,pxyz);
  end subroutine xyz_DEalloc
end module allocate_arrays

module read_input
implicit none
private :: get_init_config
private :: read_prms
public :: get_indata
contains
 subroutine get_indata
  use data_module  ! read the input data
  call get_init_config(trim(file_name_cfg));
  call read_prms(trim(file_name_prm));
  flag_integration_method=1;
 end subroutine get_indata
 
 subroutine get_init_config(nf)  ! read the config (coordinates, masses, velocities)
 use allocate_arrays
 use data_module
 character(*), intent(IN) :: nf
 character*1000 line
 integer i,j,k
 logical in_loop

 open(unit=10,file=trim(nf))   
 read(10,*) N_objects
 print*,trim(nf),N_objects
 call xyz_alloc
 read(10,*)
 
  i=0
  in_loop = .true.
  do while (in_loop)

    read(10,'(A1000)') line
   if (len(trim(line))/=0)then ! if line is not empty
   in_loop = .not.(line(1:4)=='vxyz')
    print*,trim(line)
    if (line(1:1) /= '!'.and.in_loop)  then !if is not a coment and if still in loop
      i = i + 1
      if (i < N_objects+1) read(line,*) mass(i),xyz(i,:)
    endif
!    print*,i,mass(i),xyz(i,:)
    endif !line not empty
  enddo
  if (i < N_objects) then
    print*,'ERROR: More date required in config.grv when read coordinates' ;STOP
  endif
  i=0
  in_loop = .true.
  do while (in_loop)
    read(10,'(A1000)',end=10) line
   if (len(trim(line))/=0)then ! if line is not empty
    print*,trim(line)
    if (line(1:1) /= '!')  then !if is not a coment
      i = i + 1
      if (i < N_objects+1) read(line,*) vxyz(i,:)
    endif
!    print*,i,vxyz(i,:) ! line not empty
   endif
  enddo

10 continue

  if (i < N_objects) then
    print*,'ERROR: More date required in config.grv when read velocities';STOP
  endif

 mass = mass * unit_mass ! from Earth mass to SI
 xyz = xyz * unit_length ! from au to SI
 vxyz = vxyz * 1000.0_rp/3600.0_rp  ! from km/h to SI
 call validate_initial_momentum; ! re-assign it to zero by extracting the excess momentum from heaviest object
 close (10)
  print*,'fcg read'
 end  subroutine get_init_config
 
 subroutine read_prms(nf)  ! read some other parameters like number of steps or the size of the integration step
  use data_module
  character(*), intent(IN) :: nf
  open(unit=10,file=trim(nf))
  read(10,*) time_step, N_steps
  read(10,*) Nprint
  time_step = time_step * 3600.00_rp   ! from hours (input) to seconds
  read(10,*) flag_integration_method ! 1=VV 2=RK4
  print*,'prms read'
 end  subroutine read_prms

 subroutine validate_initial_momentum  ! Make sure we start from a linear momentum zero. (not required)
 use data_module
 implicit none
 integer i,imax
 real(rp) mmax , t(3)
 mmax=mass(1);imax=1
 do i = 1, N_objects
  if (mass(i) > mmax) then
    mmax=mass(i)
   imax=i
  endif
 enddo
 do i = 1, N_objects
  t=0.0_rp
  if (i /= imax) then
   t = t + mass(i)*vxyz(i,:)
  endif
  vxyz(imax,:) = -t( /mmax
 enddo
 end  subroutine validate_initial_momentum


end module read_input

module compute
implicit none

contains
subroutine forces  ! EVALUATE forces and energies
  use data_module
  integer i,j,k,N
  real(rp) t(3), mass_i,r_ij,r_sq_ij,i_r_ij,ff,dr(3),fi(3),Uij

  close_contact=.false.
  Enpot = 0.00_rp;
  fxyz=0.00_rp
  do i = 1, N_objects-1
    mass_i = mass(i); t = xyz(i,:)
    fi( = 0.0_rp
    do j = i+1, N_objects
       dr = t-xyz(j,:)
       r_sq_ij = dot_product(dr,dr)
       r_ij = sqrt(r_sq_ij) ; i_r_ij = 1.0_rp/r_ij
      if (r_ij < CONTACT) then
        close_contact=.true.
       return
      endif
       Uij = - (KAPPA * mass(j) * i_r_ij)  * mass_i
       Enpot = Enpot + Uij
       ff = Uij * i_r_ij * i_r_ij
       fi = fi + (ff*dr)
       fxyz(j,:)  = fxyz(j,:) - (ff*dr)
     enddo
    fxyz(i,:)  = fxyz(i,:) + fi(
  enddo

  do i = 1, N_objects
  axyz(i,:) = fxyz(i,:) / mass(i)
  enddo
end subroutine forces



subroutine integrate_VV   ! the time-reversible integrator
 use data_module;
 implicit none
 integer i
 real(rp) d
 xyz = xyz + vxyz*time_step+(0.5_rp*time_step*time_step)*axyz !r(t+1)=r(t)+v(t)*dt+dt^2/2*a(t)
 Enkin=0.00_rp;linmom=0.00_rp;d=0.00_rp;
 do i = 1, ubound(mass,dim=1)
  Enkin=Enkin + 0.50_rp * mass(i) * dot_product(vxyz(i,:),vxyz(i,:))
 enddo
! linmom=linmom/d
 Entot = Enkin+Enpot;
 vxyz = vxyz + axyz *  (0.50_rp*time_step) ! v(t+1/2) = v(t) + a(t)*dt/2
 call forces  ! a(t+1)
 vxyz = vxyz + axyz * (time_step * 0.50_rp) ! v(t+1)=v(t+1/2)+a(t+1)*dt/2
 do i = 1,3 ; pxyz(:,i) = vxyz(:,i)*mass(; enddo
end subroutine integrate_VV

end module compute

module out_module
contains
subroutine write_them(i)
 use data_module;
 implicit none
 integer, intent(IN) ::  i
 integer k
 real(rp) t(3),ox
 logical, save :: first_time = .true.
 if (first_time) then
   first_time=.false.
   open(unit=10,file='eng.hist',recl=500)   
   open(unit=22,file='xyz.hist',recl=20000)
   open(unit=33,file='vxyz.hist',recl=40000)
 else
   open(unit=10,file='eng.hist',recl=500,access='append')   
   open(unit=22,file='xyz.hist',recl=20000,access='append')
   open(unit=33,file='vxyz.hist',recl=40000,access='append')
 endif

 ox = real(i,rp)*time_step/3600.00_rp/24.0_rp
 t=xyz(1,:)-xyz(2,:)
 print*,ox,Entot/1e6_rp,sqrt(dot_product((xyz(1,:)-xyz(2,:)),(xyz(1,:)-xyz(2,:)) ) )/unit_length!,&
                     !linmom,sum(fxyz);
 write(10,*) ox,Enpot/1e6_rp,Enkin/1e6,Entot/1e6_rp,linmom/unit_linmom
 write(22,*) ox ,(xyz(k,1:3)/unit_length,k=1,N_objects),&
             sqrt(dot_product(t,t))/unit_length
 write(33,*) ox ,(vxyz(k,1:3)*3600.0_rp/1000.0_rp,k=1,N_objects),&
 (sqrt(dot_product(vxyz(k,1:3),vxyz(k,1:3))) *3600.0_rp/1000.0_rp,k=1,N_objects)

close(10)
close(22)
close(33)

end subroutine write_them
end module out_module

program main_program
use data_module
use read_input
use compute
use out_module

implicit none
integer i,j,k

call get_indata;  ! go read the data
do i = 1, N_steps
 call integrate_VV   ! integrate
 if (close_contact) then;  print*,'COLISION'; STOP;  endif
 if (mod(i,Nprint)==0) call write_them(i)
enddo
end program main_program

[florin_try]

1. Din pacate au aparut acesti stupizi smiles. (ca si cum cineva i-ar folosi la ceva),

 Adi, poti scoate sau dezactiva acei "smiles" din postul de mai sus cu codul asa incit sa apara doar codul cum e el? Poate cineva va gasi util acel cod sau il poate traduce cu usurinta intr-un limbaj pentru web.(interactiv, dupa ce modifica partea IO).


2. Iata si niste rezultate cu acel cod.
 2.1 Conservarea energiei:
 

 2.2 Energia potentiala versus cea cinetica trebuie sa fie o dreapta de panta -1 si toate punctele sa fie localizate pe aceasta drapta daca F-Nabla U
   

2.3 Coordonatele y in functie de timp arata periodicitatea si diferenta intre periodicitate intre planete. Observati caci cam corespund cu cele reala.
   
2.4 Traiectoriile in spatiu ale planetelor. (Am plotat y versus x)
   

EDIT: Nu stiu de ce calitatea imaginilor e asa de proasta.

Voi reveni insa cu intrebari.

[Adi]

Felicitari pentru initiativa si ducerea ei la bun sfarsit. Eu sunt cel care am propus problema. Iar in ceea ce priveste smiley-urile, nu pot sa dezactivez in acel post, dar tu poti pune codul intr-un fisier .txt si sa il atasezi, sau cineva poate da reply la codul tau si atunci ii apare codul fara smileyuri.

[Electron]

Pentru a evita problema cu interpretarea caracterelor, trebuie folosite tagurile [ code ] si [ / code ] (accesibile si prin butonul cu "#" de pe formularul de editare al mesajelor).

Cod: [Selectează]

:)

e-

[florin_try]

 Iata cum arata traiectoriile planetelor in jurul Pamintului. In origine (0,0) e Terra.

 1. Mercur cum ar fi vazut de pe Terra:



 2. Venus in Jurul terrei:


 3. Marte in jurul Terrei:


 4. Si in fine Jupiter in Jurul Terei:


 Voi reveni cu o versiune interactiva ce va face animatie, probabil in java (caci tot urasc acest limbaj) in care parametrii de executie (numar planete, masa fiecarei planete, viteza animatie) sa fie introdusi si modificati interactiv in timpul executiei.

[Adi]

Foarte frumoase grafice.

In Java ce imi imaginez eu ca trebuie vazut este sistemul solar vazut de un observator de sus din afara sistemuliu solar. Avantajul in Java este ca poate rula pe orice browser si oricine poate schimba parametri. Asa devine o unealta educationala pentru scoala. Elevii vor putea chiar sa schimba valoarea constantei gravitationale si sa vada cum ar evolua sistemul solar atunci. Asta chiar ar fi educativ ...

[florin_try]

Cum as putea include un applet java .class ?
E posibil prin bb-coduri?

[Adi]

Cum as putea include un applet java .class ?
E posibil prin bb-coduri?

Am adaugat extensia "class" ca permisa la atasamente la posturi.

[florin_try]

 Nu cred ca mi-am dat seama, era interesant daca putea fi pus pe forum asa incit cineva deschide sa vada animatia. Probabil trebuiesc privilegii de administrator... In html appletul se acceseaza prin <applet></applet>, credeam ca se poate si cu bb-coduri, aparent nu sau nu e asa simplu.
nevermind, insa atasez o arhiva zip.

 Arhiva trebuie dezarhivata si fisierul index.html din noul folder test1 creat la de-zipare trebuie deschis intr-un browser java enabled.
 Daca totul e OK (adica java instalat si fisierul index.html deschis cu modzila (de exemplu) )  atunci o animatie ar trebui sa inceapa, cu un fel de sistem solar definit cam la intimplare.

 Drag cu mousul dreapta e rotatie, drag cu mousul stinga e mutare, si rotita de la mause e zoom dar vad ca l-am pierdut in browser.
 (daca apletul e rulat din eclipse ( http://www.eclipse.org/ - vezi IDE pentru  java )sau netbeans ( www.netbeans.org/  ) evenimeentele de la rotita mousului sunt ascultate corect - vezi http://bugs.sun.com/bugdatabase/view_bug.do?bug_id=6516675 ; nu am avut timp de workaround )

 Codul ar fi desigur o versiune de plecare ..... basic ... abc-ul .... scheletul de plecare ... o animatie + basic Mouse control.... cine stie programare stie bine la ce ma refer.... mai departe trebuie adaugat o modalitate de a modifica parametri de intrare, viteza animatie, si multe altele, nu e greu dar consuma timp, aam sa mai agaug in weekend cite ceva.

[florin_try]

 
 Codul e prea amre ca sa incapa ca mesaj, dar a fost oricum atasat anterior in fisieul ./src/test1.java

 Cind o sa am timp o sa adaug si niste comenzi de intrare (modificarea interactiva a unor parametri (viteze planete, viteza animatie) si multe altele.