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