Test pianeta per ODE (a=1, e=?)
La teoria
module testpianeta
! Revisione giugno 2008
!
! Questo modulo serve per effettuare test numerici di ODE usando
! l'equazione del moto centrale di un pianeta attorno all'origine.
! L'eccentricita' deve essere inferiore ad uno in modo che
! l'orbita sia una ellisse (o un cerchio quando l'eccentricita e' zero).
!
implicit none
!
! Modificando qui si possono fare calcoli in quadruplice precisione.
!
integer,parameter,public::q_testpianeta=selected_real_kind(15,300)
!
integer,parameter,private::q=q_testpianeta
real(kind=q),private::ecciusata=-1.0_q,periodo=0.0_q
real(kind=q),private,allocatable,dimension(:)::eatempi
private :: pianeta,eqkeplero,stampatab
public :: orbita_esempio
contains
function pianeta(eccentri,eanomali)result(yt)
!
! Il semiasse maggiore vale 1, il periodo vale 2*π
! eccentri == eccentricita' compresa tra 0 incluso e 1 escluso.
! eanomali == anomalia eccentrica, la variabile indipendente.
! Ad eanomali=t=0 il pianeta sta al perielio.
! yt(1) == ascissa del pianeta
! yt(2) == ordinata del pianeta
! yt(3) == velocita' lungo l'ascissa
! yt(4) == velocita' lungo l'ordinata
! yt(5) == tempo a cui corrisponde l'anomalia eccentrica.
!
real(kind=q),intent(in)::eccentri,eanomali
real(kind=q),dimension(5)::yt
real(kind=q) :: eacos,easin,eraqu
if(0.0_q>eccentri.or.eccentri>=1.0) then
yt(:)=0.0_q
return
end if
easin=sin(eanomali)
eacos=cos(eanomali)
eraqu=sqrt(1.0_q-eccentri*eccentri)
yt(5)=eanomali-eccentri*easin ! il tempo
yt(1)=eacos-eccentri ! x
yt(2)=eraqu*easin ! y
yt(3)=-easin/(1.0-eccentri*eacos) ! vx
yt(4)=eraqu*eacos/(1.0_q-eccentri*eacos) ! vy
return
end function pianeta
!
subroutine eqkeplero(t,eccentri,eanomali)
!
! Noto il tempo e l'eccentricita' trova l'anomalia eccentrica.
!
real(kind=q),intent(in)::t,eccentri
real(kind=q),intent(out)::eanomali
real(kind=q) :: tt
integer::j
if(6.0_q>periodo) then
periodo=8.0_q*atan(1.0_q)
allocate(eatempi(0:10001))
eatempi(:)=0.0_q
end if
if(1.0e-8_q>abs(eccentri-ecciusata)) then
ecciusata=eccentri
do j=0,10001
eanomali=(periodo*j)/10000.0_q
eatempi(j)=eanomali-eccentri*sin(eanomali)
end do
end if
tt=modulo(t,periodo)
if(abs(tt)>0.0_q) then
eanomali=tt
ea:do j=0,10001
if(eatempi(j)>tt) then
eanomali=(periodo*j)/10000.0_q
exit ea
end if
end do ea
do j=1,100
eanomali=tt+eccentri*sin(eanomali)
end do
eanomali=t-tt+eanomali
else
eanomali=t
end if
return
end subroutine eqkeplero
!
subroutine stampatab(ecce,fi,posizioni,passo)
!
! Genera i file con i risultati.
!
real(kind=q),intent(in)::ecce
character(len=*),intent(in)::fi
real(kind=q),intent(in),optional::passo
real(kind=q),dimension(0:,1:),intent(in),optional::posizioni
!
integer::j
real(kind=q)::ea
real(kind=q),dimension(5)::yt
character(len=*),parameter::a=char(60),c=char(62)
!
open(unit=10,file=fi,status="replace",action="write")
write(unit=10,fmt="(3a)")a,"html",c
write(unit=10,fmt=*)a,"head",c,a,"title",c,"Moto planetario",a,"/title",c
write(unit=10,fmt=*)a,"/head",c,a,"body style='margin-left:20'",c
write(unit=10,fmt=*)a,"h3",c,&
"Test dei solutori di sistemi di equazioni differenziali ordinarie (ODE)"
write(unit=10,fmt=*)a,"/h3",c
write(unit=10,fmt=*) &
"Si consideri il seguente sistema di quattro equazioni differenziali:"
write(unit=10,fmt=*)a,"br/",c,a,"br/",c
write(unit=10,fmt=*)"`y[1] = y[3]",a,"br/",c
write(unit=10,fmt=*)"`y[2] = y[4]",a,"br/",c
write(unit=10,fmt=*)"`y[3] = −y[1] / ( y[1]",a,"sup",c,"2",a,"/sup",c," +"
write(unit=10,fmt=*)"y[2]",a,"sup",c,"2",a,"/sup",c
write(unit=10,fmt=*)" )",a,"sup",c,"3/2",a,"/sup",c,a,"br/",c
write(unit=10,fmt=*)"`y[4] = −y[2] / ( y[1]",a,"sup",c,"2",a,"/sup",c," +"
write(unit=10,fmt=*)"y[2]",a,"sup",c,"2",a,"/sup",c
write(unit=10,fmt=*)" )",a,"sup",c,"3/2",a,"/sup",c,a,"br/",c
write(unit=10,fmt=*)a,"br/",c
write(unit=10,fmt=*)"Le condizioni iniziali al tempo t = 0, sono:"
write(unit=10,fmt=*)a,"br/",c,a,"br/",c
write(unit=10,fmt=*)" y[1] = 1 − e",a,"br/",c
write(unit=10,fmt=*)" y[2] = 0",a,"br/",c
write(unit=10,fmt=*)" y[3] = 0",a,"br/",c
write(unit=10,fmt=*)" y[4] = ((1 + e) / (1 − e))",a,"sup",c,"1/2",a,"/sup",c
write(unit=10,fmt=*)a,"br/",c,a,"br/",c
!
write(unit=10,fmt=*)"Il semiasse maggiore a vale 1 e"
write(unit=10,fmt=*)"il periodo T vale 2·π",a,"br/",c
write(unit=10,fmt=*)"Eccentricità dell'orbita: e=",ecce,a,"pre",c
if(present(passo)) then
write(unit=10,fmt=*)" "
write(unit=10,fmt=*)"Passo di integrazione numerica: ",passo
write(unit=10,fmt=*)"Vengono riportati i valori di x ed y ottenuti numericamente."
write(unit=10,fmt=*)" "
end if
!
write(unit=10,fmt=*)" n",&
" x ",&
" y "," ",&
" vx ",&
" vy "," ",&
" t "
do j=0,630
call eqkeplero(j*0.01_q,ecce,ea)
yt=pianeta(ecce,ea)
write(unit=10,fmt= "(i6,2f13.9,f15.5,f13.5,f15.9)")j,yt
if(present(posizioni))then
write(unit=10,fmt="(i6,2f13.9,f15.5,f13.5,f15.9)") &
j,posizioni(j,1),posizioni(j,2)
end if
end do
write(unit=10,fmt=*)a,"/pre",c,a,"hr/",c,a,"/body",c
write(unit=10,fmt=*)a,"/html",c
close(unit=10)
return
end subroutine stampatab
!
subroutine orbita_esempio(jec,posizioni,passo,jcoda)
integer,intent(in),optional::jec
character(len=*),intent(in),optional::jcoda
real(kind=q_testpianeta),intent(in),optional::passo
real(kind=q_testpianeta),dimension(0:,1:),intent(in),optional::posizioni
integer::j,kec,nc
real(kind=q_testpianeta)::ecce
character(len=12)::ff,coda
character(len=*),parameter::fi="eccepianeta"
if(present(jcoda)) then
coda=jcoda
else
coda=".htm"
end if
nc=len_trim(coda)
if(present(jec)) then
kec=jec
else
kec=-1
end if
ripeto:do
write(unit=*,fmt=*)" "
if(0>kec) then
if(present(jec)) then
exit ripeto
end if
write(unit=*,fmt=*)"Test pianeta: per verificare i solutori di ODE."
write(unit=*,fmt=*)"Assegna l'eccentricita' ",&
"(da 0 ma minore di 0.999999 altrimenti amen)"
read(unit=*,fmt=*) ecce
else
write(unit=*,fmt=*)"Stampa la soluzione analitica"
ecce=kec*0.001_q_testpianeta
kec=-1
end if
if(0.9999995>=ecce.and.ecce>=0) then
write(unit=ff,fmt="(1f12.9)") ecce
do j=1,12
if("!">ff(j:j)) then
ff(j:j)="_"
end if
end do
if(present(posizioni)) then
if(present(passo)) then
call stampatab(ecce,fi//ff//coda(1:nc),posizioni,passo)
else
call stampatab(ecce,fi//ff//coda(1:nc),posizioni)
end if
else
call stampatab(ecce,fi//ff//coda(1:nc))
end if
write(unit=*,fmt=*)"Guardare : ",fi,ff,coda(1:nc)
else
write(unit=*,fmt=*)"Fuori dai limiti.... termina."
exit ripeto
end if
end do ripeto
write(unit=*,fmt=*)"...un invio... "
read(unit=*,fmt=*)
return
end subroutine orbita_esempio
!
end module testpianeta
!
!
! Commentare, eventualmente, quanto segue.
!
!program t
! use testpianeta
! call orbita_esempio()
! stop
!end program t
!