1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
/***********************************************/
/**
* @file simulateKeplerOrbit.cpp
*
* @brief Compute Keplerian orbit.
*
* @author Torsten Mayer-Guerr
* @date 2009-10-31
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
This program simulates a Keplerian \file{orbit}{instrument} at a given \config{timeSeries}
starting from the given \config{integrationConstants}.
)";
/***********************************************/
#include "programs/program.h"
#include "base/kepler.h"
#include "base/equinoctial.h"
#include "files/fileInstrument.h"
#include "classes/timeSeries/timeSeries.h"
/***** CLASS ***********************************/
/** @brief Compute Keplerian orbit.
* @ingroup programsGroup */
class SimulateKeplerOrbit
{
public:
void run(Config &config, Parallel::CommunicatorPtr comm);
};
GROOPS_REGISTER_PROGRAM(SimulateKeplerOrbit, SINGLEPROCESS, "compute Keplerian orbit", Simulation, Orbit, Instrument)
/***********************************************/
void SimulateKeplerOrbit::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
try
{
FileName outName;
Double GM;
TimeSeriesPtr timeSeries;
Vector3d position, velocity;
std::string choice;
Time time0;
readConfig(config, "outputfileOrbit", outName, Config::MUSTSET, "", "");
readConfig(config, "timeSeries", timeSeries, Config::MUSTSET, "", "");
readConfig(config, "GM", GM, Config::DEFAULT, STRING_DEFAULT_GM, "Geocentric gravitational constant");
readConfigChoice(config, "integrationConstants", choice, Config::MUSTSET, "", "");
if(readConfigChoiceElement(config, "kepler", choice, ""))
{
Double a, e;
Angle i, Omega, omega, M;
readConfig(config, "majorAxis", a, Config::MUSTSET, "", "[m]");
readConfig(config, "eccentricity", e, Config::MUSTSET, "", "[-]");
readConfig(config, "inclination", i, Config::MUSTSET, "", "[degree]");
readConfig(config, "ascendingNode", Omega, Config::MUSTSET, "", "[degree]");
readConfig(config, "argumentOfPerigee", omega, Config::MUSTSET, "", "[degree]");
readConfig(config, "meanAnomaly", M, Config::MUSTSET, "", "[degree]");
readConfig(config, "time", time0, Config::MUSTSET, "", "integration constants are valid at this epoch");
if(!isCreateSchema(config))
{
Kepler kepler(time0, Omega, i, omega, a, e, M, GM);
kepler.orbit(time0, position, velocity);
}
}
if(readConfigChoiceElement(config, "positionAndVelocity", choice, ""))
{
readConfig(config, "position0x", position.x(), Config::MUSTSET, "", "[m] in CRF");
readConfig(config, "position0y", position.y(), Config::MUSTSET, "", "[m] in CRF");
readConfig(config, "position0z", position.z(), Config::MUSTSET, "", "[m] in CRF");
readConfig(config, "velocity0x", velocity.x(), Config::MUSTSET, "", "[m/s]");
readConfig(config, "velocity0y", velocity.y(), Config::MUSTSET, "", "[m/s]");
readConfig(config, "velocity0z", velocity.z(), Config::MUSTSET, "", "[m/s]");
readConfig(config, "time", time0, Config::MUSTSET, "", "integration constants are valid at this epoch");
}
endChoice(config);
if(isCreateSchema(config)) return;
std::vector<Time> times = timeSeries->times();
Equinoctial kepler(time0, position, velocity, GM);
// Info
// ----
Kepler k(kepler);
logInfo<<"Kepler elements"<<Log::endl;
logInfo<<" majorAxis a : "<<k.a/1000<<" km"<<Log::endl;
logInfo<<" eccentricity e : "<<k.e<<Log::endl;
logInfo<<" inclination i : "<<k.i*RAD2DEG<<" Degree"<<Log::endl;
logInfo<<" ascending node : "<<k.Omega*RAD2DEG<<" Degree"<<Log::endl;
logInfo<<" perigee : "<<k.omega*RAD2DEG<<" Degree"<<Log::endl;
// Computation
// -----------
logStatus<<"computing"<<Log::endl;
OrbitArc orbit;
Single::forEach(times.size(), [&](UInt i)
{
OrbitEpoch epoch;
epoch.time = times.at(i);
kepler.orbit(times.at(i), epoch.position, epoch.velocity, epoch.acceleration);
orbit.push_back(epoch);
});
// write results
// -------------
logStatus<<"write orbit data to file <"<<outName<<">"<<Log::endl;
InstrumentFile::write(outName, orbit);
}
catch(std::exception &e)
{
GROOPS_RETHROW(e)
}
}
/***********************************************/
|
