IT şi electronică > Software
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: ---:)
--- Terminare cod ---
e-
Navigare
[#] Pagina următoare
Du-te la versiunea completă