#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <stdexcept>
#include <algorithm>
#include <cmath>

int cosind(int N, int M, int n, int m) {
 return m * N - m * (m - 1) / 2 + n;
}

int sinind(int N, int M, int n, int m) {
 return cosind(N, M, n, m) - (N + 1);
}

void storecos(std::vector<double>& C, int N, int M,
 double v, int n, int m) {
 if (v == 0)
 return;
 if (n < 0 || n > N || m < 0 || m > M)
 throw std::runtime_error("Invalid coefficient");
 C[cosind(N, M, n, m)] = v;
}

void storesin(std::vector<double>& S, int N, int M,
 double v, int n, int m) {
 if (v == 0)
 return;
 if (n < 0 || n > N || m <= 0 || m > M)
 throw std::runtime_error("Invalid coefficient");
 S[sinind(N, M, n, m)] = v;
}

void storecoeff(std::ostream& str,
 const std::vector<double>& C, const std::vector<double>& S,
 int M, int N) {
 int K = (M + 1) * (2*N - M + 2) / 2;
 int ind[2] = {N, M};
 str.write(reinterpret_cast<const char*>(ind), 2 * sizeof(int));
 str.write(reinterpret_cast<const char*>(&C[0]), K * sizeof(double));
 str.write(reinterpret_cast<const char*>(&S[0]), (K-N-1) * sizeof(double));
}

int main(int argc, char* argv[]) {
 if (argc != 3) {
 std::cerr << "Usage: reformat model outfile\n";
 return 1;
 }
 std::string model = argv[1];
 std::string file = argv[2];
 try {
 // model = wmm2010
 //http://ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2010_Sph_Windows_Linux.zip
 // unpack
 // coefficients are in EMM_Sph_Windows_Linux/WMM2010.COF
 // N=739, M=718
 std::ofstream fout(file.c_str(), std::ios::binary);
 if (model == "emm2010") {
 // download
 // http://ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2010_Sph_Windows_Linux.zip
 // unpack
 // coefficients are in EMM_Sph_Windows_Linux/EMM2010.COF
 // and EMM_Sph_Windows_Linux/EMM2010SV.COF
 std::string id = "EMM2010A";
 fout.write(id.c_str(), 8);
 for (int i = 0; i < 2; ++i) {
 std::string filename = i == 0 ? "EMM_Sph_Windows_Linux/EMM2010.COF" :
 "EMM_Sph_Windows_Linux/EMM2010SV.COF";
 int
 N = i == 0 ? 739 : 16,
 M = i == 0 ? 718 : 16,
 K = (M + 1) * (2*N - M + 2) / 2;
 std::vector<double> C(K, 0.0);
 std::vector<double> S(K - (N + 1), 0.0);
 std::ifstream fin(filename.c_str());
 std::string ss;
 if (i == 0) std::getline(fin, ss); // Skip first line
 while (std::getline(fin, ss)) {
 int n, m;
 double c, s;
 std::istringstream is(ss);
 if (!(is >> n >> m >> c >> s))
 throw std::runtime_error("Short read");
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 }
 storecoeff(fout, C, S, M, N);
 }
 } else if (model == "emm2015") {
 // download
 //http://www.ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2015_Sph_Linux.zip
 // unpack
 //
 // * The only coefficients needed are EMM20{00-15}.COF and EMM2015SV.COF.
 //
 // * The other SV files can be ignored because the time dependence can be
 // obtained using linear interpolation for dates < 2015 (or
 // extrapolation for dates < 2000).
 //
 // * The time varying part of the field is of degree 15. This
 // constitutes all the coefficients in EMM20{00-14}.COF and
 // EMM2015SV.COF and a subset of the coefficients in EMM2015.COF.
 //
 // * To this should be added a time independent short wavelength field
 // which is given by the degree > 15 terms in EMM2015.COF.
 //
 // * These time independent terms are of degree 729 and order 718. There
 // are higher degree and order coefficients listed in the file, but
 // these are zero.
 //
 // * The EMM2015 coefficients as used by GeographicLib compress much
 // better than those for EMM2010 (660 kB instead of 3700 kB).
 // Presumably this is because the EMM2015 are only given with 4 decimal
 // digits.
 //
 // * The GeographicLib implementation produces the same results as listed
 // in EMM2015_TEST_VALUES.txt

 std::string id = "EMM2015A";
 fout.write(id.c_str(), 8);
 for (int i = 0; i <= 17; ++i) {
 std::string filename;
 {
 std::ostringstream os;
 os << "EMM2015_linux/EMM" << (2000 + std::min(15, i))
 << (i == 16 ? "SV" : "") << ".COF";
 filename = os.str();
 }
 int
 N = i == 17 ? 729 : 15,
 M = i == 17 ? 718 : 15,
 K = (M + 1) * (2*N - M + 2) / 2;
 std::vector<double> C(K, 0.0);
 std::vector<double> S(K - (N + 1), 0.0);
 std::ifstream fin(filename.c_str());
 std::string ss;
 if (i != 16) std::getline(fin, ss); // Skip first line
 while (std::getline(fin, ss)) {
 int n, m;
 double c, s;
 std::istringstream is(ss);
 if (!(is >> n >> m >> c >> s))
 throw std::runtime_error("Short read " + filename + ": " + ss);
 if (i == 15 && n > 15)
 continue;
 if (i == 17 && n <= 15)
 continue;
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 }
 storecoeff(fout, C, S, M, N);
 }
 } else if (model == "emm2017") {
 // download
 //https://www.ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2017_Sph_Linux.zip
 // unpack
 //
 // * The only coefficients needed are EMM20{00-17}.COF and EMM2017SV.COF.
 //
 // * The other SV files can be ignored because the time dependence can be
 // obtained using linear interpolation for dates < 2017 (or
 // extrapolation for dates < 2000).
 //
 // * The time varying part of the field is of degree 15. This
 // constitutes all the coefficients in EMM20{00-16}.COF and
 // EMM2017SV.COF and a subset of the coefficients in EMM2017.COF.
 //
 // * To this should be added a time independent short wavelength field
 // which is given by the degree > 15 terms in EMM2017.COF.
 //
 // * These time independent terms are of degree 790 and order 790.
 //
 // * The GeographicLib implementation produces the same results as listed
 // in EMM2017TestValues.txt

 std::string id = "EMM2017A";
 fout.write(id.c_str(), 8);
 int maxy = 17;
 for (int i = 0; i <= 19; ++i) {
 // i = maxy for low res components for 2000+maxy
 // i = maxy+1 for SV at 2000+maxy
 // i = maxy+2 for high res components for 2000+maxy
 std::string filename;
 {
 std::ostringstream os;
 os << "EMM2017_Linux/EMM" << (2000 + std::min(maxy, i))
 << (i == maxy+1 ? "SV" : "") << ".COF";
 filename = os.str();
 }
 int
 N = i == maxy+2 ? 790 : 15,
 M = i == maxy+2 ? 790 : 15,
 K = (M + 1) * (2*N - M + 2) / 2;
 std::vector<double> C(K, 0.0);
 std::vector<double> S(K - (N + 1), 0.0);
 std::ifstream fin(filename.c_str());
 std::string ss;
 if (i != maxy+1) std::getline(fin, ss); // Skip first line
 while (std::getline(fin, ss)) {
 int n, m;
 double c, s;
 std::istringstream is(ss);
 if (!(is >> n >> m >> c >> s))
 throw std::runtime_error("Short read " + filename + ": " + ss);
 if (i == maxy && n > 15)
 continue;
 if (i == maxy+2 && n <= 15)
 continue;
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 }
 storecoeff(fout, C, S, M, N);
 }
 } else if (model == "wmm2010" || model == "igrf11") {
 // Download
 // http://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
 // unpack
 // wmm2010 coefficients are in geomag70_linux/WMM2010.COF
 // igrf11 coefficients are in geomag70_linux/IGRF11.COF
 std::string id = model == "wmm2010" ? "WMM2010A" : "IGRF11-A";
 fout.write(id.c_str(), 8);
 std::string filename = model == "wmm2010" ? "geomag70_linux/WMM2010.COF"
 : "geomag70_linux-2010/IGRF11.COF";
 std::ifstream fin(filename.c_str());
 std::string ss;
 bool start = true;
 std::getline(fin, ss);
 std::vector<double> C;
 std::vector<double> S;
 std::vector<double> C1;
 std::vector<double> S1;
 int N, M, N1, M1;
 while (true) {
 if (ss.size() == 0)
 break;
 std::istringstream is(ss);
 int n, m;
 double c, s, c1, s1;
 if (start) {
 std::string mm;
 if (!(is >> mm >> mm >> N >> N1))
 throw std::runtime_error("Short read on header");
 M = N; M1 = N1;
 int
 K = (M + 1) * (2*N - M + 2) / 2,
 K1 = (M1 + 1) * (2*N1 - M1 + 2) / 2;
 C.resize(K);
 S.resize(K - (N + 1));
 C1.resize(K1);
 S1.resize(K1 - (N1 + 1));
 std::fill(C.begin(), C.end(), 0.0);
 std::fill(S.begin(), S.end(), 0.0);
 std::fill(C1.begin(), C1.end(), 0.0);
 std::fill(S1.begin(), S1.end(), 0.0);
 start = false;
 } else {
 if (!(is >> n >> m >> c >> s >> c1 >> s1))
 throw std::runtime_error("Short read on data");
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 storecos(C1, N1, M1, c1, n, m);
 storesin(S1, N1, M1, s1, n, m);
 }
 if (!std::getline(fin, ss))
 ss = "";
 if (ss.size() == 0 || (ss.size() >= 3 && ss.substr(0, 3) == " ")) {
 if (ss.size() > 0 && N1 != 0)
 throw std::runtime_error("Secular coeffs not last");
 if (ss.size() == 0 && N1 == 0)
 throw std::runtime_error("No secular coeffs");
 storecoeff(fout, C, S, M, N);
 if (ss.size() == 0)
 storecoeff(fout, C1, S1, M1, N1);
 start = true;
 }
 }
 } else if (model == "wmm2015" || model == "igrf12") {
 // Download WMM2015COF.zip
 // http://ngdc.noaa.gov/geomag/WMM/WMM_coeff.shtml
 // wmm2015 coefficients are in WMM2015COF/WMM.COF
 //
 // igrf12 coefficients
 // http://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
 // igrf12 coefficients are in geomag70_linux/IGRF11.COF
 std::string id = model == "wmm2015" ? "WMM2015A" : "IGRF12-A";
 fout.write(id.c_str(), 8);
 std::string filename = model == "wmm2015" ? "WMM2015COF/WMM.COF"
 : "geomag70_linux-2015/IGRF12.COF";
 std::ifstream fin(filename.c_str());
 std::string ss;
 bool start = true;
 std::getline(fin, ss);
 std::vector<double> C;
 std::vector<double> S;
 std::vector<double> C1;
 std::vector<double> S1;
 int N, M, N1, M1;
 while (true) {
 if (ss.size() == 0)
 break;
 std::istringstream is(ss);
 int n, m;
 double c, s, c1, s1;
 if (start) {
 if (model == "wmm2015") {
 N = 12; N1 = 12;
 } else {
 std::string mm;
 if (!(is >> mm >> mm >> N >> N1))
 throw std::runtime_error("Short read on header");
 }
 M = N; M1 = N1;
 int
 K = (M + 1) * (2*N - M + 2) / 2,
 K1 = (M1 + 1) * (2*N1 - M1 + 2) / 2;
 C.resize(K);
 S.resize(K - (N + 1));
 C1.resize(K1);
 S1.resize(K1 - (N1 + 1));
 std::fill(C.begin(), C.end(), 0.0);
 std::fill(S.begin(), S.end(), 0.0);
 std::fill(C1.begin(), C1.end(), 0.0);
 std::fill(S1.begin(), S1.end(), 0.0);
 start = false;
 } else {
 if (!(is >> n >> m >> c >> s >> c1 >> s1)) {
 std::cerr << ss << "\n";
 throw std::runtime_error("Short read on data");
 }
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 storecos(C1, N1, M1, c1, n, m);
 storesin(S1, N1, M1, s1, n, m);
 }
 if (!std::getline(fin, ss))
 ss = "";
 if (model == "wmm2015") {
 if (ss.size() && ss[0] == '9')
 ss = "";
 if (ss.size() == 0) {
 storecoeff(fout, C, S, M, N);
 storecoeff(fout, C1, S1, M1, N1);
 }
 } else {
 if (ss.size() == 0 || (ss.size() >= 3 && ss.substr(0, 3) == " ")) {
 if (ss.size() > 0 && N1 != 0)
 throw std::runtime_error("Secular coeffs not last");
 if (ss.size() == 0 && N1 == 0)
 throw std::runtime_error("No secular coeffs");
 storecoeff(fout, C, S, M, N);
 if (ss.size() == 0)
 storecoeff(fout, C1, S1, M1, N1);
 start = true;
 }
 }
 }
 } else if (model == "egm2008" || model == "egm96" || model == "egm84"
 || model == "wgs84") {
 std::string id = model == "egm2008" ? "EGM2008A" :
 (model == "egm96" ? "EGM1996A" :
 (model == "egm84" ? "EGM1984A" : "WGS1984A"));
 fout.write(id.c_str(), 8);
 for (int i = 0; i < 2; ++i) {
 std::string filename = "../gravity/";
 filename +=
 model == "egm2008" ?
 (i == 0 ? "EGM2008_to2190_TideFree" : "Zeta-to-N_to2160_egm2008") :
 (model == "egm96" ?
 (i == 0 ? "EGM96" : "CORRCOEF") :
 (model == "egm84" ?
 (i == 0 ? "egm180.nor" : "zeta84") :
 (i == 0 ? "wgs84.cof" : "zeta84")));
 std::ifstream fin(filename.c_str());
 if (!fin.good())
 throw std::runtime_error("File not found");
 int
 N = model == "egm2008" ? (i == 0 ? 2190 : 2160) :
 (model == "egm96" ? 360 :
 (model == "egm84" ? (i == 0 ? 180 : -1) :
 (i == 0 ? 20 : -1))),
 M = model == "egm2008" ? (i == 0 ? 2159 : 2160) :
 (model == "egm96" ? 360 :
 (model == "egm84" ? (i == 0 ? 180 : -1) :
 (i == 0 ? 0 : -1))),
 K = (M + 1) * (2*N - M + 2) / 2;
 std::vector<double> C(K, 0.0);
 std::vector<double> S(K - (N + 1), 0.0);
 std::string ss;
 while (std::getline(fin, ss)) {
 std::string::size_type p = 0;
 while (true) {
 p = ss.find_first_of('D', p);
 if (p < ss.size())
 ss[p] = 'E';
 else
 break;
 ++p;
 }
 std::istringstream is(ss);
 int n, m;
 double c, s;
 if (!(is >> n >> m >> c >> s))
 throw std::runtime_error("Short read");
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 }
 storecoeff(fout, C, S, M, N);
 }
 } else if (model == "dtm2006") {
 std::string id = "DTM2006A";
 fout.write(id.c_str(), 8);
 int N = 2190, M = N;
 int K = (M + 1) * (2*N - M + 2) / 2;
 std::string filename
 = "../gravity/Coeff_Height_and_Depth_to2190_DTM2006.0";
 std::ifstream fin(filename.c_str());
 std::vector<double> C(K, 0.0);
 std::vector<double> S(K - (N + 1), 0.0);
 std::string ss;
 while (std::getline(fin, ss)) {
 std::string::size_type p = 0;
 while (true) {
 p = ss.find_first_of('D', p);
 if (p < ss.size())
 ss[p] = 'E';
 else
 break;
 ++p;
 }
 std::istringstream is(ss);
 int n, m;
 double c, s;
 if (!(is >> n >> m >> c >> s))
 throw std::runtime_error("Short read");
 storecos(C, N, M, c, n, m);
 storesin(S, N, M, s, n, m);
 }
 storecoeff(fout, C, S, M, N);
 } else {
 std::cerr << "UNKNOWN MODEL " << model << "\n";
 return 1;
 }
 }
 catch (const std::exception& e) {
 std::cerr << "ERROR: " << e.what() << "\n";
 return 1;
 }
 return 0;
}