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
131
132
|
/***********************************************/
/**
* @file gravityfieldReplacePotentialCoefficients.cpp
*
* @brief Replace single potential coefficients in a gravity field
*
* @author Andreas Kvas
* @author Torsten Mayer-Guerr
* @date 2015-11-17
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
Replaces single potential coefficients in a gravity field.
Both \configClass{gravityfield}{gravityfieldType}
and \configClass{gravityfieldReplacement}{gravityfieldType} are evaluated
at \config{time} and converted to spherical harmonic coefficients.
Single \config{coefficients} are then replaced in \configClass{gravityfield}{gravityfieldType}
by the values from \configClass{gravityfieldReplacement}{gravityfieldType}
and the result is written to \configFile{outputfilePotentialCoefficients}{potentialCoefficients}
from \config{minDegree} to \config{maxDegree},
)";
/***********************************************/
#include "programs/program.h"
#include "files/fileSphericalHarmonics.h"
#include "classes/gravityfield/gravityfield.h"
/***** CLASS ***********************************/
/** @brief Replace single potential coefficients in a gravity field.
* @ingroup programsGroup */
class GravityfieldReplacePotentialCoefficients
{
public:
class Coeff
{
public:
Bool isCnm;
UInt degree, order;
};
void run(Config &config, Parallel::CommunicatorPtr comm);
};
GROOPS_REGISTER_PROGRAM(GravityfieldReplacePotentialCoefficients, SINGLEPROCESS, "Replace single potential coefficients in a gravity field", Gravityfield)
GROOPS_RENAMED_PROGRAM(GravityfieldReplaceC20, GravityfieldReplacePotentialCoefficients, date2time(2020, 5, 24))
/***********************************************/
template<> Bool readConfig(Config &config, const std::string &name, GravityfieldReplacePotentialCoefficients::Coeff &var, Config::Appearance mustSet, const std::string &defaultValue, const std::string &annotation)
{
std::string choice;
if(!readConfigChoice(config, name, choice, mustSet, defaultValue, annotation))
return FALSE;
if(readConfigChoiceElement(config, "cnm", choice, ""))
{
var.isCnm = TRUE;
readConfig(config, "degree", var.degree, Config::MUSTSET, "", "");
readConfig(config, "order", var.order, Config::MUSTSET, "", "");
}
if(readConfigChoiceElement(config, "snm", choice, ""))
{
var.isCnm = FALSE;
readConfig(config, "degree", var.degree, Config::MUSTSET, "", "");
readConfig(config, "order", var.order, Config::MUSTSET, "", "");
}
endChoice(config);
return TRUE;
}
/***********************************************/
void GravityfieldReplacePotentialCoefficients::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
try
{
FileName fileNameCoeff;
GravityfieldPtr gravityfield1, gravityfield2;
std::vector<Coeff> coefficients;
UInt minDegree, maxDegree = INFINITYDEGREE;
Double GM, R;
Time time;
readConfig(config, "outputfilePotentialCoefficients", fileNameCoeff, Config::MUSTSET, "", "");
readConfig(config, "gravityfield", gravityfield1, Config::MUSTSET, "", "single coefficients are replaced by the other gravityfield");
readConfig(config, "gravityfieldReplacement", gravityfield2, Config::MUSTSET, "", "contains the coefficients for replacement");
readConfig(config, "coefficients", coefficients, Config::MUSTSET, "", "");
readConfig(config, "minDegree", minDegree, Config::DEFAULT, "0", "");
readConfig(config, "maxDegree", maxDegree, Config::OPTIONAL, "", "");
readConfig(config, "GM", GM, Config::DEFAULT, STRING_DEFAULT_GM, "Geocentric gravitational constant");
readConfig(config, "R", R, Config::DEFAULT, STRING_DEFAULT_R, "reference radius");
readConfig(config, "time", time, Config::OPTIONAL, "", "at this time the gravity field will be evaluated");
if(isCreateSchema(config)) return;
logStatus<<"create spherical harmonics"<<Log::endl;
SphericalHarmonics harm1 = gravityfield1->sphericalHarmonics(time, maxDegree, minDegree, GM, R);
SphericalHarmonics harm2 = gravityfield2->sphericalHarmonics(time, maxDegree, minDegree, GM, R);
for(auto &coeff : coefficients)
{
if(coeff.isCnm)
{
harm1.cnm()(coeff.degree, coeff.order) = harm2.cnm()(coeff.degree, coeff.order);
if(harm1.sigma2cnm().size())
harm1.sigma2cnm()(coeff.degree, coeff.order) = 0;
if(harm1.sigma2cnm().size() && harm2.sigma2cnm().size())
harm1.sigma2cnm()(coeff.degree, coeff.order) = harm2.sigma2cnm()(coeff.degree, coeff.order);
}
else
{
harm1.snm()(coeff.degree, coeff.order) = harm2.snm()(coeff.degree, coeff.order);
if(harm1.sigma2snm().size())
harm1.sigma2snm()(coeff.degree, coeff.order) = 0;
if(harm1.sigma2snm().size() && harm2.sigma2snm().size())
harm1.sigma2snm()(coeff.degree, coeff.order) = harm2.sigma2snm()(coeff.degree, coeff.order);
}
}
logStatus<<"writing potential coefficients to file <"<<fileNameCoeff<<">"<<Log::endl;
writeFileSphericalHarmonics(fileNameCoeff, harm1);
}
catch(std::exception &e)
{
GROOPS_RETHROW(e)
}
}
/***********************************************/
|
