29 double improve_mean_to_eccentric_annomaly_guess(
const double e,
const double M,
const double x);
31 void calc_cartesian_for_ellipse(
double&
x,
double& y,
double & z,
double &vx,
double &vy,
double &vz,
const double a,
const double e,
const double i,
const double O,
const double w,
const double M,
const double GM);
32 void calc_keplerian_for_cartesian(
double& a,
double& e,
double& i,
double& O,
double& w,
double& M,
const double x,
const double y,
const double z,
const double vx,
const double vy,
const double vz,
const double GM);
36 double improve_mean_to_eccentric_annomaly_guess(
const double e,
const double M,
const double x)
50 dx = -f/(fp+dx*fpp/2.);
51 dx = -f/(fp+dx*fpp/2.+dx*dx*fppp/6.);
59 const int ORBEL_EHIE_NMAX = 3;
61 int nper =
static_cast<int>(M/(2.*M_PI));
63 if(M<0.) M = M + 2.*M_PI;
68 double x = pow(6.*M,1./3.);
69 for(
int i=1;i<=ORBEL_EHIE_NMAX;++i)
70 x = improve_mean_to_eccentric_annomaly_guess(e,M,x);
76 double x = pow(6.*M,1./3.);
77 for(
int i=1;i<=ORBEL_EHIE_NMAX;++i)
78 x = improve_mean_to_eccentric_annomaly_guess(e,M,x);
84 void calc_cartesian_for_ellipse(
double&
x,
double& y,
double & z,
double &vx,
double &vy,
double &vz,
const double a,
const double e,
const double i,
const double O,
const double w,
const double M,
const double GM)
88 sincos(cape,&scap,&ccap);
89 double sqe = sqrt(1.-e*e);
90 double sqgma = sqrt(GM*a);
91 double xfac1 = a*(ccap-e);
92 double xfac2 = a*sqe*scap;
93 double ri = 1./(a*(1.-e*ccap));
94 double vfac1 = -ri*sqgma * scap;
95 double vfac2 = ri*sqgma*sqe*ccap;
97 double sw, cw, so, co, si, ci;
101 double d1[] = { cw*co-sw*so*ci, cw*so+sw*co*ci, sw*si};
102 double d2[] = {-sw*co-cw*so*ci,-sw*so+cw*co*ci, cw*si};
103 x = d1[0]*xfac1+d2[0]*xfac2;
104 y = d1[1]*xfac1+d2[1]*xfac2;
105 z = d1[2]*xfac1+d2[2]*xfac2;
106 vx = d1[0]*vfac1+d2[0]*vfac2;
107 vy = d1[1]*vfac1+d2[1]*vfac2;
108 vz = d1[2]*vfac1+d2[2]*vfac2;
110 void calc_keplerian_for_cartesian(
double& a,
double& e,
double& i,
double& O,
double& w,
double& M,
const double x,
const double y,
const double z,
const double vx,
const double vy,
const double vz,
const double GM)
112 const double TINY = 1.e-8;
114 double h[] = {y*vz-z*vy, z*vx-x*vz, x*vy-y*vx};
115 double h2 = h[0]*h[0]+h[1]*h[1]+h[2]*h[2];
116 double hh = sqrt(h2);
118 double fac = sqrt(h[0]*h[0]+h[1]*h[1])/hh;
124 if(fabs(i-M_PI)<10.*TINY) u = -u;
128 O = atan2(h[0],-h[1]);
129 u = atan2(z/sin(i),x*cos(O)+y*sin(O));
131 if(O<0.) O += 2.*M_PI;
132 if(u<0.) u += 2.*M_PI;
133 double r = sqrt(x*x+y*y+z*z);
134 double energy = (vx*vx+vy*vy+vz*vz)*0.5-GM/r;
136 if(fabs(energy*r/GM)<sqrt(TINY))
140 double ww = acos(2.*a/r-1.);
141 if(vx*x+vy*y+vz*z<0.) w = 2.*M_PI-w;
142 double tmpf = tan(0.5*w);
143 M = tmpf*(1.+tmpf*tmpf/3.);
145 if(w<0.) w+= 2.*M_PI;
146 w -= round(w/(2.*M_PI))*2.*M_PI;
156 double face = (a-r)/(a*e);
157 if(face>1.) cape = 0.;
158 else if (face>-1.) cape = acos(face);
161 if(vx*x+vy*y+vz*z<0.) cape = 2.*M_PI-cape;
162 double cw = (cos(cape)-e)/(1.-e*cos(cape));
163 double sw = sqrt(1.-e*e)*sin(cape)/(1.-e*cos(cape));
165 if(ww<0.) ww += 2.*M_PI;
173 M = cape - e*sin(cape);
175 if(w<0.) w += 2.*M_PI;
176 w -= round(w/(2.*M_PI))*2.*M_PI;
186 double tmpf = (a+r)/(a*e);
187 capf = log(tmpf+sqrt(tmpf*tmpf-1.));
188 if(vx*x+vy*y+vz*z<0.) capf = -capf;
189 double cw = (e-cosh(capf))/(e*cosh(capf)-1.);
190 double sw = sqrt(e*e-1.)*sinh(capf)/(e*cosh(capf)-1.);
192 if(ww<0.) ww += 2.*M_PI;
197 double tmpf = 0.5*h2/GM;
198 ww = acos(2.*tmpf/r-1.);
199 if(vx*x+vy*y+vz*z<0.) ww = 2.*M_PI-ww;
201 capf = log(tmpf+sqrt(tmpf*tmpf-1.));
203 M = e*sinh(capf)-capf;
206 w -= round(w/(2.*M_PI))*2.*M_PI;