Код:
#restart;read("P:\\src\\math\\earth\\s1.m12");
with(ScientificConstants):GetConstants():
with(plots):
with(MTM):
GM:=evalf(Constant(G) * Constant(M[Sun])):
entry_km_day_kg:=proc(title,km,day,kg)
local x;
global GM;
x:=struct('R',km*1000,'m',kg,'T',day*24*3600,'M',0,'a',0,'w',0);
x:-w:=evalf(2*Pi)/x:-T; # radian/s
x:-M:=x:-m*(x:-R)^2*x:-w; # moment
x:-a:= GM*x:-m;
return x;
end proc:
entry_ae2_day_kg:=proc(title,aep,aea,day,kg) # perigelii, afelii, ? big poluos
return entry_km_day_kg(title,149597871*(aep+aea)/2,day,kg);
end proc:
Sun:=entry_km_day_kg("Sun",2.955036102,27.35, 1.988e30):
A:=[]:
xadd:=proc(x) global A; A:=[op(A),x]; end proc:
xadd(entry_km_day_kg("Merk",58e6, 88, 3.3e23)):
xadd(entry_km_day_kg("Venera",108e6, 225, 4.9e24)):
xadd(entry_km_day_kg("Earth",150e6, 365, 6e24)): # 3
xadd(entry_km_day_kg("Mars",228e6, 687, 6.4e23)):
J:=entry_km_day_kg("Jupiter",778e6, 12*365, 1.9e27): xadd(J): # 5
S:=entry_km_day_kg("Saturn",1429e6, 29*365, 5.7e26): xadd(S): # 6
xadd(entry_km_day_kg("Uran",2871e6, 84*365, 8.7e25)):
xadd(entry_km_day_kg("Neptun",4504e6, 165*365, 1e26)):
xadd(entry_ae2_day_kg("Pluton",29.6,49.3,90613.3, 1.3e22)):
#10
xadd(entry_ae2_day_kg("Pallada",2.132, 3.144, 1686.643, 2.06e20)): # http://ru.wikipedia.org/wiki/(2)_%D0%9F%D0%B0%D0%BB%D0%BB%D0%B0%D0%B4%D0%B0
xadd(entry_ae2_day_kg("Vesta",2.151, 2.571, 1326.081, 2.59e20)): # http://ru.wikipedia.org/wiki/(4)_%D0%92%D0%B5%D1%81%D1%82%D0%B0
xadd(entry_ae2_day_kg("Сerera",2.5465, 2.9842, 1680.5, 9.43e20)): # http://ru.wikipedia.org/wiki/%D0%A6%D0%B5%D1%80%D0%B5%D1%80%D0%B0
xadd(entry_ae2_day_kg("Gigey",2.772, 3.506, 2031.353, 9.03e20)): # http://ru.wikipedia.org/wiki/(10)_%D0%93%D0%B8%D0%B3%D0%B5%D1%8F
ug:=proc(i,j,r)
local t;
global A;
t := abs(A[j]:-R - r);
return -(Constant(G)*A[i]:-m*A[j]:-m)/t;
end proc:
U:=proc(i,r) # LL1,15.2
local u,t,m;
global A,Sun,J,S;
u := (A[i]:-M)^2/(2*A[i]:-m*r^2); # 1/r^2
u := u + (Constant(G)*A[i]:-m*Sun:-m)/(-r);
u := u + ug(i,5,r); # +Uran
u := u + ug(i,6,r);
m := 2e34;
return piecewise(u>m,m,u<-m,-m,u);
end proc:
U1:=proc(i,r)
local u,m;
u := U(i,r);
m := 2e29;
return piecewise(u>m,m,u<-m,-m,u); # nice plot
end proc:
a:=0.5e11:
b:=1.5e12:
plot([U(2,r),U(3,r),U(4,r)],r=a..b,color = ["green","blue","red","orange"]);
plot([U1(10,r),U1(11,r),U1(12,r),U1(13,r)],r=a..b,color = ["green","blue","red","orange"]);
#plot([U(2,r),U(3,r),U(4,r),U1(10,r),U1(11,r),U1(12,r),U1(13,r)],r=a..b,color = ["green","blue","red", "orange","green","blue","red","orange"]);