00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_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>
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
00086 if (status == st) return;
00087
00088
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
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 void AstorbFile::Read() {
00199
00200
00201
00202
00203
00204 Open();
00205
00206 if (status != OPEN_R) {
00207 ORSA_ERROR("Status error!");
00208 return;
00209 }
00210
00211
00212
00213 db->clear();
00214
00215
00216
00217
00218
00219
00220
00221
00222 char line[300];
00223
00224 double a,e,i,omega_node,omega_pericenter,M;
00225
00226 string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
00227 string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
00228
00229
00230 string year,month,day;
00231 int y,m,d;
00232
00233
00234
00235
00236 Asteroid ast;
00237
00238
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
00248
00249 if (strlen(line) < 100) continue;
00250
00251
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
00260 sleep(1);
00261 read_progress(local_index,bool_pause,bool_stop);
00262 }
00263
00264
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
00282
00283
00284
00285 ast.n = atoi(number.c_str());
00286
00287 ast.name = name;
00288 remove_leading_trailing_spaces(ast.name);
00289
00290
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
00319 ast.orb.mu = GetG()*GetMSun();
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
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
00362
00363
00364
00365
00366
00367 AstorbFile::AstorbFile() : AsteroidDatabaseFile() {
00368
00369
00370 db = new AsteroidDatabase();
00371 }
00372
00373 AstorbFile::~AstorbFile() {
00374 delete db;
00375 db = 0;
00376 }
00377
00378
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
00400
00401
00402
00403 char line[300];
00404
00405 double a,e,i,omega_node,omega_pericenter,M;
00406
00407 string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
00408 string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
00409
00410
00411 string year,month,day;
00412 int y=0,m=0,d=0;
00413
00414
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
00424 Date tmp_date;
00425
00426
00427
00428
00429
00430
00431
00432
00433 while ( (GETS_FILE(line,300,file)) != 0 ) {
00434
00435
00436
00437 ++local_index;
00438
00439 if (strlen(line) < 201) continue;
00440
00441
00442
00443
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
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
00476
00477 {
00478
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
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
00515
00516
00517 year.assign(epoch,0,1);
00518 ch = (int)year.c_str()[0];
00519
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
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
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
00583
00584 tmp_date.SetGregor(y,m,d,TDT);
00585 ast.orb.epoch.SetDate(tmp_date);
00586
00587
00588 ast.orb.mu = GetG()*GetMSun();
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
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
00624
00625 db->push_back(ast);
00626
00627 }
00628
00629 read_finished();
00630 }
00631
00632
00633
00634 MPCCometFile::MPCCometFile() : AsteroidDatabaseFile() {
00635 db = new AsteroidDatabase();
00636
00637 }
00638
00639 MPCCometFile::~MPCCometFile() {
00640 delete db;
00641 db = 0;
00642 }
00643
00644
00645
00646
00647
00648
00649
00650 void MPCCometFile::Read() {
00651
00652
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
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
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;
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
00701 sleep(1);
00702 read_progress(local_index,bool_pause,bool_stop);
00703 }
00704
00705
00706 number.assign(line,0,4);
00707 type.assign(line,4,1);
00708
00709 name.assign(line,102,strlen(line)-102-1);
00710
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
00721
00722 ast.name = name;
00723 remove_leading_trailing_spaces(ast.name);
00724
00725 ast.n = 0;
00726
00727 ast.mag = atof(absolute_magnitude.c_str());
00728
00729
00730 e = atof(eccentricity.c_str());
00731
00732
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
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
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
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
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
00837
00838 MPCObsFile::MPCObsFile() { }
00839
00840 void MPCObsFile::Read() {
00841
00842
00843
00844
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
00858
00859 Observation dummy_obs;
00860
00861
00862
00863 double gradi,primi,secondi;
00864
00865
00866
00867 char line[256];
00868
00869
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;
00879
00880
00881
00882
00883
00884
00885
00886 number.assign(line,0,5);
00887
00888
00889 designation.assign(line,5,7);
00890 remove_leading_trailing_spaces(designation);
00891
00892 discovery_asterisk.assign(line,12,1);
00893
00894
00895 note1.assign(line,13,1);
00896
00897 note2.assign(line,14,1);
00898
00899 date.assign(line,15,17);
00900
00901
00902 ra.assign(line,32,12);
00903
00904 dec.assign(line,44,12);
00905
00906 magnitude.assign(line,65,6);
00907
00908 observatory_code.assign(line,77,3);
00909 remove_leading_trailing_spaces(observatory_code);
00910
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
00921
00922
00923 double y=0.0, m=0.0, d=0.0;
00924 sscanf(date.c_str(),"%lf %lf %lf",&y,&m,&d);
00925
00926
00927
00928
00929
00930 dummy_obs.date.SetGregor((int)y,(int)m,d,UTC);
00931
00932
00933
00934
00935 sscanf(ra.c_str(),"%lf %lf %lf",&gradi,&primi,&secondi);
00936
00937
00938 dummy_obs.ra.SetHMS(gradi,primi,secondi);
00939
00940
00941
00942 sscanf(dec.c_str(),"%lf %lf %lf",&gradi,&primi,&secondi);
00943
00944
00945 dummy_obs.dec.SetDPS(gradi,primi,secondi);
00946
00947
00948
00949
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
00969
00970
00971
00972
00973
00974
00975 }
00976
00977 }
00978
00979 bool MPCObsFile::ReadNominalOrbit(OrbitWithEpoch &) {
00980
00981 return false;
00982 }
00983
00984
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
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
01010 if (strlen(line) < 130) continue;
01011 if (line[0] == '!') continue;
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
01027
01028 dummy_obs.obscode = observatory_code;
01029
01030
01031
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
01045 if ((designation != "") &&
01046 (observatory_code != "") ) {
01047 obs.push_back(dummy_obs);
01048 }
01049
01050 }
01051
01052 }
01053
01054
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
01074
01075 REWIND_FILE(file);
01076
01077 string stag, spre="<pre>", spreend="</pre>";
01078 char tag[256];
01079
01080
01081 while (GETS_FILE(line,300,file) != 0) {
01082 sscanf(line,"%s",tag);
01083 stag = tag;
01084 if (stag == spre) {
01085
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
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);
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
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
01157 SWIFTFile::SWIFTFile(OrbitStream &osin) : ReadFile() {
01158 os = &osin;
01159
01160 }
01161
01162
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
01169
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
01191 file_time = FromUnits(file_time,YEAR);
01192
01193 return (good);
01194 }
01195
01196 int SWIFTFile::AsteroidsInFile() {
01197
01198
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
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
01232
01233 OrbitWithEpoch fo;
01234
01235
01236
01237 double time_old = 0, timestep;
01238
01239 int jump = 0, i_jump = 0;
01240
01241
01242
01243 ost.clear();
01244 ost.timestep = 0.0;
01245 const int asteroid_number = ost.asteroid_number;
01246
01247 char label[10];
01248 sprintf(label,"%04i",ost.asteroid_number);
01249 ost.label = label;
01250
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
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;
01307
01308 ost.push_back(fo);
01309
01310
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
01325 if (ost.size() == 2) {
01326 ost.timestep = timestep;
01327 }
01328
01329 }
01330
01331 } while (good != 0);
01332 }
01333
01334
01335
01336 OrsaConfigFile::OrsaConfigFile() : ReadWriteFile() {
01337
01338
01339
01340
01341
01342 char path[1024], command[1024];
01343
01344
01345
01346
01347
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
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
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
01405
01406
01407 Open(OPEN_R);
01408
01409
01410
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
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
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
01460 Close();
01461
01462
01463
01464 Open(OPEN_W);
01465
01466 if (status != OPEN_W) {
01467 ORSA_ERROR("Status error!");
01468 return;
01469 }
01470
01471
01472
01473
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
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
01572 byte_order = BYTEORDER;
01573 Write(&byte_order);
01574
01575
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);
01642
01643
01644
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
01661 Write((*e)->GetIntegrator());
01662
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
01680 if ((*e)->size() > 0) Write(&((*(*e))[0]));
01681
01682
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
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
01742 Frame f;
01743
01744 Read(&last_ofdt_read);
01745
01746
01747
01748
01749
01750 if (last_ofdt_read == OFDT_FRAME) {
01751
01752 Read(&f);
01753 (*e)->push_back(f);
01754 }
01755
01756 Read(&last_ofdt_read);
01757
01758
01759
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
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
01907 void OrsaFile::Write(const Integrator * i) {
01908
01909
01910 IntegratorType it = i->GetType();
01911 Write(&it);
01912
01913
01914 UniverseTypeAwareTimeStep ts = i->timestep;
01915 Write(&ts);
01916
01917
01918 double a = i->accuracy;
01919 Write(&a);
01920
01921
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
01944
01945
01946
01947
01948
01949
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
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
02008
02009
02010 free(name);
02011 {
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
02027
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
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092 void OrsaFile::Write(IntegratorType *it) {
02093
02094 unsigned int i = *it;
02095 Write(&i);
02096 }
02097
02098 void OrsaFile::Read(IntegratorType *it) {
02099
02100 unsigned int i;
02101 Read(&i);
02102 convert(*it,i);
02103 }
02104
02105 void OrsaFile::Write(InteractionType *it) {
02106
02107 unsigned int i = *it;
02108 Write(&i);
02109 }
02110
02111 void OrsaFile::Read(InteractionType *it) {
02112
02113 unsigned int i;
02114 Read(&i);
02115 convert(*it,i);
02116 }
02117
02118 void OrsaFile::Write(time_unit *tu) {
02119
02120 unsigned int i = *tu;
02121 Write(&i);
02122 }
02123
02124 void OrsaFile::Read(time_unit *tu) {
02125
02126 unsigned int i;
02127 Read(&i);
02128 convert(*tu,i);
02129 }
02130
02131 void OrsaFile::Write(length_unit *lu) {
02132
02133 unsigned int i = *lu;
02134 Write(&i);
02135 }
02136
02137 void OrsaFile::Read(length_unit *lu) {
02138
02139 unsigned int i;
02140 Read(&i);
02141 convert(*lu,i);
02142 }
02143
02144 void OrsaFile::Write(mass_unit *mu) {
02145
02146 unsigned int i = *mu;
02147 Write(&i);
02148 }
02149
02150 void OrsaFile::Read(mass_unit *mu) {
02151
02152 unsigned int i;
02153 Read(&i);
02154 convert(*mu,i);
02155 }
02156
02157 void OrsaFile::Write(ReferenceSystem *rs) {
02158
02159 unsigned int i = *rs;
02160 Write(&i);
02161 }
02162
02163 void OrsaFile::Read(ReferenceSystem *rs) {
02164
02165 unsigned int i;
02166 Read(&i);
02167 convert(*rs,i);
02168 }
02169
02170 void OrsaFile::Write(UniverseType *ut) {
02171
02172 unsigned int i = *ut;
02173 Write(&i);
02174 }
02175
02176 void OrsaFile::Read(UniverseType *ut) {
02177
02178 unsigned int i;
02179 Read(&i);
02180 convert(*ut,i);
02181 }
02182
02183 void OrsaFile::Write(TimeScale *ts) {
02184
02185 unsigned int i = *ts;
02186 Write(&i);
02187 }
02188
02189 void OrsaFile::Read(TimeScale *ts) {
02190
02191 unsigned int i;
02192 Read(&i);
02193 convert(*ts,i);
02194 }
02195
02196 void OrsaFile::Write(OrsaFileDataType *ofdt) {
02197
02198 unsigned int i = *ofdt;
02199 Write(&i);
02200 }
02201
02202
02203
02204
02205
02206
02207
02208
02209 void OrsaFile::Read(OrsaFileDataType *ofdt) {
02210
02211
02212 unsigned int i;
02213 const int val = read_swap(&i,sizeof(unsigned int));
02214
02215
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
02225 unsigned int i = *jp;
02226 Write(&i);
02227 }
02228
02229 void OrsaFile::Read(JPL_planets *jp) {
02230
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
02258 bool OrsaFile::GoodFile(const string &filename) {
02259
02260
02261
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
02288 RadauModIntegrationFile::RadauModIntegrationFile(OrbitStream &osin) : ReadFile() {
02289 os = &osin;
02290 }
02291
02292 void RadauModIntegrationFile::Read() {
02293
02294
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
02314
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);
02326
02327 }
02328
02329 fo.epoch = FromUnits(time,YEAR);
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
02341 if (fo.e >= 1.0) {
02342 ORSA_ERROR("reading eccentricity > 1.0, returning.");
02343 return;
02344 }
02345
02346 }
02347
02348 }
02349
02350
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
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;
02397
02398 if (line[0] == ' ') continue;
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;
02449
02450 if (line[0] != ' ') break;
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
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
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
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
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
02593 string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
02594 string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
02595
02596
02597 string year,month,day;
02598
02599
02600 Asteroid ast;
02601
02602
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;
02612
02613 if (line[0]=='-') continue;
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
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
02652 if ((ast.n==0) || (a==0.0)) {
02653
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
02665
02666
02667
02668
02669
02670
02671
02672 tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
02673 ast.orb.epoch.SetDate(tmp_date);
02674
02675 ast.orb.mu = GetG()*GetMSun();
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
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
02741 string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
02742 string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
02743
02744
02745 string year,month,day;
02746
02747
02748 Asteroid ast;
02749
02750
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;
02760
02761 if (line[0]=='-') continue;
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
02774
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
02800 if ((a==0.0)) {
02801
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
02813
02814
02815
02816
02817
02818
02819
02820 tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
02821 ast.orb.epoch.SetDate(tmp_date);
02822
02823 ast.orb.mu = GetG()*GetMSun();
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837
02838
02839
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
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
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
02912
02913 if (line[0]=='-') continue;
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
02922 sleep(1);
02923 read_progress(local_index,bool_pause,bool_stop);
02924 }
02925
02926
02927
02928
02929
02930 name.assign(line,0,37);
02931 epoch.assign(line,38,5);
02932
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
02942
02943 ast.name = name;
02944 remove_leading_trailing_spaces(ast.name);
02945
02946
02947 ast.n = 0;
02948
02949
02950
02951
02952 e = atof(eccentricity.c_str());
02953
02954
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
02963 if ((q==0.0)) {
02964
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
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
03002
03003
03004
03005
03006
03007
03008
03009
03010
03011
03012
03013
03014
03015
03016
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
03043
03044
03045 NEODYSCAT::NEODYSCAT() : AsteroidDatabaseFile() {
03046
03047 db = new AsteroidDatabase();
03048 }
03049
03050 NEODYSCAT::~NEODYSCAT() {
03051 delete db;
03052 db = 0;
03053 }
03054
03055 void NEODYSCAT::Read() {
03056
03057
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
03072 string number,name,orbit_computer,absolute_magnitude,arc,numobs,epoch;
03073 string mean_anomaly,pericenter,node,inclination,eccentricity,semimajor_axis;
03074
03075
03076 string year,month,day;
03077
03078
03079 Asteroid ast;
03080
03081
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;
03091
03092 if (line[0]=='!') continue;
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
03105
03106 name.assign(line,0,12);
03107
03108 {
03109
03110 string::size_type pos;
03111 while ((pos=name.find("'",0)) != string::npos) {
03112
03113 name.erase(pos,1);
03114 }
03115
03116 }
03117
03118
03119
03120
03121
03122
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
03135
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
03159
03160
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
03177
03178
03179
03180
03181
03182
03183
03184 tmp_date.SetJulian(2400000.5+atof(epoch.c_str()),TDT);
03185 ast.orb.epoch.SetDate(tmp_date);
03186
03187 ast.orb.mu = GetG()*GetMSun();
03188
03189
03190
03191
03192
03193
03194
03195
03196
03197
03198
03199
03200
03201
03202
03203
03204
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
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
03284
03285
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;
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;
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));
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);
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
03400
03401 have_one = have_two = false;
03402 }
03403 }
03404 }
03405
03406 }