1 // Write the coefficient files needed for approximating the normal gravity
2 // field with a GravityModel.  WARNING: this creates files, wgs84.egm.cof and
3 // grs80.egm.cof, in the current directory.
4 
5 #include <cmath>
6 #include <fstream>
7 #include <iostream>
8 #include <GeographicLib/NormalGravity.hpp>
9 #include <GeographicLib/Utility.hpp>
10 
11 using namespace std;
12 using namespace GeographicLib;
13 
main()14 int main() {
15   try {
16     Utility::set_digits();
17     const char* filenames[] = {"wgs84.egm.cof", "grs80.egm.cof"};
18     const char* ids[] = {"WGS1984A", "GRS1980A"};
19     for (int grs80 = 0; grs80 < 2; ++grs80) {
20       ofstream file(filenames[grs80], ios::binary);
21       Utility::writearray<char, char, false>(file, ids[grs80], 8);
22       const int N = 20, M = 0,
23         cnum = (M + 1) * (2 * N - M + 2) / 2; // cnum = N + 1
24       vector<int> num(2);
25       num[0] = N; num[1] = M;
26       Utility::writearray<int, int, false>(file, num);
27       vector<Math::real> c(cnum, 0);
28       const NormalGravity& earth(grs80 ? NormalGravity::GRS80() :
29                                  NormalGravity::WGS84());
30       for (int n = 2; n <= N; n += 2)
31         c[n] = - earth.DynamicalFormFactor(n) / sqrt(Math::real(2*n + 1));
32       Utility::writearray<double, Math::real, false>(file, c);
33       num[0] = num[1] = -1;
34       Utility::writearray<int, int, false>(file, num);
35     }
36   }
37   catch (const exception& e) {
38     cerr << "Caught exception: " << e.what() << "\n";
39     return 1;
40   }
41   catch (...) {
42     cerr << "Caught unknown exception\n";
43     return 1;
44   }
45 }
46