ВНИМАНИЕ! На форуме начался конкурс - астрофотография месяца СЕНТЯБРЬ!
0 Пользователей и 1 Гость просматривают эту тему.
unit modDetOr;interfaceuses Math;// в систеие СИ, углы - радианыconst GM_EA=3.9860043289694E+14; type tVector=record X,Y,Z:Double;end;procedure deterOrbit(vP,vV:tVector;var a,e,i,node,w,M:Double); // e<1implementation//========function modVector(v:tVector):Double; // = |v|begin Result:=Sqrt(v.X*v.X+v.Y*v.Y+v.Z*v.Z);end;function proScal(v1,v2:tVector):Double; // =v1*v2 скалярноеbegin Result:=v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;end;function to2PI(x:Double):Double; // угол x -> 0..2пиbegin if x<0 then Result:=x+2*PI else if x>2*PI then Result:=x-2*PI else Result:=x;end;procedure deterOrbit(vP,vV:tVector;var a,e,i,node,w,M:Double);var v0,r0,v2c,v0v2c,cosFi,e2,v_p,p,nu,EE:Double;begin v0:=modVector(vV); r0:=modVector(vP); v2c:=GM_EA/r0; v0v2c:=v0*v0/v2c; a:=r0/(2-v0v2c); v_p:=proScal(vP,vV); cosFi:=v_p/(r0*v0); e2:=(v0v2c-1)*(v0v2c-1)+(r0/a)*v0v2c*cosFi*cosFi; e:=Sqrt(e2); p:=a*(1-e2); nu:=ArcTan2(v_p*Sqrt(p/GM_EA),p-r0); i:=ArcCos((vP.X*vV.Y-vP.Y*vV.X)/Sqrt(p*GM_EA)); node:=to2PI(ArcTan2(vP.Y*vV.Z-vP.Z*vV.Y,vP.X*vV.Z-vP.Z*vV.X)); w:=to2PI(ArcTan2(vP.Z*Cosecant(i),vP.X*Cos(node)+vP.Y*Sin(node))-nu); EE:=2*ArcTan(Tan(nu/2)*Sqrt((1-e)/(1+e))); M:=to2PI(EE-e*Sin(EE));end;end.
procedure deterOrbitAll(vK,vV:tVector;var a,e,i,node,w,M:Double);var r0,V0,r0V0,t,E0,pp,eSinE:Double;begin r0:=ModVector(vK); V0:=ModVector(vV); t:=2/r0-V0*V0/GM_EA; a:=1/t; r0V0:=proScal(vK,vV); eSinE:=r0V0/(Sqrt(GM_EA*Abs(a))); t:=1-r0/a; if a>0 then begin // эллипс E0:=ArcTan2(eSinE,t); e:=t/Cos(E0); M:=E0-e*Sin(E0); if M<0 then M:=M+2*PI; end else begin // гипербола E0:=ArcTanH(eSinE/t); e:=t/CosH(E0); M:=-E0+e*SinH(E0); end; pp:=a*(1-e*e); t:=(vK.X*vV.Y-vK.Y*vV.X)/Sqrt(pp*GM_EA); i:=ArcCos(t); Node:=ArcTan2(vK.Y*vV.Z-vK.Z*vV.Y,vK.X*vV.Z-vK.Z*vV.X); if Node<0 then Node:=Node+2*PI; w:=ArcTan2(vK.Z/Sin(i),vK.X*Cos(Node)+vK.Y*Sin(Node))-ArcTan2(Sqrt(pp)*r0V0,Sqrt(GM_EA)*(pp-r0)); if w<0 then w:=w+2*PI; if w>(2*PI) then w:=w-2*PI;end;
// a - большая полуось// e - эксцентриситет// i - наклонение орбиты// node - долгота восходящего узла// w - аргумент перицентра// M - средняя аномалияvoid allObdit(const vec3d &P, const vec3d &V, double GM, double &a, double &e, double &i, double &node, double &w, double &M){ static const double _2PI = 6.283185307179586476925286766559; double v0 = V.length(); double r0 = P.length(); double t = 2 / r0 - v0*v0 / GM; a = 1 / t; double r0v0 = P * V; double eSinE = r0v0 / sqrt(GM*abs(a)); t = 1 - r0 / a; double E0; if (a > 0) // эллипс { E0 = atan2(eSinE, t); e = t / cos(E0); M = E0 - e*sin(E0); if (M < 0) M = M + _2PI; } else // гипербола { E0 = atanh(eSinE / t); e = t / cosh(E0); M = -E0 + e*sinh(E0); } double pp = a*(1 - e*e); t = (P.x()*V.y() - P.y()*V.x()) / sqrt(pp*GM); i = acos(t); node = atan2(P.y()*V.z() - P.z()*V.y(), P.x()*V.z() - P.z()*V.x()); if (node<0) node = node + _2PI; w = atan2(P.z() / sin(i), P.x()*cos(node) + P.y()*sin(node)) - atan2(sqrt(pp)*r0v0, sqrt(GM)*(pp - r0)); if (w<0) w = w + _2PI; if (w>_2PI) w = w - _2PI;}
Паскаль... Кровь из глаз
Вот, если надо , код
Я до сих пор с содроганием впоминаю конструкции С++ вида а+++b
Это будет работать если, объект будет перемещаться по орбите и вектор скорости так же будет меняться?
Код: [Выделить]procedure deterOrbitAll(vK,vV:tVector;var a,e,i,node,w,M:Double);
procedure deterOrbitAll(vK,vV:tVector;var a,e,i,node,w,M:Double);
Это вектор положения относительно центра притяжения или центра эллипса?