orsa_file.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_file.h"
00026 
00027 #include <algorithm>
00028 #include <iostream>
00029 
00030 #include <ctype.h>
00031 
00032 #include "orsa_common.h"
00033 #include "orsa_units.h"
00034 #include "orsa_error.h"
00035 #include "orsa_version.h"
00036 
00037 #include "sdncal.h"
00038 
00039 #include <cstdio>
00040 #include "jpleph.h"
00041 #include "jpl_int.h"
00042 
00043 #ifndef _WIN32
00044 #include <unistd.h> // sleep()
00045 #else
00046 #include <shlobj.h>
00047 #include <direct.h>
00048 #endif
00049 
00050 #include "support.h"
00051 
00052 #ifdef HAVE_CONFIG_H
00053 #include "config.h"
00054 #endif // HAVE_CONFIG_H
00055 
00056 using namespace std;
00057 
00058 namespace orsa {
00059   
00060   void ReadFile::Open() {
00061     if (status != CLOSE) return;
00062     
00063     file = OPEN_FILE(filename.c_str(),OPEN_READ);
00064     
00065     if (file == 0) { 
00066       ORSA_ERROR("Can't open file %s",filename.c_str());
00067     } else {
00068       status = OPEN_R;
00069     }
00070   }
00071   
00072   void WriteFile::Open() {
00073     if (status != CLOSE) return;
00074     
00075     file = OPEN_FILE(filename.c_str(),OPEN_WRITE);
00076     
00077     if (file == 0) { 
00078       ORSA_ERROR("Can't open file %s",filename.c_str());
00079     } else {
00080       status = OPEN_W;
00081     }
00082   }
00083   
00084   void ReadWriteFile::Open(const FILE_STATUS st) {
00085     
00086     // already in the right status
00087     if (status == st) return;
00088     
00089     // anomalous...
00090     if (st == CLOSE) {
00091       Close();
00092       return;
00093     }
00094     
00095     Close();
00096     
00097     if ((st == OPEN_R) && ((file = OPEN_FILE(filename.c_str(),OPEN_READ)) != 0)) {
00098       status = OPEN_R;
00099       return;
00100     }
00101     
00102     if ((st == OPEN_W) && ((file = OPEN_FILE(filename.c_str(),OPEN_WRITE)) != 0)) {
00103       status = OPEN_W;
00104       return;
00105     }
00106     
00107     if (file == 0) {
00108       ORSA_ERROR("Can't open file %s",filename.c_str());
00109     }
00110     
00111     status = CLOSE;
00112   }
00113     
00114     
00115   void File::Close() {
00116     if (status != CLOSE) {
00117       CLOSE_FILE(file);
00118       status = CLOSE;
00119     }
00120   }
00121   
00122   // AsteroidDatabaseFile
00123   
00124   // a copy of ReadFile::Open()
00125   /* 
00126      void AsteroidDatabaseFile::Open() {
00127      if (status != CLOSE) return;
00128      
00129      file = OPEN_FILE(filename.c_str(),OPEN_READ);
00130      
00131      if (file == 0) { 
00132      cerr << "Can't open file " << filename << endl;
00133      } else {
00134      status = OPEN_R;
00135      }
00136      }
00137   */
00138   
00139   // AstorbFile
00140   
00141   /* 
00142      int AstorbFile::GetSize() {
00143      
00144      if (status == CLOSE) Open();
00145      
00146      char line[300];
00147      
00148      REWIND_FILE(file);
00149      
00150      int lines = 0;
00151      // while ( (fgets(line,300,file)) != 0 ) {
00152      while ((GETS_FILE(line,300,file)) != 0) {
00153      lines++;
00154      }
00155      
00156      REWIND_FILE(file);
00157      
00158      return lines;
00159      }
00160   */
00161  
00162   /* 
00163      Date AstorbFile::GetEpoch() {
00164      
00165      if (status == CLOSE) Open();
00166      
00167      char line[1024];
00168      
00169      string epoch;
00170      
00171      REWIND_FILE(file);
00172      
00173      GETS_FILE(line,1024,file);
00174      while (strlen(line) < 100) {
00175      GETS_FILE(line,1024,file);
00176      }      
00177      
00178      epoch.assign(line,105,8);
00179      
00180      Date date(TDT);
00181      
00182      string year,month,day;
00183      int    y,m,d;
00184      
00185      year.assign(epoch,0,4);
00186      month.assign(epoch,4,2);
00187      day.assign(epoch,6,2);
00188      
00189      y = atoi(year.c_str());
00190      m = atoi(month.c_str());
00191      d = atoi(day.c_str());
00192      
00193      date.SetGregor(y,m,d);
00194      
00195      return (date);
00196      }
00197   */
00198   
00199   void AstorbFile::Read() {
00200     
00201     // if (db_updated) return;
00202     
00203     // if (status == CLOSE) Open();
00204     
00205     Open();
00206     
00207     if (status != OPEN_R) {
00208       ORSA_ERROR("Status error!");
00209       return;
00210     }
00211     
00212     // reset database 
00213     // cerr << "pre-resize..." << endl;
00214     db->clear();
00215     
00216     // db->reserve(GetSize()); // a little too time consuming...
00217     
00218     // cerr << "reading pass (1)" << endl;
00219     
00220     // int print_step = 0;
00221     // int print_step = 10000;
00222     
00223     char line[300];
00224     
00225     double a,e,i,omega_node,omega_pericenter,M;
00226     // int    n;
00227     string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
00228     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
00229     // string ceu;
00230     
00231     string year,month,day;
00232     int    y,m,d;
00233     
00234     // Body orb_ref_body("",GetMSun());
00235     
00236     // AstorbDataEntry ast;
00237     Asteroid ast;
00238     
00239     // Date tmp_date(TDT);
00240     Date tmp_date;
00241     
00242     unsigned int local_index = 0;
00243     bool bool_stop=false;
00244     bool bool_pause=false;
00245     REWIND_FILE(file);
00246     while ((GETS_FILE(line,300,file)) != 0) {
00247       
00248       // cerr << "AstorbFile::Read() inside the while() loop..." << endl;
00249       
00250       if (strlen(line) < 100) continue; // not a good line, maybe a comment or a white line...
00251       
00252       // cerr << "inside while loop..." << endl;
00253       
00254       local_index++;
00255       read_progress(local_index,bool_pause,bool_stop);
00256       
00257       if (bool_stop) break;
00258       
00259       while (bool_pause) {
00260         // cerr << "AstorbFile::Read() sleeping..." << endl;
00261         sleep(1);
00262         read_progress(local_index,bool_pause,bool_stop);
00263       }
00264       
00265       // uncomment the ones used
00266       number.assign(line,0,5);
00267       name.assign(line,6,18);
00268       orbit_computer.assign(line,25,15);
00269       absolute_magnitude.assign(line,41,5);
00270       //
00271       arc.assign(line,94,5);
00272       numobs.assign(line,99,5);
00273       //
00274       epoch.assign(line,105,8);
00275       //
00276       mean_anomaly.assign(line,114,10);
00277       pericenter.assign(line,125,10);
00278       node.assign(line,136,10);
00279       inclination.assign(line,146,10);
00280       eccentricity.assign(line,157,10);
00281       semimajor_axis.assign(line,168,12);
00282       // ceu.assign(line,190,7);
00283       
00284       //////////////
00285       
00286       ast.n = atoi(number.c_str());
00287       
00288       ast.name = name;
00289       remove_leading_trailing_spaces(ast.name);
00290       
00291       // ast.ceu  = atof(ceu.c_str());
00292       
00293       ast.mag  = atof(absolute_magnitude.c_str());
00294       
00295       a                = atof(semimajor_axis.c_str());
00296       e                = atof(eccentricity.c_str());
00297       i                = (pi/180)*atof(inclination.c_str());
00298       omega_node       = (pi/180)*atof(node.c_str());
00299       omega_pericenter = (pi/180)*atof(pericenter.c_str());
00300       M                = (pi/180)*atof(mean_anomaly.c_str());
00301       
00302       ast.orb.a                = FromUnits(a,AU);
00303       ast.orb.e                = e;
00304       ast.orb.i                = i;
00305       ast.orb.omega_node       = omega_node;
00306       ast.orb.omega_pericenter = omega_pericenter;
00307       ast.orb.M                = M;
00308       
00309       year.assign(epoch,0,4);
00310       month.assign(epoch,4,2);
00311       day.assign(epoch,6,2);
00312       
00313       y = atoi(year.c_str());
00314       m = atoi(month.c_str());
00315       d = atoi(day.c_str());
00316       
00317       tmp_date.SetGregor(y,m,d,TDT);
00318       ast.orb.epoch.SetDate(tmp_date);
00319       // ast.orb.T = sqrt(4*pisq/(GetG()*GetMSun())*pow(FromUnits(ast.orb.a,AU),3));
00320       ast.orb.mu = GetG()*GetMSun();
00321       // ast.orb.ref_body = orb_ref_body;
00322       
00323       /* 
00324          switch (universe->GetReferenceSystem()) {
00325          case ECLIPTIC: break;
00326          case EQUATORIAL:
00327          { 
00328          // cerr << "Rotating astorb orbit..." << endl;
00329          const double obleq_rad = obleq(tmp_date).GetRad();
00330          Vector position,velocity;
00331          ast.orb.RelativePosVel(position,velocity);
00332          position.rotate(0.0,obleq_rad,0.0);
00333          velocity.rotate(0.0,obleq_rad,0.0);
00334          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
00335          }
00336          break;
00337          }
00338       */
00339       
00340       switch (universe->GetReferenceSystem()) {
00341       case ECLIPTIC: break;
00342       case EQUATORIAL:
00343         { 
00344           Vector position,velocity;
00345           ast.orb.RelativePosVel(position,velocity);
00346           EclipticToEquatorial_J2000(position);
00347           EclipticToEquatorial_J2000(velocity);
00348           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
00349         }
00350         break;
00351       }
00352       
00353       db->push_back(ast);
00354       
00355     }
00356     
00357     read_finished();
00358     
00359   }
00360   
00361   /* 
00362      bool AstorbFile::GuessCorrectType() {
00363      return true;
00364      }  
00365   */
00366   
00367   // AstorbFile::AstorbFile() : ReadFile() {
00368   AstorbFile::AstorbFile() : AsteroidDatabaseFile() {
00369     // status = CLOSE;
00370     // db = new AstorbDatabase();
00371     db = new AsteroidDatabase();
00372   }
00373   
00374   AstorbFile::~AstorbFile() {
00375     delete db;
00376     db = 0;
00377   }
00378   
00379   // MPC orbit database
00380   MPCOrbFile::MPCOrbFile() : AsteroidDatabaseFile() {
00381     db = new AsteroidDatabase();
00382   }
00383   
00384   MPCOrbFile::~MPCOrbFile() {
00385     delete db;
00386     db = 0;
00387   }
00388   
00389   void MPCOrbFile::Read() {
00390     
00391     Open();
00392     
00393     if (status != OPEN_R) {
00394       ORSA_ERROR("Status error!");
00395       return;
00396     }
00397     
00398     db->clear();
00399     
00400     // Body orb_ref_body("",GetMSun());
00401     
00402     // unsigned int local_index;
00403     
00404     char line[300];
00405     
00406     double a,e,i,omega_node,omega_pericenter,M;
00407     // int    n;
00408     string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
00409     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
00410     // string ceu;
00411     
00412     string year,month,day;
00413     int    y=0,m=0,d=0;
00414     
00415     // AstorbDataEntry ast;
00416     Asteroid ast;
00417     
00418     unsigned int local_index = 0;
00419     bool bool_stop=false;
00420     bool bool_pause=false;
00421 
00422     REWIND_FILE(file);
00423     
00424     // Date tmp_date(TDT);
00425     Date tmp_date;
00426     
00427     /* // skip header, not needed with the check on (strlen(line) < 201)...
00428        do { 
00429        GETS_FILE(line,300,file);
00430        } while (line[0] != '-');
00431        cerr << "below the header..." << endl;
00432     */
00433     
00434     while ( (GETS_FILE(line,300,file)) != 0 ) {
00435       
00436       // cerr << "line: -->" << line << "<--" << endl; 
00437       
00438       ++local_index;
00439       
00440       if (strlen(line) < 201) continue; // not a good line, maybe a comment or a white line...
00441       
00442       // cerr << "strlen(line): " << strlen(line) << endl;
00443       
00444       // cerr << "line: -->" << line << "<--" << endl; 
00445       
00446       read_progress(local_index,bool_pause,bool_stop);
00447       
00448       if (bool_stop) break;
00449       
00450       while (bool_pause) {
00451         sleep(1);
00452         read_progress(local_index,bool_pause,bool_stop);
00453       }
00454       
00455       // uncomment the ones used
00456       number.assign(line,0,7);
00457       absolute_magnitude.assign(line,8,5);
00458       epoch.assign(line,20,5);
00459       //
00460       mean_anomaly.assign(line,26,9);
00461       pericenter.assign(line,37,9);
00462       node.assign(line,48,9);
00463       inclination.assign(line,59,9);
00464       eccentricity.assign(line,70,9);
00465       semimajor_axis.assign(line,92,11);
00466       //
00467       numobs.assign(line,117,5);
00468       orbit_computer.assign(line,150,10);
00469       
00470       if (strlen(line) > 160) {
00471         name.assign(line,166,28);
00472       } else {
00473         name = "";      
00474       }
00475       
00476       // conversions
00477       
00478       {
00479         // remove -->(NUMBER)<--
00480         
00481         string::size_type pos_open  = name.find("(",0);
00482         string::size_type pos_close = name.find(")",0);
00483         
00484         if ( (pos_open  != pos_close) &&
00485              (pos_open  != string::npos) &&
00486              (pos_close != string::npos) ) {
00487           
00488           name.erase(pos_open,pos_close-pos_open+1);
00489         }
00490       }    
00491       //
00492       ast.name = name;
00493       remove_leading_trailing_spaces(ast.name);
00494       
00495       // ast.n = 0; // arbitrary, for the moment
00496       ast.n = atoi(number.c_str());
00497       
00498       ast.mag  = atof(absolute_magnitude.c_str());
00499          
00500       a                = atof(semimajor_axis.c_str());
00501       e                = atof(eccentricity.c_str());
00502       i                = (pi/180)*atof(inclination.c_str());
00503       omega_node       = (pi/180)*atof(node.c_str());
00504       omega_pericenter = (pi/180)*atof(pericenter.c_str());
00505       M                = (pi/180)*atof(mean_anomaly.c_str());
00506       
00507       ast.orb.a                = FromUnits(a,AU);
00508       ast.orb.e                = e;
00509       ast.orb.i                = i;
00510       ast.orb.omega_node       = omega_node;
00511       ast.orb.omega_pericenter = omega_pericenter;
00512       ast.orb.M                = M;
00513 
00514       int ch;
00515       // cerr << "epoch string: " << epoch << endl;
00516       //// epoch
00517       // year
00518       year.assign(epoch,0,1);
00519       ch = (int)year.c_str()[0];
00520       // cerr << "ch: " << ch << "  " << ((int)('I')) << endl;
00521       switch (ch) {
00522       case 'I': y=1800; break; 
00523       case 'J': y=1900; break;
00524       case 'K': y=2000; break;
00525       }
00526       //
00527       year.assign(epoch,1,2);
00528       y += atoi(year.c_str());
00529       // month
00530       month.assign(epoch,3,1);
00531       ch = (int)month.c_str()[0];
00532       switch (ch) {
00533       case '1': m=1;  break;
00534       case '2': m=2;  break;
00535       case '3': m=3;  break;
00536       case '4': m=4;  break;
00537       case '5': m=5;  break;
00538       case '6': m=6;  break;
00539       case '7': m=7;  break;
00540       case '8': m=8;  break;
00541       case '9': m=9;  break;
00542       case 'A': m=10; break; 
00543       case 'B': m=11; break;
00544       case 'C': m=12; break;
00545       }
00546       // day
00547       day.assign(epoch,4,1);
00548       ch = (int)day.c_str()[0];
00549       switch (ch) {
00550       case '1': d=1;  break;
00551       case '2': d=2;  break;
00552       case '3': d=3;  break;
00553       case '4': d=4;  break;
00554       case '5': d=5;  break;
00555       case '6': d=6;  break;
00556       case '7': d=7;  break;
00557       case '8': d=8;  break;
00558       case '9': d=9;  break;
00559       case 'A': d=10; break; 
00560       case 'B': d=11; break;
00561       case 'C': d=12; break;
00562       case 'D': d=13; break;
00563       case 'E': d=14; break;
00564       case 'F': d=15; break;
00565       case 'G': d=16; break;
00566       case 'H': d=17; break;
00567       case 'I': d=18; break;
00568       case 'J': d=19; break;
00569       case 'K': d=20; break;
00570       case 'L': d=21; break;
00571       case 'M': d=22; break;
00572       case 'N': d=23; break;
00573       case 'O': d=24; break;
00574       case 'P': d=25; break;
00575       case 'Q': d=26; break;
00576       case 'R': d=27; break;
00577       case 'S': d=28; break;
00578       case 'T': d=29; break;
00579       case 'U': d=30; break;
00580       case 'V': d=31; break;
00581       }
00582       //
00583       // cerr << "MPC::Read() --> y: " << y << "  m: " << m << "  d: " << d << endl;
00584       // ast.orb.time = FromUnits(GregorianToSdn(y,m,d)-0.5,DAY);
00585       tmp_date.SetGregor(y,m,d,TDT);
00586       ast.orb.epoch.SetDate(tmp_date);
00587       
00588       // ast.orb.T = sqrt(4*pisq/(GetG()*GetMSun())*pow(FromUnits(ast.orb.a,AU),3));
00589       ast.orb.mu = GetG()*GetMSun();
00590       // ast.orb.ref_body.mass=GetMSun();
00591       // ast.orb.ref_body = Body("",GetMSun());
00592       // ast.orb.ref_body = orb_ref_body;
00593       
00594       /* 
00595          switch (universe->GetReferenceSystem()) {
00596          case ECLIPTIC: break;
00597          case EQUATORIAL:
00598          { 
00599          // cerr << "Rotating astorb orbit..." << endl;
00600          const double obleq_rad = obleq(tmp_date).GetRad();
00601          Vector position,velocity;
00602          ast.orb.RelativePosVel(position,velocity);
00603          position.rotate(0.0,obleq_rad,0.0);
00604          velocity.rotate(0.0,obleq_rad,0.0);
00605          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
00606          }
00607          break;
00608          }
00609       */
00610       
00611       switch (universe->GetReferenceSystem()) {
00612       case ECLIPTIC: break;
00613       case EQUATORIAL:
00614         { 
00615           Vector position,velocity;
00616           ast.orb.RelativePosVel(position,velocity);
00617           EclipticToEquatorial_J2000(position);
00618           EclipticToEquatorial_J2000(velocity);
00619           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
00620         }
00621         break;
00622       }
00623       
00624       // cerr << "adding object: " << ast.name << "  a: " << ast.orb.a << endl;
00625       
00626       db->push_back(ast);
00627       
00628     }
00629     
00630     read_finished();
00631   }
00632   
00633   
00634   // MPC Comets orbit database
00635   MPCCometFile::MPCCometFile() : AsteroidDatabaseFile() {
00636     db = new AsteroidDatabase();
00637     // status = CLOSE;
00638   }
00639   
00640   MPCCometFile::~MPCCometFile() {
00641     delete db;
00642     db = 0;
00643   }
00644   
00645   /* 
00646      bool MPCCometFile::GuessCorrectType() {
00647      return true;
00648      }
00649   */
00650   
00651   void MPCCometFile::Read() {
00652     
00653     // if (status == CLOSE) Open();
00654     
00655     Open();
00656     
00657     if (status != OPEN_R) {
00658       ORSA_ERROR("Status error!");
00659       return;
00660     }
00661     
00662     db->clear();
00663     
00664     char line[300];
00665     
00666     double a,e,i,omega_node,omega_pericenter,M;
00667     string number,type,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
00668     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
00669     // string ceu;
00670     string pericenter_distance,pericenter_epoch;
00671     
00672     string year,month,day;
00673     int    y=0,m=0,d=0;
00674     double frac_day;
00675     
00676     Asteroid ast;
00677     
00678     double q;
00679     
00680     REWIND_FILE(file);
00681     
00682     // Date tmp_date(TDT);
00683     Date tmp_date;
00684     
00685     bool have_perturbed_solution_epoch;
00686     
00687     unsigned int local_index = 0;
00688     bool bool_stop=false;
00689     bool bool_pause=false;
00690     
00691     while ( (GETS_FILE(line,300,file)) != 0 ) {
00692       
00693       if (strlen(line) < 90) continue; // not a good line, maybe a comment or a white line...
00694       
00695       ++local_index;
00696       read_progress(local_index,bool_pause,bool_stop);
00697       
00698       if (bool_stop) break;
00699       
00700       while (bool_pause) {
00701         // cerr << "AstorbFile::Read() sleeping..." << endl;
00702         sleep(1);
00703         read_progress(local_index,bool_pause,bool_stop);
00704       }
00705       
00706       // uncomment the ones used
00707       number.assign(line,0,4);
00708       type.assign(line,4,1);
00709       // name.assign(line,5,7);
00710       name.assign(line,102,strlen(line)-102-1); // the last -1 is set to avoid the '\n' character in the name
00711       // cerr << "comet name: " << name << endl;
00712       pericenter_epoch.assign(line,14,15);
00713       pericenter_distance.assign(line,30,9);
00714       eccentricity.assign(line,41,8);
00715       pericenter.assign(line,51,8);
00716       node.assign(line,61,8);
00717       inclination.assign(line,71,8);
00718       //
00719       epoch.assign(line,81,8);
00720       
00721       // conversions
00722       
00723       ast.name = name;
00724       remove_leading_trailing_spaces(ast.name);
00725       
00726       ast.n = 0; // arbitrary, for the moment
00727       
00728       ast.mag  = atof(absolute_magnitude.c_str());
00729       
00730       // a                = atof(semimajor_axis.c_str());
00731       e                = atof(eccentricity.c_str());
00732       
00733       // to be tested...
00734       q = atof(pericenter_distance.c_str());
00735       if (e == 1.0) {
00736         a = q;
00737       } else {
00738         a = q/fabs(1.0-e);
00739       }
00740       
00741       i                = (pi/180)*atof(inclination.c_str());
00742       omega_node       = (pi/180)*atof(node.c_str());
00743       omega_pericenter = (pi/180)*atof(pericenter.c_str());
00744       // M                = (pi/180)*atof(mean_anomaly.c_str());
00745       
00746       year.assign(epoch,0,4);
00747       month.assign(epoch,4,2);
00748       day.assign(epoch,6,2);
00749       
00750       y = atoi(year.c_str());
00751       m = atoi(month.c_str());
00752       d = atoi(day.c_str());
00753       
00754       // check on the presence of the 'perturbed solutions' epoch...
00755       if (y < 1700) {
00756         have_perturbed_solution_epoch=false;
00757       } else {
00758         have_perturbed_solution_epoch=true;
00759       } 
00760       
00761       tmp_date.SetGregor(y,m,d,TDT);
00762       
00763       year.assign(pericenter_epoch,0,4);
00764       month.assign(pericenter_epoch,5,2);
00765       day.assign(pericenter_epoch,8,7);
00766       
00767       y = atoi(year.c_str());
00768       m = atoi(month.c_str());
00769       frac_day = atof(day.c_str());
00770       
00771       Date peri_date;
00772       peri_date.SetGregor(y,m,frac_day,TDT);
00773       UniverseTypeAwareTime pericenter_passage(peri_date);
00774       
00775       if (have_perturbed_solution_epoch) {
00776         ast.orb.epoch.SetDate(tmp_date);
00777       } else {
00778         ast.orb.epoch.SetDate(peri_date);
00779       }
00780       
00781       ast.orb.mu = GetG()*GetMSun();
00782       
00783       ast.orb.a                = FromUnits(a,AU);
00784       ast.orb.e                = e;
00785       ast.orb.i                = i;
00786       ast.orb.omega_node       = omega_node;
00787       ast.orb.omega_pericenter = omega_pericenter;
00788       //
00789       if (have_perturbed_solution_epoch) {
00790         M = ((ast.orb.epoch.Time() - pericenter_passage.Time())/ast.orb.Period())*twopi;  
00791         M = fmod(10*twopi+fmod(M,twopi),twopi);
00792         //
00793         ast.orb.M = M;
00794       } else {
00795         ast.orb.M = 0.0;
00796       }
00797       
00798       // cerr << "comet: " << ast.name << "  q: " << q << "  e: " << e << "  i: " << i*(180/pi) << endl;
00799       
00800       /* 
00801          switch (universe->GetReferenceSystem()) {
00802          case ECLIPTIC: break;
00803          case EQUATORIAL:
00804          { 
00805          // cerr << "Rotating astorb orbit..." << endl;
00806          const double obleq_rad = obleq(tmp_date).GetRad();
00807          Vector position,velocity;
00808          ast.orb.RelativePosVel(position,velocity);
00809          position.rotate(0.0,obleq_rad,0.0);
00810          velocity.rotate(0.0,obleq_rad,0.0);
00811          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
00812          }
00813          break;
00814          }
00815       */
00816       
00817       switch (universe->GetReferenceSystem()) {
00818       case ECLIPTIC: break;
00819       case EQUATORIAL:
00820         { 
00821           Vector position,velocity;
00822           ast.orb.RelativePosVel(position,velocity);
00823           EclipticToEquatorial_J2000(position);
00824           EclipticToEquatorial_J2000(velocity);
00825           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
00826         }
00827         break;
00828       }
00829       
00830       db->push_back(ast);
00831       
00832     }
00833     
00834     read_finished();
00835   }
00836   
00837   //! MPC observation file
00838   
00839   MPCObsFile::MPCObsFile() { }
00840   
00841   void MPCObsFile::Read() {
00842     
00843     // if (file == 0) Open();
00844     // if (status != OPEN_R) Open();
00845     // if (status == CLOSE) Open();
00846     
00847     Open();
00848     
00849     if (status != OPEN_R) {
00850       ORSA_ERROR("Status error!");
00851       return;
00852     }
00853     
00854     obs.clear();
00855     
00856     REWIND_FILE(file);
00857     
00858     // cerr << "...inside read_MPC()...\n";
00859     
00860     Observation dummy_obs;
00861     // dummy_obs.date.SetTimeScale(UTC); // checked, is UTC
00862     
00863     // double y,m,d;
00864     double gradi,primi,secondi;
00865     
00866     // double tmp;
00867     
00868     char line[256];
00869     
00870     // int entries = 0;
00871     
00872     string number,designation,discovery_asterisk,note1,note2;
00873     string date,ra,dec;
00874     string magnitude;
00875     string observatory_code;
00876     
00877     while (GETS_FILE(line,256,file) != 0) {
00878       
00879       if (strlen(line) < 80) continue; // not a good line, maybe a comment or a white line...
00880       // if (strlen(line) < 79) continue; // not a good line, maybe a comment or a white line...
00881       // if (strlen(line) < 78) continue; // not a good line, maybe a comment or a white line...
00882       
00883       // cerr << "inside while loop..." << endl;
00884       
00885       // entries--;
00886       
00887       number.assign(line,0,5);
00888       // cerr << "...pass (1)" << endl;
00889       // cerr <<"number: " << atoi(number.c_str()) << endl;
00890       designation.assign(line,5,7); 
00891       remove_leading_trailing_spaces(designation);
00892       // cerr <<"designation: " << designation << endl;
00893       discovery_asterisk.assign(line,12,1); 
00894       // cerr << "asterisk: [" << discovery_asterisk << "]\n";
00895       // cerr << "...pass (1/B)" << endl;
00896       note1.assign(line,13,1);
00897       // cerr << "...pass (1/C)" << endl;
00898       note2.assign(line,14,1);
00899       // cerr << "...pass (1/D)" << endl;
00900       date.assign(line,15,17);
00901       // cerr << "...pass (2)" << endl;
00902       // cerr << "date: " << date << endl;
00903       ra.assign(line,32,12);
00904       // cerr << "ra: " << ra << endl;
00905       dec.assign(line,44,12);
00906       // cerr << "dec: " << dec << endl;
00907       magnitude.assign(line,65,6);
00908       // cerr << "magnitude: " << magnitude << endl;
00909       observatory_code.assign(line,77,3);
00910       remove_leading_trailing_spaces(observatory_code);
00911       // cerr << "...pass (3)" << endl;
00912       
00913       dummy_obs.designation        = designation;
00914       dummy_obs.discovery_asterisk = discovery_asterisk;
00915       dummy_obs.obscode            = observatory_code;
00916       
00917       double _tmp = 0.0;
00918       sscanf(magnitude.c_str(),"%lf",&_tmp);
00919       dummy_obs.mag = _tmp;
00920       
00921       // cerr << "DES\n";
00922       
00923       // printf("DATE: %s\n",date.c_str());
00924       double y=0.0, m=0.0, d=0.0;
00925       sscanf(date.c_str(),"%lf %lf %lf",&y,&m,&d);
00926       // printf("LETTI: %f %f %f\n",y,m,d);
00927       //
00928       
00929       // dummy_obs.julian_date = GregorianToSdn((int)y,(int)m,(int)d) + ((d - floor(d)) - 0.5);
00930       // cerr << "JD\n";
00931       dummy_obs.date.SetGregor((int)y,(int)m,d,UTC);
00932       
00933       // cerr << "...pass (5)" << endl;
00934       
00935       // cout << "designation: " << designation << " " << ra << endl;      
00936       sscanf(ra.c_str(),"%lf %lf %lf",&gradi,&primi,&secondi);
00937       // angle_rad_hms(&tmp,gradi,primi,secondi);  //   dummy_obs.alpha
00938       // dummy_obs.alpha = tmp;
00939       dummy_obs.ra.SetHMS(gradi,primi,secondi);
00940       
00941       
00942       // cout <<"designation: " << designation << " " << dec << endl;      
00943       sscanf(dec.c_str(),"%lf %lf %lf",&gradi,&primi,&secondi);
00944       // angle_rad(&tmp,gradi,primi,secondi); //   dummy_obs.delta       
00945       // dummy_obs.delta = tmp;
00946       dummy_obs.dec.SetDPS(gradi,primi,secondi);
00947       
00948       // cerr << "OBS: " << dummy_obs.date.GetJulian() << "  " << dummy_obs.ra.GetRad() << "  " << dummy_obs.dec.GetRad() << endl;
00949       
00950       //check for a good observation
00951       if ((designation != "") && 
00952           (observatory_code != "") &&
00953           (strlen(observatory_code.c_str())) == 3) {
00954         if ( (isalnum(observatory_code[0])) &&
00955              (isalnum(observatory_code[1])) &&
00956              (isalnum(observatory_code[2])) &&
00957              (isspace(line[19])) &&
00958              (isspace(line[22])) &&
00959              (isspace(line[31])) &&
00960              (isspace(line[34])) &&
00961              (isspace(line[37])) &&
00962              (isspace(line[43])) &&
00963              (isspace(line[47])) &&
00964              (isspace(line[50])) &&
00965              (isspace(line[55])) ) {
00966           obs.push_back(dummy_obs);
00967         }
00968       }
00969       // obs->push_back(dummy_obs);
00970       
00971       // cerr << "PB\n";
00972       
00973       // D_ar.Insert(dummy_obs.alpha);
00974       // D_dec.Insert(dummy_obs.delta);
00975       
00976     }
00977     
00978   }
00979   
00980   bool MPCObsFile::ReadNominalOrbit(OrbitWithEpoch &) {
00981     // true if succeed
00982     return false;
00983   }
00984   
00985   //! AstDyS observation file, usually with a .rwo extension
00986   
00987   RWOFile::RWOFile() { }
00988   
00989   void RWOFile::Read() {
00990     
00991     Open();
00992     
00993     obs.clear();
00994     
00995     REWIND_FILE(file);
00996     
00997     Observation dummy_obs;
00998     // dummy_obs.date.SetTimeScale(UTC); // checked, is UTC
00999     
01000     double y,m,d,p,s,h;
01001     
01002     char line[256];
01003     
01004     string designation;
01005     string date,ra,dec;
01006     string observatory_code;
01007     
01008     while (GETS_FILE(line,256,file) != 0) {
01009       
01010       // if (strlen(line) < 45) continue; // not a good line, maybe a comment or a white line...
01011       if (strlen(line) < 130) continue; // not a good line, maybe a comment or a white line... NOTE: radar observations are not read now, and have shorter lines
01012       if (line[0] == '!')     continue; // comment/header line
01013       
01014       designation.assign(line,1,9); 
01015       remove_leading_trailing_spaces(designation);
01016       
01017       date.assign(line,15,17);
01018       
01019       ra.assign(line,34,12);
01020       
01021       dec.assign(line,63,12);
01022       
01023       observatory_code.assign(line,114,3);
01024       remove_leading_trailing_spaces(observatory_code);
01025       
01026       dummy_obs.designation        = designation;
01027       // dummy_obs.discovery_asterisk = discovery_asterisk;
01028       // dummy_obs.obscode            = atoi(observatory_code.c_str());
01029       dummy_obs.obscode            = observatory_code;
01030       
01031       /* sscanf(magnitude.c_str(),"%lf",&tmp);
01032          dummy_obs.mag = tmp;
01033       */
01034       
01035      
01036       sscanf(date.c_str(),"%lf %lf %lf",&y,&m,&d);
01037       dummy_obs.date.SetGregor((int)y,(int)m,d,UTC);
01038       
01039       sscanf(ra.c_str(),"%lf %lf %lf",&h,&m,&s);
01040       dummy_obs.ra.SetHMS(h,m,s);
01041       
01042       sscanf(dec.c_str(),"%lf %lf %lf",&d,&p,&s);
01043       dummy_obs.dec.SetDPS(d,p,s);
01044       
01045       //check for a good observation
01046       if ((designation != "") && 
01047           (observatory_code != "") ) {
01048         obs.push_back(dummy_obs);
01049       }
01050       
01051     }
01052     
01053   }
01054   
01055   // Locations of the observatories
01056   
01057   LocationFile::LocationFile() : ReadFile() { }
01058   
01059   void LocationFile::Read() {
01060     
01061     Open();
01062     
01063     if (status != OPEN_R) {
01064       ORSA_ERROR("Status error!");
01065       return;
01066     }
01067     
01068     string obscode,longitude,pos_xy,pos_z,name;
01069     
01070     char line[300];
01071     
01072     Location dummy_loc;
01073     
01074     // dummy_loc.present = true;
01075     
01076     REWIND_FILE(file);
01077     
01078     string stag, spre="<pre>", spreend="</pre>";
01079     char tag[256];
01080     
01081     // reach <pre> and skip header
01082     while (GETS_FILE(line,300,file) != 0) {
01083       sscanf(line,"%s",tag);
01084       stag = tag;
01085       if (stag == spre) {
01086         // skip another line (header) and break;
01087         GETS_FILE(line,300,file);
01088         break;
01089       } 
01090     }
01091     
01092     while (GETS_FILE(line,300,file) != 0) {
01093       
01094       sscanf(line,"%s",tag);
01095       stag = tag;
01096       // cerr << "stag: " << stag << endl;
01097       if (stag == spreend) break;
01098       
01099       if (strlen(line) < 30) continue;
01100       
01101       obscode.assign(line,0,3);
01102       longitude.assign(line,4,8);
01103       pos_xy.assign(line,13,7);
01104       pos_z.assign(line,21,8);
01105       name.assign(line,30,strlen(line)-30-1); // the last -1 is set to avoid the '\n' character in the name
01106       
01107       remove_leading_trailing_spaces(longitude);
01108       remove_leading_trailing_spaces(pos_xy);
01109       remove_leading_trailing_spaces(pos_z);
01110       remove_leading_trailing_spaces(name);
01111       
01112       // set zero values, for locations like SOHO, HST...
01113       dummy_loc.lon = dummy_loc.pxy = dummy_loc.pz = 0.0;
01114       
01115       if (longitude.size()) dummy_loc.lon  = (pi/180.0) * atof(longitude.c_str());
01116       if (pos_xy.size())    dummy_loc.pxy  = FromUnits(atof(pos_xy.c_str()),REARTH);
01117       if (pos_z.size())     dummy_loc.pz   = FromUnits(atof(pos_z.c_str()), REARTH);
01118       //
01119       dummy_loc.name = name;
01120       
01121       locations[obscode] = dummy_loc;
01122       
01123       codes.push_back(obscode);
01124     }
01125     
01126   }
01127   
01128   Vector LocationFile::ObsPos(const string obscode, const Date &date) {
01129     
01130     std::list<std::string>::iterator it_find = find(codes.begin(),codes.end(),obscode);
01131     if (it_find == codes.end()) {
01132       ORSA_ERROR("obscode %s not found in file %s",obscode.c_str(),GetFileName().c_str());      
01133       return Vector();
01134     }
01135     
01136     const double lon = locations[obscode].lon;
01137     const double pxy = locations[obscode].pxy;
01138     const double pz  = locations[obscode].pz;
01139     
01140 #ifdef HAVE_SINCOS
01141     double s,c;
01142     ::sincos(lon,&s,&c);
01143     Vector obspos(pxy*c,pxy*s,pz);
01144 #else // HAVE_SINCOS    
01145     Vector obspos(pxy*cos(lon),pxy*sin(lon),pz);
01146 #endif // HAVE_SINCOS
01147     
01148     obspos.rotate(gmst(date).GetRad(),0,0);
01149     
01150     if (universe->GetReferenceSystem() == ECLIPTIC) {
01151       EquatorialToEcliptic_J2000(obspos);
01152     }
01153     
01154     return obspos;
01155   }
01156   
01157   // SWIFT files
01158   SWIFTFile::SWIFTFile(OrbitStream &osin) : ReadFile() {
01159     os = &osin;
01160     // status = CLOSE;
01161   }
01162   
01163   // SWIFT common data
01164   int    nast;
01165   double el[6],l_ts,w_ts,file_time;
01166   
01167   int SWIFTRawReadBinaryFile(FILE_TYPE file, const int version = 2) {
01168     
01169     // version == 1 --> not filetered data
01170     // version == 2 --> filtered data
01171     
01172     int good = 0;
01173     double dummy;
01174     
01175     if ( version == 1 ) {  
01176       READ_FILE(&dummy,        sizeof(int),    1,file);
01177       READ_FILE(&nast,         sizeof(int),    1,file);
01178       READ_FILE(&el,           sizeof(double), 6,file);
01179       READ_FILE(&file_time,    sizeof(double), 1,file);
01180       good = READ_FILE(&dummy, sizeof(int),    1,file);
01181     } else if ( version == 2 ) {
01182       READ_FILE(&dummy,        sizeof(int),    1,file);
01183       READ_FILE(&nast,         sizeof(int),    1,file);
01184       READ_FILE(&el,           sizeof(double), 6,file);
01185       READ_FILE(&l_ts,         sizeof(double), 1,file);
01186       READ_FILE(&w_ts,         sizeof(double), 1,file);
01187       READ_FILE(&file_time,    sizeof(double), 1,file);
01188       good = READ_FILE(&dummy, sizeof(int),    1,file);
01189     }
01190     
01191     // Units conversions
01192     file_time = FromUnits(file_time,YEAR);
01193     
01194     return (good);
01195   }
01196   
01197   int SWIFTFile::AsteroidsInFile() {
01198     
01199     // close and reopen to avoid odd zlib problems related to the gzseek function
01200     Close();
01201     Open();
01202     
01203     int number_of_asteroids_in_file=0;
01204     
01205     REWIND_FILE(file);
01206     
01207     int good;
01208     while ( (good = SWIFTRawReadBinaryFile(file,2)) != 0) {
01209       if (number_of_asteroids_in_file<nast) number_of_asteroids_in_file = nast;
01210       else if (number_of_asteroids_in_file!=0) break;
01211     }
01212     
01213     return (number_of_asteroids_in_file);
01214   }
01215   
01216   
01217   void SWIFTFile::Read() {
01218     
01219     // close and reopen to avoid odd zlib problems related to the gzseek function
01220     Close();
01221     Open();
01222     
01223     if (status != OPEN_R) {
01224       ORSA_ERROR("Status error!");
01225       return;
01226     }
01227     
01228     OrbitStream &ost = *os;
01229     
01230     const int version = 2;
01231     
01232     // cerr << "reading data from the input file...\n";
01233     
01234     OrbitWithEpoch fo;
01235     
01236     // double wwj=0,wj=0,cj=0;
01237     // double wtil,crit,wmix;
01238     double time_old = 0, timestep;
01239     
01240     int jump = 0, i_jump = 0;
01241     
01242     // reset!
01243     // ost.resize(0);
01244     ost.clear();
01245     ost.timestep = 0.0;
01246     const int asteroid_number =  ost.asteroid_number;
01247     // cerr << " SWIFTFile::Read() --> reading object: " << asteroid_number << endl;
01248     char label[10];
01249     sprintf(label,"%04i",ost.asteroid_number);
01250     ost.label = label;
01251     // cerr << "LABEL: [" << ost.label << "]" << endl;
01252     REWIND_FILE(file);
01253     
01254     if ( version == 1 )
01255       jump = 3*sizeof(int)+7*sizeof(double);
01256     else if ( version == 2 )
01257       jump = 3*sizeof(int)+9*sizeof(double);
01258     
01259     int  good = 1, check = 0, number_of_asteroids_in_file = 0;
01260     
01261     do {
01262       
01263       if (check == 0) {
01264         good = SWIFTRawReadBinaryFile(file,version);    
01265       } else {
01266         
01267         i_jump = (number_of_asteroids_in_file+asteroid_number-nast-1)%(number_of_asteroids_in_file);
01268         
01269         if (i_jump != 0) {
01270           if ( (SEEK_FILE(file,jump*i_jump,SEEK_CUR)) == -1) {
01271             cerr << "setting good=0 from SEEK_FILE..." << endl;
01272             good = 0; 
01273           }
01274         }
01275         
01276         if (good != 0) {
01277           good = SWIFTRawReadBinaryFile(file,version);
01278         }
01279         
01280       }
01281       
01282       if ( number_of_asteroids_in_file < nast ) {
01283         number_of_asteroids_in_file = nast;
01284       } else {
01285         check = 1;
01286       }
01287       
01288       //// asteroid number too big!
01289       if ( (check == 1) && (asteroid_number > number_of_asteroids_in_file) ) {
01290         ORSA_ERROR("asteroid number too big (%d > %d)", asteroid_number, number_of_asteroids_in_file);
01291         return;
01292       }
01293       
01294       if (nast == asteroid_number && good != 0) {
01295         
01296         if ((file_time >= time_old) && 
01297             (file_time >= ost.wp.window_start)) {
01298           
01299           fo.epoch.SetTime(file_time);
01300           fo.a                = el[4];
01301           fo.e                = el[3];
01302           fo.i                = (pi/180.0)*el[2];
01303           fo.omega_node       = (pi/180.0)*el[0];
01304           fo.omega_pericenter = (pi/180.0)*el[1];
01305           fo.M                = (pi/180.0)*el[5];
01306           //
01307           fo.libration_angle  = (pi/180.0)*l_ts; // temporary
01308           
01309           ost.push_back(fo);
01310           
01311           // QUICK AND DIRTY!
01312           if (fo.e >= 1.0) {
01313             cerr << "reading eccentricity > 1.0, returning." << endl;
01314             return;
01315           }
01316           
01317           
01318           if ( ((file_time) > (ost.wp.window_amplitude+ost.wp.window_start)) && (ost.wp.window_step == 0.0) ) {
01319             return;
01320           }
01321         }
01322         
01323         timestep = file_time - time_old;
01324         time_old = file_time;
01325         // one of all, but not the first...
01326         if (ost.size() == 2) {
01327           ost.timestep = timestep;
01328         }
01329         
01330       }
01331       
01332     } while (good != 0);
01333   }
01334   
01335   //! orsa configuration file
01336   // OrsaConfigFile::OrsaConfigFile(Config *conf_in) : ReadWriteFile() {
01337   OrsaConfigFile::OrsaConfigFile() : ReadWriteFile() {
01338     
01339     // status = CLOSE;
01340     
01341     // conf = conf_in;
01342     
01343     char path[1024], command[1024];
01344     
01345     // needed to avoid some odd segfaults...
01346     // OrsaPaths p; // make sure the constructor gets called
01347     
01348     // cerr << "OrsaPaths::work_path() = " << OrsaPaths::work_path() << endl;
01349     
01350     strcpy(path, OrsaPaths::work_path());
01351 #ifndef _WIN32    
01352     sprintf(command,"mkdir -p %s",path);
01353     system(command);
01354 #else
01355     _mkdir(path);
01356 #endif    
01357     strcat(path,"config");
01358     
01359     SetFileName(path);
01360     
01361     list_enum.push_back(JPL_EPHEM_FILE);
01362     list_enum.push_back(JPL_DASTCOM_NUM);
01363     list_enum.push_back(JPL_DASTCOM_UNNUM);
01364     list_enum.push_back(JPL_DASTCOM_COMET);
01365     list_enum.push_back(LOWELL_ASTORB);
01366     list_enum.push_back(MPC_MPCORB);
01367     list_enum.push_back(MPC_COMET);
01368     list_enum.push_back(MPC_NEA);
01369     list_enum.push_back(MPC_DAILY);
01370     list_enum.push_back(MPC_DISTANT);
01371     list_enum.push_back(MPC_PHA);
01372     list_enum.push_back(MPC_UNUSUALS);
01373     list_enum.push_back(ASTDYS_ALLNUM_CAT);
01374     list_enum.push_back(ASTDYS_ALLNUM_CTC);
01375     list_enum.push_back(ASTDYS_ALLNUM_CTM);
01376     list_enum.push_back(ASTDYS_UFITOBS_CAT);
01377     list_enum.push_back(ASTDYS_UFITOBS_CTC);
01378     list_enum.push_back(ASTDYS_UFITOBS_CTM);
01379     list_enum.push_back(NEODYS_CAT);
01380     list_enum.push_back(NEODYS_CTC);
01381     list_enum.push_back(OBSCODE);
01382     // TLE
01383     list_enum.push_back(TLE_NASA);
01384     list_enum.push_back(TLE_GEO);
01385     list_enum.push_back(TLE_GPS);
01386     list_enum.push_back(TLE_ISS);
01387     list_enum.push_back(TLE_KEPELE);
01388     list_enum.push_back(TLE_VISUAL);
01389     list_enum.push_back(TLE_WEATHER);
01390     // textures
01391     list_enum.push_back(TEXTURE_SUN);
01392     list_enum.push_back(TEXTURE_MERCURY);
01393     list_enum.push_back(TEXTURE_VENUS);
01394     list_enum.push_back(TEXTURE_EARTH);
01395     list_enum.push_back(TEXTURE_MOON);
01396     list_enum.push_back(TEXTURE_MARS);
01397     list_enum.push_back(TEXTURE_JUPITER);
01398     list_enum.push_back(TEXTURE_SATURN);
01399     list_enum.push_back(TEXTURE_URANUS);
01400     list_enum.push_back(TEXTURE_NEPTUNE);
01401   }
01402   
01403   void OrsaConfigFile::Read() {
01404     
01405     // if (file == 0) Open();
01406     // if (status == CLOSE) Open();
01407     
01408     Open(OPEN_R);
01409     
01410     // should improve this check
01411     // if (status != OPEN_R) return;
01412     
01413     if (status != OPEN_R) {
01414       ORSA_ERROR("Status error!");
01415       return;
01416     }
01417     
01418     char line[1024];    
01419     string stag, svalue;
01420     
01421     REWIND_FILE(file);
01422     
01423     while (GETS_FILE(line,1024,file) != 0) {
01424       
01425       {
01426         // the first white space is the separator between tag and value
01427         string s_line=line;
01428         string::size_type white_space_pos;
01429         white_space_pos = s_line.find(" ",0);
01430         if (white_space_pos != string::npos) {
01431           stag.assign(s_line,0,white_space_pos);
01432           svalue.assign(s_line,white_space_pos+1,s_line.size()-white_space_pos-2);
01433           remove_leading_trailing_spaces(stag);
01434           remove_leading_trailing_spaces(svalue);
01435           // 
01436           // cerr << "tag -->" << stag << "<--     value -->" << svalue << "<-- " << endl;
01437         }
01438       }
01439       
01440       if (svalue.size()>0) {
01441         
01442         list<ConfigEnum>::const_iterator it = list_enum.begin();
01443         while (it != list_enum.end()) {
01444           if (stag == config->paths[(*it)]->tag) {
01445             config->paths[(*it)]->SetValue(svalue);
01446             break;
01447           }
01448           ++it;
01449         }
01450         
01451       }
01452       
01453     }
01454     
01455     Close();
01456   }
01457   
01458   void OrsaConfigFile::Write() {
01459     
01460     // this close is necessary to avoid multiple write of the same options
01461     Close();
01462     
01463     // *** TODO: make a backup copy before to save the new one! *** 
01464     
01465     Open(OPEN_W);
01466     
01467     if (status != OPEN_W) {
01468       ORSA_ERROR("Status error!");
01469       return;
01470     }
01471     
01472     // cerr << "OrsaConfigFile::Write() ==> " << filename << endl;
01473     
01474     // rewind(file);
01475     
01476     char line[1024];
01477     
01478     list<ConfigEnum>::const_iterator it = list_enum.begin();
01479     while (it != list_enum.end()) {
01480       sprintf(line,"%s %s\n",config->paths[(*it)]->tag.c_str(),config->paths[(*it)]->GetValue().c_str());   
01481       PUTS_FILE(line,file);
01482       ++it;
01483     }
01484     
01485     FLUSH_FILE(file);
01486     
01487     Close(); 
01488   }
01489   
01490   
01491   //! OrsaFile
01492   
01493 #define SWAP_MACRO( A, B, TEMP)   { (TEMP) = (A);  (A) = (B);  (B) = (TEMP); }
01494   
01495   static void swap_4(void *ptr) {
01496     char *tptr = (char *)ptr, tchar;
01497     
01498     SWAP_MACRO( tptr[0], tptr[3], tchar);
01499     SWAP_MACRO( tptr[1], tptr[2], tchar);
01500   }
01501   
01502   static void swap_8(void *ptr) {
01503     char *tptr = (char *)ptr, tchar;
01504     
01505     SWAP_MACRO( tptr[0], tptr[7], tchar);
01506     SWAP_MACRO( tptr[1], tptr[6], tchar);
01507     SWAP_MACRO( tptr[2], tptr[5], tchar);
01508     SWAP_MACRO( tptr[3], tptr[4], tchar);
01509   }
01510   
01511   static void swap(void *ptr, unsigned int size) {
01512     switch(size) {
01513     case 4: swap_4(ptr); break;
01514     case 8: swap_8(ptr); break;
01515     default:
01516       ORSA_WARNING("called read_swap with size = %i",size);
01517       break;
01518     }
01519   }
01520   
01521   size_t OrsaFile::read_swap(void *ptr, unsigned int size) {
01522     const size_t s = READ_FILE(ptr,size,1,file);
01523     if (swap_bytes) {
01524       swap(ptr,size);
01525     }
01526     return s;
01527   }
01528   
01529   OrsaFile::OrsaFile() : ReadWriteFile() { }
01530   
01531   void OrsaFile::Write() {
01532     
01533     Open(OPEN_W);
01534     
01535     if (status != OPEN_W) {
01536       ORSA_ERROR("Status error!");
01537       return;
01538     }
01539     
01540     if (!universe) {
01541       ORSA_ERROR("cannot write a non-allocated universe!");
01542       return;
01543     }
01544     
01545     Write(&universe);
01546     
01547     FLUSH_FILE(file);
01548     
01549     Close();
01550   }
01551   
01552   void OrsaFile::Read() {
01553     
01554     Open(OPEN_R);
01555     
01556     if (status != OPEN_R) {
01557       ORSA_ERROR("Status error!");
01558       return;
01559     }
01560     
01561     Read(&universe);
01562     
01563     Close();
01564     
01565     ORSA_DEBUG("ORSA file %s [ORSA version: %s, byte order: %i, evolutions: %i, units: [%s,%s,%s]]",
01566                GetFileName().c_str(), orsa_version.c_str(), byte_order,universe->size(),
01567                LengthLabel().c_str(), MassLabel().c_str(),TimeLabel().c_str());
01568   }
01569   
01570   void OrsaFile::Write(Universe **u) {
01571     
01572     // endian issues
01573     byte_order = BYTEORDER; // from config.h
01574     Write(&byte_order);
01575     
01576     // various info...
01577     orsa_version = ORSA_VERSION;
01578     Write(&orsa_version);
01579     
01580     time_unit   tu = units->GetTimeBaseUnit();
01581     length_unit lu = units->GetLengthBaseUnit();
01582     mass_unit   mu = units->GetMassBaseUnit();
01583     //
01584     Write(&tu); Write(&lu); Write(&mu);
01585     
01586     UniverseType ut = (*u)->GetUniverseType();
01587     Write(&ut);
01588     
01589     ReferenceSystem rs = (*u)->GetReferenceSystem();
01590     Write(&rs);
01591     
01592     TimeScale ts = (*u)->GetTimeScale();
01593     Write(&ts);
01594     
01595     Write(&((*u)->name));
01596     Write(&((*u)->description));
01597     
01598     unsigned int j;
01599     for(j=0;j<(*u)->size();j++) {
01600       if ((*(*u))[j]!=0) Write(&(*(*u))[j]);
01601     }
01602   }
01603   
01604   void OrsaFile::make_new_universe(Universe **u, length_unit lu, mass_unit mu, time_unit tu, UniverseType ut, ReferenceSystem rs, TimeScale ts) {    
01605     delete (*u);    
01606     (*u) = new Universe(lu,mu,tu,ut,rs,ts);
01607   }
01608   
01609   void OrsaFile::Read(Universe **u) {
01610     
01611     swap_bytes = false;
01612     Read(&byte_order);
01613     if (byte_order != BYTEORDER) {
01614       swap_bytes = true;
01615       swap(&byte_order,sizeof(byte_order));
01616     }
01617     
01618     Read(&orsa_version);
01619     
01620     time_unit   tu;
01621     length_unit lu;
01622     mass_unit   mu;
01623     //
01624     Read(&tu);
01625     Read(&lu);
01626     Read(&mu);
01627     
01628     UniverseType ut;
01629     Read(&ut);
01630     
01631     ReferenceSystem rs;;
01632     Read(&rs);
01633     
01634     TimeScale ts;
01635     Read(&ts);
01636     
01637     make_new_universe(u,lu,mu,tu,ut,rs,ts);
01638     
01639     Read(&((*u)->name));
01640     Read(&((*u)->description));
01641     
01642     Read(&last_ofdt_read); // the others are read by the Read(evol..)
01643     /* 
01644        { // debug
01645        ORSA_ERROR("Read(Universe *u)  ofdt = %i",last_ofdt_read);
01646        }
01647     */
01648     while (last_ofdt_read == OFDT_EVOLUTION) {
01649       Evolution * e = 0;  
01650       Read(&e);
01651       (*u)->push_back(e);
01652     }
01653   }
01654   
01655   void OrsaFile::Write(Evolution **e) {
01656     
01657     OrsaFileDataType t = OFDT_EVOLUTION; Write(&t);
01658     
01659     Write(&((*e)->name));
01660     UniverseTypeAwareTimeStep sp = (*e)->GetSamplePeriod(); Write(&sp);
01661     // const Integrator * itg = (*e)->GetIntegrator(); Write(&itg);
01662     Write((*e)->GetIntegrator());
01663     // Write(&((*e)->interaction));
01664     Write((*e)->GetInteraction());
01665     
01666     unsigned int n = (*e)->start_bodies.size();
01667     Write(&n);
01668     for(unsigned int j=0;j<n;j++) {
01669       Write(&((*e)->start_bodies[j]));
01670     }
01671     //
01672     if (universe->GetUniverseType() == Real) {
01673       n = (*e)->start_JPL_bodies.size();
01674       Write(&n);
01675       for(unsigned int j=0;j<n;++j) {
01676         Write(&((*e)->start_JPL_bodies[j]));
01677       }
01678     }
01679     
01680     // the first only
01681     if ((*e)->size() > 0) Write(&((*(*e))[0]));
01682     
01683     // from the second on, write only position and velocity
01684     for(unsigned int j=1;j<(*e)->size();++j) {
01685       Write(&((*(*e))[j]),true);
01686     }
01687   }
01688   
01689   void OrsaFile::make_new_evolution(Evolution **e) {
01690     delete (*e);
01691     (*e) = new Evolution;
01692   }
01693   
01694   void OrsaFile::Read(Evolution **e) {
01695     
01696     string name;
01697     Read(&name);
01698     
01699     // double sample_period;
01700     UniverseTypeAwareTimeStep sample_period;
01701     Read(&sample_period);
01702     
01703     Integrator * integrator = 0;
01704     Read(&integrator);
01705     
01706     Interaction * interaction = 0;
01707     Read(&interaction);
01708     
01709     make_new_evolution(e);
01710     
01711     (*e)->clear();
01712     
01713     (*e)->name          = name;
01714     (*e)->SetSamplePeriod(sample_period);
01715     (*e)->SetIntegrator(integrator);
01716     (*e)->SetInteraction(interaction);
01717     
01718     delete integrator;
01719     integrator = 0;
01720     
01721     delete interaction;
01722     interaction = 0;
01723     
01724     unsigned int n;
01725     Read(&n);
01726     (*e)->start_bodies.resize(n);
01727     for(unsigned int j=0;j<n;j++) {
01728       Read(&((*e)->start_bodies[j]));
01729     }
01730     //
01731     if (universe->GetUniverseType() == Real) {
01732       Read(&n);
01733       (*e)->start_JPL_bodies.clear();
01734       
01735       JPL_planets tmp_jp;
01736       for(unsigned int j=0;j<n;j++) {
01737         Read(&tmp_jp);
01738         (*e)->start_JPL_bodies.push_back(tmp_jp);
01739       }
01740     }
01741     
01742     // we REALLY need a Frame to keep all the constant values....
01743     Frame f;
01744     
01745     Read(&last_ofdt_read);
01746     /* 
01747        { // debug
01748        ORSA_ERROR("Read(Evolution *e)  ofdt = %i",last_ofdt_read);
01749        }
01750     */
01751     if (last_ofdt_read == OFDT_FRAME) {
01752       // the first is different from the others
01753       Read(&f);
01754       (*e)->push_back(f);
01755     }
01756     
01757     Read(&last_ofdt_read);
01758     /* 
01759        { // debug
01760        ORSA_ERROR("Read(Evolution *e)  ofdt = %i",last_ofdt_read);
01761        }
01762     */
01763     while (last_ofdt_read == OFDT_FRAME) {
01764       Read(&f,true);
01765       (*e)->push_back(f);
01766       Read(&last_ofdt_read);
01767     }
01768   }
01769   
01770   void OrsaFile::Write(Frame * f, bool write_only_r_v) {
01771     
01772     OrsaFileDataType t = OFDT_FRAME; Write(&t);
01773     
01774     UniverseTypeAwareTime f_time = *f;
01775     Write(&f_time);
01776     unsigned int n = f->size();
01777     Write(&n);
01778     // unsigned int j;
01779     Vector v;
01780     if (write_only_r_v) {
01781       for(unsigned int j=0;j<n;j++) {
01782         v = (*f)[j].position(); Write(&v);
01783         v = (*f)[j].velocity(); Write(&v);
01784       }
01785     } else {
01786       for(unsigned int j=0;j<n;j++) {
01787         Write(&((*f)[j]));
01788       }
01789     }
01790   }
01791   
01792   void OrsaFile::Read(Frame * f, bool read_only_r_v) {
01793     
01794     UniverseTypeAwareTime f_time;
01795     Read(&f_time);
01796     f->SetTime(f_time);
01797     unsigned int n = f->size();
01798     Read(&n);
01799     f->resize(n);
01800     unsigned int j;
01801     Vector v;
01802     if (read_only_r_v) {
01803       for(j=0;j<n;j++) {
01804         Read(&v); (*f)[j].SetPosition(v);
01805         Read(&v); (*f)[j].SetVelocity(v);
01806       }
01807     } else {
01808       for(j=0;j<n;j++) { 
01809         Read(&((*f)[j]));
01810       }
01811     }
01812   }
01813   
01814   void OrsaFile::Write(Body *b) {
01815     string b_name   = b->name();   Write(&b_name);
01816     double b_mass   = b->mass();   Write(&b_mass);
01817     double b_radius = b->radius(); Write(&b_radius);
01818     JPL_planets b_planet = b->JPLPlanet(); Write(&b_planet);
01819     Vector v;
01820     v = b->position(); Write(&v);
01821     v = b->velocity(); Write(&v);
01822   }
01823   
01824   void OrsaFile::Read(Body *b) {
01825     string b_name;   Read(&b_name);
01826     double b_mass;   Read(&b_mass);
01827     double b_radius; Read(&b_radius);
01828     JPL_planets b_planet; Read(&b_planet);
01829     
01830     *b = Body(b_name,b_mass,b_radius,b_planet);
01831     
01832     Vector v;
01833     Read(&v); b->SetPosition(v);
01834     Read(&v); b->SetVelocity(v);
01835   }
01836   
01837   void OrsaFile::Write(BodyWithEpoch * b) {
01838     Write((Body*)b);
01839     UniverseTypeAwareTime b_epoch = b->Epoch(); Write(&b_epoch);
01840   }
01841   
01842   void OrsaFile::Read(BodyWithEpoch * b) {
01843     Read((Body*)b);
01844     UniverseTypeAwareTime b_epoch; Read(&b_epoch); b->SetEpoch(b_epoch);
01845   }
01846   
01847   void OrsaFile::Write(Date *d) {
01848     double j = d->GetJulian(); Write(&j);
01849   }
01850   
01851   void OrsaFile::Read(Date *d) {
01852     double j; Read(&j); d->SetJulian(j);
01853   }
01854   
01855   void OrsaFile::Write(UniverseTypeAwareTime * t) {
01856     switch (universe->GetUniverseType()) {
01857     case Real: {
01858       Date d = t->GetDate(); Write(&d);
01859       break;
01860     }
01861     case Simulated: {
01862       double tt = t->GetTime(); Write(&tt);
01863       break;
01864     }
01865     }
01866   }
01867   
01868   void OrsaFile::Read(UniverseTypeAwareTime *t) {
01869     switch (universe->GetUniverseType()) {
01870     case Real: {
01871       Date d; Read(&d); t->SetDate(d);
01872       break;
01873     }
01874     case Simulated: {
01875       double tt; Read(&tt); t->SetTime(tt);
01876       break;
01877     }
01878     }
01879   }
01880   
01881   void OrsaFile::Write(UniverseTypeAwareTimeStep *ts_in) {
01882     switch (universe->GetUniverseType()) {
01883     case Real: {
01884       TimeStep _ts = ts_in->GetTimeStep(); Write(&_ts);
01885       break;
01886     }
01887     case Simulated: {
01888       double tt = ts_in->GetDouble(); Write(&tt);
01889       break;
01890     }
01891     }
01892   }
01893   
01894   void OrsaFile::Read(UniverseTypeAwareTimeStep *ts_in) {
01895     switch (universe->GetUniverseType()) {
01896     case Real: {
01897       TimeStep _ts; Read(&_ts); ts_in->SetTimeStep(_ts);
01898       break;
01899     }
01900     case Simulated: {
01901       double tt; Read(&tt); ts_in->SetDouble(tt);
01902       break;
01903     }
01904     }
01905   }
01906   
01907   // void OrsaFile::Write(Integrator **i) {
01908   void OrsaFile::Write(const Integrator * i) {
01909     
01910     //  IntegratorType it = (*i)->GetType();
01911     IntegratorType it = i->GetType();
01912     Write(&it);
01913     
01914     // UniverseTypeAwareTimeStep ts = (*i)->timestep;
01915     UniverseTypeAwareTimeStep ts = i->timestep;
01916     Write(&ts);
01917     
01918     // double a = (*i)->accuracy;
01919     double a = i->accuracy;
01920     Write(&a);
01921     
01922     // unsigned int m = (*i)->m;
01923     unsigned int m = i->m;
01924     Write(&m);
01925   }
01926   
01927   void OrsaFile::Read(Integrator **i) {
01928     
01929     IntegratorType type; Read(&type);
01930     make_new_integrator(i, type);
01931     
01932     UniverseTypeAwareTimeStep ts;
01933     Read(&ts);
01934     (*i)->timestep = ts;
01935     
01936     double a;       Read(&a);
01937     unsigned int m; Read(&m);
01938     
01939     (*i)->accuracy = a;
01940     (*i)->m        = m;
01941   }
01942   
01943   /* 
01944      void OrsaFile::Write(const Interaction * i) {
01945      InteractionType type = i->GetType(); Write(&type); 
01946      }
01947      
01948      void OrsaFile::Read(Interaction **i) {
01949      InteractionType type; Read(&type);
01950      make_new_interaction(i, type);
01951      }
01952   */
01953   
01954   void OrsaFile::Write(const Interaction * i) {
01955     InteractionType type = i->GetType(); Write(&type); 
01956     bool b = i->IsSkippingJPLPlanets();  Write(&b);
01957     if (type == NEWTON) {
01958       const Newton * newton = dynamic_cast <const Newton *> (i);
01959       if (newton) {
01960         b = newton->IsIncludingMultipoleMoments();        Write(&b);
01961         b = newton->IsIncludingRelativisticEffects();     Write(&b);
01962         b = newton->IsIncludingFastRelativisticEffects(); Write(&b);
01963       } else {
01964         b = false;
01965         Write(&b);
01966         Write(&b);
01967         Write(&b);
01968       }
01969     }
01970   }
01971   
01972   void OrsaFile::Read(Interaction **i) {
01973     InteractionType type; Read(&type);
01974     make_new_interaction(i, type);
01975     bool b; Read(&b); (*i)->SkipJPLPlanets(b);
01976     if (type == NEWTON) {
01977       Newton * newton = dynamic_cast <Newton *> (*i);
01978       if (newton) {
01979         Read(&b); newton->IncludeMultipoleMoments(b);
01980         Read(&b); newton->IncludeRelativisticEffects(b);
01981         Read(&b); newton->IncludeFastRelativisticEffects(b);
01982       } else {
01983         b = false;
01984         Read(&b);
01985         Read(&b);
01986         Read(&b);
01987       }
01988     }
01989   }
01990   
01991   void OrsaFile::Write(std::string * s) {
01992     const unsigned int size = s->size();
01993     unsigned int n = 1 + size;
01994     Write(&n);
01995     char * name = (char*)malloc(n*sizeof(char));
01996     // 
01997     // strcpy(name,s->c_str());
01998     {
01999       unsigned int i;
02000       for (i=0;i<size;++i) {
02001         name[i] = (*s)[i];
02002       }
02003       name[size] = '\0';
02004     }
02005     //
02006     WRITE_FILE(name,sizeof(char),n,file);
02007     /* 
02008     ORSA_ERROR("Write(std::string *s)  n = %i   s = %s   s->size() = %i   strlen(s->c_str()) = %i",
02009       n,s->c_str(),s->size(),strlen(s->c_str()));
02010     */
02011     free(name);
02012     { // check
02013       if (strlen(s->c_str()) > n) {
02014         ORSA_ERROR("string length problem...");
02015       }
02016     }
02017   }
02018   
02019   void OrsaFile::Read(std::string * s) {
02020     unsigned int n; 
02021     Read(&n);
02022     if (n > 0) {
02023       char * name = (char*)malloc(n*sizeof(char));
02024       READ_FILE(name,sizeof(char),n,file);
02025       *s = name;
02026       /* 
02027       ORSA_ERROR("Read(std::string *s)  n = %i   s = %s   s->size() = %i   strlen(s->c_str()) = %i   name = %s",
02028         n,s->c_str(),s->size(),strlen(s->c_str()),name);
02029       */
02030       free(name);
02031     }
02032   }
02033   
02034   void OrsaFile::Write(orsa::Vector *v) {
02035     Write(&v->x); 
02036     Write(&v->y);
02037     Write(&v->z);
02038   }
02039   
02040   void OrsaFile::Read(orsa::Vector *v) {
02041     Read(&v->x); 
02042     Read(&v->y); 
02043     Read(&v->z);
02044   }
02045   
02046   void OrsaFile::Write(bool * b) {
02047     WRITE_FILE(b,sizeof(bool),1,file);
02048   }
02049   
02050   void OrsaFile::Read(bool * b) {
02051     read_swap(b,sizeof(bool));
02052   }
02053   
02054   void OrsaFile::Write(unsigned int * i) {
02055     WRITE_FILE(i,sizeof(unsigned int),1,file);
02056   }
02057   
02058   void OrsaFile::Read(unsigned int * i) {
02059     read_swap(i,sizeof(unsigned int));
02060   }
02061   
02062   void OrsaFile::Write(int *i) {
02063     WRITE_FILE(i,sizeof(int),1,file);
02064   }
02065   
02066   void OrsaFile::Read(int *i) {
02067     read_swap(i,sizeof(int));
02068   }
02069   
02070   void OrsaFile::Write(double * d) {
02071     WRITE_FILE(d,sizeof(double),1,file);
02072   }
02073   
02074   void OrsaFile::Read(double * d) {
02075     read_swap(d,sizeof(double));
02076   }
02077   
02078   /* 
02079      void OrsaFile::Write(bool *b) {
02080      // WRITE_FILE(b,sizeof(bool),1,file);
02081      unsigned int i = *b;
02082      Write(&i);
02083      }
02084      
02085      void OrsaFile::Read(bool *b) {
02086      // read_swap(b,sizeof(bool));
02087      unsigned int i;
02088      Read(&i);
02089      *b = (bool)i; 
02090      }
02091   */
02092   
02093   void OrsaFile::Write(IntegratorType *it) {
02094     // WRITE_FILE(it,sizeof(IntegratorType),1,file);
02095     unsigned int i = *it;
02096     Write(&i);
02097   }
02098   
02099   void OrsaFile::Read(IntegratorType *it) {
02100     // read_swap(it,sizeof(IntegratorType));
02101     unsigned int i;
02102     Read(&i);
02103     convert(*it,i);
02104   }
02105   
02106   void OrsaFile::Write(InteractionType *it) {
02107     // WRITE_FILE(it,sizeof(InteractionType),1,file);
02108     unsigned int i = *it;
02109     Write(&i);
02110   }
02111   
02112   void OrsaFile::Read(InteractionType *it) {
02113     // read_swap(it,sizeof(InteractionType));
02114     unsigned int i;
02115     Read(&i);
02116     convert(*it,i);
02117   }
02118   
02119   void OrsaFile::Write(time_unit *tu) {
02120     // WRITE_FILE(tu,sizeof(time_unit),1,file);
02121     unsigned int i = *tu;
02122     Write(&i);
02123   }
02124   
02125   void OrsaFile::Read(time_unit *tu) {
02126     // read_swap(tu,sizeof(time_unit));
02127     unsigned int i;
02128     Read(&i);
02129     convert(*tu,i);
02130   }
02131   
02132   void OrsaFile::Write(length_unit *lu) {
02133     // WRITE_FILE(lu,sizeof(length_unit),1,file);
02134     unsigned int i = *lu;
02135     Write(&i);
02136   }
02137   
02138   void OrsaFile::Read(length_unit *lu) {
02139     // read_swap(lu,sizeof(length_unit));
02140     unsigned int i;
02141     Read(&i); 
02142     convert(*lu,i);
02143   }
02144   
02145   void OrsaFile::Write(mass_unit *mu) {
02146     // WRITE_FILE(mu,sizeof(mass_unit),1,file);
02147     unsigned int i = *mu;
02148     Write(&i);
02149   }
02150   
02151   void OrsaFile::Read(mass_unit *mu) {
02152     // read_swap(mu,sizeof(mass_unit));
02153     unsigned int i;
02154     Read(&i);
02155     convert(*mu,i);
02156   }
02157   
02158   void OrsaFile::Write(ReferenceSystem *rs) {
02159     // WRITE_FILE(rs,sizeof(ReferenceSystem),1,file);
02160     unsigned int i = *rs;
02161     Write(&i);
02162   }
02163   
02164   void OrsaFile::Read(ReferenceSystem *rs) {
02165     // read_swap(rs,sizeof(ReferenceSystem));
02166     unsigned int i;
02167     Read(&i);
02168     convert(*rs,i);
02169   }
02170   
02171   void OrsaFile::Write(UniverseType *ut) {
02172     // WRITE_FILE(ut,sizeof(UniverseType),1,file);
02173     unsigned int i = *ut;
02174     Write(&i);
02175   }
02176   
02177   void OrsaFile::Read(UniverseType *ut) {
02178     // read_swap(ut,sizeof(UniverseType));
02179     unsigned int i;
02180     Read(&i);
02181     convert(*ut,i);
02182   }
02183   
02184   void OrsaFile::Write(TimeScale *ts) {
02185     // WRITE_FILE(ts,sizeof(TimeScale),1,file);
02186     unsigned int i = *ts;
02187     Write(&i);
02188   }
02189   
02190   void OrsaFile::Read(TimeScale *ts) {
02191     // read_swap(ts,sizeof(TimeScale));
02192     unsigned int i;
02193     Read(&i);
02194     convert(*ts,i);
02195   }
02196   
02197   void OrsaFile::Write(OrsaFileDataType *ofdt) {
02198     // WRITE_FILE(ofdt,sizeof(OrsaFileDataType),1,file);
02199     unsigned int i = *ofdt;
02200     Write(&i);
02201   }
02202   
02203   /* 
02204      void OrsaFile::Read(OrsaFileDataType *ofdt) {
02205      const int val = read_swap(ofdt,sizeof(OrsaFileDataType));
02206      if (val==0) *ofdt = OFDT_END_OF_FILE;
02207      }
02208   */
02209   
02210   void OrsaFile::Read(OrsaFileDataType *ofdt) {
02211     // const int val = read_swap(ofdt,sizeof(OrsaFileDataType));
02212     // if (val==0) *ofdt = OFDT_END_OF_FILE;
02213     unsigned int i;
02214     const int val = read_swap(&i,sizeof(unsigned int));
02215     // convert(*ofdt,i);
02216     // if (val==0) *ofdt = OFDT_END_OF_FILE;
02217     if (val==0) {
02218       *ofdt = OFDT_END_OF_FILE;
02219     } else {
02220       convert(*ofdt,i);
02221     }
02222   }
02223   
02224   void OrsaFile::Write(JPL_planets *jp) {
02225     // WRITE_FILE(jp,sizeof(JPL_planets),1,file);
02226     unsigned int i = *jp;
02227     Write(&i);
02228   }
02229   
02230   void OrsaFile::Read(JPL_planets *jp) {
02231     // read_swap(jp,sizeof(JPL_planets));
02232     unsigned int i;
02233     Read(&i);
02234     convert(*jp,i);
02235   }
02236   
02237   void OrsaFile::Write(TimeStep * ts) {
02238     unsigned int days = ts->days(); 
02239     Write(&days);
02240     //
02241     unsigned int day_fraction = ts->day_fraction();
02242     Write(&day_fraction); 
02243     //
02244     int sign = ts->sign();
02245     Write(&sign);
02246   }
02247   
02248   void OrsaFile::Read(TimeStep * ts) {
02249     unsigned int days;
02250     Read(&days);
02251     unsigned int day_fraction;
02252     Read(&day_fraction);
02253     int sign;
02254     Read(&sign);
02255     *ts = TimeStep(days,day_fraction,sign);
02256   }
02257   
02258   // very simple check, better than nothing for the moment...
02259   bool OrsaFile::GoodFile(const string &filename) {
02260     
02261     // define locally this variables, to allow
02262     // the GoodFile() method to be static
02263     unsigned int byte_order;
02264     FILE_TYPE    file;
02265     
02266     if ((file = OPEN_FILE(filename.c_str(),OPEN_READ)) == 0) return false;
02267     
02268     READ_FILE(&byte_order,sizeof(byte_order),1,file);
02269     
02270     if ( (byte_order != 1234) &&
02271          (byte_order != 4321) ) {
02272       
02273       swap(&byte_order,sizeof(byte_order));
02274       
02275       if ( (byte_order != 1234) &&
02276            (byte_order != 4321) ) {
02277         
02278         CLOSE_FILE(file);
02279         return false;
02280         
02281       }
02282     }
02283     
02284     CLOSE_FILE(file);
02285     return true;
02286   }
02287   
02288   //! Modified Radau input files
02289   RadauModIntegrationFile::RadauModIntegrationFile(OrbitStream &osin) : ReadFile() {
02290     os = &osin;
02291   }
02292   
02293   void RadauModIntegrationFile::Read() {
02294     
02295     // if (status == CLOSE) Open();
02296     
02297     Open();
02298     
02299     if (status != OPEN_R){ 
02300       ORSA_ERROR("problems encountered when opening file.");
02301       return;
02302     }
02303     
02304     os->clear();
02305     os->timestep = 0.0;
02306     OrbitWithEpoch fo;
02307     REWIND_FILE(file); 
02308     
02309     double a,e,i,omega_per,omega_nod,M;
02310     double time,time_old=0,timestep;
02311     
02312     char line[1024];
02313     
02314     /* while ( (fscanf(file,"%lf %lf %lf %lf %lf %lf %lf",
02315        &time,&a,&e,&i,&M,&omega_per,&omega_nod)) != EOF ) {
02316     */
02317     
02318     while (GETS_FILE(line,1024,file) != 0) {
02319       
02320       sscanf(line,"%lf %lf %lf %lf %lf %lf %lf",
02321              &time,&a,&e,&i,&M,&omega_per,&omega_nod);
02322       
02323       timestep  = time - time_old;
02324       time_old  = time;
02325       if (os->size() == 2) { 
02326         os->timestep = FromUnits(timestep,YEAR); // read in days, save in the current time units
02327         // cerr << "timestep set to: " << os->timestep << endl;
02328       }
02329       
02330       fo.epoch = FromUnits(time,YEAR); // read in days, save in the current time units
02331       
02332       fo.a                = FromUnits(a,AU);
02333       fo.e                = e;
02334       fo.i                = (pi/180.0)*i;
02335       fo.omega_node       = (pi/180.0)*omega_nod;
02336       fo.omega_pericenter = (pi/180.0)*omega_per;
02337       fo.M                = (pi/180.0)*M;
02338       
02339       os->push_back(fo);  
02340       
02341       // QUICK AND DIRTY!
02342       if (fo.e >= 1.0) {
02343         ORSA_ERROR("reading eccentricity > 1.0, returning.");
02344         return;
02345       }
02346       
02347     }
02348     
02349   }
02350   
02351   // AstDySMatrixFile
02352   
02353   AstDySMatrixFile::AstDySMatrixFile() : AsteroidDatabaseFile() {
02354     db = new AsteroidDatabase();
02355   }
02356   
02357   AstDySMatrixFile::~AstDySMatrixFile() {
02358     delete db;
02359     db = 0;
02360   }
02361   
02362   void AstDySMatrixFile::Read() {
02363     
02364     Open();
02365     
02366     if (status != OPEN_R) {
02367       ORSA_ERROR("Status error!");
02368       return;
02369     }
02370     
02371     REWIND_FILE(file);
02372     
02373     char line[1024],tag[10];
02374     string stag;
02375     const string end_of_header="END_OF_HEADER";
02376     
02377     // skip header
02378     while (GETS_FILE(line,1024,file)!=0) {
02379       sscanf(line,"%s",tag);
02380       stag = tag;
02381       if (stag == end_of_header) {
02382         break;
02383       }
02384     }
02385     
02386     Asteroid ast;
02387     
02388     double cov_content[21];
02389     int cov_line;
02390     
02391     unsigned int local_index = 0;
02392     bool bool_stop=false;
02393     bool bool_pause=false;
02394     
02395     while (GETS_FILE(line,1024,file)!=0) {
02396       
02397       if (line[0] == '!') continue; // comment/header line
02398       
02399       if (line[0] == ' ') continue; // not the line number
02400       
02401       ++local_index;
02402       read_progress(local_index,bool_pause,bool_stop);
02403       
02404       if (bool_stop) break;
02405       
02406       while (bool_pause) {
02407         sleep(1);
02408         read_progress(local_index,bool_pause,bool_stop);
02409       }
02410       
02411       sscanf(line,"%s",tag);
02412       stag = tag;
02413       remove_leading_trailing_spaces(stag);     
02414       
02415       ast.name = stag;
02416       
02417       {
02418         const string name=stag;
02419         char c;
02420         unsigned int ck;
02421         bool is_only_digit=true;
02422         for (ck=0;ck<name.size();++ck) {
02423           c = name[ck];
02424           if (isalpha(c)) { 
02425             is_only_digit=false;
02426             break;
02427           }
02428         }
02429         
02430         if (is_only_digit) {
02431           ast.n = atoi(name.c_str());
02432         } else {
02433           ast.n = 0;
02434         }
02435       }
02436       
02437 #ifdef HAVE_GSL
02438       OrbitWithCovarianceMatrixGSL orbit;
02439 #else
02440       OrbitWithEpoch orbit;
02441 #endif
02442       
02443       orbit.mu = GetG()*GetMSun();
02444       
02445       cov_line = 0;
02446       
02447       while (GETS_FILE(line,1024,file)!=0) {
02448         
02449         if (line[0] == '!') continue; // comment/header line
02450         
02451         if (line[0] != ' ') break; // new asteroid
02452         
02453         sscanf(line,"%s",tag);
02454         stag = tag;
02455         
02456         if (stag == "EQU") {
02457           
02458           double x[6];
02459           
02460           sscanf(line,"%s %lg %lg %lg %lg %lg %lg",tag,&x[0],&x[1],&x[2],&x[3],&x[4],&x[5]);
02461           
02462           const double omega_tilde = secure_atan2(x[1],x[2]);
02463           
02464           orbit.a                  = x[0];
02465           orbit.e                  = sqrt(x[1]*x[1]+x[2]*x[2]);
02466           orbit.i                  = 2.0*atan(sqrt(x[3]*x[3]+x[4]*x[4]));
02467           orbit.omega_node         = fmod(10.0*twopi+secure_atan2(x[3],x[4]),twopi);
02468           orbit.omega_pericenter   = fmod(10.0*twopi+omega_tilde-orbit.omega_node,twopi);
02469           orbit.M                  = fmod(10.0*twopi+(pi/180)*x[5]-omega_tilde,twopi);
02470           
02471         } else if (stag == "MJD") {
02472           
02473           double epoch_MJD;
02474           sscanf(line,"%s %lg",tag,&epoch_MJD);
02475           Date  epoch;
02476           epoch.SetJulian(epoch_MJD+2400000.5,TDT);
02477           
02478           orbit.epoch.SetDate(epoch);
02479           
02480         } else if (stag == "MAG") {
02481           
02482           double mag, m2;
02483           sscanf(line,"%s %lg %lg",tag,&mag,&m2);
02484           
02485           ast.mag = mag;
02486           
02487         } else if (stag == "COV") {
02488           
02489           double c0,c1,c2;          
02490           
02491           sscanf(line,"%s %lg %lg %lg",tag,&c0,&c1,&c2);
02492           
02493           cov_content[cov_line*3]   = c0;
02494           cov_content[cov_line*3+1] = c1;
02495           cov_content[cov_line*3+2] = c2;
02496           
02497           cov_line++;
02498         }
02499         
02500         if (cov_line==7) break;
02501         
02502       }
02503       
02504 #ifdef HAVE_GSL
02505       if (cov_line==7) {
02506         double covm[6][6];
02507         int i,k;
02508         int ik=0;
02509         for (i=0;i<6;++i) {
02510           for (k=i;k<6;++k) {
02511             // IMPORTANT: units correction
02512             if (i==5) cov_content[ik] *= (pi/180);
02513             if (k==5) cov_content[ik] *= (pi/180);
02514             //
02515             covm[i][k] = cov_content[ik];
02516             covm[k][i] = covm[i][k];
02517             //
02518             ++ik;
02519           }
02520         }       
02521         orbit.SetCovarianceMatrix((double**)covm,Equinoctal);
02522       } else {
02523         ORSA_ERROR("Cannot set covariance matrix for object %s from file %s.",ast.name.c_str(),filename.c_str());
02524       }      
02525 #endif // HAVE_GSL
02526       
02527       ast.orb = orbit;
02528       
02529       // CHECK THIS CODE!
02530       // NB: DON'T KNOW HOW TO 'ROTATE' THE COVARIANCE MATRIX
02531       /* 
02532          switch (universe->GetReferenceSystem()) {
02533          case ECLIPTIC: break;
02534          case EQUATORIAL:
02535          { 
02536          Date tmp_date(TDT);
02537          // cerr << "Rotating astorb orbit..." << endl;
02538          const double obleq_rad = obleq(tmp_date).GetRad();
02539          Vector position,velocity;
02540          ast.orb.RelativePosVel(position,velocity);
02541          position.rotate(0.0,obleq_rad,0.0);
02542          velocity.rotate(0.0,obleq_rad,0.0);
02543          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
02544          }
02545          break;
02546          }
02547       */
02548       //
02549       switch (universe->GetReferenceSystem()) {
02550       case ECLIPTIC: break;
02551       case EQUATORIAL:
02552         { 
02553           Vector position,velocity;
02554           ast.orb.RelativePosVel(position,velocity);
02555           EclipticToEquatorial_J2000(position);
02556           EclipticToEquatorial_J2000(velocity);
02557           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
02558         }
02559         break;
02560       }
02561       
02562       db->push_back(ast);
02563     }
02564     
02565     read_finished();
02566   }
02567   
02568   // JPL DASTCOM files
02569   
02570   JPLDastcomNumFile::JPLDastcomNumFile() : AsteroidDatabaseFile() {
02571     db = new AsteroidDatabase();
02572   }
02573   
02574   JPLDastcomNumFile::~JPLDastcomNumFile() {
02575     delete db;
02576     db = 0;
02577   }
02578   
02579   void JPLDastcomNumFile::Read() {
02580     
02581     Open();
02582     
02583     if (status != OPEN_R) {
02584       ORSA_ERROR("Status error!");
02585       return;
02586     }
02587     
02588     db->clear();
02589     
02590     char line[300];
02591     
02592     double a,e,i,omega_node,omega_pericenter,M;
02593     // int    n;
02594     string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
02595     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
02596     // string ceu;
02597     
02598     string year,month,day;
02599     // int    y,m,d;
02600     
02601     Asteroid ast;
02602     
02603     // Date tmp_date(TDT);
02604     Date tmp_date;
02605     
02606     unsigned int local_index = 0;
02607     bool bool_stop=false;
02608     bool bool_pause=false;
02609     REWIND_FILE(file);
02610     while ((GETS_FILE(line,300,file)) != 0) {
02611       
02612       if (strlen(line) < 100) continue; // not a good line, maybe a comment or a white line...
02613       
02614       if (line[0]=='-') continue; // comment
02615       
02616       local_index++;
02617       read_progress(local_index,bool_pause,bool_stop);
02618       
02619       if (bool_stop) break;
02620       
02621       while (bool_pause) {
02622         sleep(1);
02623         read_progress(local_index,bool_pause,bool_stop);
02624       }
02625       
02626       // uncomment the ones used
02627       number.assign(line,0,5);
02628       name.assign(line,6,17); 
02629       
02630       epoch.assign(line,24,5);
02631       //
02632       semimajor_axis.assign(line,30,10);
02633       eccentricity.assign(line,41,10);
02634       inclination.assign(line,52,9);
02635       pericenter.assign(line,62,9);
02636       node.assign(line,72,9);
02637       mean_anomaly.assign(line,82,11);
02638       //////////////
02639       
02640       ast.n = atoi(number.c_str());
02641       
02642       ast.name = name;
02643       remove_leading_trailing_spaces(ast.name);
02644       
02645       a                = atof(semimajor_axis.c_str());
02646       e                = atof(eccentricity.c_str());
02647       i                = (pi/180)*atof(inclination.c_str());
02648       omega_node       = (pi/180)*atof(node.c_str());
02649       omega_pericenter = (pi/180)*atof(pericenter.c_str());
02650       M                = (pi/180)*atof(mean_anomaly.c_str());
02651       
02652       // checks
02653       if ((ast.n==0) || (a==0.0)) {
02654         // bad line...
02655         continue;
02656       }
02657       
02658       ast.orb.a                = FromUnits(a,AU);
02659       ast.orb.e                = e;
02660       ast.orb.i                = i;
02661       ast.orb.omega_node       = omega_node;
02662       ast.orb.omega_pericenter = omega_pericenter;
02663       ast.orb.M                = M;
02664       
02665       // year.assign(epoch,0,4);
02666       // month.assign(epoch,4,2);
02667       // day.assign(epoch,6,2);
02668       
02669       // y = atoi(year.c_str());
02670       // m = atoi(month.c_str());
02671       // d = atoi(day.c_str());
02672       
02673       tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
02674       ast.orb.epoch.SetDate(tmp_date);
02675       // ast.orb.T = sqrt(4*pisq/(GetG()*GetMSun())*pow(FromUnits(ast.orb.a,AU),3));
02676       ast.orb.mu = GetG()*GetMSun();
02677       // ast.orb.ref_body = orb_ref_body;
02678       
02679       /* 
02680          switch (universe->GetReferenceSystem()) {
02681          case ECLIPTIC: break;
02682          case EQUATORIAL:
02683          { 
02684          // cerr << "Rotating astorb orbit..." << endl;
02685          const double obleq_rad = obleq(tmp_date).GetRad();
02686          Vector position,velocity;
02687          ast.orb.RelativePosVel(position,velocity);
02688          position.rotate(0.0,obleq_rad,0.0);
02689          velocity.rotate(0.0,obleq_rad,0.0);
02690          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
02691          }
02692          break;
02693          }
02694       */
02695       
02696       switch (universe->GetReferenceSystem()) {
02697       case ECLIPTIC: break;
02698       case EQUATORIAL:
02699         { 
02700           Vector position,velocity;
02701           ast.orb.RelativePosVel(position,velocity);
02702           EclipticToEquatorial_J2000(position);
02703           EclipticToEquatorial_J2000(velocity);
02704           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
02705         }
02706         
02707         break;
02708       }
02709       
02710       db->push_back(ast);
02711     }
02712     
02713     read_finished();
02714   }
02715   
02716   ////
02717 
02718   JPLDastcomUnnumFile::JPLDastcomUnnumFile() : AsteroidDatabaseFile() {
02719     db = new AsteroidDatabase();
02720   }
02721   
02722   JPLDastcomUnnumFile::~JPLDastcomUnnumFile() {
02723     delete db;
02724     db = 0;
02725   }
02726   
02727   void JPLDastcomUnnumFile::Read() {
02728     
02729     Open();
02730     
02731     if (status != OPEN_R) {
02732       ORSA_ERROR("Status error!");
02733       return;
02734     }
02735     
02736     db->clear();
02737     
02738     char line[300];
02739     
02740     double a,e,i,omega_node,omega_pericenter,M;
02741     // int    n;
02742     string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
02743     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
02744     // string ceu;
02745     
02746     string year,month,day;
02747     // int    y,m,d;
02748     
02749     Asteroid ast;
02750     
02751     // Date tmp_date(TDT);
02752     Date tmp_date;
02753     
02754     unsigned int local_index = 0;
02755     bool bool_stop=false;
02756     bool bool_pause=false;
02757     REWIND_FILE(file);
02758     while ((GETS_FILE(line,300,file)) != 0) {
02759       
02760       if (strlen(line) < 100) continue; // not a good line, maybe a comment or a white line...
02761       
02762       if (line[0]=='-') continue; // comment
02763       
02764       local_index++;
02765       read_progress(local_index,bool_pause,bool_stop);
02766       
02767       if (bool_stop) break;
02768       
02769       while (bool_pause) {
02770         sleep(1);
02771         read_progress(local_index,bool_pause,bool_stop);
02772       }
02773       
02774       // uncomment the ones used
02775       // number.assign(line,0,5);
02776       name.assign(line,0,11); 
02777       
02778       epoch.assign(line,12,5);
02779       //
02780       semimajor_axis.assign(line,18,11);
02781       eccentricity.assign(line,30,10);
02782       inclination.assign(line,41,9);
02783       pericenter.assign(line,51,9);
02784       node.assign(line,61,9);
02785       mean_anomaly.assign(line,71,11);
02786       //////////////
02787       
02788       ast.n = atoi(number.c_str());
02789       
02790       ast.name = name;
02791       remove_leading_trailing_spaces(ast.name);
02792       
02793       a                = atof(semimajor_axis.c_str());
02794       e                = atof(eccentricity.c_str());
02795       i                = (pi/180)*atof(inclination.c_str());
02796       omega_node       = (pi/180)*atof(node.c_str());
02797       omega_pericenter = (pi/180)*atof(pericenter.c_str());
02798       M                = (pi/180)*atof(mean_anomaly.c_str());
02799       
02800       // checks
02801       if ((a==0.0)) {
02802         // bad line...
02803         continue;
02804       }
02805       
02806       ast.orb.a                = FromUnits(a,AU);
02807       ast.orb.e                = e;
02808       ast.orb.i                = i;
02809       ast.orb.omega_node       = omega_node;
02810       ast.orb.omega_pericenter = omega_pericenter;
02811       ast.orb.M                = M;
02812       
02813       // year.assign(epoch,0,4);
02814       // month.assign(epoch,4,2);
02815       // day.assign(epoch,6,2);
02816       
02817       // y = atoi(year.c_str());
02818       // m = atoi(month.c_str());
02819       // d = atoi(day.c_str());
02820       
02821       tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
02822       ast.orb.epoch.SetDate(tmp_date);
02823       // ast.orb.T = sqrt(4*pisq/(GetG()*GetMSun())*pow(FromUnits(ast.orb.a,AU),3));
02824       ast.orb.mu = GetG()*GetMSun();
02825       // ast.orb.ref_body = orb_ref_body;
02826       
02827       /* 
02828          switch (universe->GetReferenceSystem()) {
02829          case ECLIPTIC: break;
02830          case EQUATORIAL:
02831          { 
02832          // cerr << "Rotating astorb orbit..." << endl;
02833          const double obleq_rad = obleq(tmp_date).GetRad();
02834          Vector position,velocity;
02835          ast.orb.RelativePosVel(position,velocity);
02836          position.rotate(0.0,obleq_rad,0.0);
02837          velocity.rotate(0.0,obleq_rad,0.0);
02838          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
02839          }
02840          break;
02841          }
02842       */
02843       
02844       switch (universe->GetReferenceSystem()) {
02845       case ECLIPTIC: break;
02846       case EQUATORIAL:
02847         { 
02848           Vector position,velocity;
02849           ast.orb.RelativePosVel(position,velocity);
02850           EclipticToEquatorial_J2000(position);
02851           EclipticToEquatorial_J2000(velocity);
02852           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
02853         }
02854         break;
02855       }
02856       
02857       db->push_back(ast);
02858     }
02859     
02860     read_finished();
02861   }
02862   
02863   ////
02864   
02865   JPLDastcomCometFile::JPLDastcomCometFile() : AsteroidDatabaseFile() {
02866     db = new AsteroidDatabase();
02867   }
02868   
02869   JPLDastcomCometFile::~JPLDastcomCometFile() {
02870     delete db;
02871     db = 0;
02872   }
02873   
02874   void JPLDastcomCometFile::Read() {
02875       
02876     Open();
02877     
02878     if (status != OPEN_R) {
02879       ORSA_ERROR("Status error!");
02880       return;
02881     }
02882     
02883     db->clear();
02884     
02885     char line[300];
02886     
02887     double a,e,i,omega_node,omega_pericenter,M;
02888     string number,type,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
02889     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
02890     // string ceu;
02891     string pericenter_distance,pericenter_epoch;
02892     
02893     string year,month,day;
02894     int    y=0,m=0;
02895     double frac_day;
02896     
02897     Asteroid ast;
02898     
02899     double q;
02900     
02901     REWIND_FILE(file);
02902     
02903     // Date tmp_date(TDT);
02904     Date tmp_date;
02905     
02906     unsigned int local_index = 0;
02907     bool bool_stop=false;
02908     bool bool_pause=false;
02909     
02910     while ( (GETS_FILE(line,300,file)) != 0 ) {
02911       
02912       // if (strlen(line) < 90) continue; // not a good line, maybe a comment or a white line...
02913       
02914       if (line[0]=='-') continue; // comment
02915       
02916       ++local_index;
02917       read_progress(local_index,bool_pause,bool_stop);
02918       
02919       if (bool_stop) break;
02920       
02921       while (bool_pause) {
02922         // cerr << "AstorbFile::Read() sleeping..." << endl;
02923         sleep(1);
02924         read_progress(local_index,bool_pause,bool_stop);
02925       }
02926       
02927       // uncomment the ones used
02928       // number.assign(line,0,4);
02929       // type.assign(line,4,1);
02930       // name.assign(line,5,7);
02931       name.assign(line,0,37);
02932       epoch.assign(line,38,5);
02933       // cerr << "comet name: " << name << endl;
02934       pericenter_distance.assign(line,44,10);
02935       eccentricity.assign(line,55,10);
02936       inclination.assign(line,66,9);
02937       pericenter.assign(line,76,9);
02938       node.assign(line,86,9);
02939       //
02940       pericenter_epoch.assign(line,96,14);
02941       
02942       // conversions
02943       
02944       ast.name = name;
02945       remove_leading_trailing_spaces(ast.name);
02946       
02947       // ast.n = 0; // arbitrary, for the moment
02948       ast.n = 0;
02949       
02950       // ast.mag  = atof(absolute_magnitude.c_str());
02951       
02952       // a                = atof(semimajor_axis.c_str());
02953       e  = atof(eccentricity.c_str());
02954       
02955       // to be tested...
02956       q = atof(pericenter_distance.c_str());
02957       if (e == 1.0) {
02958         a = q;
02959       } else {
02960         a = q/fabs(1.0-e);
02961       }
02962       
02963       // checks
02964       if ((q==0.0)) {
02965         // bad line...
02966         continue;
02967       }
02968       
02969       i                = (pi/180)*atof(inclination.c_str());
02970       omega_node       = (pi/180)*atof(node.c_str());
02971       omega_pericenter = (pi/180)*atof(pericenter.c_str());
02972       // M                = (pi/180)*atof(mean_anomaly.c_str());
02973       
02974       tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
02975       ast.orb.epoch.SetDate(tmp_date);
02976       
02977       year.assign(pericenter_epoch,0,4);
02978       month.assign(pericenter_epoch,4,2);
02979       day.assign(pericenter_epoch,6,8);
02980       
02981       y = atoi(year.c_str());
02982       m = atoi(month.c_str());
02983       frac_day = atof(day.c_str());
02984       
02985       Date peri_date;
02986       peri_date.SetGregor(y,m,frac_day,TDT);
02987       UniverseTypeAwareTime pericenter_passage(peri_date);
02988       
02989       ast.orb.mu = GetG()*GetMSun();
02990       
02991       ast.orb.a                = FromUnits(a,AU);
02992       ast.orb.e                = e;
02993       ast.orb.i                = i;
02994       ast.orb.omega_node       = omega_node;
02995       ast.orb.omega_pericenter = omega_pericenter;
02996       //
02997       M = ((ast.orb.epoch.Time() - pericenter_passage.Time())/ast.orb.Period())*twopi;  
02998       M = fmod(10*twopi+fmod(M,twopi),twopi);
02999       //
03000       ast.orb.M                = M;
03001       
03002       // cerr << "comet: " << ast.name << "  q: " << q << "  e: " << e << "  i: " << i*(180/pi) << endl;
03003       
03004       /* 
03005          switch (universe->GetReferenceSystem()) {
03006          case ECLIPTIC: break;
03007          case EQUATORIAL:
03008          { 
03009          // cerr << "Rotating astorb orbit..." << endl;
03010          const double obleq_rad = obleq(tmp_date).GetRad();
03011          Vector position,velocity;
03012          ast.orb.RelativePosVel(position,velocity);
03013          position.rotate(0.0,obleq_rad,0.0);
03014          velocity.rotate(0.0,obleq_rad,0.0);
03015          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
03016          }
03017          break;
03018          }
03019       */
03020       
03021       switch (universe->GetReferenceSystem()) {
03022       case ECLIPTIC: break;
03023       case EQUATORIAL:
03024         { 
03025           Vector position,velocity;
03026           ast.orb.RelativePosVel(position,velocity);
03027           EclipticToEquatorial_J2000(position);
03028           EclipticToEquatorial_J2000(velocity);
03029           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
03030         }
03031         break;
03032       }
03033       
03034       db->push_back(ast);
03035       
03036     }
03037     
03038     read_finished();
03039   }
03040   
03041   ////
03042   
03043   // NEODyS .cat file
03044   
03045   // NEODYSCAT::NEODYSCAT() : ReadFile() {
03046   NEODYSCAT::NEODYSCAT() : AsteroidDatabaseFile() {
03047     // status = CLOSE;
03048     db = new AsteroidDatabase();
03049   }
03050   
03051   NEODYSCAT::~NEODYSCAT() {
03052     delete db;
03053     db = 0;
03054   }
03055   
03056   void NEODYSCAT::Read() {
03057     
03058     // if (status == CLOSE) Open();
03059     
03060     Open();
03061     
03062     if (status != OPEN_R) {
03063       ORSA_ERROR("Status error!");
03064       return;
03065     }
03066     
03067     db->clear();
03068     
03069     char line[300];
03070     
03071     double a,e,i,omega_node,omega_pericenter,M;
03072     // int    n;
03073     string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
03074     string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
03075     // string ceu;
03076     
03077     string year,month,day;
03078     // int    y,m,d;
03079     
03080     Asteroid ast;
03081     
03082     // Date tmp_date(TDT);
03083     Date tmp_date;
03084     
03085     unsigned int local_index = 0;
03086     bool bool_stop=false;
03087     bool bool_pause=false;
03088     REWIND_FILE(file);
03089     while ((GETS_FILE(line,300,file)) != 0) {
03090       
03091       if (strlen(line) < 100) continue; // not a good line, maybe a comment or a white line...
03092       
03093       if (line[0]=='!') continue; // comment
03094       
03095       local_index++;
03096       read_progress(local_index,bool_pause,bool_stop);
03097       
03098       if (bool_stop) break;
03099       
03100       while (bool_pause) {
03101         sleep(1);
03102         read_progress(local_index,bool_pause,bool_stop);
03103       }
03104       
03105       // uncomment the ones used
03106       // number.assign(line,0,5);
03107       name.assign(line,0,12); 
03108       //
03109       {
03110         // remove -->'<--
03111         string::size_type pos;
03112         while ((pos=name.find("'",0)) != string::npos) {
03113           // cerr << "name: " << name << "  pos: " << pos << endl;
03114           name.erase(pos,1);
03115         }
03116         // cerr << "final name: " << name << endl;
03117       }
03118       
03119       // orbit_computer.assign(line,25,15);
03120       // absolute_magnitude.assign(line,41,5);
03121       //
03122       // arc.assign(line,94,5);
03123       // numobs.assign(line,99,5);
03124       //
03125       epoch.assign(line,13,15);
03126       //
03127       semimajor_axis.assign(line,29,25);
03128       eccentricity.assign(line,54,25);
03129       inclination.assign(line,79,25);
03130       node.assign(line,104,25);
03131       pericenter.assign(line,129,25);
03132       mean_anomaly.assign(line,154,25);
03133       //////////////
03134       
03135       // ast.n = atoi(number.c_str());
03136       // ast.n = 0;
03137       {
03138         char c;
03139         unsigned int ck;
03140         bool is_only_digit=true;
03141         for (ck=0;ck<name.size();++ck) {
03142           c = name[ck];
03143           if (isalpha(c)) { 
03144             is_only_digit=false;
03145             break;
03146           }
03147         }
03148         
03149         if (is_only_digit) {
03150           ast.n = atoi(name.c_str());
03151         } else {
03152           ast.n = 0;
03153         }
03154       }
03155       
03156       ast.name = name;
03157       remove_leading_trailing_spaces(ast.name);
03158       
03159       // ast.ceu  = atof(ceu.c_str());
03160       
03161       // ast.mag  = atof(absolute_magnitude.c_str());
03162       
03163       a                = atof(semimajor_axis.c_str());
03164       e                = atof(eccentricity.c_str());
03165       i                = (pi/180)*atof(inclination.c_str());
03166       omega_node       = (pi/180)*atof(node.c_str());
03167       omega_pericenter = (pi/180)*atof(pericenter.c_str());
03168       M                = (pi/180)*atof(mean_anomaly.c_str());
03169       
03170       ast.orb.a                = FromUnits(a,AU);
03171       ast.orb.e                = e;
03172       ast.orb.i                = i;
03173       ast.orb.omega_node       = omega_node;
03174       ast.orb.omega_pericenter = omega_pericenter;
03175       ast.orb.M                = M;
03176       
03177       // year.assign(epoch,0,4);
03178       // month.assign(epoch,4,2);
03179       // day.assign(epoch,6,2);
03180       
03181       // y = atoi(year.c_str());
03182       // m = atoi(month.c_str());
03183       // d = atoi(day.c_str());
03184       
03185       tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
03186       ast.orb.epoch.SetDate(tmp_date);
03187       // ast.orb.T = sqrt(4*pisq/(GetG()*GetMSun())*pow(FromUnits(ast.orb.a,AU),3));
03188       ast.orb.mu = GetG()*GetMSun();
03189       // ast.orb.ref_body = orb_ref_body;
03190       
03191       /* 
03192          switch (universe->GetReferenceSystem()) {
03193          case ECLIPTIC: break;
03194          case EQUATORIAL:
03195          { 
03196          // cerr << "Rotating astorb orbit..." << endl;
03197          const double obleq_rad = obleq(tmp_date).GetRad();
03198          Vector position,velocity;
03199          ast.orb.RelativePosVel(position,velocity);
03200          position.rotate(0.0,obleq_rad,0.0);
03201          velocity.rotate(0.0,obleq_rad,0.0);
03202          ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
03203          }
03204          
03205          break;
03206          }
03207       */
03208       
03209       switch (universe->GetReferenceSystem()) {
03210       case ECLIPTIC: break;
03211       case EQUATORIAL:
03212         { 
03213           Vector position,velocity;
03214           ast.orb.RelativePosVel(position,velocity);
03215           EclipticToEquatorial_J2000(position);
03216           EclipticToEquatorial_J2000(velocity);
03217           ast.orb.Compute(position,velocity,ast.orb.mu,ast.orb.epoch);
03218         }
03219         break;
03220       }
03221       
03222       db->push_back(ast);
03223     }
03224     
03225     read_finished();
03226   }
03227   
03228   /// OrsaPaths
03229   
03230   char * OrsaPaths::path=0;
03231   char   OrsaPaths::_path_separator=0;
03232   
03233   OrsaPaths *orsa_paths=0;
03234   
03235   OrsaPaths::OrsaPaths() {
03236     if (orsa_paths) {
03237       ORSA_ERROR("cannot create two OrsaPaths from the same session");
03238       exit(0);
03239     }
03240     
03241     set_path_separator();
03242     set_path();
03243     
03244     orsa_paths = this;
03245   }
03246   
03247   OrsaPaths::OrsaPaths(const string &config_path) {
03248     if (orsa_paths) {
03249       ORSA_ERROR("cannot create two OrsaPaths from the same session");
03250       exit(0);
03251     }
03252     
03253     set_path_separator();
03254     path = strdup(config_path.c_str());
03255     
03256     orsa_paths = this;
03257   }
03258   
03259   void OrsaPaths::set_path_separator() {
03260     if (_path_separator!=0) return;
03261 #ifdef _WIN32
03262     _path_separator = '\\';
03263 #else
03264     _path_separator = '/';
03265 #endif
03266   }
03267   
03268   void OrsaPaths::set_path() {
03269     char p[2048];
03270     char * eh = getenv("HOME");
03271     p[0] = 0;
03272     if (eh != 0) strcpy(p, eh);
03273 #ifdef _WIN32
03274     char home[2048];
03275     if (SUCCEEDED(SHGetFolderPathA(NULL, CSIDL_PERSONAL | CSIDL_FLAG_CREATE, NULL, 0, home)))
03276       strcpy(p, home);
03277     strcat(p, "\\WINORSA\\");
03278 #else
03279     strcat(p, "/.orsa/");
03280 #endif
03281     path = strdup(p);
03282   }
03283   
03284   // static OrsaPaths p; // make sure the constructor gets called
03285   
03286   // TLE
03287   
03288   TLEFile::TLEFile() : ReadFile() {
03289     
03290   }
03291   
03292   void TLEFile::Read() {
03293     Open();
03294     if (status != OPEN_R) {
03295       ORSA_ERROR("Status error!");
03296       return;
03297     }
03298     sat.clear();
03299     string name;
03300     string s_tmp;
03301     int year=0;
03302     double days=0.0;
03303     double inclination=0.0,node=0.0,eccentricity=0.0,peri=0.0,M=0.0,period=0.0;
03304     bool have_one=false;
03305     bool have_two=false;
03306     char line[1024];
03307     unsigned int local_index = 0;
03308     while (GETS_FILE(line,1024,file) != 0) {
03309       
03310       if (line[0] == '1') {
03311         
03312         if (strlen(line) < 69) continue;
03313         
03314         if (isalpha(line[6])) continue; // test for single chars...
03315         
03316         s_tmp.assign(line,18,2);
03317         year = atoi(s_tmp.c_str());
03318         if (year > 70)
03319           year += 1900;
03320         else 
03321           year += 2000;
03322         
03323         s_tmp.assign(line,20,12);
03324         days = atof(s_tmp.c_str());
03325         
03326         have_one = true;
03327         have_two = false;
03328         
03329       } else if (line[0] == '2') {
03330         
03331         if (strlen(line) < 69) continue;
03332         
03333         if (!have_one) continue;
03334         
03335         if (isalpha(line[6])) continue; // test for single chars...
03336         
03337         s_tmp.assign(line,8,8);
03338         inclination = (pi/180.0)*atof(s_tmp.c_str());
03339         
03340         s_tmp.assign(line,17,8);
03341         node = (pi/180.0)*atof(s_tmp.c_str());
03342         
03343         s_tmp.assign(line,26,7);
03344         eccentricity = 1.0e-7*atof(s_tmp.c_str());      
03345         
03346         s_tmp.assign(line,34,8);
03347         peri = (pi/180.0)*atof(s_tmp.c_str());  
03348         
03349         s_tmp.assign(line,43,8);
03350         M = (pi/180.0)*atof(s_tmp.c_str());     
03351         
03352         s_tmp.assign(line,52,11);
03353         period = FromUnits(1.0/atof(s_tmp.c_str()),DAY);        
03354         
03355         have_two = true;
03356         
03357       } else {
03358         name.assign(line,0,MIN(24,strlen(line)-1)); // the last -1 is set to avoid the '\n' character in the name
03359         remove_leading_trailing_spaces(name);
03360         have_one = false;
03361         have_two = false;
03362       }
03363       
03364       if (have_one && have_two) {
03365         
03366         Date epoch;
03367         epoch.SetGregor(year,1,1,UTC); // UTC?
03368         double jd = epoch.GetJulian(UTC);
03369         jd += days-1.0;
03370         epoch.SetJulian(jd,UTC);
03371         
03372         JPLBody Earth(EARTH,epoch);
03373         
03374         Orbit orbit;
03375         orbit.mu = GetG()*Earth.mass();
03376         orbit.a = cbrt(period*period*orbit.mu/(4*pisq));
03377         orbit.e = eccentricity;
03378         orbit.i = inclination;
03379         orbit.omega_node       = node;
03380         orbit.omega_pericenter = peri;
03381         orbit.M                = M;
03382         
03383         Vector position,velocity;
03384         orbit.RelativePosVel(position,velocity);
03385         
03386         if (universe->GetReferenceSystem() == ECLIPTIC) {
03387           Angle obl = obleq_J2000();
03388           position.rotate(0.0,-obl.GetRad(),0.0);
03389           velocity.rotate(0.0,-obl.GetRad(),0.0);
03390         }
03391         
03392         position += Earth.position();
03393         velocity += Earth.velocity();
03394         
03395         sat.push_back(BodyWithEpoch(name,0.0,position,velocity,epoch));
03396         
03397         ++local_index;
03398         read_progress(local_index);
03399         
03400         // cerr << name << " period[DAYS]: " << FromUnits(period,DAY,-1) << "  a[ER]: " << FromUnits(orbit.a,ER,-1) << endl;
03401         
03402         have_one = have_two = false;
03403       }
03404     }
03405   }
03406   
03407 } // namespace orsa

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