!Orbita di pianeta con semiasse maggiore 1

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 !