00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_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
00039 Units * units = 0;
00040
00041
00042 Universe * universe = 0;
00043
00044
00045 Config * config = 0;
00046
00047
00048 JPLFile * jpl_file = 0;
00049
00050
00051 JPLCache * jpl_cache = 0;
00052
00053
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
00106
00107
00108
00109
00110
00111
00112
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
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
00131
00132
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
00148
00149
00150
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
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
00230
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
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
00254 int sign;
00255
00256
00257 if ((time_stop - f_stop) > 0) {
00258
00259 integrator->timestep = integrator->timestep.absolute();
00260 sign = +1;
00261 } else {
00262
00263 integrator->timestep = -integrator->timestep.absolute();
00264 sign = -1;
00265 }
00266
00267
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
00284
00285
00286
00287
00288
00289
00290
00291
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
00330
00331 if ( (f_start - (*this)[size()-1] + integrator->timestep).absolute() > sample_period) {
00332 if (integrator->timestep > 0.0) last_genuine_timestep = integrator->timestep;
00333
00334 integrator->timestep = (*this)[size()-1] - f_start + sign*sample_period;
00335
00336
00337
00338
00339
00340
00341
00342
00343 }
00344
00345
00346
00347 if (integrator->timestep > 0) {
00348
00349
00350 if ( (f_start+integrator->timestep) > time_stop) {
00351
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
00357
00358 if ( (f_start+integrator->timestep) < time_stop) {
00359
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
00366 }
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
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
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
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
00445
00446
00447
00448
00449
00450
00451
00452
00453 switch (universe->GetUniverseType()) {
00454 case Real:
00455
00456
00457 if ((f_stop - time_stop).IsZero()) {
00458
00459 continue_integration = false;
00460 if (save_last_anyway) save_frame = true;
00461 } else {
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
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
00522
00523
00524
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
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
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
00584 Frame StartFrame(const vector<BodyWithEpoch> & f, vector<JPL_planets> & jf, const Interaction * itr, const Integrator * itg, const UniverseTypeAwareTime & t) {
00585
00586
00587
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
00598
00599 running_evolution.SetIntegrator(itg);
00600 running_evolution.SetInteraction(itr);
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
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
00636 running_frame.SetDate(f[j].GetEpoch().GetDate());
00637
00638 running_frame.push_back(f[j]);
00639 ++j;
00640
00641
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
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
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
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
00681
00682
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 }