00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_orbit.h"
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif // HAVE_CONFIG_H
00030 #include "support.h"
00031
00032 #include "orsa_units.h"
00033 #include "orsa_common.h"
00034 #include "orsa_universe.h"
00035 #include "orsa_file.h"
00036 #include "orsa_frame.h"
00037 #include "orsa_secure_math.h"
00038
00039 #ifndef _WIN32
00040 #include <sys/time.h>
00041 #endif
00042 #include <time.h>
00043
00044 #include <cstring>
00045 #include <limits>
00046 #include <algorithm>
00047
00048 using namespace std;
00049
00050 namespace orsa {
00051
00052 double Orbit::GetE() const {
00053
00054 if (e >= 1.0) {
00055 ORSA_WARNING("orsa::Orbit::GetE() called with eccentricity = %g; returning M.",e);
00056
00057 return M;
00058 }
00059
00060 double E = 0.0;
00061 if (e < 0.8) {
00062
00063 #ifdef HAVE_SINCOS
00064 double s,c;
00065
00066 ::sincos(M,&s,&c);
00067 const double sm = s;
00068 const double cm = c;
00069 #else // HAVE_SINCOS
00070 const double sm = sin(M);
00071 const double cm = cos(M);
00072 #endif // HAVE_SINCOS
00073
00074
00075 double x = M + e*sm*( 1.0 + e*( cm + e*( 1.0 -1.5*sm*sm)));
00076
00077 double sx,cx;
00078 E = x;
00079 double old_E;
00080 double es,ec,f,fp,fpp,fppp,dx;
00081
00082 unsigned int count = 0;
00083 const unsigned int max_count = 128;
00084 do {
00085
00086 #ifdef HAVE_SINCOS
00087 ::sincos(x,&sx,&cx);
00088 #else // HAVE_SINCOS
00089 sx = sin(x);
00090 cx = cos(x);
00091 #endif // HAVE_SINCOS
00092
00093 es = e*sx;
00094 ec = e*cx;
00095 f = x - es - M;
00096 fp = 1.0 - ec;
00097 fpp = es;
00098 fppp = ec;
00099 dx = -f/fp;
00100 dx = -f/(fp + dx*fpp/2.0);
00101 dx = -f/(fp + dx*fpp/2.0 + dx*dx*fppp/6.0);
00102
00103 old_E = E;
00104 E = x + dx;
00105 ++count;
00106
00107 x = E;
00108
00109 } while ((fabs(E-old_E) > 10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon()) && (count < max_count));
00110
00111 if (count >= max_count) {
00112 ORSA_ERROR("Orbit::GetE(): max count reached, e = %g E = %g fabs(E-old_E) = %g 10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon() = %g",e,E,fabs(E-old_E),10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon());
00113 }
00114
00115 } else {
00116
00117 double m = fmod(10*twopi+fmod(M,twopi),twopi);
00118 bool iflag = false;
00119 if (m > pi) {
00120 m = twopi - m;
00121 iflag = true;
00122 }
00123
00124
00125 double x = secure_pow(6.0*m,1.0/3.0) - m;
00126 E = x;
00127 double old_E;
00128 double sa,ca,esa,eca,f,fp,dx;
00129
00130 unsigned int count = 0;
00131 const unsigned int max_count = 128;
00132 do {
00133
00134 #ifdef HAVE_SINCOS
00135 ::sincos(x+m,&sa,&ca);
00136 #else // HAVE_SINCOS
00137 sa = sin(x+m);
00138 ca = cos(x+m);
00139 #endif // HAVE_SINCOS
00140
00141 esa = e*sa;
00142 eca = e*ca;
00143 f = x - esa;
00144 fp = 1.0 - eca;
00145 dx = -f/fp;
00146 dx = -f/(fp + 0.5*dx*esa);
00147 dx = -f/(fp + 0.5*dx*(esa+1.0/3.0*eca*dx));
00148 x += dx;
00149
00150 old_E = E;
00151 E = x + m;
00152 ++count;
00153
00154 } while ((fabs(E-old_E) > 10*(fabs(E)+fabs(M)+twopi)*std::numeric_limits<double>::epsilon()) && (count < max_count));
00155
00156 if (iflag) {
00157 E = twopi - E;
00158 old_E = twopi - old_E;
00159 }
00160
00161 if (count >= max_count) {
00162 ORSA_WARNING("Orbit::GetE(): max count reached, e = %g E = %g fabs(E-old_E) = %g 10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon() = %g",e,E,fabs(E-old_E),10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon());
00163 }
00164 }
00165
00166 if (E == 0.0) {
00167 ORSA_LOGIC_ERROR("E==0.0 in orsa::Orbit::GetE(); this should never happen.");
00168 }
00169 return E;
00170 }
00171
00172
00173 void Orbit::RelativePosVel(Vector &relative_position, Vector &relative_velocity) const {
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 #ifdef HAVE_SINCOS
00185 double s,c;
00186
00187 ::sincos(omega_pericenter,&s,&c);
00188 const double sp = s;
00189 const double cp = c;
00190
00191 ::sincos(omega_node,&s,&c);
00192 const double so = s;
00193 const double co = c;
00194
00195 ::sincos(i,&s,&c);
00196 const double si = s;
00197 const double ci = c;
00198 #else // HAVE_SINCOS
00199 const double sp = sin(omega_pericenter);
00200 const double cp = cos(omega_pericenter);
00201
00202 const double so = sin(omega_node);
00203 const double co = cos(omega_node);
00204
00205 const double si = sin(i);
00206 const double ci = cos(i);
00207 #endif // HAVE_SINCOS
00208
00209 const Vector d1(cp*co - sp*so*ci,
00210 cp*so + sp*co*ci,
00211 sp*si);
00212
00213 const Vector d2(-sp*co - cp*so*ci,
00214 -sp*so + cp*co*ci,
00215 cp*si);
00216
00217
00218
00219
00220 double cape,tmp;
00221 double xfac1,xfac2,vfac1,vfac2;
00222
00223 double capf,shcap,chcap;
00224 double zpara;
00225
00226 if (e < 1.0) {
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 cape = GetE();
00239
00240
00241
00242 #ifdef HAVE_SINCOS
00243 ::sincos(cape,&s,&c);
00244 const double scap = s;
00245 const double ccap = c;
00246 #else // HAVE_SINCOS
00247 const double scap = sin(cape);
00248 const double ccap = cos(cape);
00249 #endif // HAVE_SINCOS
00250 const double sqe = secure_sqrt(1.0 - e*e);
00251 const double sqgma = secure_sqrt(mu*a);
00252
00253 xfac1 = a*(ccap - e);
00254 xfac2 = a*sqe*scap;
00255
00256 const double ri = 1.0/(a*(1.0 - e*ccap));
00257 vfac1 = -ri * sqgma * scap;
00258 vfac2 = ri * sqgma * ccap * sqe;
00259
00260 } else if (e > 1.0) {
00261
00262 double x,shx,chx,esh,ech;
00263 double f,fp,fpp,fppp,dx;
00264
00265
00266 double local_M = M;
00267 if (fabs(local_M-twopi) < fabs(local_M)) {
00268 local_M -= twopi;
00269 }
00270
00271
00272 if (local_M < 0.0) {
00273 tmp = -2.0*local_M/e + 1.8;
00274 x = -secure_log(tmp);
00275 } else {
00276 tmp = +2.0*local_M/e + 1.8;
00277 x = secure_log(tmp);
00278 }
00279
00280 capf = x;
00281
00282 int count = 0;
00283 do {
00284 x = capf;
00285 shx = sinh(x);
00286 chx = cosh(x);
00287 esh = e*shx;
00288 ech = e*chx;
00289 f = esh-x-local_M;
00290 fp = ech - 1.0;
00291 fpp = esh;
00292 fppp = ech;
00293 dx = -f/fp;
00294 dx = -f/(fp + dx*fpp/2.0);
00295 dx = -f/(fp + dx*fpp/2.0 + dx*dx*fppp/6.0);
00296 capf = x + dx;
00297 ++count;
00298 } while ((fabs(dx) > 1.0e-14) && (count < 100));
00299
00300 shcap = sinh(capf);
00301 chcap = cosh(capf);
00302
00303 const double sqe = secure_sqrt(e*e - 1.0);
00304 const double sqgma = secure_sqrt(mu*a);
00305 xfac1 = a*(e - chcap);
00306 xfac2 = a*sqe*shcap;
00307 const double ri = 1.0/(a*(e*chcap - 1.0));
00308 vfac1 = -ri * sqgma * shcap;
00309 vfac2 = ri * sqgma * chcap * sqe;
00310
00311 } else {
00312
00313 double q = M;
00314 if (q < 1.0e-3) {
00315 zpara = q*(1.0 - (q*q/3.0)*(1.0-q*q));
00316 } else {
00317 double x = 0.5*(3.0*q+secure_sqrt(9.0*(q*q)+4.0));
00318
00319 double tmp = cbrt(x);
00320 zpara = tmp - 1.0/tmp;
00321 }
00322
00323 const double sqgma = secure_sqrt(2.0*mu*a);
00324 xfac1 = a*(1.0 - zpara*zpara);
00325 xfac2 = 2.0*a*zpara;
00326 const double ri = 1.0/(a*(1.0 + zpara*zpara));
00327 vfac1 = -ri * sqgma * zpara;
00328 vfac2 = ri * sqgma;
00329
00330 }
00331
00332 relative_position = d1*xfac1 + d2*xfac2;
00333 relative_velocity = d1*vfac1 + d2*vfac2;
00334
00335 }
00336
00337 void Orbit::Compute(const Body &b, const Body &ref_b) {
00338
00339 const Vector dr = b.position() - ref_b.position();
00340 const Vector dv = b.velocity() - ref_b.velocity();
00341
00342 const double mu = GetG()*(b.mass()+ref_b.mass());
00343
00344 Compute(dr,dv,mu);
00345 }
00346
00347 void Orbit::Compute(const Vector &relative_position, const Vector &relative_velocity, const double mu_in) {
00348
00349
00350
00351
00352
00353
00354 mu = mu_in;
00355
00356 const double tiny = 1.0e-100;
00357
00358
00359 double face,cape,capf,tmpf,cw,sw,w,u;
00360 int ialpha = 0;
00361
00362
00363 Vector h = ExternalProduct(relative_position,relative_velocity);
00364
00365 double h2 = h.LengthSquared();
00366 double hh = h.Length();
00367
00368
00369 i = secure_acos(h.z/hh);
00370
00371
00372
00373 double fac = secure_sqrt(h.x*h.x+h.y*h.y)/h2;
00374
00375 if (fac < tiny) {
00376 omega_node = 0.0;
00377 u = secure_atan2(relative_position.y, relative_position.x);
00378 if ( fabs(i-pi) < 10.0 * tiny ) u = -u;
00379 } else {
00380 omega_node = secure_atan2(h.x,-h.y);
00381 u = secure_atan2(relative_position.z/sin(i), relative_position.x*cos(omega_node)+relative_position.y*sin(omega_node));
00382 }
00383
00384 if (omega_node < 0.0) omega_node += twopi;
00385 if (u < 0.0) u += twopi;
00386
00387
00388
00389 double r = relative_position.Length();
00390 double v2 = relative_velocity.LengthSquared();
00391
00392 double vdotr = relative_position*relative_velocity;
00393
00394 double energy = v2/2.0 - mu/r;
00395
00396
00397 if (fabs(energy*r/mu) < tiny) {
00398 ialpha = 0;
00399 } else {
00400 if (energy < 0) ialpha = -1;
00401 if (energy > 0) ialpha = +1;
00402 }
00403
00404
00405
00406
00407 if (ialpha == -1) {
00408
00409 a = -mu/(2.0*energy);
00410
00411 fac = 1.0 - h2/(mu*a);
00412
00413 if (fac > tiny) {
00414 e = secure_sqrt(fac);
00415 face = (a-r)/(a*e);
00416
00417 if (face > 1.0) {
00418 cape = 0.0;
00419 } else {
00420 if (face > -1.0)
00421 cape = secure_acos(face);
00422 else
00423 cape = pi;
00424 }
00425
00426 if (vdotr < 0.0) cape = twopi - cape;
00427 cw = (cos(cape)-e)/(1.0-e*cos(cape));
00428 sw = secure_sqrt(1.0-e*e)*sin(cape)/(1.0-e*cos(cape));
00429 w = secure_atan2(sw,cw);
00430 if (w < 0.0) w += twopi;
00431 } else {
00432 e = 0.0;
00433 w = u;
00434 cape = u;
00435 }
00436
00437 M = cape - e*sin(cape);
00438 omega_pericenter = u - w;
00439 if (omega_pericenter < 0) omega_pericenter += twopi;
00440 omega_pericenter = fmod(omega_pericenter,twopi);
00441
00442 }
00443
00444
00445 if (ialpha == 1) {
00446 a = mu/(2.0*energy);
00447 fac = h2/(mu*a);
00448
00449 if (fac > tiny) {
00450 e = secure_sqrt(1.0+fac);
00451 tmpf = (a+r)/(a*e);
00452 if (tmpf < 1.0) tmpf = 1.0;
00453
00454 capf = secure_log(tmpf+secure_sqrt(tmpf*tmpf-1.0));
00455
00456 if (vdotr < 0.0) capf = -capf;
00457
00458 cw = (e-cosh(capf))/(e*cosh(capf)-1.0);
00459 sw = secure_sqrt(e*e-1.0)*sinh(capf)/(e*cosh(capf)-1.0);
00460 w = secure_atan2(sw,cw);
00461 if (w < 0.0) w += twopi;
00462 } else {
00463
00464
00465 e = 1.0;
00466 tmpf = h2/(2.0*mu);
00467 w = secure_acos(2.0*tmpf/r - 1.0);
00468 if (vdotr < 0.0) w = twopi - w;
00469 tmpf = (a+r)/(a*e);
00470 capf = secure_log(tmpf+secure_sqrt(tmpf*tmpf-1.0));
00471 }
00472
00473 M = e * sinh(capf) - capf;
00474 omega_pericenter = u - w;
00475 if (omega_pericenter < 0) omega_pericenter += twopi;
00476 omega_pericenter = fmod(omega_pericenter,twopi);
00477 }
00478
00479
00480
00481 if (ialpha == 0) {
00482 a = 0.5*h2/mu;
00483 e = 1.0;
00484 w = secure_acos(2.0*a/r -1.0);
00485 if (vdotr < 0.0) w = twopi - w;
00486 tmpf = tan(0.5*w);
00487
00488 M = tmpf*(1.0+tmpf*tmpf/3.0);
00489 omega_pericenter = u - w;
00490 if (omega_pericenter < 0) omega_pericenter += twopi;
00491 omega_pericenter = fmod(omega_pericenter,twopi);
00492 }
00493
00494 }
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620 double RMS_residuals(const vector<Observation> &obs, const OrbitWithEpoch &orbit) {
00621
00622 double sum_delta_arcsec_squared=0.0;
00623 unsigned int k=0;
00624 Sky sky;
00625 OptimizedOrbitPositions opt(orbit);
00626 while (k<obs.size()) {
00627
00628 sky = opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode);
00629
00630
00631 const double delta_arcsec = sky.delta_arcsec(obs[k]);
00632 sum_delta_arcsec_squared += delta_arcsec*delta_arcsec;
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646 ++k;
00647 }
00648
00649 return (secure_sqrt(sum_delta_arcsec_squared/obs.size()));
00650 }
00651
00652 double residual(const Observation & obs, const OrbitWithEpoch & orbit) {
00653 OptimizedOrbitPositions opt(orbit);
00654 Sky sky = opt.PropagatedSky_J2000(obs.date,obs.obscode);
00655 return fabs(sky.delta_arcsec(obs));
00656 }
00657
00658
00659
00660
00661
00662 Sky PropagatedSky_J2000(const OrbitWithEpoch & orbit, const UniverseTypeAwareTime & final_time, const std::string & obscode, const bool integrate, const bool light_time_corrections) {
00663
00664 if (!integrate) {
00665
00666 list<JPL_planets> l;
00667
00668 l.push_back(SUN);
00669 l.push_back(EARTH);
00670
00671 Frame frame;
00672
00673 SetupSolarSystem(frame,l,final_time);
00674
00675
00676
00677
00678
00679
00680
00681
00682 Vector geo = frame[1].position();
00683
00684 geo += location_file->ObsPos(obscode,final_time.GetDate());
00685
00686 Vector r,v;
00687 orbit.RelativePosVel(r,v,final_time);
00688
00689 r += frame[0].position();
00690 v += frame[0].velocity();
00691
00692 const Vector relative_position = r - geo;
00693
00694 Sky sky;
00695
00696 if (light_time_corrections) {
00697 sky.Compute_J2000(relative_position-v*relative_position.Length()/GetC());
00698 } else {
00699 sky.Compute_J2000(relative_position);
00700 }
00701
00702 return sky;
00703 }
00704
00705 Frame frame;
00706
00707 list<JPL_planets> l;
00708
00709 l.push_back(SUN);
00710 l.push_back(MERCURY);
00711 l.push_back(VENUS);
00712 l.push_back(EARTH_AND_MOON);
00713 l.push_back(MARS);
00714 l.push_back(JUPITER);
00715 l.push_back(SATURN);
00716 l.push_back(URANUS);
00717 l.push_back(NEPTUNE);
00718
00719 SetupSolarSystem(frame,l,orbit.epoch);
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735 {
00736 Body b("object",0.0);
00737
00738 Vector r,v;
00739 orbit.RelativePosVel(r,v);
00740
00741
00742
00743
00744
00745 r += frame[0].position();
00746 v += frame[0].velocity();
00747
00748 b.SetPosition(r);
00749 b.SetVelocity(v);
00750 frame.push_back(b);
00751
00752
00753
00754 }
00755
00756 const Frame orbit_epoch_frame = frame;
00757
00758 UniverseTypeAwareTime obs_utat;
00759
00760 Radau15 itg;
00761
00762 itg.accuracy = 1.0e-12;
00763 itg.timestep = FromUnits(1.0,DAY);
00764
00765 Newton newton;
00766
00767
00768
00769
00770
00771 Evolution evol;
00772
00773
00774 evol.SetInteraction(&newton);
00775
00776
00777 evol.SetIntegrator(&itg);
00778 evol.push_back(orbit_epoch_frame);
00779 evol.SetMaxUnsavedSubSteps(100000);
00780 obs_utat.SetDate(final_time.GetDate());
00781
00782
00783
00784 evol.SetSamplePeriod(orbit_epoch_frame - obs_utat);
00785 evol.Integrate(obs_utat.GetTime(),true);
00786 Frame last_frame = evol[evol.size()-1];
00787
00788
00789
00790
00791
00792
00793
00794 Vector geo = last_frame[3].position()-last_frame[0].position();
00795
00796
00797
00798 geo += location_file->ObsPos(obscode,final_time.GetDate());
00799
00800
00801
00802
00803
00804
00805 const Vector relative_position = last_frame[last_frame.size()-1].position() - last_frame[0].position() - geo;
00806 const Vector v = last_frame[last_frame.size()-1].velocity();
00807
00808
00809
00810
00811
00812 Sky sky;
00813
00814
00815
00816
00817 if (light_time_corrections) {
00818 sky.Compute_J2000(relative_position-v*relative_position.Length()/GetC());
00819 } else {
00820 sky.Compute_J2000(relative_position);
00821 }
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 return sky;
00865 }
00866
00867
00868
00869 OptimizedOrbitPositions::OptimizedOrbitPositions(const OrbitWithEpoch &orbit) : _orbit(orbit) {
00870
00871 l.push_back(SUN);
00872 l.push_back(MERCURY);
00873 l.push_back(VENUS);
00874 l.push_back(EARTH_AND_MOON);
00875 l.push_back(MARS);
00876 l.push_back(JUPITER);
00877 l.push_back(SATURN);
00878 l.push_back(URANUS);
00879 l.push_back(NEPTUNE);
00880
00881 frames.clear();
00882 }
00883
00884 OrbitWithEpoch OptimizedOrbitPositions::PropagatedOrbit(const UniverseTypeAwareTime & final_time, const bool integrate) {
00885
00886 std::sort(frames.begin(),frames.end());
00887
00888
00889
00890
00891
00892
00893
00894 if (!integrate) {
00895 OrbitWithEpoch oe(_orbit);
00896 oe.M += twopi*(final_time.Time() - _orbit.epoch.Time())/(_orbit.Period());
00897 oe.M = fmod(10*twopi+fmod(oe.M,twopi),twopi);
00898 oe.epoch = final_time;
00899 return (oe);
00900 }
00901
00902 Frame start_frame;
00903
00904 if (frames.size() > 0) {
00905
00906 if (final_time <= frames[0]) {
00907 start_frame = frames[0];
00908
00909 } else if (final_time >= frames[frames.size()-1]) {
00910 start_frame = frames[frames.size()-1];
00911
00912 } else {
00913 unsigned int k=1;
00914 while (k < frames.size()) {
00915
00916 if ( (final_time >= frames[k-1]) && (final_time <= frames[k]) ) {
00917 if ((final_time-frames[k-1]).absolute() < (frames[k]-final_time).absolute()) {
00918 start_frame = frames[k-1];
00919
00920 break;
00921 } else {
00922 start_frame = frames[k];
00923
00924 break;
00925 }
00926 }
00927
00928 ++k;
00929 }
00930 }
00931
00932 } else {
00933
00934
00935
00936 SetupSolarSystem(start_frame,l,_orbit.epoch);
00937
00938
00939 const JPLBody Sun(SUN,_orbit.epoch);
00940
00941 Body b("object");
00942 Vector r, v;
00943 _orbit.RelativePosVel(r,v);
00944 r += Sun.position();
00945 v += Sun.velocity();
00946
00947 b.SetPosition(r);
00948 b.SetVelocity(v);
00949
00950 start_frame.push_back(b);
00951
00952 frames.push_back(start_frame);
00953 }
00954
00955 Radau15 itg;
00956 itg.accuracy = 1.0e-12;
00957 itg.timestep = FromUnits(1.0,DAY);
00958
00959 Evolution evol;
00960 evol.SetInteraction(NEWTON);
00961 evol.SetIntegrator(&itg);
00962 evol.push_back(start_frame);
00963 evol.SetMaxUnsavedSubSteps(100000);
00964
00965
00966
00967
00968
00969
00970 evol.SetSamplePeriod(final_time - start_frame);
00971 evol.Integrate(final_time,true);
00972
00973 Frame last_frame = evol[evol.size()-1];
00974
00975
00976
00977
00978 if (evol.size() > 1) frames.push_back(last_frame);
00979
00980 const Vector relative_position = last_frame[last_frame.size()-1].position() - last_frame[0].position();
00981 const Vector relative_velocity = last_frame[last_frame.size()-1].velocity() - last_frame[0].velocity();
00982
00983 OrbitWithEpoch oe;
00984 oe.Compute(relative_position,relative_velocity,GetG()*GetMSun(),final_time);
00985
00986 return (oe);
00987 }
00988
00989 Sky OptimizedOrbitPositions::PropagatedSky_J2000(const UniverseTypeAwareTime &final_time, const string &obscode, const bool integrate, const bool light_time_corrections) {
00990
00991 OrbitWithEpoch oe = PropagatedOrbit(final_time,integrate);
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010 JPLBody Sun(SUN,final_time.GetDate());
01011 JPLBody Earth(EARTH,final_time.GetDate());
01012
01013
01014 Vector geo = Earth.position();
01015
01016 geo += location_file->ObsPos(obscode,final_time.GetDate());
01017
01018 Vector r,v;
01019 oe.RelativePosVel(r,v,final_time);
01020 r += Sun.position();
01021 v += Sun.velocity();
01022
01023 Vector relative_position = r - geo;
01024
01025 Sky sky;
01026
01027
01028
01029 if (light_time_corrections) {
01030 sky.Compute_J2000(relative_position-v*relative_position.Length()/GetC());
01031 } else {
01032 sky.Compute_J2000(relative_position);
01033 }
01034
01035 return sky;
01036 }
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176 void OptimizedOrbitPositions::PropagatedPosVel(const UniverseTypeAwareTime &final_time, Vector &position, Vector &velocity, bool integrate) {
01177 OrbitWithEpoch oe = PropagatedOrbit(final_time,integrate);
01178 oe.RelativePosVel(position,velocity,final_time);
01179 }
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313 void Sky::Compute_J2000(const Vector &relative_position) {
01314
01315 switch (universe->GetReferenceSystem()) {
01316 case ECLIPTIC: {
01317 Vector r = relative_position;
01318
01319
01320 EclipticToEquatorial_J2000(r);
01321 _ra.SetRad(fmod(secure_atan2(r.y,r.x)+twopi,twopi));
01322 _dec.SetRad(pi/2 - secure_acos(r.z/r.Length()));
01323 break;
01324 }
01325 case EQUATORIAL: {
01326 _ra.SetRad(fmod(secure_atan2(relative_position.y,relative_position.x)+twopi,twopi));
01327 _dec.SetRad(pi/2 - secure_acos(relative_position.z/relative_position.Length()));
01328 break;
01329 }
01330 }
01331 }
01332
01333 static Vector unit_vector(const Angle & ra, const Angle & dec) {
01334 double sin_ra, cos_ra;
01335 sincos(ra,sin_ra,cos_ra);
01336
01337 double sin_dec, cos_dec;
01338 sincos(dec,sin_dec,cos_dec);
01339
01340 return Vector(cos_dec*sin_ra,
01341 cos_dec*cos_ra,
01342 sin_dec);
01343 }
01344
01345 double Sky::delta_arcsec(const Observation & obs) const {
01346 return ((180/pi)*3600.0*acos(unit_vector(obs.ra,obs.dec)*unit_vector(_ra,_dec)));
01347 }
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373 }