#!/usr/bin/gnuplot # vsop82.plt # Analyticka teorie pohybu Zeme VSOP82 (podle Wolf aj. 1992) # a porovnani s presnou numerickou integraci. # Miroslav Broz (miroslav.broz@email.cz), Feb 7th 2007 # stredni delky planet l1(T) = 4.4026 + 2608.7903*T # pro Merkur l2(T) = 3.1761 + 1021.3286*T l3(T) = 1.7535 + 628.3076*T # <-- tj. priblizne 100 * 2pi l4(T) = 6.2035 + 334.0612*T l5(T) = 0.5995 + 52.9691*T # pro Jupiter # STREDNI nesingularni elementy pro teziste soustavy Zeme-Mesic a0(T) = 1.0000010 L0(T) = 1.7534703 + 628.3075849*T - 0.0000001*T**2 # tj. jen presnejsi verze l3(T) k0(T) = -0.0037408 - 0.0000823*T + 0.0000003*T**2 h0(T) = 0.0162845 - 0.0000620*T - 0.0000003*T**2 q0(T) = -0.0001135*T + 0.0000001*T**2 p0(T) = 0.0000102*T + 0.0000005*T**2 # PERTURBACE (pouze NEKOLIK NEJVETSICH clenu!!!) a(T) = a0(T) + 1e-7 * (\ \ + 112*cos(2*(l3(T)-l5(T))) \ + 76*cos(l2(T)-l3(T)) \ - 41*cos(2*(l2(T)-l3(T))) \ - 25*cos(3*(l2(T)-l3(T))) \ + 15*sin(2*l2(T)-3*l3(T)) \ + 11*cos(l3(T)-l5(T)) \ ) L(T) = L0(T) + 1e-7 * (\ + 322*cos(4*l3(T)-8*l4(T)+3*l5(T)) \ - 206*sin(2*(l3(T)-l5(T))) \ + 166*sin(l2(T)-l3(T)) \ ) k(T) = k0(T) + 1e-7 * (\ - 199*cos(2*l2(T)-3*l3(T)) \ + 186*cos(l3(T)-2*l5(T)) \ - 150*cos(l5(T)) \ ) h(T) = h0(T) + 1e-7 * (\ + 199*sin(2*l2(T)-3*l3(T)) \ - 186*sin(l3(T)-2*l5(T)) \ - 151*sin(l5(T)) \ ) p(T) = p0(T) q(T) = q0(T) # Keplerovy elementy e(T) = sqrt(k(T)**2+h(T)**2) varpi(T) = atan2(h(T),k(T)) I(T) = 2*asin(sqrt(q(T)**2+p(T)**2)) Omega(T) = atan2(p(T),q(T)) M(T) = L(T)-varpi(T) ######################################################################## set xl "cas T / Julianska stoleti" set yl "velka poloosa a / AU" set y2l "excentricita e" set xr [0:0.1] set xr [0:1] set yr [0.999975:1.00009] set y2r [0.01645:0.01677] #set xtics 0.01 set ytics nomirror set y2tics set samples 1000 set style line 1 lt 1 lw 1 set style line 2 lt 3 lw 1 set style line 3 lt 7 lw 2 set style line 4 lt 5 lw 2 p "follow.out_-4" u ($2*1e4):3 t "a(T) integrace" w l ls 3,\ "follow.out_-4" u ($2*1e4):4 ax x1y2 t "e(T) integrace" w l ls 4,\ a(x) t "a(T) VSOP82" ls 1,\ e(x) ax x1y2 t "e(T) VSOP82" ls 2 pa -1 set term png small set out "vsop82.png" rep q