orsa_universe.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_universe.h"
00026 #include "orsa_file.h"
00027 #include "orsa_error.h"
00028 
00029 #include <iostream>
00030 #include <string>
00031 #include <limits>
00032 #include <algorithm>
00033 
00034 using namespace std;
00035 
00036 namespace orsa {
00037   
00038   // global units
00039   Units * units = 0;
00040   
00041   // active universe
00042   Universe * universe = 0;
00043   
00044   // active configuration
00045   Config * config = 0;
00046   
00047   // global JPL file
00048   JPLFile * jpl_file = 0;
00049   
00050   // JPL cache global class
00051   JPLCache * jpl_cache = 0;
00052   
00053   // world visible LocationFile
00054   LocationFile * location_file = 0;
00055   
00056   Universe::Universe() : std::vector<Evolution*>(), type(Simulated), sys(ECLIPTIC), timescale(ET) {
00057     common_init(AU,MSUN,YEAR);
00058   }
00059   
00060   Universe::Universe(length_unit lu, mass_unit mu, time_unit tu, UniverseType ut, ReferenceSystem rs, TimeScale ts) : std::vector<Evolution*>(), type(ut), sys(rs), timescale(ts) {
00061     common_init(lu,mu,tu);
00062   }
00063   
00064   void Universe::common_init(const length_unit lu, const mass_unit mu, const time_unit tu) {
00065     if (universe) delete universe; 
00066     universe = 0;
00067     
00068     if (!orsa_paths) orsa_paths = new OrsaPaths;
00069     Debug::construct();
00070     //
00071     if (!config) {
00072       config = new Config;
00073     }
00074     config->read_from_file();
00075     //
00076     if (!units) {
00077       units = new Units;
00078     }
00079     //
00080     units->SetSystem(tu,lu,mu);
00081     //
00082     if (!jpl_file) {
00083       jpl_file = new JPLFile(config->paths[JPL_EPHEM_FILE]->GetValue().c_str());
00084     }
00085     //
00086     if (!jpl_cache) {
00087       jpl_cache = new JPLCache;
00088     }
00089     //
00090     if (!location_file) {
00091       location_file = new LocationFile; 
00092       location_file->SetFileName(config->paths[OBSCODE]->GetValue().c_str());
00093       location_file->Open();
00094       location_file->Read();
00095       location_file->Close();
00096     }
00097     //
00098     modified = true;
00099     default_Date_timescale = timescale;
00100     universe = this;
00101   }
00102   
00103   Universe::~Universe() {
00104     /* 
00105        int k;
00106        k = size();
00107        while (k>0) {
00108        --k;
00109        if ((*this)[k]) (*this)[k]->clear();
00110        }
00111     */
00112     // NOTE: keep these two loops separate!
00113     int k = size();
00114     while (k>0) {
00115       --k;
00116       delete (*this)[k];
00117       (*this)[k] = 0;
00118     }
00119     
00120     universe = 0;
00121   }
00122   
00123   // declared in orsa_units.h -- to put elsewhere in future
00124   double GetG()        { return (orsa::units->GetG()); }
00125   double GetG_MKS()    { return (orsa::units->GetG_MKS()); }
00126   double GetMSun()     { return (orsa::units->GetMSun()); }
00127   double GetC()        { return (orsa::units->GetC()); }
00128   
00129   /* 
00130      double FromUnits(double x, time_unit   u, double power) { return orsa::units->FromUnits(x,u,power); };
00131      double FromUnits(double x, length_unit u, double power) { return orsa::units->FromUnits(x,u,power); };
00132      double FromUnits(double x, mass_unit   u, double power) { return orsa::units->FromUnits(x,u,power); };
00133   */
00134   
00135   double FromUnits(const double x, const time_unit   u, const int power) { return orsa::units->FromUnits(x,u,power); };
00136   double FromUnits(const double x, const length_unit u, const int power) { return orsa::units->FromUnits(x,u,power); };
00137   double FromUnits(const double x, const mass_unit   u, const int power) { return orsa::units->FromUnits(x,u,power); };
00138   
00139   string TimeLabel()   { return orsa::units->TimeLabel();   };
00140   string LengthLabel() { return orsa::units->LengthLabel(); };
00141   string MassLabel()   { return orsa::units->MassLabel();   };
00142   
00143   unsigned int Evolution::used_evolution_id = 0;
00144   
00145   Evolution::Evolution() : std::vector<Frame>(), id(used_evolution_id++) {
00146     /* 
00147        if (universe->GetUniverseType() == Real) {
00148        sample_period = FromUnits(1,DAY);
00149        } else {
00150        sample_period = FromUnits(0.01,YEAR);
00151        }
00152     */
00153     sample_period = FromUnits(0.01,YEAR);
00154     integrator    = 0;
00155     interaction   = 0;
00156     max_unsaved_substeps_active = false;
00157     _integrating=false;
00158   }
00159   
00160   //! This copy constructor copies the parameters, not the frames!
00161   Evolution::Evolution(const Evolution & e) : std::vector<Frame>(), id(used_evolution_id++) {
00162     sample_period = e.sample_period;
00163     integrator    = e.integrator->clone();
00164     interaction   = e.interaction->clone();
00165     //
00166     max_unsaved_substeps_active = false;
00167     _integrating=false;
00168   }
00169   
00170   Evolution::~Evolution() {
00171     delete integrator;
00172     integrator = 0;
00173     delete interaction;
00174     interaction = 0;
00175   }
00176   
00177   void Evolution::SetIntegrator(const IntegratorType type) {
00178     make_new_integrator(&integrator,type);
00179   }
00180   
00181   void Evolution::SetIntegrator(const Integrator * itg) {
00182     delete integrator;
00183     integrator = itg->clone();
00184   }
00185   
00186   void Evolution::SetIntegratorTimeStep(const UniverseTypeAwareTimeStep ts) {
00187     integrator->timestep = ts;
00188   }
00189   
00190   const UniverseTypeAwareTimeStep & Evolution::GetIntegratorTimeStep() const {
00191     return integrator->timestep;
00192   }
00193   
00194   void Evolution::SetIntegratorAccuracy(const double a) {
00195     integrator->accuracy = a;
00196   }
00197   
00198   double Evolution::GetIntegratorAccuracy() const {
00199     return integrator->accuracy;
00200   }
00201   
00202   void Evolution::SetInteraction(const InteractionType type) {
00203     make_new_interaction(&interaction,type);
00204   }
00205   
00206   void Evolution::SetInteraction(const Interaction * itr) {
00207     delete interaction;
00208     interaction = itr->clone();
00209   }
00210   
00211   void Evolution::SetSamplePeriod(const UniverseTypeAwareTimeStep & sp) {
00212     sample_period = sp;
00213   }
00214   
00215   void Evolution::Integrate(const UniverseTypeAwareTime & time_stop, const bool save_last_anyway) {
00216     
00217     if (integrator==0) {
00218       ORSA_WARNING("Integrator not initialized");
00219       return;
00220     }
00221     
00222     if (interaction==0) {
00223       ORSA_WARNING("Interaction not initialized!");
00224       return;
00225     }
00226     
00227     integration_started();
00228     
00229     // FIXME: if the timestep is not set?
00230     // if (sample_period == 0) sample_period = integrator->timestep;
00231     if (sample_period.IsZero()) sample_period = integrator->timestep;
00232     sample_period = sample_period.absolute();
00233     
00234     if (size() == 0) {
00235       ORSA_ERROR("no starting frame in integration.");
00236       return;
00237     }
00238     
00239     Frame f_start = (*this)[size()-1];
00240     Frame f_stop  = f_start;
00241     
00242     if (f_start == time_stop) {
00243       // no integration needed...
00244       return;
00245     }
00246     
00247     const double time_start = f_start.Time();
00248     
00249     bool save_frame = false;
00250     
00251     bool continue_integration = true;
00252     
00253     // double sign;
00254     int sign;
00255     
00256     // if ((time_stop.Time()-f_stop.Time()) > 0) {
00257     if ((time_stop - f_stop) > 0) {
00258       // integrator->timestep = fabs(integrator->timestep);
00259       integrator->timestep = integrator->timestep.absolute();
00260       sign = +1;
00261     } else {
00262       // integrator->timestep = -fabs(integrator->timestep);
00263       integrator->timestep = -integrator->timestep.absolute();
00264       sign = -1;
00265     }
00266     
00267     // important init value!
00268     UniverseTypeAwareTimeStep last_genuine_timestep;
00269     if (integrator->timestep.absolute() > 0.0) {
00270       last_genuine_timestep = integrator->timestep;
00271     } else {
00272       last_genuine_timestep = time_stop - f_start;
00273     }
00274     
00275     unsigned int unsaved_substeps = 0;
00276     unsigned int total_substeps = 0;
00277     
00278     while ( (continue_integration) && (((time_stop.Time()-f_stop.Time())/(time_stop.Time()-time_start)) > 0) ) {
00279       
00280       f_start = f_stop;
00281       
00282       /* 
00283          fprintf(stderr,
00284          "==========================         \n"
00285          "f_start.Time()           : %20.12e \n"
00286          "integrator->timestep     : %20.12e \n"
00287          "(*this)[size()-1].Time() : %20.12e \n"
00288          "sample_period            : %20.12e \n"
00289          "frames                   : %zi     \n"
00290          "==========================         \n",
00291          f_start.Time(),integrator->timestep,(*this)[size()-1].Time(),sample_period,size()
00292          );
00293       */
00294       
00295       if (0) { 
00296         ORSA_DEBUG(
00297                 "\n"
00298                 "==========================         \n"
00299                 "f_start.Time()           : %20.12e \n"
00300                 "time_stop.Time()         : %20.12e \n"
00301                 "integrator->timestep     : %20.12e \n"
00302                 "(*this)[size()-1].Time() : %20.12e \n"
00303                 "sample_period            : %20.12e \n"
00304                 "frames                   : %zi     \n"
00305                 "==========================         \n",
00306                 f_start.Time(),time_stop.Time(),integrator->timestep.GetDouble(),
00307                 (*this)[size()-1].Time(),sample_period.GetDouble(),size());
00308         //
00309         if (universe->GetUniverseType() == Real) {
00310           int y,m,d;
00311           double frac;
00312           // 
00313           f_start.GetDate().GetGregor(y,m,d);
00314           frac = f_start.GetDate().GetDayFraction();
00315           ORSA_ERROR(
00316                   "f_start date:    %i/%02i/%013.10f",
00317                   y,m,d+frac);
00318           //
00319           time_stop.GetDate().GetGregor(y,m,d);
00320           frac = time_stop.GetDate().GetDayFraction();
00321           ORSA_ERROR(
00322                   "time_stop date:  %i/%02i/%013.10f",
00323                   y,m,d+frac);
00324         }
00325       }
00326       
00327       if (integrator->timestep.IsZero()) integrator->timestep = last_genuine_timestep;
00328       
00329       // if ( fabs(f_start.Time() + integrator->timestep - (*this)[size()-1].Time()) > sample_period)  {
00330       // if ( (f_start + integrator->timestep - (*this)[size()-1]).absolute() > sample_period) {
00331       if ( (f_start - (*this)[size()-1] + integrator->timestep).absolute() > sample_period) {
00332         if (integrator->timestep > 0.0) last_genuine_timestep = integrator->timestep;
00333         // integrator->timestep = (*this)[size()-1].Time() + sign*sample_period - f_start.Time();
00334         integrator->timestep = (*this)[size()-1] - f_start + sign*sample_period;
00335         /* 
00336            fprintf(stderr,
00337            "timestep changed: from %.20e to %.20e      val: %.20e sp: %.20e\n",
00338            last_genuine_timestep.GetDouble(),
00339            integrator->timestep.GetDouble(),
00340            (f_start + integrator->timestep - (*this)[size()-1]).absolute().GetDouble(),
00341            sample_period.GetDouble());
00342         */
00343       }
00344       
00345       // cerr << "timestep (tmp):    " << integrator->timestep.GetDouble() << endl;
00346       
00347       if (integrator->timestep > 0) {
00348         // cerr << "timestep > 0 !!!" << endl;
00349         // if ( (f_start.Time()+integrator->timestep) > (time_stop.Time()) ) {
00350         if ( (f_start+integrator->timestep) > time_stop) {
00351           // integrator->timestep = time_stop.Time() - f_start.Time();
00352           if (integrator->timestep > 0.0) last_genuine_timestep = integrator->timestep;
00353           integrator->timestep = time_stop - f_start;
00354         }
00355       } else if (integrator->timestep < 0) {
00356         // cerr << "timestep < 0 !!!" << endl;
00357         // if ( (f_start.Time()+integrator->timestep) < (time_stop.Time()) ) {
00358         if ( (f_start+integrator->timestep) < time_stop) {
00359           // integrator->timestep = time_stop.Time() - f_start.Time();
00360           if (integrator->timestep > 0.0) last_genuine_timestep = integrator->timestep;
00361           integrator->timestep = time_stop - f_start;
00362         }
00363       } else {
00364         ORSA_WARNING("Timestep equal to zero!");
00365         // timestep == 0...
00366       }
00367       
00368       // cerr << "timestep (tmp-II): " << integrator->timestep.GetDouble() << endl;
00369       
00370       /* 
00371          if (f_start == time_stop) {
00372          return;
00373          }
00374       */
00375       
00376       // debug
00377       /* 
00378          if (integrator->timestep.absolute().GetDouble() < FromUnits(1.0,SECOND)) {
00379          ORSA_ERROR("timestep smaller than 1.0 seconds!");
00380          // sleep(1.0);
00381          }
00382       */
00383       
00384       // step
00385       if (integrator->timestep.IsZero()) {
00386         ORSA_ERROR("zero timestep, no integrator->Step() call...");
00387         f_stop = f_start;
00388       } else {
00389         integrator->Step(f_start,f_stop,interaction);
00390         if (f_start.Time() != f_stop.Time()) {
00391           ++total_substeps;
00392         }
00393       }
00394       
00395       if (0) {
00396         if (universe->GetUniverseType() == Real) {
00397           int y,m,d;
00398           double frac;
00399           // 
00400           f_stop.GetDate().GetGregor(y,m,d);
00401           frac = f_stop.GetDate().GetDayFraction();
00402           ORSA_ERROR(
00403                   "f_stop date:     %i/%02i/%013.10f",
00404                   y,m,d+frac);
00405         }
00406       }
00407       
00408       /* 
00409          {
00410          // debug
00411          fprintf(stderr,"..after step: %g\n",((f_stop - (*this)[size()-1]).absolute() - sample_period).absolute().GetDouble());
00412          fprintf(stderr,"....timestep: %g\n",integrator->timestep.GetDouble());
00413          fprintf(stderr,"....time....: %g\n",f_stop.Time());
00414          }
00415       */
00416       
00417       /* 
00418          { 
00419          // debug
00420          fprintf(stderr,
00421          "save frame test: %.20e\n"
00422          "sign*sample_period: %.20e      (*this)[size()-1].GetTime(): %.9f\n",
00423          // ((f_stop-(*this)[size()-1]).absolute()-sample_period).absolute().GetDouble());
00424          ( (*this)[size()-1] - f_stop + sign*sample_period).GetDouble(),
00425          sign*sample_period.GetDouble(),
00426          (*this)[size()-1].GetTime());
00427          }
00428       */
00429       
00430       // save frame
00431       switch (universe->GetUniverseType()) {
00432       case Real:
00433         if ( ( (*this)[size()-1] - f_stop + sign*sample_period).IsZero() ) { 
00434           save_frame = true;
00435         }
00436         break;
00437       case Simulated:
00438         if ( fabs( fabs(f_stop.Time() - (*this)[size()-1].Time()) - sample_period.GetDouble()) < (fabs(f_stop.Time())+fabs((*this)[size()-1].Time())+fabs(sample_period.GetDouble()))*total_substeps*std::numeric_limits<double>::epsilon()) {
00439           save_frame = true;
00440         }
00441         break;
00442       }
00443       
00444       // end of integration
00445       /* if ( fabs(f_stop.Time() - time_stop.Time()) < fabs(time_stop.Time())*std::numeric_limits<double>::epsilon()) {
00446          continue_integration=false;
00447          if (save_last_anyway) save_frame=true;
00448          }
00449       */
00450       //
00451       // cerr << "end of integration checks..." << endl;
00452       // end of integration
00453       switch (universe->GetUniverseType()) {
00454       case Real:
00455         // if (f_stop == time_stop) {
00456         // if (f_stop.GetDate() == time_stop.GetDate()) {
00457         if ((f_stop - time_stop).IsZero()) { // IMPORTANT! This one works properly, while "(f_stop == time_stop)" or "(f_stop.GetDate() == time_stop.GetDate())" don't... WHY??
00458           // cerr << "(f_stop == time_stop)" << endl;
00459           continue_integration = false;
00460           if (save_last_anyway) save_frame = true;
00461         } else {
00462           /* 
00463              cerr << "(f_stop != time_stop)" << endl;
00464              //
00465              const Date fs_date = f_stop.GetDate();
00466              const Date ts_date = time_stop.GetDate();
00467              int y,m,d;
00468              double frac;
00469              //
00470              fs_date.GetGregor(y,m,d);
00471              frac = fs_date.GetDayFraction();
00472              ORSA_ERROR(
00473              "     fs_date:    %i/%02i/%013.10f",
00474              y,m,d+frac);
00475              //
00476              ts_date.GetGregor(y,m,d);
00477              frac = ts_date.GetDayFraction();
00478              ORSA_ERROR(
00479              "     ts_date:    %i/%02i/%013.10f",
00480              y,m,d+frac);
00481              //
00482              ORSA_ERROR(
00483              "   ts-fs: %013.10f seconds",
00484              FromUnits(ts_date.GetTime()-fs_date.GetTime(),SECOND,-1));
00485              //
00486              UniverseTypeAwareTime fsw = fs_date;
00487              UniverseTypeAwareTime tsw = ts_date;
00488              ORSA_ERROR(
00489              "2- ts-fs: %013.10f seconds",
00490              FromUnits((tsw-fsw).GetDouble(),SECOND,-1));
00491           */
00492         }
00493         break;
00494       case Simulated:
00495         if ( fabs(f_stop.Time() - time_stop.Time()) < fabs(time_stop.Time())*total_substeps*std::numeric_limits<double>::epsilon()) {
00496           continue_integration = false;
00497           if (save_last_anyway) save_frame = true;
00498         }
00499         break;
00500       }
00501       
00502       if (save_frame) {
00503         step_done(time_start,time_stop,last_genuine_timestep,f_stop,continue_integration);
00504         if (interaction->IsSkippingJPLPlanets()) {
00505           f_stop.ForceJPLEphemerisData();
00506         }
00507         push_back(f_stop);
00508         save_frame = false;
00509         integrator->timestep = last_genuine_timestep;
00510         unsaved_substeps = 0;
00511       } else {
00512         ++unsaved_substeps;
00513         if (max_unsaved_substeps_active) {
00514           if (unsaved_substeps > max_unsaved_substeps) {
00515             continue_integration=false;
00516             ORSA_WARNING("max unsaved substeps hit!!! stopping integration!");
00517           }
00518         }
00519       }
00520       
00521       // test
00522       /* 
00523          if (integrator->timestep.IsZero()) {
00524          continue_integration=false;
00525          }
00526       */
00527     }
00528     
00529     integration_finished();
00530   }
00531   
00532   void Evolution::Integrate(const Frame & f, const UniverseTypeAwareTime & utat_start, const UniverseTypeAwareTime& utat_stop) {
00533     
00534     // clear Evolution
00535     clear();
00536     
00537     UniverseTypeAwareTime t_start = utat_start;
00538     UniverseTypeAwareTime t_stop  = utat_stop;
00539     
00540     if (utat_start > utat_stop) {
00541       t_start = utat_stop;
00542       t_stop  = utat_start;
00543     }
00544     
00545     const UniverseTypeAwareTimeStep saved_sample_period = sample_period;
00546     const UniverseTypeAwareTimeStep huge_sample_period  = FromUnits(1e3,YEAR);
00547     
00548     Frame running_frame;
00549     
00550     if (t_start < f) {
00551       if (t_stop < f) {
00552         push_back(f);
00553         sample_period = huge_sample_period;
00554         Integrate(t_stop,true);
00555         running_frame = (*this)[size()-1];
00556         clear();
00557         push_back(running_frame);
00558         sample_period = saved_sample_period;
00559         Integrate(t_start);
00560       } else {
00561         push_back(f);
00562         Integrate(t_start);
00563         // swap
00564         (*this)[0] = (*this)[size()-1];
00565         (*this)[size()-1] = f;
00566         //
00567         Integrate(t_stop);
00568       }
00569     } else {
00570       push_back(f);
00571       sample_period = huge_sample_period;
00572       Integrate(t_start,true);
00573       running_frame = (*this)[size()-1];
00574       clear();
00575       push_back(running_frame);
00576       sample_period = saved_sample_period;
00577       Integrate(t_stop);
00578     }
00579     
00580     std::sort(begin(),end());
00581   }
00582   
00583   // Frame StartFrame(const vector<BodyWithEpoch> &f, vector<JPLBody> &jf, Interaction *itr, Integrator *itg, const UniverseTypeAwareTime &t) {
00584   Frame StartFrame(const vector<BodyWithEpoch> & f, vector<JPL_planets> & jf, const Interaction * itr, const Integrator * itg, const UniverseTypeAwareTime & t) {
00585   // Frame StartFrame(const vector<BodyWithEpoch> & f, vector<JPL_planets> & jf, InteractionType interaction_type, IntegratorType integrator_type, const UniverseTypeAwareTime & t) {
00586     
00587     // const double original_timestep = itg->timestep;
00588     const UniverseTypeAwareTimeStep original_timestep = itg->timestep;
00589     
00590     Frame frame_start;
00591     
00592     switch (universe->GetUniverseType()) {
00593     case Real: { 
00594       frame_start.SetDate(t.GetDate());
00595       Frame          running_frame;
00596       Evolution      running_evolution;
00597       // running_evolution.interaction = itr->clone();
00598       // running_evolution.integrator  = itg->clone();
00599       running_evolution.SetIntegrator(itg);
00600       running_evolution.SetInteraction(itr);
00601       
00602       // put the JPL bodies into the frame_start
00603       /* 
00604          if (jf.size()) {
00605          Date ddf = t.GetDate();
00606          unsigned int k;
00607          for (k=0;k<jf.size();++k) {
00608          jf[k].SetEpoch(ddf);
00609          Body b(jf[k].name(),
00610          jf[k].mass(),
00611          jf[k].radius(),
00612          jf[k].position(),
00613          jf[k].velocity());
00614          frame_start.push_back(b);
00615          }
00616          }
00617       */
00618       //
00619       if (jf.size()) {
00620         unsigned int k;
00621         for (k=0;k<jf.size();++k) {
00622           frame_start.push_back(JPLBody(jf[k],t.GetDate()));
00623         }
00624       }
00625       
00626       if (f.size()) {
00627         
00628         Frame last_frame;
00629         unsigned int base_bodies, l;
00630         
00631         unsigned int j=0;
00632         while (j < f.size()) {
00633           
00634           running_frame.clear();
00635           // running_frame.SetDate(f[j].epoch.GetDate());
00636           running_frame.SetDate(f[j].GetEpoch().GetDate());
00637           
00638           running_frame.push_back(f[j]);
00639           ++j;
00640           
00641           // while ((j<f.size()) && (f[j].epoch.GetDate() == running_frame.GetDate())) {
00642           while ((j<f.size()) && (f[j].Epoch().GetDate() == running_frame.GetDate())) {
00643             running_frame.push_back(f[j]);
00644             ++j;
00645           }
00646           
00647           base_bodies = running_frame.size();
00648           
00649           // JPL 
00650           /* 
00651              Date ddt = running_frame.GetDate();
00652              unsigned int k;
00653              for (k=0;k<jf.size();++k) {
00654              jf[k].SetEpoch(ddt);
00655              Body b(jf[k].name(),
00656              jf[k].mass(),
00657              jf[k].radius(),
00658              jf[k].position(),
00659              jf[k].velocity());
00660              running_frame.push_back(b);
00661              }
00662           */
00663           {
00664             unsigned int k;
00665             for (k=0;k<jf.size();++k) {
00666               running_frame.push_back(JPLBody(jf[k],running_frame.GetDate()));
00667             }
00668           }
00669           
00670           running_evolution.clear();
00671           running_evolution.push_back(running_frame);
00672           // running_evolution.sample_period = fabs(running_frame.GetTime() - frame_start.GetTime());
00673           running_evolution.SetSamplePeriod((running_frame - frame_start).absolute());
00674           running_evolution.SetIntegratorTimeStep(original_timestep);
00675           running_evolution.Integrate(t);
00676           
00677           last_frame = running_evolution[running_evolution.size()-1];
00678           if (fabs(last_frame.Time()-frame_start.Time()) > FromUnits(1.0e-3,SECOND)) {
00679             /* 
00680                fprintf(stderr,"!!! last_frame.Time() != frame_start.Time() --->> %30.20f != %30.20f\n",last_frame.Time(),frame_start.Time());
00681                cerr << "|dT| = " << FromUnits(fabs(last_frame.Time()-frame_start.Time()),SECOND,-1) << " seconds" << endl;
00682                print(last_frame);
00683             */
00684             ORSA_WARNING(
00685                     "last_frame.Time() != frame_start.Time() --->> %30.20f != %30.20f\n"
00686                     "|dT| = %g seconds\n",
00687                     last_frame.Time(),frame_start.Time(),FromUnits(fabs(last_frame.Time()-frame_start.Time()),SECOND,-1));
00688             ORSA_WARNING("  objects in frame: ");
00689             unsigned int k;
00690             for (k=0;k<last_frame.size();++k) {
00691               ORSA_WARNING("    %s", last_frame[k].name().c_str());
00692             }
00693           }
00694           
00695           for (l=0;l<base_bodies;++l) {
00696             frame_start.push_back(last_frame[l]);
00697           }
00698         }
00699       }
00700       break;
00701     }
00702     case Simulated: {
00703       frame_start.SetTime(0.0);
00704       unsigned int k;  
00705       for (k=0;k<f.size();++k) {
00706         frame_start.push_back(f[k]);
00707       }
00708       break;
00709     }
00710     }
00711     
00712     return frame_start;
00713   }
00714   
00715 } // namespace orsa

Generated on Wed May 30 13:04:53 2007 for liborsa by  doxygen 1.5.2