coursepointer 0.3.4

#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 degree or order");
  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 degree or order");
  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 = 0, M = 0, N1 = 0, M1 = 0;
      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 == "wmm2015v2" ||
               model == "wmm2020" ||
               model == "wmm2025" || model == "wmmhr2025" ||
               // defer treatment of igfr1[23] tp igrf14 below
               model == "igrf12xx" || model == "igrf13xx") {
      // Download WMM2015COF.zip
      // http://ngdc.noaa.gov/geomag/WMM/WMM_coeff.shtml
      // wmm2015 coefficients are in WMM2015COF/WMM.COF
      //
      // Download WMM2015v2COF.zip
      // https://www.ngdc.noaa.gov/geomag/WMM/data/WMM2015/WMM2015v2COF.zip
      // wmm2015v2 coefficients are in WMM2015v2COF/WMM.COF
      //
      // Download WMM2020COF.zip
      // http://ngdc.noaa.gov/geomag/WMM/WMM_coeff.shtml
      // wmm2020 coefficients are in WMM2020COF/WMM.COF
      //
      // Download WMM2025COF.zip
      // https://www.ncei.noaa.gov/products/world-magnetic-model
      // wmm2025 coefficients are in WMM2025COF/WMM.COF
      //
      // Download WMMHR2025COF.zip
      // https://www.ncei.noaa.gov/products/world-magnetic-model-high-resolution
      // wmmhr2025 coefficients are in WMMHR2025COF/WMMHR.COF
      //
      // igrf12 coefficients
      // http://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
      // igrf12 coefficients are in geomag70_linux/IGRF12.COF
      //
      // igrf13 coefficients
      // https://www.ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
      // igrf13 coefficients are in geomag70_linux/IGRF13.COF
      std::string id = model == "wmm2015" ? "WMM2015A" :
        (model == "wmm2015v2" ?  "WMM2015B" :
         (model == "wmm2020" ? "WMM2020A" :
          (model == "wmm2025" ? "WMM2025A" :
           (model == "wmmhr2025" ? "WMM2025H" :
            (model == "igrf12" ? "IGRF12-A" : "IGRF13-A")))));
      fout.write(id.c_str(), 8);
      std::string filename = model == "wmm2015" ? "WMM2015COF/WMM.COF"
        : (model == "wmm2015v2" ? "WMM2015v2COF/WMM.COF"
           : (model == "wmm2020" ? "WMM2020COF/WMM.COF"
              : (model == "wmm2025" ? "WMM2025COF/WMM.COF"
                 : (model == "wmmhr2025" ? "WMMHR2025COF/WMMHR.COF"
                    : (model == "igrf12" ? "geomag70_linux-2015/IGRF12.COF"
                       : "geomag70_linux-2020/IGRF13.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 = 0, M = 0, N1 = 0, M1 = 0;
      while (true) {
        if (ss.size() == 0)
          break;
        std::istringstream is(ss);
        int n, m;
        double c, s, c1, s1;
        if (start) {
          if (model == "wmm2015" || model == "wmm2015v2" ||
              model == "wmm2020" || model == "wmm2025" ) {
            N = 12; N1 = 12;
          } else if (model == "wmmhr2025") {
            N = 133; N1 = 133;
          } 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" || model == "wmm2015v2" ||
            model == "wmm2020" ||
            model == "wmm2025" || model == "wmmhr2025") {
          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 == "igrf12" || model == "igrf13" || model == "igrf14") {
      std::string id = "IGRFnn-A";
      id[4] = model[4]; id[5] = model[5];
      fout.write(id.c_str(), 8);
      // https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf12coeffs.txt
      // https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf13coeffs.txt
      // https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt
      std::string filename = model + "coeffs.txt";
      std::ifstream fin(filename.c_str());
      // format for igrf14coeffs.txt (gen = 14)
      // comments start with #
      // 2 lines of headers
      // data = c/sflag n m coeffsx(gen+12) sv
      // coeffs are for 1990+5*i for i in [0, gen+12)
      // sv covers period 1900+5*(gen+12) + [-5,0]
      // order for first 20 columns 1900-1995 is 10
      // order for rest 13
      // order of last sv column is 8
      int L = 10*(model[4]-'0') + (model[5]-'0') + 13, // gen + 13
        L0 = 20;
      std::vector<std::vector<double>> C(L);
      std::vector<std::vector<double>> S(L);
      std::vector<int> N(L);
      std::vector<int> M(L);
      for (int l = 0; l < L; ++l) {
        N[l] = l == L-1 ? 8 : (l < L0 ? 10 : 13);
        M[l] = N[l];
        int K = (M[l] + 1) * (2*N[l] - M[l] + 2) / 2;
        C[l].resize(K, 0.0);
        S[l].resize(K - (N[l] + 1), 0.0);
      }
      std::string ss;
      int nhead = 0;
      while (nhead < 2) {
        if (!std::getline(fin, ss) || ss.empty())
          throw std::runtime_error("Short read on file");
        if (ss[0] == '#') continue;
        ++nhead;
      }
      char flag;
      int n, m;
      double v;
      while (std::getline(fin, ss)) {
        std::istringstream is(ss);
        if (!(is >> flag >> n >> m))
          throw std::runtime_error("Short read on data");
        for (int l = 0; l < L; ++l) {
          if (!(is >> v))
            throw std::runtime_error("Short read on coeffs");
          if (flag == 'g')
            storecos(C[l], N[l], M[l], v, n, m);
          else if (flag == 'h')
            storesin(S[l], N[l], M[l], v, n, m);
          else
            throw std::runtime_error("Bad flag value");
        }
      }
      for (int l = 0; l < L; ++l)
        storecoeff(fout, C[l], S[l], M[l], N[l]);
    } 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;
}