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