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

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