#!/usr/bin/gnuplot # slapy1.plt # Disipace energie pri slapovem pusobeni v soustave Zeme-Mesic. # (viz Bertotti et al. (2003), str. 488) # Miroslav Broz (miroslav.broz@email.cz), Nov 13th 2007 G = 6.7e-11 # N/kg^2*m^2; Newtonova gravitacni konstanta M_p = 6.e24 # kg; hmotnost Zeme M_s = 7.4e22 # kg; hmotnost Mesice I_p = 8.e37 # kg*m^2; moment setrvacnosti Zeme (Stacey 1977) omega_p = 2.*pi/((23.+56./60.)*3600.) # rad/s; uhlova rychlost rotace Zeme omega_s = 2.*pi/(27.3*24.*3600.) # totez pro Mesic R_p = 6378.e3 # m; polomer Zeme R_s = 1737.e3 # m; polomer Mesice r = 384000.e3 # m; stredni vzdalenost Mesice od Zeme I_p_koule = 2./5.*M_p*R_p**2. # kg*m^2; moment setrvacnosti Zeme jako homogenni koule I_s_koule = 2./5.*M_s*R_s**2. # kg*m^2; totez pro Mesic L_rot_p = I_p*omega_p # rotacni moment hybnosti Zeme L_rot_s = I_s_koule*omega_s # totez pro Mesic l(r) = sqrt(G)*M_p*M_s/sqrt(M_p+M_s) * sqrt(abs(r)) * sgn(r) # orbitalni moment hybnosti r(l) = (abs(l) / (sqrt(G)*M_p*M_s/sqrt(M_p+M_s)))**2. * sgn(l) # inverzni vztah omega(r) = (L-l(r))/I_p # uhlova rychlost Zeme (vypoctena z L) E(r) = 0.5*I_p*omega(r)**2. - G*M_p*M_s/(2.*abs(r)) # celkova energie n(r) = sqrt(G*(M_p+M_s)/r**3.) # stredni pohyb Mesice L = L_rot_p + l(r) # kg*m^2/s ; soucasna hodnota CELKOVEHO momentu hybnosti soustavy Zeme-Mesic ######################################################################## print "soucasna uhlova rychlost rotace Zeme a jeji moment setrvacnosti:" print "omega_p = ", omega_p, " rad/s" print "I_p = ", I_p, " kg*m^2" print "I_p_koule = ", I_p_koule, " kg*m^2" print "omega_s = ", omega_s, " rad/s" print "I_s_koule = ", I_s_koule, " kg*m^2" print "\nrotacni a orbitalni moment hybnosti Zeme, resp. soustavy Zeme-Mesic:" print "L_rot_p = ", L_rot_p, " kg*m^2/s" print "L_rot_s = ", L_rot_s, " kg*m^2/s = ", L_rot_s/L_rot_p, " L_rot_p" print "l = ", l(r), " kg*m^2/s = ", l(r)/L_rot_p, " L_rot_p" print "L_celkovy = ", L, " kg*m^2/s" ######################################################################## set xl "vzdalenost Zeme-Mesic r / km" set x2l "vzdalenost Zeme-Mesic r / R_Zeme" set yl "orbitalni moment hybnosti l / kg*m^2/s" set samples 1000 km = 1.e3 tmp = 1.67*r x1=-tmp/km x2=tmp/km set xr [x1:x2] set x2r [-tmp/R_p:tmp/R_p] set xtics nomirror set x2tics y2=1.2*L y1=-y2 set yr [y1:y2] lw=3 set style line 1 lt 1 lw lw set style line 2 lt 3 lw lw set style line 3 lt 7 lw lw set style line 4 lt 4 set style line 5 lt 5 set label "{/Helvetica-Oblique r} = 60 {/Helvetica-Oblique R}_{Z}" at r/km, graph 1.03 set label "{/Helvetica-Oblique L}" at graph 1.02, first L p l(x*km) not ls 1,\ "vline.dat" u ($1*r/km):(y1+(y2-y1)*$2) w l not ls 4,\ L not ls 5 #pa -1 call "png.plt" "slapy0.png" ######################################################################## set yl "rotacni a orbitalni frekvence omega, n / rad/s" set y2l "celkova energie E / J" y1=-1.e-3 y2=1.e-3 set yr [y1:y2] set ytics nomirror y21=-6.e30 y22=3.e31 set y2r [y21:y22] set y2tics set nolabel 2 set key bottom left width -18 p omega(x*km) t "{/Symbol w}" ls 1,\ n(x*km) t "{/Helvetica-Oblique n}" ls 2,\ "vline.dat" u ($1*r/km):(y1+(y2-y1)*$2) w l not ls 4,\ E(x*km) ax x1y2 t "{/Helvetica-Oblique E}" ls 3 #pa -1 call "png.plt" "slapy1.png" ######################################################################## set encoding iso_8859_2 lw=3 set style line 1 lt 2 lw lw set style line 2 lt 4 lw lw set style line 3 lt 1 lw lw set style line 4 lt 1 lw 1 set style line 5 lt 1 lw lw set style line 1 lt 7 lw 1 set style line 2 lt 8 lw 4 set style line 3 lt 7 lw 4 set style line 4 lt 0 lw 1 set style line 5 lt 1 lw 5 set lmargin 10.0 set rmargin 11.0 set tmargin 3.5 set bmargin 3.0 set nokey set xl "vzdálenost Země-Měsíc {/Helvetica-Oblique r} / km" set x2l "vzdálenost Země-Měsíc {/Helvetica-Oblique r} / R_{Země}" set yl "rotační nebo orbitální frekvence {/Symbol w}, {/Helvetica-Oblique n} / rad/s" offset +1,0 set y2l "celková energie {/Helvetica-Oblique E} / J" offset -1,0 print "R_p = ", R_p, " m" x3=87.5*R_p/km y3=-4e-5 x4=x3 y4=8.e29 set arrow from x3,y3 to x4,second y4 heads size 0.02*(x2-x1),30 filled ls 5 set label " konečný stav" at x4-0.03*(x2-x1),second y4 rotate x5=60.*R_p/km y5=2.e31 set label " dnešní stav" at x5-0.015*(x2-x1),second y5 rotate font "Helvetica,10" x6=2.28*R_p/km y6=4.27e30+0.01*(y22-y21) x7=x6+0.04*(x2-x1) y7=y6+0.08*(y22-y21) set arrow from x7,second y7 to x6,second y6 ls 1 set label " geostacionární" at x7,second y7+0.05*(y22-y21) font "Helvetica,10" center set label " dráha" at x7,second y7+0.02*(y22-y21) font "Helvetica,10" center x8=-60*R_p/km y8=2.4e31 y8=2e31 set label "retrográdní pohyb" at x8,second y8 rotate by -38.5 font "Helvetica,10" call "eps.plt" "slapy1.eps" ######################################################################## set xl "orbitalni moment hybnosti l / kg*m^2/s" set x2l "orbitalni moment hybnosti l / L" x1=-L x2=L set xr [x1:x2] set x2r [x1/L:x2/L] set nolabel set label "l" at first l(r), graph 1.03 center p omega(r(x)) ls 1,\ n(r(x)) ls 2,\ E(r(x)) ax x1y2 ls 3,\ "vline.dat" u ($1*l(r)):(y1+(y2-y1)*$2) w l not ls 5 #pa -1 call "png.plt" "slapy2.png" q