#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#if defined(_OPENMP)
#define HAVE_OPENMP 1
#else
#define HAVE_OPENMP 0
#endif
#if HAVE_OPENMP
# include <omp.h>
#endif
#include <GeographicLib/GravityModel.hpp>
#include <GeographicLib/GravityCircle.hpp>
#include <GeographicLib/Utility.hpp>
using namespace std;
using namespace GeographicLib;
int main(int argc, char* argv[]) {
if (argc != 4) {
std::cerr << "Usage: " << argv[0] << " model ninterval output\n";
return 1;
}
try {
typedef Math::real real;
using std::pow;
enum format {
PGM = 0,
PGM4 = 1,
GTX = 2,
};
format mode = PGM;
int ndigits = Utility::set_digits();
std::string model(argv[1]);
int ndeg = Utility::val<int>(std::string(argv[2]));
std::string filename(argv[3]);
GravityModel g(model);
int
nlat = 180 * ndeg + 1,
nlon = 360 * ndeg;
real
offset = mode == PGM ? -108 : -128,
scale = mode == PGM ? 3/real(1000) : pow(1/real(2), 24);
ofstream file(filename.c_str(), ios::binary);
switch (mode) {
case PGM:
file << "P5\n"
<< "# Geoid file in PGM format for the GeographicLib::Geoid class\n"
<< "# Description WGS84 " << model
<< ", " << 60/real(ndeg) << "-minute grid\n"
<< "# Offset " << offset << "\n"
<< "# Scale " << scale << "\n"
<< "# Origin 90N 0E\n"
<< "# AREA_OR_POINT Point\n"
<< "# Vertical_Datum WGS84\n"
<< nlon << " " << nlat << "\n" << (1<<16)-1 << "\n";
break;
case PGM4:
file << std::setprecision(17)
<< "P5\n"
<< "# Geoid file in PGM format for the GeographicLib::Geoid class\n"
<< "# Description WGS84 " << model
<< ", " << 60/real(ndeg) << "-minute grid\n"
<< "# Offset " << offset << "\n"
<< "# Scale " << scale << "\n"
<< "# Origin 90N 0E\n"
<< "# AREA_OR_POINT Point\n"
<< "# Vertical_Datum WGS84\n"
<< nlon << " " << nlat << "\n" << 0xffffffffu << "\n";
break;
case GTX:
break;
}
const int nbatch = 64;
vector< vector<float> > Nf(mode == GTX ? nbatch : 0, vector<float>(nlon));
vector< vector<unsigned short> >
Ns(mode == PGM ? nbatch : 0, vector<unsigned short>(nlon));
vector< vector<unsigned> >
Nu(mode == PGM4 ? nbatch : 0, vector<unsigned>(nlon));
for (int ilat0 = 0; ilat0 < nlat; ilat0 += nbatch) { int nlat0 = min(nlat, ilat0 + nbatch);
#if HAVE_OPENMP
# pragma omp parallel for
#endif
for (int ilat = ilat0; ilat < nlat0; ++ilat) { Utility::set_digits(ndigits); real lat = (90 * ndeg - ilat) / real(ndeg), h = 0, half = 1/real(2);
GravityCircle c(g.Circle(lat, h, GravityModel::GEOID_HEIGHT));
for (int ilon = 0; ilon < nlon; ++ilon) { real lon = (2*ilon < nlon ? ilon : ilon - nlon) / real(ndeg),
gh = c.GeoidHeight(lon);
switch (mode) {
case PGM:
Ns[ilat - ilat0][ilon] = (unsigned short)
(std::max(0, std::min((1<<16)-1,
int((gh - offset)/scale + half))));
break;
case PGM4:
Nu[ilat - ilat0][ilon] = unsigned((gh - offset)/scale + half);
break;
case GTX:
Nf[ilat - ilat0][ilon] = float(gh);
break;
}
} } for (int ilat = ilat0; ilat < nlat0; ++ilat)
switch (mode) {
case PGM:
Utility::writearray<short unsigned, short unsigned, true>
(file, Ns[ilat - ilat0]);
break;
case PGM4:
Utility::writearray<unsigned, unsigned, true>
(file, Nu[ilat - ilat0]);
break;
case GTX:
file.write(reinterpret_cast<char *>(Nf[ilat - ilat0].data()),
nlon * sizeof(float));
break;
}
} }
catch (const std::exception& e) {
std::cerr << "Caught exception: " << e.what() << "\n";
return 1;
}
catch (...) {
std::cerr << "Caught unknown exception\n";
return 1;
}
return 0;
}