orsa_orbit.cc

Go to the documentation of this file.
00001 /* 
00002    ORSA - Orbit Reconstruction, Simulation and Analysis
00003    Copyright (C) 2002-2004 Pasquale Tricarico
00004    
00005    This program is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU General Public License
00007    as published by the Free Software Foundation; either version 2
00008    of the License, or (at your option) any later version.
00009    
00010    As a special exception, Pasquale Tricarico gives permission to
00011    link this program with Qt commercial edition, and distribute the
00012    resulting executable, without including the source code for the Qt
00013    commercial edition in the source distribution.
00014    
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019    
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the Free Software
00022    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
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       // begin with a guess accurate to order e^3 
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         // update x, ready for the next iteration
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       // Make a first guess that works well for e near 1.  
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     // This alghoritm is taken from the swift package, 
00177     // ORBEL_EL2XV.F file written by M. Duncan.
00178     /////////////////////////////////////////////////////
00179     
00180     // Generate rotation matrices (on p. 42 of Fitzpatrick)
00181     
00182     // double sp,cp,so,co,si,ci;
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     // Get the other quantities depending on orbit type
00218     
00219     // double  cape,scap,ccap,sqe,sqgma,tmp;
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          int count = 0;
00230          cape = M;
00231          do {
00232          tmp  = cape;
00233          cape = e * sin(cape) + M;
00234          ++count;
00235          } while ( (fabs((cape-tmp)/cape) > 1.0e-15) && (count < 100) );
00236       */
00237       
00238       cape = GetE();
00239       
00240       // cerr << "tmp: " << tmp << "  cape: " << cape << "  fabs(cape - tmp)" << fabs(cape - tmp) << endl;
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       // ri = 1/r
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       // use the 'right' value for M -- NEEDED!
00266       double local_M = M;
00267       if (fabs(local_M-twopi) < fabs(local_M)) {
00268         local_M -= twopi;
00269       }
00270       
00271       // begin with a guess proposed by Danby   
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 { // e = 1.0 within roundoff errors
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         // double tmp = secure_pow(x,(1.0/3.0));
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     // This alghoritm is taken from the swift package, 
00351     // ORBEL_XV2EL.F file written by M. Duncan.
00352     /////////////////////////////////////////////////////
00353     
00354     mu = mu_in;
00355     
00356     const double tiny = 1.0e-100; // about 4.0e-15
00357     
00358     // internals
00359     double  face,cape,capf,tmpf,cw,sw,w,u;
00360     int ialpha = 0;
00361     
00362     // angular momentum
00363     Vector h = ExternalProduct(relative_position,relative_velocity);
00364     
00365     double h2 = h.LengthSquared();
00366     double hh = h.Length();
00367     
00368     // inclination
00369     i = secure_acos(h.z/hh);
00370     
00371     // Compute longitude of ascending node omega_node and the argument of latitude u
00372     // double fac = secure_sqrt(secure_pow(h.x,2)+secure_pow(h.y,2))/h2;
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     //  Compute the radius r and velocity squared v2, and the dot
00388     //  product rdotv, the energy per unit mass energy 
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     // Determine type of conic section and label it via ialpha
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     // Depending on the conic type, determine the remaining elements 
00405     
00406     // ellipse 
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     // hyperbola
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         // we only get here if a hyperbola is essentially a parabola 
00464         // so we calculate e and w accordingly to avoid singularities
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     // parabola
00480     //  NOTE - in this case we use "a" to mean pericentric distance
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   /* //! The RMS of the residuals in arcsec units
00497      double RMS_residuals(vector<Observation> &obs, OrbitWithEpoch &orbit) {
00498      
00499      Frame frame;
00500      
00501      list<JPL_planets> l;
00502      l.push_back(SUN);
00503      l.push_back(MERCURY);
00504      l.push_back(VENUS);
00505      l.push_back(EARTH_AND_MOON);
00506      l.push_back(MARS);
00507      l.push_back(JUPITER);
00508      l.push_back(SATURN);
00509      l.push_back(URANUS);
00510      l.push_back(NEPTUNE);
00511      
00512      SetupSolarSystem(frame,l,orbit.epoch);
00513      
00514      Config conf;
00515      OrsaConfigFile ocf(&conf);
00516      ocf.Read();
00517      ocf.Close();
00518      // vector<Location> loc;
00519      // LocationFile lf(&loc);
00520      LocationFile lf;
00521      lf.SetFileName(conf.paths[OBSCODE]->GetValue().c_str());
00522      lf.Read();
00523      lf.Close();
00524      
00525      {
00526      Body b("object",0.0);
00527      // RelativePosVel(b.position,b.velocity);
00528      Vector r,v;
00529      orbit.RelativePosVel(r,v);
00530      // cerr << "b.position: " << b.position.x << "  " << b.position.y << "  " << b.position.z << endl;
00531      // add Sun postion and velocity
00532      r += frame[0].position();
00533      v += frame[0].velocity();
00534      // cerr << "b.position: " << b.position.x << "  " << b.position.y << "  " << b.position.z << endl;
00535      b.SetPosition(r);
00536      b.SetVelocity(v);
00537      frame.push_back(b);
00538      }
00539      
00540      const Frame orbit_epoch_frame = frame;
00541      
00542      // cerr << "orbit_epoch_frame" << endl;
00543      // print(orbit_epoch_frame);
00544      
00545      UniverseTypeAwareTime obs_utat;
00546      Radau15 itg; 
00547      itg.accuracy = 1.0e-10;
00548      itg.timestep = FromUnits(1,HOUR);
00549      //
00550      Newton newton;
00551      // evol_base.integrator = itg;
00552      // evol_base.sample_period = FromUnits(3,HOUR);
00553      // evol_base.SetMaxUnsavedSubSteps(1000);
00554      Frame  last_frame;
00555      Vector relative_position;
00556      Sky    sky;
00557      double sum_delta_arcsec_squared=0.0, tmp_delta_arcsec;
00558      
00559      unsigned int k=0;
00560      while (k<obs.size()) {
00561      
00562      // Evolution *evol = new Evolution(evol_base);
00563      // Evolution *evol = new Evolution;
00564      Evolution evol;
00565      // push back only when debugging..
00566      // universe->push_back(evol);
00567      // evol->interaction = new Newton();
00568      evol.interaction = &newton;
00569      evol.integrator  = &itg;
00570      evol.push_back(orbit_epoch_frame);
00571      evol.SetMaxUnsavedSubSteps(100000);
00572      obs_utat.SetDate(obs[k].date); 
00573      // cerr << "integration period: " << FromUnits(orbit_epoch_frame.Time() - obs_utat.GetTime(),DAY,-1) << " days" << endl;
00574      evol.sample_period = orbit_epoch_frame.Time() - obs_utat.Time(); // important, otherwise the final frame is not saved
00575      evol.Integrate(obs_utat.GetTime());
00576      last_frame = evol[evol.size()-1];
00577      
00578      Vector geo = last_frame[3].position(); // earth center position
00579      
00580      // cerr << "last_frame[3].name() = " << last_frame[3].name() << endl;
00581      
00582      geo += lf.ObsPos(obs[k].obscode,obs[k].date);
00583      
00584      if (0) {
00585      unsigned int k=0;
00586      while (k<last_frame.size()) {
00587      cerr << "name: " << last_frame[k].name() << "    position: " << last_frame[k].position().x << "  " << last_frame[k].position().y << "  " << last_frame[k].position().z << endl;
00588      ++k;
00589      }
00590      }
00591      
00592      // relative_position = last_frame[last_frame.size()-1].position - last_frame[3].position; // object - earth
00593      relative_position = last_frame[last_frame.size()-1].position() - geo; // object - (earth + observer)
00594      
00595      // cerr << "relative_position: " << relative_position.x << "  " << relative_position.y << "  " << relative_position.z << "  " << endl;
00596      
00597      sky.Compute(relative_position,obs_utat);
00598      tmp_delta_arcsec = sky.delta_arcsec(obs[k]);
00599      // cerr << "obs. " << k << " delta_arcsec: " << tmp_delta_arcsec << endl;
00600      
00601      // cerr << "OBS:  R.A. = " << obs[k].ra.GetRad() << "   DEC. = " << obs[k].dec.GetRad() << endl;
00602      // cerr << "OBJ:  R.A. = " << sky.ra().GetRad()  << "   DEC. = " << sky.dec().GetRad()  << endl;
00603      
00604      // cerr << obs[k].ra.GetRad() << "  " << obs[k].dec.GetRad() << "  " << sky.ra().GetRad() << "  " << sky.dec().GetRad() << endl;
00605      
00606      sum_delta_arcsec_squared += secure_pow(tmp_delta_arcsec,2);
00607      
00608      // cerr << "RMS: adding " << tmp_delta_arcsec << endl;
00609      
00610      ++k;
00611      }
00612      
00613      return secure_sqrt(sum_delta_arcsec_squared/obs.size()); 
00614      }
00615   */
00616   
00617   //
00618   
00619   //! The RMS of the residuals in arcsec units
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       // 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   }
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   // Sky PropagatedSky(const OrbitWithEpoch &orbit, const UniverseTypeAwareTime &final_time, const std::string &obscode, bool integrate) {
00661   // Sky PropagatedSky_J2000(const OrbitWithEpoch &orbit, const UniverseTypeAwareTime &final_time, const std::string &obscode, bool integrate) {
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          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   }
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        for (unsigned int k=0;k<frames.size();++k) {
00890        ORSA_ERROR("OptimizedOrbitPositions::PropagatedOrbit(...)  ==> frame[%i] time: %f   dt: %f",k,frames[k].GetTime(),(frames[k]-final_time).absolute().GetDouble());
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         // ORSA_ERROR("OptimizedOrbitPositions::PropagatedOrbit(...)  ==> starting from a saved frame...");
00909       } else if (final_time >= frames[frames.size()-1]) {
00910         start_frame = frames[frames.size()-1];  
00911         // ORSA_ERROR("OptimizedOrbitPositions::PropagatedOrbit(...)  ==> starting from a saved frame...");
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               // ORSA_ERROR("OptimizedOrbitPositions::PropagatedOrbit(...)  ==> starting from a saved frame...");
00920               break;
00921             } else {
00922               start_frame = frames[k];
00923               // ORSA_ERROR("OptimizedOrbitPositions::PropagatedOrbit(...)  ==> starting from a saved frame...");
00924               break;
00925             }
00926           }
00927           
00928           ++k;
00929         }
00930       }
00931       
00932     } else {
00933       
00934       // ORSA_ERROR("OptimizedOrbitPositions::PropagatedOrbit(...)  ==> first call...");
00935       
00936       SetupSolarSystem(start_frame,l,_orbit.epoch);
00937       // start_frame.SetDate(_orbit.epoch);
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       // cerr << "b.position: " << b.position.x << "  " << b.position.y << "  " << b.position.z << endl;
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     // obs_utat.SetDate(final_time.GetDate()); 
00965     // evol.sample_period = orbit_epoch_frame.Time() - obs_utat.Time(); // important, otherwise the final frame is not saved
00966     // evol.sample_period = (orbit_epoch_frame.Time() - obs_utat.Time())*(1.0-1.0e-8); // important, otherwise the final frame is not saved
00967     // evol.sample_period = (start_frame.Time() - obs_utat.Time())*(1.0-1.0e-8); // important, otherwise the final frame is not saved
00968     // evol.sample_period = (start_frame.Time() - obs_utat.Time());
00969     // evol.sample_period = start_frame - final_time;
00970     evol.SetSamplePeriod(final_time - start_frame);
00971     evol.Integrate(final_time,true);
00972     
00973     Frame last_frame = evol[evol.size()-1];
00974     
00975     // don't add the frame if the size of evol is 1,
00976     // because this means that the integration was not needed,
00977     // and the 'frames' already contains it.
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     // Frame frame;
00994     // SetupSolarSystem(frame,l,final_time);
00995     
00996     /* 
00997        Config conf;
00998        OrsaConfigFile ocf(&conf);
00999        ocf.Read();
01000        ocf.Close();
01001     */
01002     
01003     /* 
01004        LocationFile lf;
01005        lf.SetFileName(config->paths[OBSCODE]->GetValue().c_str());
01006        lf.Read();
01007        lf.Close();
01008     */
01009     
01010     JPLBody Sun(SUN,final_time.GetDate());
01011     JPLBody Earth(EARTH,final_time.GetDate());
01012     
01013     // Vector geo = frame[3].position()-frame[0].position();
01014     Vector geo = Earth.position();
01015     // geo += lf.ObsPos(obscode,final_time.GetDate());
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     // sky.Compute(relative_position,final_time);
01027     // sky.Compute_J2000(relative_position);
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      Sky OptimizedOrbitPositions::PropagatedSky(const UniverseTypeAwareTime &final_time, const string &obscode, bool integrate) {
01040      
01041      sort(frames.begin(),frames.end());
01042      
01043      if (!integrate) {
01044      
01045      Vector r,v;
01046      _orbit.RelativePosVel(r,v,final_time);
01047      
01048      Frame frame;
01049      
01050      SetupSolarSystem(frame,l,final_time);
01051      
01052      Config conf;
01053      OrsaConfigFile ocf(&conf);
01054      ocf.Read();
01055      ocf.Close();
01056      
01057      LocationFile lf;
01058      lf.SetFileName(conf.paths[OBSCODE]->GetValue().c_str());
01059      lf.Read();
01060      lf.Close();
01061      
01062      Vector geo = frame[3].position()-frame[0].position();
01063      geo += lf.ObsPos(obscode,final_time.GetDate());
01064      
01065      Vector relative_position = r - geo;
01066      
01067      Sky sky;
01068      sky.Compute(relative_position,final_time);
01069      
01070      return sky;
01071      }
01072      
01073      Frame start_frame;
01074      
01075      if (frames.size()>1) {
01076      
01077      if (final_time.GetTime() < frames[0].GetTime()) {
01078      start_frame = frames[0];
01079      // cerr << "starting from a saved frame..." << endl;
01080      } else if (final_time.GetTime() > frames[frames.size()-1].GetTime()) {
01081      start_frame = frames[frames.size()-1];     
01082      // cerr << "starting from a saved frame..." << endl;
01083      } else {
01084      unsigned int k=1;
01085      while (k<frames.size()) {
01086      if ( (final_time.GetTime() > frames[k-1].GetTime()) &&
01087      (final_time.GetTime() < frames[k].GetTime()) ) {
01088      if ((final_time.GetTime()-frames[k-1].GetTime()) < (frames[k].GetTime()-final_time.GetTime())) {
01089      start_frame = frames[k-1];
01090      // cerr << "starting from a saved frame..." << endl;
01091      break;
01092      } else {
01093      start_frame = frames[k];
01094      // cerr << "starting from a saved frame..." << endl;
01095      break;
01096      }
01097      }
01098      ++k;
01099      }
01100      }    
01101      
01102      } else {
01103      
01104      // cerr << "don't using a saved frame..." << endl;
01105      
01106      SetupSolarSystem(start_frame,l,_orbit.epoch);
01107      
01108      Body b("object",0.0);
01109      Vector r,v;
01110      _orbit.RelativePosVel(r,v);
01111      r += start_frame[0].position();
01112      v += start_frame[0].velocity();
01113      // cerr << "b.position: " << b.position.x << "  " << b.position.y << "  " << b.position.z << endl;
01114      b.SetPosition(r);
01115      b.SetVelocity(v);
01116      start_frame.push_back(b);
01117      
01118      frames.push_back(start_frame);
01119      }
01120      
01121      UniverseTypeAwareTime obs_utat;
01122      
01123      Radau15 itg; 
01124      itg.accuracy = 1.0e-12;
01125      itg.timestep = FromUnits(1,MINUTE);
01126      //
01127      Newton newton;
01128      
01129      Frame  last_frame;
01130      Vector relative_position;
01131      
01132      Evolution evol;
01133      
01134      evol.interaction = &newton;
01135      // evol.interaction = &itr;
01136      evol.integrator  = &itg;
01137      evol.push_back(start_frame);
01138      evol.SetMaxUnsavedSubSteps(100000);
01139      obs_utat.SetDate(final_time.GetDate()); 
01140      // evol.sample_period = orbit_epoch_frame.Time() - obs_utat.Time(); // important, otherwise the final frame is not saved
01141      // evol.sample_period = (orbit_epoch_frame.Time() - obs_utat.Time())*(1.0-1.0e-8); // important, otherwise the final frame is not saved
01142      evol.sample_period = (start_frame.Time() - obs_utat.Time())*(1.0-1.0e-8); // important, otherwise the final frame is not saved
01143      evol.Integrate(obs_utat.GetTime());
01144      last_frame = evol[evol.size()-1];
01145      
01146      frames.push_back(last_frame);
01147      
01148      Vector geo = last_frame[3].position()-last_frame[0].position(); // earth center position
01149      // Vector geo = frame[3].position; // earth center position
01150      
01151      Config conf;
01152      OrsaConfigFile ocf(&conf);
01153      ocf.Read();
01154      ocf.Close();
01155      
01156      LocationFile lf;
01157      lf.SetFileName(conf.paths[OBSCODE]->GetValue().c_str());
01158      lf.Read();
01159      lf.Close();
01160      
01161      // IMPORTANT!!!
01162      geo += lf.ObsPos(obscode,final_time.GetDate());
01163      
01164      relative_position = last_frame[last_frame.size()-1].position() - last_frame[0].position() - geo; // object - (earth + observer)
01165      
01166      Sky sky;
01167      // sky.Compute(relative_position,obs_utat);
01168      sky.Compute(relative_position,last_frame);
01169      
01170      return sky;
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      // OrbitWithCovarianceMatrix
01185      
01186      OrbitWithCovarianceMatrix::OrbitWithCovarianceMatrix(Orbit &nominal_orbit, double **covariance_matrix, CovarianceMatrixElements base) : OrbitWithEpoch(nominal_orbit), cov_base(base) {
01187      memcpy((void*)covm, (const void*)covariance_matrix, 6*6*sizeof(double)); // works?
01188      }
01189   
01190      void OrbitWithCovarianceMatrix::generate(vector<Orbit> &list, const int num) {
01191      
01192      list.clear();
01193      
01194      if (1) {
01195      // correlation matrix
01196      int i,k;
01197      for (i=0;i<6;++i) {
01198      for (k=0;k<6;++k) {
01199      // fprintf(stderr,"correlation[%i][%i]= %g\n",i,k,covm[i][k]/secure_sqrt(covm[i][i]*covm[k][k]));
01200      fprintf(stderr,"correlation[%i][%i]= %20.14f\n",i,k,covm[i][k]/secure_sqrt(covm[i][i]*covm[k][k]));
01201      }
01202      }
01203      }
01204      
01205      const size_t n=6;
01206      
01207      double meanv[n],parm[n*(n+3)/2+1],x[n],work[n][n];
01208      
01209      double *meanv_p = &meanv[0];
01210      double *covm_p  =  &covm[0][0];
01211      double *parm_p  =  &parm[0];
01212      double *x_p     =     &x[0];
01213      double *work_p  =  &work[0][0];
01214      
01215      switch (cov_base) {
01216      case Osculating: {
01217      meanv[0] = a;
01218      meanv[1] = e;
01219      meanv[2] = i;
01220      meanv[3] = omega_node;
01221      meanv[4] = omega_pericenter;
01222      meanv[5] = M;
01223      break;
01224      }
01225      case Equinoctal: {
01226      meanv[0] = a;
01227      meanv[1] = e*sin(omega_node+omega_pericenter);
01228      meanv[2] = e*cos(omega_node+omega_pericenter);
01229      meanv[3] = tan(i/2)*sin(omega_node);
01230      meanv[4] = tan(i/2)*cos(omega_node);
01231      meanv[5] = fmod(omega_node+omega_pericenter+M,twopi);
01232      break;
01233      }
01234      }
01235      
01236      // init
01237      setgmn(meanv_p,covm_p,n,parm_p);
01238      
01239      Orbit gen_orb = *this;
01240      int generated=0;
01241      
01242     switch (cov_base) {
01243     case Osculating: {
01244     while (1) {
01245     genmn(parm_p,x_p,work_p);
01246     ++generated;
01247     
01248     gen_orb.a                 = x[0];
01249     gen_orb.e                 = x[1];
01250     gen_orb.i                 = x[2];
01251     gen_orb.omega_node        = x[3];
01252     gen_orb.omega_pericenter  = x[4];
01253     gen_orb.M                 = x[5];
01254     
01255     list.push_back(gen_orb);
01256     
01257     if (generated==num) break;
01258     }
01259     break;
01260     }
01261     case Equinoctal: {
01262     double omega_tilde;
01263     while (1) {
01264     genmn(parm_p,x_p,work_p);
01265     ++generated;
01266     
01267     omega_tilde               = secure_atan2(x[1],x[2]);
01268     
01269     gen_orb.a                 = x[0];
01270     gen_orb.e                 = secure_sqrt(x[1]*x[1]+x[2]*x[2]);
01271     gen_orb.i                 = 2*atan(secure_sqrt(x[3]*x[3]+x[4]*x[4]));
01272     gen_orb.omega_node        = fmod(10*twopi+secure_atan2(x[3],x[4]),twopi); // works?
01273     gen_orb.omega_pericenter  = fmod(10*twopi+omega_tilde-omega_node,twopi);
01274     gen_orb.M                 = fmod(10*twopi+x[5]-omega_tilde,twopi);
01275     
01276     list.push_back(gen_orb);
01277     
01278     if (generated==num) break;
01279     }
01280     break;
01281     }
01282     }
01283     }
01284   */
01285   
01286   // Sky
01287   
01288   /* 
01289      void Sky::Compute(const Vector &relative_position, const UniverseTypeAwareTime& epoch) {
01290      // NOTE: relative_position = object - observer
01291      
01292      // ra.SetRad(fmod(secure_atan2(relative_position.y,relative_position.x)+twopi,twopi));
01293      // dec.SetRad(pi/2 - secure_acos(relative_position.z/Length(relative_position)));
01294      
01295      switch (universe->GetReferenceSystem()) {
01296      case ECLIPTIC: {
01297      Vector r = relative_position;
01298      Angle obl = obleq(epoch);
01299      r.rotate(0,obl.GetRad(),0);
01300      _ra.SetRad(fmod(secure_atan2(r.y,r.x)+twopi,twopi));
01301      _dec.SetRad(pi/2 - secure_acos(r.z/r.Length()));
01302      break;
01303      }
01304      case EQUATORIAL: {
01305      _ra.SetRad(fmod(secure_atan2(relative_position.y,relative_position.x)+twopi,twopi));
01306      _dec.SetRad(pi/2 - secure_acos(relative_position.z/relative_position.Length()));
01307      break;
01308      }
01309      }
01310      }
01311   */
01312   
01313   void Sky::Compute_J2000(const Vector &relative_position) {
01314     // NOTE: relative_position = object - observer
01315     switch (universe->GetReferenceSystem()) {
01316     case ECLIPTIC: {
01317       Vector r = relative_position;
01318       // Angle obl = obleq_J2000();
01319       // r.rotate(0,obl.GetRad(),0);
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      double Sky::delta_arcsec(const Observation & obs) const {
01351      
01352      // return (180/pi)*3600*(secure_sqrt(secure_pow(((obs.ra.GetRad()-_ra.GetRad())*cos(obs.dec.GetRad())),2)+secure_pow((obs.dec.GetRad()-_dec.GetRad()),2)));
01353      
01354      // fprintf(stderr,"SKY ra: %16.10f    dec: %16.10f\n",ra().GetRad(),dec().GetRad());
01355      // fprintf(stderr,"OBS ra: %16.10f    dec: %16.10f\n",obs.ra.GetRad(),obs.dec.GetRad());
01356      
01357      double d_ra = fabs(obs.ra.GetRad()-_ra.GetRad());
01358      if (fabs(obs.ra.GetRad()-_ra.GetRad()+twopi) < d_ra) d_ra = fabs(obs.ra.GetRad()-_ra.GetRad()+twopi);
01359      if (fabs(obs.ra.GetRad()-_ra.GetRad()-twopi) < d_ra) d_ra = fabs(obs.ra.GetRad()-_ra.GetRad()-twopi);
01360      
01361      double d_dec = fabs(obs.dec.GetRad()-_dec.GetRad());
01362      
01363      // const double delta = (180/pi)*3600.0*secure_sqrt(secure_pow(d_ra*cos(obs.dec.GetRad()),2)+secure_pow(d_dec,2));
01364      const double delta = (180/pi)*3600.0*secure_sqrt((d_ra*cos(obs.dec.GetRad()))*(d_ra*cos(obs.dec.GetRad()))+d_dec*d_dec);
01365      
01366      // cerr << "delta_arcsec: " << delta << endl;
01367      
01368      // return (180/pi)*3600.0*secure_sqrt(secure_pow(d_ra*cos(obs.dec.GetRad()),2)+secure_pow(d_dec,2));
01369      return delta;
01370      }
01371   */
01372   
01373 } // namespace orsa

Generated on Fri Nov 3 20:37:42 2006 for liborsa by  doxygen 1.4.7