File: grailCdr2Orbit.cpp

package info (click to toggle)
groops 0%2Bgit20250907%2Bds-1
  • links: PTS, VCS
  • area: non-free
  • in suites: forky, sid
  • size: 11,140 kB
  • sloc: cpp: 135,607; fortran: 1,603; makefile: 20
file content (148 lines) | stat: -rw-r--r-- 3,677 bytes parent folder | download | duplicates (2)
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
/***********************************************/
/**
* @file grailCdr2Orbit.cpp
*
* @brief read CDR GRAIL data.
*
* @author Beate Klinger
* @date 2013-01-22
*
*/
/***********************************************/

// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
This program converts the orbit from the GRAIL SDS format into  \file{instrument file (ORBIT)}{instrument}.
)";

/***********************************************/

#include "programs/program.h"
#include "base/string.h"
#include "files/fileInstrument.h"

/***** CLASS ***********************************/

/** @brief read CDR GRAIL data.
* @ingroup programsConversionGroup */
class GrailCdr2Orbit
{
  void readFile(const FileName fileName, OrbitArc &arc);

public:
  void run(Config &config, Parallel::CommunicatorPtr comm);
};

GROOPS_REGISTER_PROGRAM(GrailCdr2Orbit, SINGLEPROCESS, "read CDR GRAIL data", Conversion, Orbit, Instrument)

/***********************************************/

void GrailCdr2Orbit::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
  try
  {
    FileName outName;
    std::vector<FileName> fileNames;

    readConfig(config, "outputfileOrbit", outName,   Config::MUSTSET,  "", "");
    readConfig(config, "inputfile",       fileNames, Config::MUSTSET,  "", "");
    if(isCreateSchema(config)) return;

    logStatus<<"read input files"<<Log::endl;
    OrbitArc orbit;
    for(auto &fileName : fileNames)
    {
      logStatus<<"read file <"<<fileName<<">"<<Log::endl;
      readFile(fileName, orbit);
    }

    logStatus<<"write data to <"<<outName<<">"<<Log::endl;
    InstrumentFile::write(outName, orbit);
    Arc::printStatistics(orbit);
  }
  catch(std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}

/***********************************************/

void GrailCdr2Orbit::readFile(const FileName fileName, OrbitArc &arc)
{
  try
  {
    InFile file(fileName);
    if(!file.good())
    {
      logWarning<<"cannot open file: "<<fileName.str()<<", continue..."<<Log::endl;
      return;
    }
    file.exceptions(std::ios::badbit|std::ios::failbit);

    // Header
    UInt numberOfRecords = 0;
    std::string line;
    for(;;)
    {
      getline(file, line);
      if(line.find("NUMBER OF DATA RECORDS")==0)
        numberOfRecords = String::toInt(line.substr(31, 10));
      if(line.find("END OF HEADER")==0)
        break;
    }

    // Data
    for(UInt i=0; i<numberOfRecords; i++)
    {
      std::string line;
      try
      {
        std::getline(file, line);
      }
      catch(std::exception &/*e*/)
      {
        //logWarning<<std::endl<<e.what()<<" continue..."<<Log::endl;
        break;
      }
      std::stringstream ss(line);

      // time (TDB) >> seconds past 12:00:00 noon 01-Jan-2000
      UInt seconds;
      ss>>seconds;
      OrbitEpoch epoch;
      epoch.time = timeTT2GPS(mjd2time(51544.5)+seconds2time(seconds));

      // satellite ID, coordinate reference frame
      Int32 AB;
      Byte  IE; //Byte
      ss>>AB>>IE;

      // position + formal error in position
      Double x_error, y_error, z_error;
      ss>>epoch.position.x()>>epoch.position.y()>>epoch.position.z();
      ss>>x_error>>y_error>>z_error;

      // velocity + formalerror in velocity
      Double dx_error, dy_error, dz_error;
      ss>>epoch.velocity.x()>>epoch.velocity.y()>>epoch.velocity.z();
      ss>>dx_error>>dy_error>>dz_error;

      // data quality flags >> data quality flags not defined!
      //Int flag;
      //ss>>flag;
      //if(flag!=0)
      //  continue;

      arc.push_back(epoch);
    }
  }

  catch(std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}

/***********************************************/