orsa_orbit.cc File Reference

#include "orsa_orbit.h"
#include "support.h"
#include "orsa_units.h"
#include "orsa_common.h"
#include "orsa_universe.h"
#include "orsa_file.h"
#include "orsa_frame.h"
#include "orsa_secure_math.h"
#include <sys/time.h>
#include <time.h>
#include <cstring>
#include <limits>
#include <algorithm>

Include dependency graph for orsa_orbit.cc:

Go to the source code of this file.

Namespaces

namespace  orsa

Functions

double RMS_residuals (const vector< Observation > &obs, const OrbitWithEpoch &orbit)
 The RMS of the residuals in arcsec units.
double residual (const Observation &obs, const OrbitWithEpoch &orbit)
 The RMS of the residuals in arcsec units.
Sky PropagatedSky_J2000 (const OrbitWithEpoch &orbit, const UniverseTypeAwareTime &final_time, const std::string &obscode, const bool integrate, const bool light_time_corrections)
static Vector unit_vector (const Angle &ra, const Angle &dec)


Function Documentation

Sky PropagatedSky_J2000 const OrbitWithEpoch &  orbit,
const UniverseTypeAwareTime &  final_time,
const std::string &  obscode,
const bool  integrate = true,
const bool  light_time_corrections = true
 

Definition at line 662 of file orsa_orbit.cc.

References Integrator::accuracy, Sky::Compute_J2000(), orsa::DAY, Radau15::e, orsa::EARTH, orsa::EARTH_AND_MOON, orsa::FromUnits(), orsa::GetC(), UniverseTypeAwareTime::GetTime(), Evolution::Integrate(), orsa::JUPITER, Vector::Length(), orsa::location_file, orsa::MARS, orsa::MERCURY, orsa::NEPTUNE, LocationFile::ObsPos(), Evolution::push_back(), orsa::SATURN, UniverseTypeAwareTime::SetDate(), Evolution::SetIntegrator(), Evolution::SetInteraction(), Evolution::SetMaxUnsavedSubSteps(), Body::SetPosition(), Evolution::SetSamplePeriod(), orsa::SetupSolarSystem(), Body::SetVelocity(), Frame::size(), Evolution::size(), orsa::SUN, Integrator::timestep, orsa::URANUS, and orsa::VENUS.

Referenced by orsa::gauss_v_diff_f(), and orsa::least_sq_diff_f().

00662                                                                                                                                                                                       {
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          LocationFile lf;
00677          lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str());
00678          lf.Read();
00679          lf.Close();
00680       */
00681       
00682       Vector geo = frame[1].position();
00683       // geo += lf.ObsPos(obscode,final_time.GetDate());
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        Config conf;
00723        OrsaConfigFile ocf(&conf);
00724        ocf.Read();
00725        ocf.Close();
00726     */
00727     
00728     /* 
00729        LocationFile lf;
00730        lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str());
00731        lf.Read();
00732        lf.Close();
00733     */
00734     
00735     {
00736       Body b("object",0.0);
00737       // RelativePosVel(b.position,b.velocity);
00738       Vector r,v;
00739       orbit.RelativePosVel(r,v);
00740       //
00741       // cerr << "PropagatedSky ---> RelativePosVel ---> r: " << Length(r) << endl;
00742       //
00743       // cerr << "b.position: " << b.position.x << "  " << b.position.y << "  " << b.position.z << endl;
00744       // add Sun postion and velocity
00745       r += frame[0].position();
00746       v += frame[0].velocity();
00747       // cerr << "b.position: " << b.position.x << "  " << b.position.y << "  " << b.position.z << endl;
00748       b.SetPosition(r);
00749       b.SetVelocity(v);
00750       frame.push_back(b);
00751       
00752       // debug
00753       // cerr << "object r: " << Length(r) << endl; 
00754     }
00755     
00756     const Frame orbit_epoch_frame = frame;
00757     
00758     UniverseTypeAwareTime obs_utat;
00759     
00760     Radau15 itg; 
00761     // itg.accuracy = 1.0e-12;
00762     itg.accuracy = 1.0e-12;
00763     itg.timestep = FromUnits(1.0,DAY);
00764     //
00765     Newton newton;
00766     
00767     // Frame  last_frame;
00768     // Vector relative_position;
00769     // Vector v;
00770     
00771     Evolution evol;
00772     
00773     // evol.interaction = &newton;
00774     evol.SetInteraction(&newton);
00775     // evol.interaction = &itr;
00776     // evol.integrator  = &itg;
00777     evol.SetIntegrator(&itg);
00778     evol.push_back(orbit_epoch_frame);
00779     evol.SetMaxUnsavedSubSteps(100000);
00780     obs_utat.SetDate(final_time.GetDate()); 
00781     // evol.sample_period = orbit_epoch_frame.Time() - obs_utat.Time(); // important, otherwise the final frame is not saved
00782     // evol.sample_period = (orbit_epoch_frame.Time() - obs_utat.Time())*(1.0-1.0e-8); // important, otherwise the final frame is not saved
00783     // evol.sample_period = (orbit_epoch_frame.Time() - obs_utat.Time());
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     // cerr << "sample_period: " << evol.sample_period << endl;
00789     // fprintf(stderr,"final time JD: %20.8f   last frame JD: %20.8f\n",final_time.GetDate().GetJulian(),last_frame.GetDate().GetJulian());
00790     
00791     // cerr << "..is Sun? -----> " << frame[0].name() << endl;
00792     // cerr << "..is Earth?  --> " << frame[3].name() << endl;
00793     //
00794     Vector geo = last_frame[3].position()-last_frame[0].position(); // earth center position
00795     // Vector geo = frame[3].position; // earth center position
00796     
00797     // IMPORTANT!!!
00798     geo += location_file->ObsPos(obscode,final_time.GetDate());
00799     
00800     // cerr << "obs. pos. length: " << Length(lf.ObsPos(obscode,final_time.GetDate())) << "  " << LengthLabel() << endl;
00801     // cerr << "obs. pos. length: " << FromUnits(Length(lf.ObsPos(obscode,final_time.GetDate())),REARTH,-1) << "  " << units.label(REARTH) << endl;
00802     
00803     // cerr << "Length(geo): " << Length(geo) << endl; 
00804     
00805     const Vector relative_position = last_frame[last_frame.size()-1].position() - last_frame[0].position() - geo; // object - (earth + observer)
00806     const Vector v                 = last_frame[last_frame.size()-1].velocity();
00807     
00808     // cerr << "relative_position: " << relative_position.x << "  " << relative_position.y << "  " << relative_position.z << "  " << endl;
00809     
00810     // cerr << "relative distance: " << Length(relative_position) << "  " << LengthLabel() << endl;
00811     
00812     Sky sky;
00813     // sky.Compute(relative_position,obs_utat);
00814     // sky.Compute(relative_position,last_frame);
00815     // sky.Compute_J2000(relative_position);
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     // fprintf(stderr,"propagated sky:   ra = %16.12f    dec = %16.12f   orbit JD: %14.5f   final time JD: %14.5f\n",sky.ra().GetRad(),sky.dec().GetRad(),orbit.epoch.GetDate().GetJulian(),final_time.GetDate().GetJulian());
00824     
00825     // debug test
00826     /* 
00827        {
00828        
00829        // random number generator
00830        // const int random_seed = 124323;
00831        struct timeval  tv;
00832        struct timezone tz;
00833        gettimeofday(&tv,&tz); 
00834        const int random_seed = tv.tv_usec;
00835        gsl_rng *rnd;
00836        rnd = gsl_rng_alloc(gsl_rng_gfsr4);
00837        gsl_rng_set(rnd,random_seed);
00838        
00839        Vector geo,rand;
00840        int oc;
00841        Sky sky;
00842        // double h,m,s,d,p;
00843        for (oc=1;oc<400;++oc) {
00844        // geo = last_frame[3].position + lf.ObsPos(oc,final_time.GetDate());
00845        rand.Set(gsl_rng_uniform(rnd)*2-1,
00846        gsl_rng_uniform(rnd)*2-1,
00847        gsl_rng_uniform(rnd)*2-1);
00848        rand /= Length(rand);
00849        rand *= FromUnits(1.0,REARTH);
00850        geo = last_frame[3].position+rand;
00851        relative_position = last_frame[last_frame.size()-1].position - geo;
00852        sky.Compute(relative_position,obs_utat);
00853        // sky.ra().GetHMS(h,m,s);
00854        // printf("R.A.:   %2f %2f %10.4f    ",h,m,s);
00855        // sky.dec().GetDPS(d,p,s);
00856        // printf("dec.:   %2f %2f %10.4f\n",d,p,s);
00857        printf("%20.12f  %20.12f\n",sky.ra().GetRad(),sky.dec().GetRad());
00858        }
00859        
00860        gsl_rng_free(rnd);      
00861        }
00862     */
00863     
00864     return sky;
00865   }

Here is the call graph for this function:

double residual const Observation &  obs,
const OrbitWithEpoch &  orbit
 

The RMS of the residuals in arcsec units.

Definition at line 652 of file orsa_orbit.cc.

References Sky::delta_arcsec(), and OptimizedOrbitPositions::PropagatedSky_J2000().

00652                                                                          {
00653     OptimizedOrbitPositions opt(orbit);
00654     Sky sky = opt.PropagatedSky_J2000(obs.date,obs.obscode);
00655     return fabs(sky.delta_arcsec(obs));
00656   }

Here is the call graph for this function:

double RMS_residuals const vector< Observation > &  obs,
const OrbitWithEpoch &  orbit
 

The RMS of the residuals in arcsec units.

Definition at line 620 of file orsa_orbit.cc.

References Sky::delta_arcsec(), OptimizedOrbitPositions::PropagatedSky_J2000(), and orsa::secure_sqrt().

Referenced by orsa::Compute_Gauss(), PreliminaryOrbit::ComputeRMS(), and orsa::E1().

00620                                                                                     {
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       // sky = PropagatedSky(orbit,obs[k].date,obs[k].obscode);
00628       sky = opt.PropagatedSky_J2000(obs[k].date,obs[k].obscode);
00629       // sum_delta_arcsec_squared += secure_pow(PropagatedSky(orbit,obs[k].date,obs[k].obscode).delta_arcsec(obs[k]),2);
00630       // sum_delta_arcsec_squared += secure_pow(sky.delta_arcsec(obs[k]),2);
00631       const double delta_arcsec = sky.delta_arcsec(obs[k]);
00632       sum_delta_arcsec_squared += delta_arcsec*delta_arcsec;
00633       //
00634       /* if (k>0) {
00635          fprintf(stderr,"T-%i-%i: %g days\n",k,k-1,obs[k].date.GetJulian()-obs[k-1].date.GetJulian());
00636          }
00637       */
00638       
00639       // debug
00640       /*
00641         printf("%16.12f %16.12f\n",obs[k].ra.GetRad(),obs[k].dec.GetRad());
00642         printf("%16.12f %16.12f\n",sky.ra().GetRad(),sky.dec().GetRad());
00643       */
00644       //
00645       
00646       ++k;
00647     }
00648     
00649     return (secure_sqrt(sum_delta_arcsec_squared/obs.size())); 
00650   }

Here is the call graph for this function:

static Vector unit_vector const Angle &  ra,
const Angle &  dec
[static]
 

Definition at line 1333 of file orsa_orbit.cc.

References orsa::sincos().

Referenced by Sky::delta_arcsec().

01333                                                                  {
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   }

Here is the call graph for this function:


Generated on Tue Jan 11 15:27:11 2005 for liborsa by  doxygen 1.4.0