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
123
124
125
126
127
128
129
130
|
/***********************************************/
/**
* @file orbit2Kepler.cpp
*
* @brief Keplerian elements from orbit position and velocity at each epoch.
*
* @author Torsten Mayer-Guerr
* @date 2010-07-20
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
This program computes the osculating Keplerian elements from position and velocity
of a given \configFile{inputfileOrbit}{instrument}.
The \configFile{inputfileOrbit}{instrument} must contain positions and velocities (see \program{OrbitAddVelocityAndAcceleration}).
The \config{outputfileKepler} is an \file{instrument file}{instrument} (MISCVALUES)
with the Keplerian elements at each epoch in the following order
\begin{itemize}
\item Ascending Node $\Omega$ [degree]
\item Inclination $i$ [degree]
\item Argument of perigee $\omega$ [degree]
\item major axis $a$ [m]
\item eccentricity $e$
\item mean anomaly $M$ [degree]
\item transit time of perigee $\tau$ [mjd]
\end{itemize}
The data of \configFile{inputfileInstrument}{instrument} are appended as values to each epoch.
)";
/***********************************************/
#include "programs/program.h"
#include "base/kepler.h"
#include "files/fileInstrument.h"
/***** CLASS ***********************************/
/** @brief Keplerian elements from orbit position and velocity at each epoch.
* @ingroup programsGroup */
class Orbit2Kepler
{
public:
void run(Config &config, Parallel::CommunicatorPtr comm);
};
GROOPS_REGISTER_PROGRAM(Orbit2Kepler, PARALLEL, "keplerian elements from orbit position and velocity at each epoch", Orbit, Instrument, TimeSeries)
/***********************************************/
void Orbit2Kepler::run(Config &config, Parallel::CommunicatorPtr comm)
{
try
{
FileName fileNameOut, fileNameOrbit;
std::vector<FileName> fileNamesInstrument;
Double GM;
readConfig(config, "outputfileKepler", fileNameOut, Config::MUSTSET, "", "instrument file (MISCVALUES: Omega, i, omega [degree], a [m], e, M [degree], tau [mjd], ...)");
readConfig(config, "inputfileOrbit", fileNameOrbit, Config::MUSTSET, "", "position and velocity at each epoch define the kepler orbit");
readConfig(config, "inputfileInstrument", fileNamesInstrument, Config::OPTIONAL, "", "data is appended");
readConfig(config, "GM", GM, Config::DEFAULT, STRING_DEFAULT_GM, "Geocentric gravitational constant");
if(isCreateSchema(config)) return;
// =======================
logStatus<<"calculate Keplerian elements"<<Log::endl;
InstrumentFile orbitFile(fileNameOrbit);
UInt dataCount = 8;
std::vector<InstrumentFilePtr> instrumentFile;
for(auto &fileName : fileNamesInstrument)
{
instrumentFile.push_back(InstrumentFile::newFile(fileName));
InstrumentFile::checkArcCount({orbitFile, *instrumentFile.back()});
dataCount += instrumentFile.back()->dataCount(TRUE/*mustDefined*/);
}
std::vector<Arc> arcList(orbitFile.arcCount());
Parallel::forEach(arcList, [&] (UInt arcNo)
{
const OrbitArc orbit = orbitFile.readArc(arcNo);
Matrix A(orbit.size(), dataCount);
for(UInt i=0; i<orbit.size(); i++)
{
if(orbit.at(i).velocity.r()==0)
throw(Exception("no velocity given"));
const Kepler kepler(orbit.at(i).time, orbit.at(i).position, orbit.at(i).velocity, GM);
const Double M = (orbit.at(i).time-kepler.tau).seconds() * std::sqrt(GM/std::pow(kepler.a,3));
A(i, 1) = std::fmod(kepler.Omega+2*PI, 2*PI) * RAD2DEG;
A(i, 2) = kepler.i * RAD2DEG;
A(i, 3) = std::fmod(kepler.omega+2*PI, 2*PI) * RAD2DEG;
A(i, 4) = kepler.a;
A(i, 5) = kepler.e;
A(i, 6) = std::fmod(M+2*PI, 2*PI) * RAD2DEG;
A(i, 7) = kepler.tau.mjd();
}
UInt idx = 8;
for(auto &file: instrumentFile)
{
Arc arc = file->readArc(arcNo);
Arc::checkSynchronized({orbit, arc});
Matrix B = arc.matrix();
copy(B.column(1, B.columns()-1), A.column(idx, B.columns()-1));
idx += B.columns()-1;
}
return Arc(orbit.times(), A);
}, comm);
// write results
// -------------
if(Parallel::isMaster(comm))
{
logStatus<<"write time series of Keplerian elements to file <"<<fileNameOut<<">"<<Log::endl;
InstrumentFile::write(fileNameOut, arcList);
Arc::printStatistics(arcList);
}
}
catch(std::exception &e)
{
GROOPS_RETHROW(e)
}
}
/***********************************************/
|
