La teoria
module cashkarp
implicit none
integer,private,parameter::q=selected_real_kind(14,300)
real(kind=q),private,parameter::n11776=11776,n16632=16632,n6237=6237,&
n24948=24948,n1540=1540,n3262=3262,n37800=37800,n4600=4600, &
n44275=44275,n6831=6831,n8855=8855,n9361=9361,n38500=38500,&
n20125=20125,n27648=27648,n10240=10240,n39550=39550,&
n148600=148600,n94675=94675,n7479=7479,n96768=96768,&
n2530=2530,n10=10
real(kind=q),private,parameter::a21=n11776*n16632
real(kind=q),private,parameter::a31=n11776*n6237,a32=3*n11776*n6237
real(kind=q),private,parameter::a41=n11776*n24948,a42=-3*n11776*n24948,&
a43=4*n11776*n24948
real(kind=q),private,parameter::a51=-11*n11776*n1540,a52=135*n11776*n1540,&
a53=-140*n11776*n1540,a54=70*n11776*n1540
real(kind=q),private,parameter::a61=n8855*n3262,a62=n8855*n37800,&
a63=n8855*n4600,a64=n8855*n44275,a65=n8855*n6831
real(kind=q),private,parameter::b51=n10240*n9361,b53=n10240*n38500,&
b54=n10240*n20125,b56=n10240*n27648
real(kind=q),private,parameter::b41=n2530*n39550,b43=n2530*n148600,&
b44=n2530*n94675,b45=n2530*n7479,b46=n2530*n96768
real(kind=q),private,parameter::b31=344565760,b33=-362700800,b34=997427200
real(kind=q),private,parameter::b21=-1468938240,b22=n10*244823040
real(kind=q),private,parameter::b11=979292160
real(kind=q),private,parameter::c2=195858432,c3=293787648,c4=587575296
real(kind=q),private,parameter::c5=979292160,c6=856880640
!
real(kind=q),dimension(:),private,allocatable::yp,ytn ! corrisponde a k1.
real(kind=q),dimension(:),private,allocatable::k2,k3,k4,k5,k6
integer,private::nv=0
private:: esempio_cashkarp
public:: program_test,xcontrollo,cash_karp,hermite,xscashkarp
!
contains
subroutine xcontrollo(tot)
! Se il valore di tot e' zero non ci sono errori evidenti.
real(kind=q),intent(out)::tot
tot=abs(a21-c2)
tot=tot+abs(a31+a32-c3)
tot=tot+abs(a41+a42+a43-c4)
tot=tot+abs(a51+a52+a53+a54-c5)
tot=tot+abs(a61+a62+a63+a64+a65-c6)
tot=tot+abs(b51+b53+b54+b56-b11)
tot=tot+abs(b41+b43+b44+b45+b46-b11)
tot=tot+abs(b31+b33+b34-b11)
tot=tot+abs(b21+b22-b11)
return
end subroutine xcontrollo
!
subroutine xscashkarp(y,nm,not_ok)
real(kind=q),intent(in),dimension(1:)::y
integer,optional,intent(in)::nm
logical,intent(out)::not_ok
not_ok=.false.
if(nv==0) then
nv=size(y,1)
if(present(nm)) then
nv=max(nv,nm)
end if
allocate(ytn(nv),k2(nv),k3(nv),k4(nv))
allocate(k5(nv),k6(nv),yp(nv))
ytn(:)=0.0_q
k2(:)=0.0_q
k3(:)=0.0_q
k4(:)=0.0_q
k5(:)=0.0_q
k6(:)=0.0_q
yp(:)=0.0_q
else if(size(y,1)>size(ytn,1)) then
not_ok=.true.
end if
return
end subroutine xscashkarp
!
subroutine cash_karp(yn,nv,y,hv,f,yt)
!
! nv rappresenta il numero di variabili da integrare e deve
! valere almeno 2 perche' la prima variabile e' quella
! indipendente.
! hv rappresenta l'ampiezza del passo di integrazione
! y(:) rappresenta il vettore delle variabili a inizio passo.
! yn(:) rappresenta il vettore delle variabili a fine passo.
! yt(:) rappresenta il vettore delle derivate prime a inizio passo
! ma e' opzionale perche', se non viene dato, viene
! calcolato internamente.
! ytn(:) rappresenta il vettore delle derivate prime a fine passo, che
! viene calcolato solo se yt(:) e' passato alla subroutine.
! f(z,w) e' la subroutine che calcola il vettore delle derivate prime
! quando è noto il vettore delle variabili.
! yp(:) rappresenta il vettore delle variabili a fine passo ma
! calcolato con una precisione minore. Se yp(:) differisce
! eccessivamente da yn(:) bisogna rifare il passo riducendo
! il valore di hv. Si noti che yp non è in argomento poiche'
! viene data come memoria del modulo.
!
integer,intent(in)::nv
real(kind=q),intent(in),dimension(1:)::y
real(kind=q),intent(in),dimension(1:),optional::yt
real(kind=q),intent(in)::hv
real(kind=q),intent(out),dimension(1:)::yn
interface
subroutine f(z,w)
integer,parameter::q=selected_real_kind(15,300)
real(kind=q),intent(in),dimension(1:)::w
real(kind=q),intent(out),dimension(1:)::z
end subroutine f
end interface
real(kind=q)::h
h=hv/b11
if(present(yt)) then
ytn(1:nv)=yt(1:nv)
else
call f(ytn,y)
end if
yp(1:nv)=y(1:nv)+(h*a21)*ytn(1:nv)
call f(k2,yp)
yp(1:nv)=y(1:nv)+(h*a31)*ytn(1:nv)+(h*a32)*k2(1:nv)
call f(k3,yp)
yp(1:nv)=y(1:nv)+(h*a41)*ytn(1:nv)+(h*a42)*k2(1:nv)+&
(h*a43)*k3(1:nv)
call f(k4,yp)
yp(1:nv)=y(1:nv)+(h*a51)*ytn(1:nv)+(h*a52)*k2(1:nv)+&
(h*a53)*k3(1:nv)+(h*a54)*k4(1:nv)
call f(k5,yp)
yp(1:nv)=y(1:nv)+(h*a61)*ytn(1:nv)+(h*a62)*k2(1:nv)+&
(h*a63)*k3(1:nv)+(h*a64)*k4(1:nv)+(h*a65)*k5(1:nv)
call f(k6,yp)
yn=y ! per trascrivere eventuali dati oltre i primi nv.
yn(1:nv)=yn(1:nv)+(h*b51)*ytn(1:nv)+(h*b53)*k3(1:nv)+&
(h*b54)*k4(1:nv)+(h*b56)*k6(1:nv)
! Valore meno preciso usabile solo per controllo:
yp(1:nv)=y(1:nv)+(h*b41)*ytn(1:nv)+(h*b43)*k3(1:nv)+&
(h*b44)*k4(1:nv)+(h*b45)*k5(1:nv)+(h*b46)*k6(1:nv)
if(present(yt)) then
call f(ytn,yn)
end if
return
end subroutine cash_karp
!
subroutine hermite(yora,y,yt,yn,ytn,theta,n)
!
! interpolazione cubica (di Hermite)
!
integer,intent(in)::n
real(kind=q),intent(in)::theta
real(kind=q),intent(in),dimension(1:)::y,yt,yn,ytn
real(kind=q),intent(out),dimension(1:)::yora
real(kind=q):: h,d1,d2,d3,d4,ds,dsq,tq
!
! La variabile indipendente deve essere la prima.
!
h=yn(1)-y(1)
tq=theta*theta
ds=theta-1.0_q
dsq=ds*ds
d1=dsq*(2.0_q*theta+1.0_q)
d2=dsq*theta*h
d3=tq*(3.0_q-2.0_q*theta)
d4=tq*ds*h
yora(1:n)=d1*y(1:n)+d2*yt(1:n)+d3*yn(1:n)+d4*ytn(1:n)
return
end subroutine hermite
!
subroutine esempio_cashkarp(z,w)
!
! Il tempo e' ottenuto risolvendo l'equazione differenziale:
! `t = 1 e il valore del tempo e' memorizzato nella prima
! componente di w. Il vettore z rappresenta il valore
! delle derivate prime corrispondenti ai valori delle
! variabili memorizzate nel vettore w.
!
real(kind=q),intent(in),dimension(1:)::w
real(kind=q),intent(out),dimension(1:)::z
real(kind=q)::r32
!
if(w(1)>=0.0_q) then
z(1)=1 ! `t
z(2)=w(4) ! `x
z(3)=w(5) ! `y
r32=w(2)*w(2)+w(3)*w(3)
r32=1.0_q/(r32*sqrt(r32))
z(4)=-w(2)*r32 ! `vx
z(5)=-w(3)*r32 ! `vy
else
!
! Per poter inizializzare il calcolo , w(1) deve essere
! negativo e il vettore z(:) rappresenta allora il vettore
! delle variabili e non il vettore delle loro derivate prime.
! Al tempo w(1).lt.0 il valore assoluto di w(1) indica
! l'eccentricita' che deve essere maggiore di 0 e inferiore
! ad 1.
! Non e' ammessa una eccentricita' esattamente nulla perche'
! allora w(1) non sarebbe negativo ma qualsiasi valore pur
! piccolissimo va bene...
!
r32=min(0.999999_q,abs(w(1)))
z(1)=0.0_q ! t=0
z(2)=1-r32 ! 1-e
z(3)=0.0_q
z(4)=0.0_q
z(5)=sqrt((1.0_q+r32)/(1.0_q-r32)) ! sqrt((1+e)/(1-e))
write(unit=*,fmt="(1a,4f18.10)") " ",z(2:)
end if
return
!
! Questo esempio e' il problema del moto di un pianeta attorno al sole
! che occupa l'origine delle coordinate. Il semiasse maggiore e' sempre
! uno, se il moto e' ellittico ossia se l'eccentricita', mai negativa,
! risulta essere inferiore ad 1. Se valesse esattamente 1 il moto
! sarebbe parabolico. Se e indica l'eccentricita', allora a t=w(1)=0
! bisogna imporre queste condizioni iniziali: z(1)=0, z(2)=1-e, z(3)=0,
! z(4)=0, z(5)=sqrt((1+e)/(1-e)) ossia il moto inizia al perielio.
!
end subroutine esempio_cashkarp
!
subroutine program_test(ece,passie,girie,posizioni)
real(kind=q),intent(in),optional::ece
integer,intent(in),optional::passie,girie
real(kind=q),dimension(0:,1:),intent(out),optional::posizioni
integer::passi,giri,j,kp
real(kind=q)::ec,hv,t0,s,x
real(kind=q),dimension(5)::y,yn,yt,yora
logical::boh
!----
write(unit=*,fmt=*)"Esegue alcune prove di funzionamento",&
" di questo modulo (module cashkarp)"
write(unit=*,fmt=*)" "
call xcontrollo(x)
write(unit=*,fmt=*)"1) Se vale 0.0 ",&
"i coefficienti sono forse giusti : ",x
if(present(ece)) then
ec=ece
else
write(unit=*,fmt="(a)",advance="no") &
" Assegna l'eccentricita' dell'orbita "
read(unit=*,fmt=*) ec
ec=max(1.0e-30_q,min(0.999999_q,ec))
end if
if(present(passie)) then
passi=passie
else
write(unit=*,fmt="(a)",advance="no") &
" Assegna il numero di passi"
read(unit=*,fmt=*) passi
end if
if(present(girie)) then
giri=girie
else
write(unit=*,fmt="(a)",advance="no") &
" Assegna il numero di giri"
read(unit=*,fmt=*) giri
end if
!
if(present(posizioni)) then
posizioni=0.0_q
end if
!
call xscashkarp(y,5,boh)
if(boh) then
write(unit=*,fmt=*)"Fallisce l'inizializzazione di cashkarp: invia e stop!"
read(unit=*,fmt=*)
stop
end if
! hv*passi deve corrispondere ad un dodicesimo di giro ossia π/6
hv=2.0_q*atan(1.0_q)/(3.0_q*max(1,passi))
yn(:)=0
yn(1)=-ec
call esempio_cashkarp(y,yn)
! Ogni j=passi dovrebbe aver fatto un dodicesimo di giro.
do kp=1,giri*12-12
do j=1,passi
call cash_karp(yn,5,y,hv,esempio_cashkarp)
y=yn
end do
end do
do kp=1,12
do j=1,passi
call cash_karp(yn,5,y,hv,esempio_cashkarp)
y=yn
end do
write(unit=*,fmt="(i5,4f18.10)") kp,y(2:5)
end do
write(unit=*,fmt=*)"Ora il tempo vale: ",y(1)
t0=giri*atan(1.0_q)*8.0_q
write(unit=*,fmt=*)"Il tempo voluto sarebbe:",t0
read(unit=*,fmt=*)
!
if(.not.present(posizioni)) then
return
end if
!
! Esegue questi calcoli solo se deve fare una tabella con
! le posizioni distanziate un centesimo di tempo ovvero
! un centesimo di radiante, dato che fa un giro in 6.28318....
!
call esempio_cashkarp(yt,y)
call cash_karp(yn,5,y,hv,esempio_cashkarp,yt)
do kp=0,630
procedo:do
if(yn(1)>t0) then
exit procedo
end if
y=yn
yt(1:5)=ytn(1:5)
call cash_karp(yn,5,y,hv,esempio_cashkarp,yt)
end do procedo
s=(t0-y(1))/(yn(1)-y(1))
call hermite(yora,y,yt,yn,ytn,s,5)
posizioni(kp,1)=yora(2)
posizioni(kp,2)=yora(3)
! read(unit=*,fmt=*)
t0=t0+0.01_q
end do
return
end subroutine program_test
!
end module cashkarp
!
!Per verificare scommenta il program seguente:
!
!program xcashkarp
! use cashkarp
! implicit none
! call program_test()
! write(unit=*,fmt="(1x,1a)",advance="no") "Un invio per finire: "
! read(unit=*,fmt=*)
! stop
!end program xcashkarp
!