A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
geolocate_oci.cpp
Go to the documentation of this file.
1 #include <fstream>
2 #include <iostream>
3 #include <sstream>
4 #include <string>
5 
6 #include <gsl/gsl_blas.h>
7 #include <gsl/gsl_linalg.h>
8 #include <gsl/gsl_matrix_double.h>
9 #include <boost/algorithm/string/trim.hpp>
10 #include <boost/date_time/posix_time/posix_time.hpp>
11 #include <boost/numeric/ublas/matrix.hpp>
12 
13 #include "geolocate_oci.h"
14 
15 #include <netcdf>
16 
17 #include "nc4utils.h"
18 #include "timeutils.h"
19 #include <clo.h>
20 #include "global_attrs.h"
21 #include <terrain.h>
22 #include <genutils.h>
23 #define ERROR_EXIT(dim) {if(dim==0){printf("--Error--: the dimension %s is zero. Exiting. See line %d in file %s",#dim,__LINE__, __FILE__);exit(EXIT_FAILURE);} }
24 
25 // Modification history:
26 // Programmer Organization Date Ver Description of change
27 // ---------- ------------ ---- --- ---------------------
28 // Joel Gales SAIC 09/13/21 0.01 Original development
29 // Joel Gales SAIC 06/14/22 0.03 Add support for tilt change
30 // Gwyn Fireman SAIC 07/10/23 0.07 Read global metadata from json file
31 // Joel Gales SAIC 08/08/23 0.08 Check scan times
32 // Joel Gales SAIC 08/18/23 0.80 Add earth_sun distance
33 // correction attribute, support
34 // for reflectance output
35 // Joel Gales SAIC 08/29/23 0.81 Change fill values
36 // Joel Gales SAIC 10/03/23 0.85 Change tilt output name to
37 // tilt_angle
38 // Change to callable routine
39 // rather than standalone program
40 
41 using namespace std;
42 
51 void read_def_eph(string ephfile, vector<double>& otime,
52  vector<vector<double>>& posj, vector<vector<double>>& velj, double l1aEpoch) {
53 
54  ifstream inStream(ephfile, ios::in);
55  if (!inStream.is_open()) {
56  cerr << "-E- Unable to open definitive ephemeris file\n";
57  exit(EXIT_FAILURE);
58  }
59 
60  string tmpBuf;
61 
62  //size_t numRecords = 1921; // Number of records per file (32 hours)
63  vector<double> posVec(3);
64  vector<double> velVec(3);
65  string ephStr;
66  vector<string> tokens; // Datetime, 3d position, 3d velocity.
67  vector<double> secondsVec(1921);
68  double oldTime = 0;
69 
70  while (inStream.peek() != EOF) {
71  getline(inStream, ephStr, '\n');
72 
73  istringstream lineStream(ephStr);
74  string token; // Each individual token in this line
75 
76  tokens.clear();
77  while (lineStream >> token) {
78  tokens.push_back(token);
79  }
80 
81  if(tokens.size() != 7 || tokens[0].find("T") == string::npos) {
82  continue;
83  }
84 
85  // skip duplicate records
86  double theTime = isodate2unix(tokens[0].c_str());
87  if(oldTime == theTime)
88  continue;
89  oldTime = theTime;
90 
91  posVec.clear();
92  posVec.push_back(stod(tokens[1]));
93  posVec.push_back(stod(tokens[2]));
94  posVec.push_back(stod(tokens[3]));
95 
96  velVec.clear();
97  velVec.push_back(stod(tokens[4]));
98  velVec.push_back(stod(tokens[5]));
99  velVec.push_back(stod(tokens[6]));
100 
101  posj.push_back(posVec);
102  velj.push_back(velVec);
103  otime.push_back(theTime - l1aEpoch);
104 
105  ephStr.clear();
106  lineStream.clear();
107  }
108  inStream.close();
109 
110  if (otime.size() == 0) {
111  cerr << "-E- No records found in definitive ephemeris file\n";
112  exit(EXIT_FAILURE);
113  }
114 
115 }
116 
117 int geolocate_oci( string l1a_filename, string geo_lut_filename,
118  geo_struct& geoLUT,
119  string l1b_filename, string dem_file, bool radiance,
120  string doi, const string ephFile, string pversion) { // add ephemeris file
121 
122  ofstream tempOut;
123 
124  // Open and read data from L1Afile
125  NcFile *l1afile = new NcFile( l1a_filename, NcFile::read);
126 
127  NcGroup ncGrps[6];
128 
129  ncGrps[0] = l1afile->getGroup( "scan_line_attributes");
130  ncGrps[1] = l1afile->getGroup( "spatial_spectral_modes");
131  ncGrps[2] = l1afile->getGroup( "engineering_data");
132  ncGrps[3] = l1afile->getGroup( "navigation_data");
133  ncGrps[4] = l1afile->getGroup( "onboard_calibration_data");
134  ncGrps[5] = l1afile->getGroup( "science_data");
135 
136  NcGroupAtt att;
137  NcVar var;
138 
139  // Append call sequence to existing history
140  string history = get_history(l1afile);
141  // history.append(call_sequence(argc, argv));
142 
143  // Get date (change this when year and day are added to time field)
144  string tstart, tend;
145  att = l1afile->getAtt("time_coverage_start");
146  att.getValues(tstart);
147 
148  att = l1afile->getAtt("time_coverage_end");
149  att.getValues(tend);
150 
151  size_t pos = tstart.find("T");
152  if(pos == string::npos) {
153  cout << "-E- L1A file time_coverage_start, bad format\n";
154  exit(EXIT_FAILURE);
155  }
156  string granule_date = tstart.substr(0,pos);
157 
158  // Vector of the positions of OCI for 32 hours starting at tstart with a resolution of 1 minute
159  vector<vector<double>> ephPosVec;
160  // Vector of the velocities of OCI for 32 hours starting at tstart with a resolution of 1 minute
161  vector<vector<double>> ephVelVec;
162  vector<double> ephOTime;
163 
164  if (!ephFile.empty()) {
165  double l1a_start = isodate2unix(tstart.c_str());
166  double l1a_end = isodate2unix(tend.c_str());
167 
168  int16_t yr,mo,dy;
169  double secs;
170  unix2ymds(l1a_start, &yr, &mo, &dy, &secs);
171  double l1a_epoch = ymds2unix(yr, mo, dy, 0.0);
172 
173  read_def_eph(ephFile, ephOTime, ephPosVec, ephVelVec, l1a_epoch);
174 
175  // Ensure ephem file covers L1A file time
176  if((l1a_start < (ephOTime.front() + l1a_epoch)) ||
177  (l1a_end > (ephOTime.back() + l1a_epoch))) {
178  cerr << "-E- Definitive ephemeris file does not cover L1A file time" << endl;
179  exit(EXIT_FAILURE);
180  }
181  }
182 
183  // Scan time, spin ID and HAM side
184  NcDim nscan_dim = l1afile->getDim("number_of_scans");
185  uint32_t nscan = nscan_dim.getSize();
186 
187  // Number of swir bands
188 
189  uint32_t nswband = l1afile->getDim("SWIR_bands").getSize();
190 
191  NcDim mcescan_dim = l1afile->getDim("number_of_mce_scans");
192  uint32_t nmcescan = mcescan_dim.getSize();
193 
194  double *sstime = new double[nscan];
195  var = ncGrps[0].getVar( "scan_start_time");
196  var.getVar( sstime);
197 
198  int32_t *spin = new int32_t[nscan];
199  var = ncGrps[0].getVar( "spin_ID");
200  var.getVar( spin);
201 
202  uint8_t *hside = new uint8_t[nscan];
203  var = ncGrps[0].getVar( "HAM_side");
204  var.getVar( hside);
205 
206  uint32_t sdim=0;
207  for (size_t i=0; i<nscan; i++) {
208  if ( spin[i] > 0) {
209  sstime[sdim] = sstime[i];
210  hside[sdim] = hside[i];
211  sdim++;
212  }
213  }
214 
215  // Check for and fill in missing scan times
216  short *sfl = new short[sdim]();
217  check_scan_times( sdim, sstime, sfl);
218 
219  NcDim natt_dim = l1afile->getDim("att_records");
220  uint32_t n_att_rec = natt_dim.getSize();
221  ERROR_EXIT(n_att_rec);
222 
223  double *att_time = new double[n_att_rec];
224  var = ncGrps[3].getVar( "att_time");
225  var.getVar( att_time);
226 
227  NcDim nquat_dim = l1afile->getDim("quaternion_elements");
228  uint32_t n_quat_elems = nquat_dim.getSize();
229  ERROR_EXIT(n_quat_elems);
230 
231  float **att_quat = new float *[n_att_rec];
232  att_quat[0] = new float[n_quat_elems*n_att_rec];
233  for (size_t i=1; i<n_att_rec; i++)
234  att_quat[i] = att_quat[i-1] + n_quat_elems;
235 
236  var = ncGrps[3].getVar( "att_quat");
237  var.getVar( &att_quat[0][0]);
238 
239  NcDim norb_dim = l1afile->getDim("orb_records");
240  uint32_t n_orb_rec = norb_dim.getSize();
241  ERROR_EXIT(n_orb_rec);
242 
243  vector<double> orb_time(n_orb_rec);
244  if (ephFile.empty()) {
245  var = ncGrps[3].getVar( "orb_time");
246  var.getVar( orb_time.data());
247  }
248 
249  float **orb_pos = new float *[n_orb_rec];
250  orb_pos[0] = new float[3*n_orb_rec];
251  for (size_t i=1; i<n_orb_rec; i++) orb_pos[i] = orb_pos[i-1] + 3;
252  var = ncGrps[3].getVar( "orb_pos");
253  var.getVar( &orb_pos[0][0]);
254 
255  float **orb_vel = new float *[n_orb_rec];
256  orb_vel[0] = new float[3*n_orb_rec];
257  for (size_t i=1; i<n_orb_rec; i++) orb_vel[i] = orb_vel[i-1] + 3;
258  var = ncGrps[3].getVar( "orb_vel");
259  var.getVar( &orb_vel[0][0]);
260 
261  NcDim ntilt_dim = l1afile->getDim("tilt_samples");
262  uint32_t n_tilts = ntilt_dim.getSize();
263  ERROR_EXIT(n_tilts);
264 
265  float *tiltin = new float[n_tilts];
266  var = ncGrps[3].getVar( "tilt");
267  var.getVar( tiltin);
268 
269  double *ttime = new double[n_tilts];
270  var = ncGrps[3].getVar( "tilt_time");
271  var.getVar( ttime);
272 
273 
274  // **************************** //
275  // *** Read geolocation LUT *** //
276  // **************************** //
277 
278  NcFile *geoLUTfile = new NcFile( geo_lut_filename, NcFile::read);
279  // geo_struct geoLUT;
280 
281  NcGroup gidTime, gidCT, gidRTA_HAM, gidPlanarity;
282 
283  gidTime = geoLUTfile->getGroup( "time_params");
284  var = gidTime.getVar( "master_clock");
285  var.getVar( &geoLUT.master_clock);
286  var = gidTime.getVar( "MCE_clock");
287  var.getVar( &geoLUT.MCE_clock);
288 
289  gidCT = geoLUTfile->getGroup( "coord_trans");
290  var = gidCT.getVar( "sc_to_tilt");
291  var.getVar( &geoLUT.sc_to_tilt);
292  var = gidCT.getVar( "tilt_axis");
293  var.getVar( &geoLUT.tilt_axis);
294  var = gidCT.getVar( "tilt_angles");
295  var.getVar( &geoLUT.tilt_angles);
296  var = gidCT.getVar( "tilt_home");
297  var.getVar( &geoLUT.tilt_home);
298  var = gidCT.getVar( "tilt_to_oci_mech");
299  var.getVar( &geoLUT.tilt_to_oci_mech);
300  var = gidCT.getVar( "oci_mech_to_oci_opt");
301  var.getVar( &geoLUT.oci_mech_to_oci_opt);
302 
303  gidRTA_HAM = geoLUTfile->getGroup( "RTA_HAM_params");
304  var = gidRTA_HAM.getVar( "RTA_axis");
305  var.getVar( &geoLUT.rta_axis);
306  var = gidRTA_HAM.getVar( "HAM_axis");
307  var.getVar( &geoLUT.ham_axis);
308  var = gidRTA_HAM.getVar( "HAM_AT_angles");
309  var.getVar( &geoLUT.ham_at_angles);
310  var = gidRTA_HAM.getVar( "HAM_CT_angles");
311  var.getVar( &geoLUT.ham_ct_angles);
312  var = gidRTA_HAM.getVar( "RTA_enc_scale");
313  var.getVar( &geoLUT.rta_enc_scale);
314  var = gidRTA_HAM.getVar( "HAM_enc_scale");
315  var.getVar( &geoLUT.ham_enc_scale);
316  var = gidRTA_HAM.getVar( "RTA_nadir");
317  var.getVar( geoLUT.rta_nadir);
318 
319  gidPlanarity = geoLUTfile->getGroup( "planarity");
320  var = gidPlanarity.getVar( "along_scan_planarity");
321  var.getVar( &geoLUT.as_planarity);
322  var = gidPlanarity.getVar( "along_track_planarity");
323  var.getVar( &geoLUT.at_planarity);
324 
325  geoLUTfile->close();
326 
327  // MCE telemetry
328  int32_t ppr_off;
329  double revpsec, secpline;
330 
331  int32_t *mspin = new int32_t[nmcescan];
332  int32_t *ot_10us = new int32_t[nmcescan];
333  uint8_t *enc_count = new uint8_t[nmcescan];
334  int16_t board_id;
335 
336  NcDim nenc_dim = l1afile->getDim("encoder_samples");
337  uint32_t nenc = nenc_dim.getSize();
338 
339  float **hamenc = new float *[nmcescan];
340  hamenc[0] = new float[nenc*nmcescan];
341  for (size_t i=1; i<nmcescan; i++) hamenc[i] = hamenc[i-1] + nenc;
342 
343  float **rtaenc = new float *[nmcescan];
344  rtaenc[0] = new float[nenc*nmcescan];
345  for (size_t i=1; i<nmcescan; i++) rtaenc[i] = rtaenc[i-1] + nenc;
346 
347  int16_t iret;
348  read_mce_tlm( l1afile, geoLUT, ncGrps[2], nmcescan, nenc, ppr_off, revpsec,
349  secpline, board_id, mspin, ot_10us, enc_count, hamenc, rtaenc,
350  iret);
351 
352  int32_t iyr, iday, msec;
353  isodate2ydmsec((char *) tstart.c_str(), &iyr, &iday, &msec);
354 
355  double ecmat[3][3];
356  quat_array *quatr = new quat_array[n_att_rec];
357  orb_array *posr; // Array of 1x3 vectors
358  orb_array *velr; // Array of 1x3 vectors
359 
360  // Transform orbit (if needed) and attitude from J2000 to ECR
361  if (!ephFile.empty()) {
362 
363  size_t ephOTimeLBIndex = -1;
364  size_t ephOTimeUBIndex = -1;
365 
366  for (size_t i = 0; i < ephOTime.size(); i++) {
367  double value = ephOTime[i];
368 
369  if (value <= sstime[0]) {
370  ephOTimeLBIndex = i;
371  }
372  if (value >= sstime[nscan - 1]) {
373  ephOTimeUBIndex = i;
374  break;
375  }
376  }
377 
378 
379  double omega_e = 7.29211585494e-5;
380  size_t num_eph_recs = ephOTimeUBIndex - ephOTimeLBIndex + 1;
381 
382  n_orb_rec = num_eph_recs;
383  posr = new orb_array[n_orb_rec]; // Array of 1x3 vectors
384  velr = new orb_array[n_orb_rec]; // Array of 1x3 vectors
385 
386  bzero(posr, n_orb_rec * sizeof(orb_array));
387  bzero(velr, n_orb_rec * sizeof(orb_array));
388 
389  orb_time.resize(n_orb_rec);
390 
391  for (size_t idxRec = 0; idxRec < num_eph_recs; idxRec++) {
392  size_t ephIndex = idxRec + ephOTimeLBIndex;
393  orb_time[idxRec] = ephOTime[ephIndex];
394 
395  j2000_to_ecr(iyr, iday, ephOTime[ephIndex], ecmat);
396 
397  for (size_t idxLine = 0; idxLine < 3; idxLine++) {
398  for (size_t idxValue = 0; idxValue < 3; idxValue++) {
399  posr[idxRec][idxLine] += ecmat[idxLine][idxValue] * ephPosVec[ephIndex][idxValue];
400  velr[idxRec][idxLine] += ecmat[idxLine][idxValue] * ephVelVec[ephIndex][idxValue];
401  }
402  }
403 
404  velr[idxRec][0] = velr[idxRec][0] + posr[idxRec][1] * omega_e;
405  velr[idxRec][1] = velr[idxRec][1] - posr[idxRec][0] * omega_e;
406  }
407  } else {
408  posr = new orb_array[n_orb_rec]; // Array of 1x3 vectors
409  velr = new orb_array[n_orb_rec]; // Array of 1x3 vectors
410  // Orbit
411  for (size_t i = 0; i < n_orb_rec; i++) {
412  // j2000_to_ecr(iyr, iday, orb_time[i], ecmat);
413 
414  for (size_t j = 0; j < 3; j++) {
415  posr[i][j] = orb_pos[i][j] / 1000;
416  velr[i][j] = orb_vel[i][j] / 1000;
417  }
418  } // i loop
419  }
420 
421 
422  // Attitude
423  for (size_t i = 0; i < n_att_rec; i++) {
424  double ecq[4], qt2[4];
425  float qt1[4];
426  j2000_to_ecr(iyr, iday, att_time[i], ecmat);
427  mtoq(ecmat, ecq);
428 
429  memcpy(qt1, &att_quat[i][0], 3 * sizeof(float));
430  qt1[3] = att_quat[i][3];
431 
432  qprod(ecq, qt1, qt2);
433  for (size_t j = 0; j < 4; j++) quatr[i][j] = qt2[j];
434  } // i loop
435 
436  // ******************************************** //
437  // *** Get spatial and spectral aggregation *** //
438  // ******************************************** //
439  NcDim spatzn_dim = l1afile->getDim("spatial_zones");
440  uint32_t spatzn = spatzn_dim.getSize();
441 
442  int16_t *dtype = new int16_t[spatzn];
443  var = ncGrps[1].getVar( "spatial_zone_data_type");
444  var.getVar( dtype);
445 
446  int16_t *lines = new int16_t[spatzn];
447  var = ncGrps[1].getVar( "spatial_zone_lines");
448  var.getVar( lines);
449 
450  int16_t *iagg = new int16_t[spatzn];
451  var = ncGrps[1].getVar( "spatial_aggregation");
452  var.getVar( iagg);
453 
454  float clines[32400], slines[4050];
455  uint16_t pcdim, psdim;
456  double ev_toff, deltc[32400], delts[4050];
457  bool dark = false;
458  get_ev( secpline, dtype, lines, iagg, pcdim, psdim, ev_toff,
459  clines, slines, deltc, delts, dark, iret);
460  if ( iret < 0) {
461  cout << "No Earth view in file: " << l1a_filename << endl;
462  l1afile->close();
463  return 1;
464  }
465 
466  double *evtime = new double[sdim];
467  for (size_t i = 0; i < sdim; i++)
468  if (sstime[i] == -32767) evtime[i] = sstime[i];
469  else
470  evtime[i] = sstime[i] + ev_toff;
471 
472  // Interpolate orbit, attitude and tilt to scan times
473  orb_array *posi = new orb_array[sdim]();
474  orb_array *veli = new orb_array[sdim]();
475  orb_array *rpy = new orb_array[sdim]();
476  orb_interp(n_orb_rec, sdim, orb_time.data(), posr, velr, evtime, posi, veli);
477 
478  quat_array *quati = new quat_array[sdim]();
479  q_interp(n_att_rec, sdim, att_time, quatr, evtime, quati);
480 
481  float *tilt = new float[sdim];
482  tilt_interp(n_tilts, sdim, ttime, tiltin, evtime, tilt);
483  for (size_t i = 0; i < sdim; i++)
484  sfl[i] |= 1; // set tilt flag for all scans
485  for (size_t i = 0; i < sdim; i++) {
486  tilt[i] += geoLUT.tilt_home; // Add tilt home position to angles
487  tilt[i] = tilt[i] < geoLUT.tilt_angles[0] ? geoLUT.tilt_angles[0] : tilt[i];
488  tilt[i] = tilt[i] > geoLUT.tilt_angles[1] ? geoLUT.tilt_angles[1] : tilt[i];
489  // Fred's tilt detection
490  // for any 2 adjecent and = tilts clear tilt flag (ie, bit AND with 110b)
491  if( ( i > 0 ) && ( tilt[i-1] == tilt[i] ) ) {
492  sfl[i-1] &= 6;
493  sfl[i] &= 6;
494  }
495  }
496  float *xlon = new float[sdim * pcdim]();
497  float *xlat = new float[sdim * pcdim]();
498  short *solz = new short[sdim * pcdim]();
499  short *sola = new short[sdim * pcdim]();
500  short *senz = new short[sdim * pcdim]();
501  short *sena = new short[sdim * pcdim]();
502 
503  uint8_t *qfl = new uint8_t[sdim * pcdim]();
504  // short *range = new short[sdim * psdim]();
505  short *height = new short[sdim * pcdim]();
506 
507  // Initialize output arrays
508  for (size_t i = 0; i < sdim * pcdim; i++) {
509  xlon[i] = -32767;
510  xlat[i] = -32767;
511 
512  solz[i] = -32767;
513  sola[i] = -32767;
514  senz[i] = -32767;
515  sena[i] = -32767;
516  }
517 
518  // Generate pointing vector array
519  float **pview = new float *[pcdim];
520  pview[0] = new float[3*pcdim];
521  for (size_t i=1; i<pcdim; i++) pview[i] = pview[i-1] + 3;
522 
523  // Get Sun vectors
524  orb_array *sunr = new orb_array[sdim]();
525 
526  double *rs = new double[sdim];
527  l_sun(sdim, iyr, iday, evtime, sunr, rs);
528  double distcorr = pow(rs[sdim/2],2);
529 
530  // S/C matrices
531  gsl_matrix *sc2tiltp = gsl_matrix_alloc(3, 3);
532  gsl_matrix *sc2ocim = gsl_matrix_alloc(3, 3);
533  gsl_matrix *sc_to_oci = gsl_matrix_alloc(3, 3);
534  gsl_matrix *smat = gsl_matrix_alloc(3, 3);
535 
536  double *thetap = new double[pcdim]();
537 
539  // Geolocate each scan line //
541  float lonwest = +180;
542  float loneast = -180;
543  float latsouth = +90;
544  float latnorth = -90;
545  for (size_t iscn = 0; iscn < sdim; iscn++) {
546 
547  if (sstime[iscn] != -999.0) {
548 
549  gsl_matrix_view A;
550  gsl_matrix_view B;
551  double *ptr_C;
552 
553  // Model tilt rotation using a quaternion
554  float qt[4];
555  qt[0] = geoLUT.tilt_axis[0]*sin(tilt[iscn]/2/RADEG);
556  qt[1] = geoLUT.tilt_axis[1]*sin(tilt[iscn]/2/RADEG);
557  qt[2] = geoLUT.tilt_axis[2]*sin(tilt[iscn]/2/RADEG);
558  qt[3] = cos(tilt[iscn]/2/RADEG);
559 
560  double tiltm[3][3];
561  qtom(qt, tiltm);
562 
563  // Combine tilt and alignments
564 
565  // sc2tiltp = tiltm#geo_lut.sc_to_tilt
566  A = gsl_matrix_view_array( (double *) geoLUT.sc_to_tilt, 3, 3);
567  B = gsl_matrix_view_array( (double *) tiltm, 3, 3);
568  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
569  &A.matrix, &B.matrix, 0.0, sc2tiltp);
570  ptr_C = gsl_matrix_ptr(sc2tiltp, 0, 0);
571 
572  // sc2ocim = geo_lut.tilt_to_oci_mech#sc2tiltp
573  B = gsl_matrix_view_array( (double *) geoLUT.tilt_to_oci_mech, 3, 3);
574  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
575  sc2tiltp, &B.matrix, 0.0, sc2ocim);
576  ptr_C = gsl_matrix_ptr(sc2ocim, 0, 0);
577 
578  // sc_to_oci = geo_lut.oci_mech_to_oci_opt#sc2ocim
579  B = gsl_matrix_view_array( (double *) geoLUT.oci_mech_to_oci_opt, 3, 3);
580  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
581  sc2ocim, &B.matrix, 0.0, sc_to_oci);
582  ptr_C = gsl_matrix_ptr(sc_to_oci, 0, 0);
583 
584 
585  // Convert quaternion to matrix
586  double qmat[3][3];
587  qtom(quati[iscn], qmat);
588 
589  // smat = sc_to_oci#qmat
590  B = gsl_matrix_view_array(&qmat[0][0], 3, 3);
591  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0,
592  sc_to_oci, &B.matrix, 0.0, smat);
593 
594  // Compute attitude angles (informational only)
595  mat2rpy(posi[iscn], veli[iscn], qmat, rpy[iscn]);
596 
597  // Get scan ellipse coefficients
598  ptr_C = gsl_matrix_ptr(smat, 0, 0);
599  double coef[10];
600  scan_ell(posi[iscn], (double(*)[3])ptr_C, coef);
601 
602  // Generate pointing vector and relative time arrays in instrument frame
603  get_oci_vecs( nscan, pcdim, geoLUT.as_planarity, geoLUT.at_planarity,
604  geoLUT.rta_nadir, geoLUT.ham_ct_angles,
605  ev_toff, spin[iscn], hside[iscn], clines, deltc,
606  revpsec, ppr_off, board_id, nmcescan, mspin, enc_count,
607  &hamenc[0], &rtaenc[0], pview, thetap, iret);
608 
609  sfl[iscn] |= 4*iret;
610 
611  // Geolocate pixels
612  size_t ip = iscn * pcdim;
613  oci_geonav(dem_file.c_str(), posi[iscn], veli[iscn], (double(*)[3])ptr_C,
614  coef, sunr[iscn], pview, pcdim, delts, &xlat[ip], &xlon[ip],
615  &solz[ip], &sola[ip], &senz[ip], &sena[ip], &height[ip],
616  lonwest, loneast, latsouth, latnorth);
617  } else {
618  qfl[iscn] |= 4;
619  }
620  } // scan loop
621  gsl_matrix_free(smat);
622  gsl_matrix_free(sc2tiltp);
623  gsl_matrix_free(sc2ocim);
624  gsl_matrix_free(sc_to_oci);
625 
626  int geo_ncid;
627  int geo_gid[10];
628 
629  // Get number of bands
630  NcDim ntaps_dim = l1afile->getDim("number_of_taps");
631  uint32_t ntaps = ntaps_dim.getSize();
632 
633  int16_t *bagg = new int16_t[ntaps];
634  var = ncGrps[1].getVar( "blue_spectral_mode");
635  var.getVar( bagg);
636 
637  int16_t *ragg = new int16_t[ntaps];
638  var = ncGrps[1].getVar( "red_spectral_mode");
639  var.getVar( ragg);
640 
641  size_t *ia = new size_t[ntaps];
642  uint32_t ntb[16];
643  uint32_t bib, rib;
644  uint32_t bbb, rbb;
645  get_nib_nbb( ntaps, ia, ntb, bagg, bib, bbb);
646  get_nib_nbb( ntaps, ia, ntb, ragg, rib, rbb);
647 
648  static geoFile outfile;
649 
650  outfile.createFile(l1b_filename.c_str(),
651  // "$OCDATAROOT/oci/OCI_GEO_Data_Structure.cdl",
652  "$OCDATAROOT/oci/OCI_GEO_Level-1B_Data_Structure.cdl",
653  sdim, &geo_ncid, geo_gid, bbb, rbb,
654  pcdim, psdim, nswband, geoLUT.rta_nadir, radiance);
655 
656  string varname;
657 
658  /*
659  char buf[32];
660  strcpy(buf, unix2isodate(now(), 'G'));
661  nc_put_att_text(geo_ncid, NC_GLOBAL, "date_created", strlen(buf), buf);
662 
663  varname.assign("scan_time");
664  status = nc_inq_varid(geo_gid[scan_attr_geo], varname.c_str(), &varid);
665  status = nc_put_var_double(geo_gid[scan_attr_geo], varid, stime);
666  check_err(status, __LINE__, __FILE__);
667 
668  status = nc_put_att_text(geo_ncid, NC_GLOBAL,
669  "time_coverage_start", strlen(tstart), tstart);
670 
671  status = nc_put_att_text(geo_ncid, NC_GLOBAL,
672  "time_coverage_end", strlen(tend), tend);
673  */
674 
675 
676  vector<size_t> start, count;
677 
678  start.clear();
679  start.push_back(0);
680  start.push_back(0);
681  start.push_back(0);
682 
683  count.push_back(sdim);
684 
685  // Need to increment ncGrps index by 1 for new geo_l1b cdl file
686 
687  var = outfile.ncGrps[1].getVar( "time");
688  var.putVar( start, count, evtime);
689  var.putAtt("units", "seconds since " + granule_date);
690 
691  var = outfile.ncGrps[1].getVar( "HAM_side");
692  var.putVar( start, count, hside);
693 
694  var = outfile.ncGrps[1].getVar( "scan_quality_flags");
695  var.putVar( start, count, sfl);
696 
697  var = outfile.ncGrps[3].getVar( "tilt_angle");
698  var.putVar( start, count, tilt);
699 
700  count.push_back(4);
701 
702  var = outfile.ncGrps[3].getVar( "att_quat");
703  var.putVar( start, count, quati);
704 
705  count.pop_back();
706  count.push_back(3);
707 
708  var = outfile.ncGrps[3].getVar( "att_ang");
709  var.putVar( start, count, rpy);
710 
711  var = outfile.ncGrps[3].getVar( "orb_pos");
712  var.putVar( start, count, posi);
713 
714  var = outfile.ncGrps[3].getVar( "orb_vel");
715  var.putVar( start, count, veli);
716 
717  var = outfile.ncGrps[3].getVar( "sun_ref");
718  var.putVar( start, count, sunr);
719 
720  count.pop_back();
721  count.push_back(pcdim);
722 
723  var = outfile.ncGrps[2].getVar( "latitude");
724  var.putVar( start, count, xlat);
725 
726  var = outfile.ncGrps[2].getVar( "longitude");
727  var.putVar( start, count, xlon);
728 
729  var = outfile.ncGrps[2].getVar( "quality_flag");
730  var.putVar( start, count, qfl);
731 
732  var = outfile.ncGrps[2].getVar( "sensor_azimuth");
733  var.putVar( start, count, sena);
734 
735  var = outfile.ncGrps[2].getVar( "sensor_zenith");
736  var.putVar( start, count, senz);
737 
738  var = outfile.ncGrps[2].getVar( "solar_azimuth");
739  var.putVar( start, count, sola);
740 
741  var = outfile.ncGrps[2].getVar( "solar_zenith");
742  var.putVar( start, count, solz);
743 
744  // var = outfile.ncGrps[2].getVar( "range");
745  //var.putVar( start, count, range);
746 
747  var = outfile.ncGrps[2].getVar( "height");
748  var.putVar( start, count, height);
749 
750  // write global attributes, including history and date_created
751  set_global_attrs(outfile.geofile, history, doi, pversion);
752 
753  outfile.geofile->putAtt("earth_sun_distance_correction",
754  NC_DOUBLE, 1, &distcorr);
755 
756  outfile.geofile->putAtt( "cdm_data_type", "swath");
757 
758  outfile.geofile->putAtt("geospatial_lat_min", NC_FLOAT, 1, &latsouth);
759  outfile.geofile->putAtt("geospatial_lat_max", NC_FLOAT, 1, &latnorth);
760  outfile.geofile->putAtt("geospatial_lon_min", NC_FLOAT, 1, &lonwest);
761  outfile.geofile->putAtt("geospatial_lon_max", NC_FLOAT, 1, &loneast);
762 
763  outfile.geofile->putAtt( "geospatial_bounds_crs", "EPSG:4326");
764  /*
765 geospatial_bounds = "POLYGON ((85.129562 -78.832314,81.489693 49.646988...))"
766  */
767  outfile.close();
768 
769  delete[](att_time);
770  // delete[](orb_time);
771  delete[](quatr);
772  delete[](posr);
773  delete[](velr);
774  delete[](posi);
775  delete[](veli);
776  delete[](rpy);
777  delete[](quati);
778  delete[](tilt);
779  delete[](xlon);
780  delete[](xlat);
781  delete[](solz);
782  delete[](sola);
783  delete[](senz);
784  delete[](sena);
785  // delete[](range);
786  delete[](height);
787  delete[](tiltin);
788  delete[](ttime);
789  delete[](qfl);
790  delete[]pview[0]; delete[](pview);
791  delete[](sunr);
792 
793  delete [] hamenc[0]; delete [] hamenc;
794  delete [] rtaenc[0]; delete [] rtaenc;
795 
796  delete [](sstime);
797  delete [](evtime);
798  delete [](thetap);
799  delete [](spin);
800  delete [](mspin);
801  delete [](ot_10us);
802  delete []att_quat[0]; delete [](att_quat);
803  delete []orb_vel[0]; delete [](orb_vel);
804  delete []orb_pos[0]; delete[](orb_pos);
805  delete [](sfl);
806  delete [](hside);
807  delete [](enc_count);
808  delete [](ragg);
809  delete [](bagg);
810  delete [](iagg);
811  delete [](dtype);
812  delete [](lines);
813  delete(l1afile);
814  delete(geoLUTfile);
815  // clo_deleteList(optionList);
816  return 0;
817 }
818 
819 int j2000_to_ecr(int32_t iyr, int32_t idy, double sec, double ecmat[3][3]) {
820  // Get J2000 to ECEF transformation matrix
821 
822  // Arguments:
823 
824  // Name Type I/O Description
825  // --------------------------------------------------------
826  // iyr I I Year, four digits
827  // idy I I Day of year
828  // sec R*8 I Seconds of day
829  // ecmat(3,3) R O J2000 to ECEF matrix
830 
831  // Get transformation from J2000 to mean-of-date inertial
832  double j2mod[3][3];
833  j2000_to_mod(iyr, idy, sec, j2mod);
834 
835  // Get nutation and UT1-UTC (once per run)
836  double xnut[3][3], ut1utc;
837  get_nut(iyr, idy, xnut);
838  get_ut1(iyr, idy, ut1utc);
839 
840  // Compute Greenwich hour angle for time of day
841  double day = idy + (sec + ut1utc) / 86400;
842  double gha, gham[3][3];
843  gha2000(iyr, day, gha);
844 
845  gham[0][0] = cos(gha / RADEG);
846  gham[1][1] = cos(gha / RADEG);
847  gham[2][2] = 1;
848  gham[0][1] = sin(gha / RADEG);
849  gham[1][0] = -sin(gha / RADEG);
850 
851  gham[0][2] = 0;
852  gham[2][0] = 0;
853  gham[1][2] = 0;
854  gham[2][1] = 0;
855 
856  // Combine all transformations
857  gsl_matrix_view A = gsl_matrix_view_array(&gham[0][0], 3, 3); // gham
858  gsl_matrix_view B = gsl_matrix_view_array(&xnut[0][0], 3, 3); // xnut
859  gsl_matrix *C = gsl_matrix_alloc(3, 3);
860 
861  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, &A.matrix, &B.matrix, 0.0, C);
862 
863  gsl_matrix_view D = gsl_matrix_view_array(&j2mod[0][0], 3, 3); // j2mod
864  gsl_matrix *E = gsl_matrix_alloc(3, 3);
865  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, C, &D.matrix, 0.0, E);
866  double *ptr_E = gsl_matrix_ptr(E, 0, 0);
867 
868  memcpy(ecmat, ptr_E, 9 * sizeof(double));
869 
870  gsl_matrix_free(C);
871  gsl_matrix_free(E);
872 
873  return 0;
874 }
875 
876 int j2000_to_mod(int32_t iyr, int32_t idy, double sec, double j2mod[3][3]) {
877  // Get J2000 to MOD (precession) transformation
878 
879  // Arguments:
880 
881  // Name Type I/O Description
882  // --------------------------------------------------------
883  // iyr I I Year, four digits
884  // idy I I Day of year
885  // sec R*8 I Seconds of day
886  // j2mod(3,3) R O J2000 to MOD matrix
887 
888  int16_t iyear = iyr;
889  int16_t iday = idy;
890 
891  double t = jday(iyear, 1, iday) - (double)2451545.5 + sec / 86400;
892  t /= 36525;
893 
894  double zeta0 = t * (2306.2181 + 0.302 * t + 0.018 * t * t) / RADEG / 3600;
895  double thetap = t * (2004.3109 - 0.4266 * t - 0.04160 * t * t) / RADEG / 3600;
896  double xip = t * (2306.2181 + 1.095 * t + 0.018 * t * t) / RADEG / 3600;
897 
898  j2mod[0][0] = -sin(zeta0) * sin(xip) + cos(zeta0) * cos(xip) * cos(thetap);
899  j2mod[0][1] = -cos(zeta0) * sin(xip) - sin(zeta0) * cos(xip) * cos(thetap);
900  j2mod[0][2] = -cos(xip) * sin(thetap);
901  j2mod[1][0] = sin(zeta0) * cos(xip) + cos(zeta0) * sin(xip) * cos(thetap);
902  j2mod[1][1] = cos(zeta0) * cos(xip) - sin(zeta0) * sin(xip) * cos(thetap);
903  j2mod[1][2] = -sin(xip) * sin(thetap);
904  j2mod[2][0] = cos(zeta0) * sin(thetap);
905  j2mod[2][1] = -sin(zeta0) * sin(thetap);
906  j2mod[2][2] = cos(thetap);
907 
908  return 0;
909 }
910 
911 int get_nut(int32_t iyr, int32_t idy, double xnut[3][3]) {
912  int16_t iyear = iyr;
913  int16_t iday = idy;
914 
915  double t = jday(iyear, 1, iday) - (double)2451545.5;
916 
917  double xls, gs, xlm, omega;
918  double dpsi, eps, epsm;
919  ephparms(t, xls, gs, xlm, omega);
920  nutate(t, xls, gs, xlm, omega, dpsi, eps, epsm);
921 
922  xnut[0][0] = cos(dpsi / RADEG);
923  xnut[1][0] = -sin(dpsi / RADEG) * cos(epsm / RADEG);
924  xnut[2][0] = -sin(dpsi / RADEG) * sin(epsm / RADEG);
925  xnut[0][1] = sin(dpsi / RADEG) * cos(eps / RADEG);
926  xnut[1][1] = cos(dpsi / RADEG) * cos(eps / RADEG) * cos(epsm / RADEG) +
927  sin(eps / RADEG) * sin(epsm / RADEG);
928  xnut[2][1] = cos(dpsi / RADEG) * cos(eps / RADEG) * sin(epsm / RADEG) -
929  sin(eps / RADEG) * cos(epsm / RADEG);
930  xnut[0][2] = sin(dpsi / RADEG) * sin(eps / RADEG);
931  xnut[1][2] = cos(dpsi / RADEG) * sin(eps / RADEG) * cos(epsm / RADEG) -
932  cos(eps / RADEG) * sin(epsm / RADEG);
933  xnut[2][2] = cos(dpsi / RADEG) * sin(eps / RADEG) * sin(epsm / RADEG) +
934  cos(eps / RADEG) * cos(epsm / RADEG);
935 
936  return 0;
937 }
938 
939 int ephparms(double t, double &xls, double &gs, double &xlm, double &omega) {
940  // This subroutine computes ephemeris parameters used by other Mission
941  // Operations routines: the solar mean longitude and mean anomaly, and
942  // the lunar mean longitude and mean ascending node. It uses the model
943  // referenced in The Astronomical Almanac for 1984, Section S
944  // (Supplement) and documented for the SeaWiFS Project in "Constants
945  // and Parameters for SeaWiFS Mission Operations", in TBD. These
946  // parameters are used to compute the solar longitude and the nutation
947  // in longitude and obliquity.
948 
949  // Sun Mean Longitude
950  xls = (double)280.46592 + ((double)0.9856473516) * t;
951  xls = fmod(xls, (double)360);
952 
953  // Sun Mean Anomaly
954  gs = (double)357.52772 + ((double)0.9856002831) * t;
955  gs = fmod(gs, (double)360);
956 
957  // Moon Mean Longitude
958  xlm = (double)218.31643 + ((double)13.17639648) * t;
959  xlm = fmod(xlm, (double)360);
960 
961  // Ascending Node of Moon's Mean Orbit
962  omega = (double)125.04452 - ((double)0.0529537648) * t;
963  omega = fmod(omega, (double)360);
964 
965  return 0;
966 }
967 
968 int nutate(double t, double xls, double gs, double xlm, double omega,
969  double &dpsi, double &eps, double &epsm) {
970  // This subroutine computes the nutation in longitude and the obliquity
971  // of the ecliptic corrected for nutation. It uses the model referenced
972  // in The Astronomical Almanac for 1984, Section S (Supplement) and
973  // documented for the SeaWiFS Project in "Constants and Parameters for
974  // SeaWiFS Mission Operations", in TBD. These parameters are used to
975  // compute the apparent time correction to the Greenwich Hour Angle and
976  // for the calculation of the geocentric Sun vector. The input
977  // ephemeris parameters are computed using subroutine ephparms. Terms
978  // are included to 0.1 arcsecond.
979 
980  // Nutation in Longitude
981  dpsi = -17.1996 * sin(omega / RADEG) + 0.2062 * sin(2. * omega / RADEG) -
982  1.3187 * sin(2. * xls / RADEG) + 0.1426 * sin(gs / RADEG) -
983  0.2274 * sin(2. * xlm / RADEG);
984 
985  // Mean Obliquity of the Ecliptic
986  epsm = (double)23.439291 - ((double)3.560e-7) * t;
987 
988  // Nutation in Obliquity
989  double deps = 9.2025 * cos(omega / RADEG) + 0.5736 * cos(2. * xls / RADEG);
990 
991  // True Obliquity of the Ecliptic
992  eps = epsm + deps / 3600;
993 
994  dpsi = dpsi / 3600;
995 
996  return 0;
997 }
998 
999 int get_ut1(int32_t iyr, int32_t idy, double &ut1utc) {
1000  int16_t iyear = iyr;
1001  int16_t iday = idy;
1002 
1003  static int32_t ijd[25000];
1004  static double ut1[25000];
1005  string utcpole = "$OCVARROOT/modis/utcpole.dat";
1006  static bool first = true;
1007  int k = 0;
1008  if (first) {
1009  string line;
1010  expandEnvVar(&utcpole);
1011  istringstream istr;
1012 
1013  ifstream utcpfile(utcpole.c_str());
1014  if (utcpfile.is_open()) {
1015  getline(utcpfile, line);
1016  getline(utcpfile, line);
1017  while (getline(utcpfile, line)) {
1018  istr.clear();
1019  istr.str(line.substr(0, 5));
1020  istr >> ijd[k];
1021  istr.clear();
1022  istr.str(line.substr(44, 9));
1023  istr >> ut1[k];
1024  k++;
1025  }
1026  ijd[k] = 0;
1027  utcpfile.close();
1028  first = false;
1029  } else {
1030  cout << utcpole.c_str() << " not found" << endl;
1031  exit(1);
1032  }
1033  } // if (first)
1034 
1035  k = 0;
1036  int mjd = jday(iyear, 1, iday) - 2400000;
1037  while (ijd[k] > 0) {
1038  if (mjd == ijd[k]) {
1039  ut1utc = ut1[k];
1040  break;
1041  }
1042  k++;
1043  }
1044 
1045  return 0;
1046 }
1047 
1048 int gha2000(int32_t iyr, double day, double &gha) {
1049  // This subroutine computes the Greenwich hour angle in degrees for the
1050  // input time. It uses the model referenced in The Astronomical Almanac
1051  // for 1984, Section S (Supplement) and documented for the SeaWiFS
1052  // Project in "Constants and Parameters for SeaWiFS Mission Operations",
1053  // in TBD. It includes the correction to mean sidereal time for nutation
1054  // as well as precession.
1055  //
1056 
1057  // Compute days since J2000
1058  int16_t iday = day;
1059  double fday = day - iday;
1060  int jd = jday(iyr, 1, iday);
1061  double t = jd - (double)2451545.5 + fday;
1062 
1063  // Compute Greenwich Mean Sidereal Time (degrees)
1064  double gmst = (double)100.4606184 + (double)0.9856473663 * t +
1065  (double)2.908e-13 * t * t;
1066 
1067  // Check if need to compute nutation correction for this day
1068  double xls, gs, xlm, omega;
1069  double dpsi, eps, epsm;
1070  ephparms(t, xls, gs, xlm, omega);
1071  nutate(t, xls, gs, xlm, omega, dpsi, eps, epsm);
1072 
1073  // Include apparent time correction and time-of-day
1074  gha = gmst + dpsi * cos(eps / RADEG) + fday * 360;
1075  gha = fmod(gha, (double)360);
1076 
1077  return 0;
1078 }
1079 
1080 int mtoq(double rm[3][3], double q[4]) {
1081  // Convert direction cosine matrix to equivalent quaternion
1082 
1083  double e[3];
1084 
1085  // Compute Euler angle
1086  double phi;
1087  double cphi = (rm[0][0] + rm[1][1] + rm[2][2] - 1) / 2;
1088  if (fabs(cphi) < 0.98) {
1089  phi = acos(cphi);
1090  } else {
1091  double ssphi = (pow(rm[0][1] - rm[1][0], 2) +
1092  pow(rm[2][0] - rm[0][2], 2) +
1093  pow(rm[1][2] - rm[2][1], 2)) /
1094  4;
1095  phi = asin(sqrt(ssphi));
1096  if (cphi < 0) phi = PI - phi;
1097  }
1098 
1099  // Compute Euler axis
1100  e[0] = (rm[2][1] - rm[1][2]) / (sin(phi) * 2);
1101  e[1] = (rm[0][2] - rm[2][0]) / (sin(phi) * 2);
1102  e[2] = (rm[1][0] - rm[0][1]) / (sin(phi) * 2);
1103  double norm = sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);
1104  e[0] = e[0] / norm;
1105  e[1] = e[1] / norm;
1106  e[2] = e[2] / norm;
1107 
1108  // Compute quaternion
1109  q[0] = e[0] * sin(phi / 2);
1110  q[1] = e[1] * sin(phi / 2);
1111  q[2] = e[2] * sin(phi / 2);
1112  q[3] = cos(phi / 2);
1113 
1114  return 0;
1115 }
1116 
1117 int qprod(double q1[4], float q2[4], double q3[4]) {
1118  // Compute the product of two quaternions q3 = q1*q2
1119 
1120  q3[0] = q1[0] * q2[3] + q1[1] * q2[2] - q1[2] * q2[1] + q1[3] * q2[0];
1121  q3[1] = -q1[0] * q2[2] + q1[1] * q2[3] + q1[2] * q2[0] + q1[3] * q2[1];
1122  q3[2] = q1[0] * q2[1] - q1[1] * q2[0] + q1[2] * q2[3] + q1[3] * q2[2];
1123  q3[3] = -q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] + q1[3] * q2[3];
1124 
1125  return 0;
1126 }
1127 
1128 int qprod(float q1[4], float q2[4], float q3[4]) {
1129  // Compute the product of two quaternions q3 = q1*q2
1130 
1131  q3[0] = q1[0] * q2[3] + q1[1] * q2[2] - q1[2] * q2[1] + q1[3] * q2[0];
1132  q3[1] = -q1[0] * q2[2] + q1[1] * q2[3] + q1[2] * q2[0] + q1[3] * q2[1];
1133  q3[2] = q1[0] * q2[1] - q1[1] * q2[0] + q1[2] * q2[3] + q1[3] * q2[2];
1134  q3[3] = -q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] + q1[3] * q2[3];
1135 
1136  return 0;
1137 }
1138 
1139 int orb_interp(size_t n_orb_rec, size_t sdim, double *torb, orb_array *p,
1140  orb_array *v, double *time, orb_array *posi, orb_array *veli) {
1141  // Purpose: Interpolate orbit position and velocity vectors to scan line
1142  // times
1143  //
1144  //
1145  // Calling Arguments:
1146  //
1147  // Name Type I/O Description
1148  // -------- ---- --- -----------
1149  // torb(*) double I Array of orbit vector times in seconds of day
1150  // p(3,*) float I Array of orbit position vectors for
1151  // each time torb
1152  // v(3,*) float I Array of orbit velocity vectors for
1153  // each time torb
1154  // time(*) double I Array of time in seconds of day
1155  // for every scan line
1156  // pi(3,*) float O Array of interpolated positions
1157  // vi(3,*) float O Array of interpolated velocities
1158  //
1159  //
1160  // By: Frederick S. Patt, GSC, August 10, 1993
1161  //
1162  // Notes: Method uses cubic polynomial to match positions and velocities
1163  // at input data points.
1164  //
1165  // Modification History:
1166  //
1167  // Created IDL version from FORTRAN code. F.S. Patt, SAIC, November 29, 2006
1168  //
1169 
1170  double a0[3], a1[3], a2[3], a3[3];
1171 
1172  // initialize the arrays
1173  for (int i=0; i<3; i++) {
1174  a0[i] = 0.0;
1175  a1[i] = 0.0;
1176  a2[i] = 0.0;
1177  a3[i] = 0.0;
1178  }
1179 
1180  /*
1181  // Make sure that first orbit vector precedes first scan line
1182  k = where (time lt torb(0))
1183  if (k(0) ne -1) then begin
1184  posi(*,k) = 0.0
1185  veli(*,k) = 0.0
1186  orbfl(k) = 1
1187  print, 'Scan line times before available orbit data'
1188  i1 = max(k) + 1
1189  endif
1190  */
1191 
1192  // orbit time
1193  double startOrbTime = torb[0];
1194  double endOrbTime = torb[n_orb_rec-1];
1195  const double NAV_FILL = -9999999.0;
1196 
1197  size_t ind = 0;
1198  for (size_t i = 0; i < sdim; i++) {
1199 
1200  // scantime outside valid orbit time, use fill value
1201  if (time[i] < startOrbTime || time[i] > endOrbTime) {
1202  for (int vector=0; vector<3; vector++) {
1203  posi[i][vector] = NAV_FILL;
1204  veli[i][vector] = BAD_FLT;
1205  }
1206  // skip calulations below
1207  continue;
1208  }
1209 
1210  // Find input orbit vectors bracketing scan
1211  for (size_t j = ind; j < n_orb_rec; j++) {
1212  if (time[i] > torb[j] && time[i] <= torb[j + 1]) {
1213  ind = j;
1214  break;
1215  }
1216  }
1217 
1218  // Set up cubic interpolation
1219  double dt = torb[ind + 1] - torb[ind];
1220  for (size_t j = 0; j < 3; j++) {
1221  a0[j] = p[ind][j];
1222  a1[j] = v[ind][j] * dt;
1223  if (dt >= 3.5) {
1224  a2[j] = 3 * p[ind + 1][j] - 3 * p[ind][j] - 2 * v[ind][j] * dt -
1225  v[ind + 1][j] * dt;
1226  a3[j] = 2 * p[ind][j] - 2 * p[ind + 1][j] + v[ind][j] * dt +
1227  v[ind + 1][j] * dt;
1228  } else {
1229  a2[j] = (v[ind + 1][j] - v[ind][j]) / 2;
1230  a3[j] = 0.0;
1231  // a2(*) = (v(*,ind+1)-v(*,ind))*dt/2
1232  // a3(*) = 0.d0
1233  }
1234  }
1235 
1236  // Interpolate orbit position and velocity components to scan line time
1237  double x = (time[i] - torb[ind]) / dt;
1238  double x2 = x * x;
1239  double x3 = x2 * x;
1240  for (size_t j = 0; j < 3; j++) {
1241  posi[i][j] = a0[j] + a1[j] * x + a2[j] * x2 + a3[j] * x3;
1242  veli[i][j] = (a1[j] + 2 * a2[j] * x + 3 * a3[j] * x2) / dt;
1243  }
1244  } // i-loop
1245 
1246  return 0;
1247 }
1248 
1249 int q_interp(size_t n_att_rec, size_t sdim, double *tq, quat_array *q,
1250  double *time, quat_array *qi) {
1251  // Purpose: Interpolate quaternions to scan line times
1252  //
1253  //
1254 
1255  // Attitude
1256  double startAttTime = tq[0];
1257  double endAttTime = tq[n_att_rec-1];
1258  const double NAV_FILL = -32767.0;
1259  size_t ind = 0;
1260  for (size_t i = 0; i < sdim; i++) {
1261 
1262  // if scan time is outside the valid att time, use fill value
1263  if (time[i] < startAttTime || time[i] > endAttTime) {
1264  for (int element=0; element<4; element++) {
1265  qi[i][element] = NAV_FILL;
1266 
1267  }
1268  continue;
1269  }
1270 
1271  // Find input attitude vectors bracketing scan
1272  for (size_t j = ind; j < n_att_rec; j++) {
1273  if (time[i] > tq[j] && time[i] <= tq[j + 1]) {
1274  ind = j;
1275  break;
1276  }
1277  }
1278 
1279  // Set up quaternion interpolation
1280  double dt = tq[ind + 1] - tq[ind];
1281  double qin[4];
1282  qin[0] = -q[ind][0];
1283  qin[1] = -q[ind][1];
1284  qin[2] = -q[ind][2];
1285  qin[3] = q[ind][3];
1286 
1287  double e[3], qr[4];
1288  qprod(qin, q[ind + 1], qr);
1289  memcpy(e, qr, 3 * sizeof(double));
1290  double sto2 = sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);
1291  for (size_t j = 0; j < 3; j++) e[j] /= sto2;
1292 
1293  // Interpolate quaternion to scan times
1294  double x = (time[i] - tq[ind]) / dt;
1295  float qri[4], qp[4];
1296  for (size_t j = 0; j < 3; j++) qri[j] = e[j] * sto2 * x;
1297  qri[3] = 1.0;
1298  qprod(q[ind], qri, qp);
1299  memcpy(qi[i], qp, 4 * sizeof(float));
1300  }
1301 
1302  return 0;
1303 }
1304 
1305 
1306 int tilt_interp(size_t n_tilts, size_t sdim, double *ttilt, float *tiltin,
1307  double *time, float *tilt) {
1308 
1309  double startTiltTime = ttilt[0];
1310  double endTiltTime = ttilt[n_tilts-1];
1311  const double NAV_FILL = -32767.0;
1312 
1313  size_t ind = 0;
1314  for (size_t i = 0; i < sdim; i++) {
1315 
1316  if (time[i] < startTiltTime || time[i] > endTiltTime) {
1317  tilt[i] = NAV_FILL;
1318  continue;
1319  }
1320 
1321  // Find input tilt vectors bracketing scan
1322  for (size_t j = ind; j < n_tilts; j++) {
1323  if (time[i] > ttilt[j] && time[i] <= ttilt[j + 1]) {
1324  ind = j;
1325  break;
1326  }
1327  }
1328 
1329  double x = (time[i] - ttilt[ind])/(ttilt[ind+1] - ttilt[ind]);
1330  tilt[i] = (1-x)*tiltin[ind] + x*tiltin[ind+1];
1331  }
1332 
1333  return 0;
1334 }
1335 
1336 
1337 int l_sun(size_t sdim, int32_t iyr, int32_t iday, double *sec,
1338  orb_array *sunr, double *rs) {
1339  // Computes unit Sun vector in geocentric rotating coordinates, using
1340  // subprograms to compute inertial Sun vector and Greenwich hour angle
1341 
1342  // Get unit Sun vector in geocentric inertial coordinates
1343  sun2000(sdim, iyr, iday, sec, sunr, rs);
1344 
1345  // Get Greenwich mean sidereal angle
1346  for (size_t i = 0; i < sdim; i++) {
1347  double day = iday + sec[i] / 86400;
1348  double gha;
1349  gha2000(iyr, day, gha);
1350  double ghar = gha / RADEG;
1351 
1352  // Transform Sun vector into geocentric rotating frame
1353  float temp0 = sunr[i][0] * cos(ghar) + sunr[i][1] * sin(ghar);
1354  float temp1 = sunr[i][1] * cos(ghar) - sunr[i][0] * sin(ghar);
1355  sunr[i][0] = temp0;
1356  sunr[i][1] = temp1;
1357  }
1358 
1359  return 0;
1360 }
1361 
1362 int sun2000(size_t sdim, int32_t iyr, int32_t idy, double *sec,
1363  orb_array *sun, double *rs) {
1364  // This subroutine computes the Sun vector in geocentric inertial
1365  // (equatorial) coordinates. It uses the model referenced in The
1366  // Astronomical Almanac for 1984, Section S (Supplement) and documented
1367  // for the SeaWiFS Project in "Constants and Parameters for SeaWiFS
1368  // Mission Operations", in TBD. The accuracy of the Sun vector is
1369  // approximately 0.1 arcminute.
1370 
1371  float xk = 0.0056932; // Constant of aberration
1372 
1373  // Compute floating point days since Jan 1.5, 2000
1374  // Note that the Julian day starts at noon on the specified date
1375  int16_t iyear = iyr;
1376  int16_t iday = idy;
1377 
1378  for (size_t i = 0; i < sdim; i++) {
1379  double t =
1380  jday(iyear, 1, iday) - (double)2451545.0 + (sec[i] - 43200) / 86400;
1381 
1382  double xls, gs, xlm, omega;
1383  double dpsi, eps, epsm;
1384  // Compute solar ephemeris parameters
1385  ephparms(t, xls, gs, xlm, omega);
1386 
1387  // Compute nutation corrections for this day
1388  nutate(t, xls, gs, xlm, omega, dpsi, eps, epsm);
1389 
1390  // Compute planet mean anomalies
1391  // Venus Mean Anomaly
1392  double g2 = 50.40828 + 1.60213022 * t;
1393  g2 = fmod(g2, (double)360);
1394 
1395  // Mars Mean Anomaly
1396  double g4 = 19.38816 + 0.52402078 * t;
1397  g4 = fmod(g4, (double)360);
1398 
1399  // Jupiter Mean Anomaly
1400  double g5 = 20.35116 + 0.08309121 * t;
1401  g5 = fmod(g5, (double)360);
1402 
1403  // Compute solar distance (AU)
1404  rs[i]
1405  = 1.00014 - 0.01671 * cos(gs / RADEG) - 0.00014 * cos(2. * gs / RADEG);
1406 
1407  // Compute Geometric Solar Longitude
1408  double dls = (6893. - 4.6543463e-4 * t) * sin(gs / RADEG) +
1409  72. * sin(2. * gs / RADEG) - 7. * cos((gs - g5) / RADEG) +
1410  6. * sin((xlm - xls) / RADEG) +
1411  5. * sin((4. * gs - 8. * g4 + 3. * g5) / RADEG) -
1412  5. * cos((2. * gs - 2. * g2) / RADEG) -
1413  4. * sin((gs - g2) / RADEG) +
1414  4. * cos((4. * gs - 8. * g4 + 3. * g5) / RADEG) +
1415  3. * sin((2. * gs - 2. * g2) / RADEG) - 3. * sin(g5 / RADEG) -
1416  3. * sin((2. * gs - 2. * g5) / RADEG);
1417 
1418  double xlsg = xls + dls / 3600;
1419 
1420  // Compute Apparent Solar Longitude// includes corrections for nutation
1421  // in longitude and velocity aberration
1422  double xlsa = xlsg + dpsi - xk / rs[i];
1423 
1424  // Compute unit Sun vector
1425  sun[i][0] = cos(xlsa / RADEG);
1426  sun[i][1] = sin(xlsa / RADEG) * cos(eps / RADEG);
1427  sun[i][2] = sin(xlsa / RADEG) * sin(eps / RADEG);
1428  } // i-loop
1429 
1430  return 0;
1431 }
1432 
1433 int qtom(float q[4], double rm[3][3]) {
1434  // Convert quaternion to equivalent direction cosine matrix
1435 
1436  rm[0][0] = q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
1437  rm[0][1] = 2 * (q[0] * q[1] + q[2] * q[3]);
1438  rm[0][2] = 2 * (q[0] * q[2] - q[1] * q[3]);
1439  rm[1][0] = 2 * (q[0] * q[1] - q[2] * q[3]);
1440  rm[1][1] = -q[0] * q[0] + q[1] * q[1] - q[2] * q[2] + q[3] * q[3];
1441  rm[1][2] = 2 * (q[1] * q[2] + q[0] * q[3]);
1442  rm[2][0] = 2 * (q[0] * q[2] + q[1] * q[3]);
1443  rm[2][1] = 2 * (q[1] * q[2] - q[0] * q[3]);
1444  rm[2][2] = -q[0] * q[0] - q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
1445 
1446  return 0;
1447 }
1448 
1449 int scan_ell(float p[3], double sm[3][3], double coef[10]) {
1450  // This program calculates the coefficients which
1451  // represent the Earth scan track in the sensor frame.
1452 
1453  // The reference ellipsoid uses an equatorial radius of 6378.137 km and
1454  // a flattening factor of 1/298.257 (WGS 1984).
1455 
1456  // Calling Arguments
1457  //
1458  // Name Type I/O Description
1459  //
1460  // pos(3) R*4 I ECR Orbit Position Vector (km)
1461  // smat(3,3) R*4 I Sensor Orientation Matrix
1462  // coef(10) R*4 O Scan path coefficients
1463 
1464  double re = 6378.137;
1465  double f = 1 / 298.257;
1466  double omf2 = (1 - f) * (1 - f);
1467 
1468  // Compute constants for navigation model using Earth radius values
1469  double rd = 1 / omf2;
1470 
1471  // Compute coefficients of intersection ellipse in scan plane
1472  coef[0] = 1 + (rd - 1) * sm[0][2] * sm[0][2];
1473  coef[1] = 1 + (rd - 1) * sm[1][2] * sm[1][2];
1474  coef[2] = 1 + (rd - 1) * sm[2][2] * sm[2][2];
1475  coef[3] = (rd - 1) * sm[0][2] * sm[1][2] * 2;
1476  coef[4] = (rd - 1) * sm[0][2] * sm[2][2] * 2;
1477  coef[5] = (rd - 1) * sm[1][2] * sm[2][2] * 2;
1478  coef[6] = (sm[0][0] * p[0] + sm[0][1] * p[1] + sm[0][2] * p[2] * rd) * 2;
1479  coef[7] = (sm[1][0] * p[0] + sm[1][1] * p[1] + sm[1][2] * p[2] * rd) * 2;
1480  coef[8] = (sm[2][0] * p[0] + sm[2][1] * p[1] + sm[2][2] * p[2] * rd) * 2;
1481  coef[9] = p[0] * p[0] + p[1] * p[1] + p[2] * p[2] * rd - re * re;
1482 
1483  return 0;
1484 }
1485 
1486 int oci_geonav(const char* dem_file, float pos[3], float vel[3],
1487  double smat[3][3], double coef[10], float sunr[3],
1488  float **pview, size_t npix, double *delt, float *xlat,
1489  float *xlon, short *solz, short *sola,
1490  short *senz, short *sena, short *height,
1491  float& lonwest, float& loneast, float& latsouth, float& latnorth)
1492  {
1493  // This subroutine performs navigation of a scanning sensor on the
1494  // surface of an ellipsoid based on an input orbit position vector and
1495  // spacecraft orientation matrix. It uses a closed-form algorithm for
1496  // determining points on the ellipsoidal surface which involves
1497  // determining the intersection of the scan plan with the ellipsoid.
1498  // The sensor view vectors in the sensor frame are passed in as a 3xN array.
1499 
1500  // The reference ellipsoid is set according to the scan
1501  // intersection coefficients in the calling sequence// an equatorial
1502  // radius of 6378.137 km. and a flattening factor of 1/298.257 are
1503  // used by both the Geodetic Reference System (GRS) 1980 and the
1504  // World Geodetic System (WGS) 1984.
1505 
1506  // It then computes geometric parameters using the pixel locations on
1507  // the Earth, the spacecraft position vector and the unit Sun vector in
1508  // the geocentric rotating reference frame. The outputs are arrays of
1509  // geodetic latitude and longitude, solar zenith and azimuth and sensor
1510  // zenith and azimuth. The azimuth angles are measured from local
1511  // North toward East. Flag values of 999. are returned for any pixels
1512  // whose scan angle is past the Earth's horizon.
1513 
1514  // Reference: "Exact closed-form geolocation geolocation algorithm for
1515  // Earth survey sensors", F. S. Patt and W. W. Gregg, IJRS, Vol. 15
1516  // No. 18, 1994.
1517 
1518  // Calling Arguments
1519 
1520  // Name Type I/O Description
1521  // dem_file char* I Digital elevation map file
1522  // pos(3) R*4 I ECR Orbit Position Vector (km) at scan
1523  // mid-time
1524  // vel(3) R*4 I ECR Orbit Velocity Vector (km/sec)
1525  // smat(3,3) R*4 I Sensor Orientation Matrix
1526  // coef(10) R*4 I Scan path coefficients
1527  // sunr(3) R*4 I Sun unit vector in geocentric rotating
1528  // reference frame
1529  // pview(3,*) R*4 I Array of sensor view vectors
1530  // npix R*4 I Number of pixels to geolocate
1531  // xlat(*) R*4 O Pixel geodetic latitudes
1532  // xlon(*) R*4 O Pixel geodetic longitudes
1533  // solz(*) I*2 O Pixel solar zenith angles
1534  // sola(*) I*2 O Pixel solar azimuth angles
1535  // senz(*) I*2 O Pixel sensor zenith angles
1536  // sena(*) I*2 O Pixel sensor azimuth angles
1537  // height(*) I*2 O Pixel terrain height
1538  //
1539 
1540  // Program written by: Frederick S. Patt
1541  // General Sciences Corporation
1542  // October 20, 1992
1543  //
1544  // Modification History:
1545 
1546  // Created universal version of the SeaWiFS geolocation algorithm
1547  // by specifying view vectors as an input. F. S. Patt, SAIC, 11/17/09
1548 
1549  // Earth ellipsoid parameters
1550  float f = 1 / 298.257;
1551  float omf2 = (1 - f) * (1 - f);
1552  gsl_vector *C = gsl_vector_alloc(3);
1553 
1554  for (size_t i = 0; i < npix; i++) {
1555  // Compute sensor-to-surface vectors for all scan angles
1556  // Compute terms for quadratic equation
1557  double o =
1558  coef[0] * pview[i][0] * pview[i][0] +
1559  coef[1] * pview[i][1] * pview[i][1] +
1560  coef[2] * pview[i][2] * pview[i][2] +
1561  coef[3] * pview[i][0] * pview[i][1] +
1562  coef[4] * pview[i][0] * pview[i][2] +
1563  coef[5] * pview[i][1] * pview[i][2];
1564 
1565  double p =
1566  coef[6] * pview[i][0] + coef[7] * pview[i][1] + coef[8] * pview[i][2];
1567 
1568  double q = coef[9];
1569 
1570  double r = p * p - 4 * q * o;
1571 
1572  xlat[i] = -32767;
1573  xlon[i] = -32767;
1574 
1575  solz[i] = -32767;
1576  sola[i] = -32767;
1577  senz[i] = -32767;
1578  sena[i] = -32767;
1579  height[i] = -32767;
1580  // Check for scan past edge of Earth
1581  if (r >= 0) {
1582 
1583  // Solve for magnitude of sensor-to-pixel vector and compute components
1584  double d = (-p - sqrt(r)) / (2 * o);
1585  double x1[3];
1586  for (size_t j = 0; j < 3; j++) x1[j] = d * pview[i][j];
1587 
1588  // Convert velocity vector to ground speed
1589  float re = 6378.137;
1590  // double omegae = 7.29211585494e-05;
1591  double pm = sqrt(pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2]);
1592  double clatg = sqrt(pos[0]*pos[0] + pos[1]*pos[1]) / pm;
1593  double rg = re*(1.-f)/sqrt(1.-(2.-f)*f*clatg*clatg);
1594  double v[3];
1595  // v[0] = (vel[0] - pos[1]*omegae) * rg / pm;
1596  //v[1] = (vel[1] - pos[0]*omegae) * rg / pm;
1597  v[0] = vel[0] * rg / pm;
1598  v[1] = vel[1] * rg / pm;
1599  v[2] = vel[2] * rg / pm;
1600 
1601  // Transform vector from sensor to geocentric frame
1602  gsl_matrix_view A = gsl_matrix_view_array((double *)smat, 3, 3);
1603  gsl_vector_view B = gsl_vector_view_array(x1, 3);
1604 
1605  gsl_blas_dgemv(CblasTrans, 1.0, &A.matrix, &B.vector, 0.0, C);
1606 
1607  float rh[3], geovec[3];
1608  double *ptr_C = gsl_vector_ptr(C, 0);
1609  for (size_t j = 0; j < 3; j++) {
1610  rh[j] = ptr_C[j];
1611  geovec[j] = pos[j] + rh[j] + v[j]*delt[i];
1612  }
1613 
1614  // Compute the local vertical, East and North unit vectors
1615  float uxy = geovec[0] * geovec[0] + geovec[1] * geovec[1];
1616  float temp = sqrt(geovec[2] * geovec[2] + omf2 * omf2 * uxy);
1617 
1618  float up[3];
1619  up[0] = omf2 * geovec[0] / temp;
1620  up[1] = omf2 * geovec[1] / temp;
1621  up[2] = geovec[2] / temp;
1622  float upxy = sqrt(up[0] * up[0] + up[1] * up[1]);
1623 
1624  float ea[3];
1625  ea[0] = -up[1] / upxy;
1626  ea[1] = up[0] / upxy;
1627  ea[2] = 0.0;
1628 
1629  // no = crossp(up,ea)
1630  float no[3];
1631  no[0] = -up[2] * ea[1];
1632  no[1] = up[2] * ea[0];
1633  no[2] = up[0] * ea[1] - up[1] * ea[0];
1634 
1635  // Compute geodetic latitude and longitude
1636  float xlat_ = RADEG * asin(up[2]);
1637  float xlon_ = RADEG * atan2(up[1], up[0]);
1638  //check if it is in the range
1639  if( std::abs(xlat_) > 90.0e0 || std::abs(xlon_) > 180.0e0){
1640  continue;
1641  }
1642  xlat[i] = xlat_;
1643  xlon[i] = xlon_;
1644  // Transform the pixel-to-spacecraft and Sun vectors into local frame
1645  float rl[3], sl[3];
1646  rl[0] = -ea[0] * rh[0] - ea[1] * rh[1] - ea[2] * rh[2];
1647  rl[1] = -no[0] * rh[0] - no[1] * rh[1] - no[2] * rh[2];
1648  rl[2] = -up[0] * rh[0] - up[1] * rh[1] - up[2] * rh[2];
1649 
1650  sl[0] = sunr[0] * ea[0] + sunr[1] * ea[1] + sunr[2] * ea[2];
1651  sl[1] = sunr[0] * no[0] + sunr[1] * no[1] + sunr[2] * no[2];
1652  sl[2] = sunr[0] * up[0] + sunr[1] * up[1] + sunr[2] * up[2];
1653 
1654  // Compute the solar zenith and azimuth
1655  solz[i] = (short)(100 * RADEG *
1656  atan2(sqrt(sl[0] * sl[0] + sl[1] * sl[1]), sl[2]));
1657 
1658  // Check for zenith close to zero
1659  if (solz[i] > 0.01) sola[i] = (short)(100 * RADEG * atan2(sl[0], sl[1]));
1660 
1661  // Compute the sensor zenith and azimuth
1662  senz[i] = (short)(100 * RADEG *
1663  atan2(sqrt(rl[0] * rl[0] + rl[1] * rl[1]), rl[2]));
1664 
1665  // Check for zenith close to zero
1666  if (senz[i] > 0.01) sena[i] = (short)(100 * RADEG * atan2(rl[0], rl[1]));
1667 
1668  // terrain correct the pixel
1669  float tmp_sena = sena[i] / 100.0;
1670  float tmp_senz = senz[i] / 100.0;
1671  float tmp_height = 0;
1672  get_nc_height(dem_file, &xlon[i], &xlat[i],
1673  &tmp_senz, &tmp_sena, &tmp_height);
1674  // there is also a bug in terraing correction
1675  if( std::abs(xlat[i]) > 90.0e0 || std::abs(xlon[i]) > 180.0e0){
1676  xlat[i] = -32767;
1677  xlon[i] = -32767;
1678  solz[i] = -32767;
1679  sola[i] = -32767;
1680  senz[i] = -32767;
1681  sena[i] = -32767;
1682  continue;
1683  }
1684 
1685  if (xlon[i] > loneast) loneast = xlon[i];
1686  if (xlon[i] < lonwest) lonwest = xlon[i];
1687  if (xlat[i] < latsouth) latsouth = xlat[i];
1688  if (xlat[i] > latnorth) latnorth = xlat[i];
1689 
1690  sena[i] = (short)(tmp_sena * 100.0);
1691  senz[i] = (short)(tmp_senz * 100.0);
1692  height[i] = (short)tmp_height;
1693  } // if on-earth
1694  } // pixel loop
1695 
1696  gsl_vector_free(C);
1697 
1698  return 0;
1699 }
1700 
1701 int mat2rpy(float pos[3], float vel[3], double smat[3][3], float rpy[3]) {
1702  // This program calculates the attitude angles from the ECEF orbit vector and
1703  // attitude matrix. The rotation order is (1,2,3).
1704 
1705  // The reference ellipsoid uses an equatorial radius of 6378.137 km and
1706  // a flattening factor of 1/298.257 (WGS 1984).
1707 
1708  // Calling Arguments
1709 
1710  // Name Type I/O Description
1711  //
1712  // pos(3) R*4 I Orbit position vector (ECEF)
1713  // vel(3) R*4 I Orbit velocity vector (ECEF)
1714  // smat(3,3) R*8 I Sensor attitude matrix (ECEF to sensor)
1715  // rpy(3) R*4 O Attitude angles (roll, pitch, yaw)
1716 
1717  double rem = 6371;
1718  double f = 1 / (double)298.257;
1719  double omegae = 7.29211585494e-5;
1720  double omf2 = (1 - f) * (1 - f);
1721 
1722  // Determine local vertical reference axes
1723  double p[3], v[3];
1724  for (size_t j = 0; j < 3; j++) {
1725  p[j] = (double)pos[j];
1726  v[j] = (double)vel[j];
1727  }
1728  v[0] -= p[1] * omegae;
1729  v[1] += p[0] * omegae;
1730 
1731  // Compute Z axis as local nadir vector
1732  double pm = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
1733  double omf2p = (omf2 * rem + pm - rem) / pm;
1734  double pxy = p[0] * p[0] + p[1] * p[1];
1735  double temp = sqrt(p[2] * p[2] + omf2p * omf2p * pxy);
1736 
1737  double z[3];
1738  z[0] = -omf2p * p[0] / temp;
1739  z[1] = -omf2p * p[1] / temp;
1740  z[2] = -p[2] / temp;
1741 
1742  // Compute Y axis along negative orbit normal
1743  double on[3];
1744  on[0] = v[1] * z[2] - v[2] * z[1];
1745  on[1] = v[2] * z[0] - v[0] * z[2];
1746  on[2] = v[0] * z[1] - v[1] * z[0];
1747 
1748  double onm = sqrt(on[0] * on[0] + on[1] * on[1] + on[2] * on[2]);
1749 
1750  double y[3];
1751  for (size_t j = 0; j < 3; j++) y[j] = -on[j] / onm;
1752 
1753  // Compute X axis to complete orthonormal triad (velocity direction)
1754  double x[3];
1755  x[0] = y[1] * z[2] - y[2] * z[1];
1756  x[1] = y[2] * z[0] - y[0] * z[2];
1757  x[2] = y[0] * z[1] - y[1] * z[0];
1758 
1759  // Store local vertical reference vectors in matrix
1760  double om[3][3];
1761  memcpy(&om[0][0], &x, 3 * sizeof(double));
1762  memcpy(&om[1][0], &y, 3 * sizeof(double));
1763  memcpy(&om[2][0], &z, 3 * sizeof(double));
1764 
1765  // Compute orbital-to-spacecraft matrix
1766  double rm[3][3];
1767  gsl_matrix_view rm_mat = gsl_matrix_view_array(&rm[0][0], 3, 3);
1768 
1769  int s;
1770 
1771  gsl_permutation *perm = gsl_permutation_alloc(3);
1772 
1773  // Compute the LU decomposition of this matrix
1774  gsl_matrix_view B = gsl_matrix_view_array(&om[0][0], 3, 3);
1775  gsl_linalg_LU_decomp(&B.matrix, perm, &s);
1776 
1777  // Compute the inverse of the LU decomposition
1778  double inv[3][3];
1779  gsl_matrix_view inv_mat = gsl_matrix_view_array(&inv[0][0], 3, 3);
1780 
1781  gsl_linalg_LU_invert(&B.matrix, perm, &inv_mat.matrix);
1782 
1783  gsl_matrix_view A = gsl_matrix_view_array(&smat[0][0], 3, 3);
1784 
1785  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &A.matrix, &inv_mat.matrix,
1786  0.0, &rm_mat.matrix);
1787 
1788  gsl_permutation_free(perm);
1789 
1790  // Compute attitude angles
1791  rpy[0] = RADEG * atan(-rm[2][1] / rm[2][2]);
1792  double cosp = sqrt(rm[2][1] * rm[2][1] + rm[2][2] * rm[2][2]);
1793  if (rm[2][2] < 0) cosp *= -1;
1794  rpy[1] = RADEG * atan2(rm[2][0], cosp);
1795  rpy[2] = RADEG * atan(-rm[1][0] / rm[0][0]);
1796 
1797  return 0;
1798 }
1799 
1800 
1801 /*----------------------------------------------------------------- */
1802 /* Create a Generic NETCDF4 file */
1803 /* ---------------------------------------------------------------- */
1804 int geoFile::createFile( const char* filename, const char* cdlfile,
1805  size_t sdim, int *ncid, int *gid,
1806  uint32_t bbb, uint32_t rbb,
1807  int16_t pcdim, int16_t psdim, size_t nswband, int32_t *rta_nadir,
1808  bool radiance) {
1809 
1810  try {
1811  geofile = new NcFile( filename, NcFile::replace);
1812  }
1813  catch ( NcException& e) {
1814  e.what();
1815  cerr << "Failure creating OCI GEO file: " << filename << endl;
1816  exit(1);
1817  }
1818 
1819  ifstream output_data_structure;
1820  string line;
1821  string dataStructureFile;
1822 
1823  dataStructureFile.assign( cdlfile);
1824  expandEnvVar( &dataStructureFile);
1825 
1826  output_data_structure.open( dataStructureFile.c_str(), ifstream::in);
1827  if ( output_data_structure.fail() == true) {
1828  cout << "\"" << dataStructureFile.c_str() << "\" not found" << endl;
1829  exit(1);
1830  }
1831 
1832  // Find "dimensions" section of CDL file
1833  while(1) {
1834  getline( output_data_structure, line);
1835  size_t pos = line.find("dimensions:");
1836  if ( pos == 0) break;
1837  }
1838 
1839  // Define dimensions from "dimensions" section of CDL file
1840  ndims = 0;
1841  // int dimid[1000];
1842  while(1) {
1843  getline( output_data_structure, line);
1844  boost::trim(line);
1845  if (line.substr(0,2) == "//") continue;
1846 
1847  size_t pos = line.find(" = ");
1848  if ( pos == string::npos) break;
1849 
1850  uint32_t dimSize;
1851  istringstream iss(line.substr(pos+2, string::npos));
1852  iss >> dimSize;
1853 
1854  iss.clear();
1855  iss.str( line);
1856  iss >> skipws >> line;
1857 
1858  if (line.compare("number_of_scans") == 0) {
1859  dimSize = sdim;
1860  }
1861 
1862  if (line.compare("blue_bands") == 0) {
1863  dimSize = bbb;
1864  }
1865 
1866  if (line.compare("red_bands") == 0) {
1867  dimSize = rbb;
1868  }
1869 
1870  if (line.compare("ccd_pixels") == 0) {
1871  dimSize = pcdim;
1872  }
1873 
1874  if (line.compare("SWIR_pixels") == 0) {
1875  dimSize = psdim;
1876  }
1877 
1878  if (line.compare("SWIR_bands") == 0) {
1879  dimSize = nswband;
1880  }
1881  try {
1882  ncDims[ndims++] = geofile->addDim( line, dimSize);
1883  }
1884  catch ( NcException& e) {
1885  e.what();
1886  cerr << "Failure creating dimension: " << line.c_str() << endl;
1887  exit(1);
1888  }
1889 
1890  } // while loop
1891 
1892  // Read global attributes (string attributes only)
1893  while(1) {
1894  getline( output_data_structure, line);
1895  boost::trim(line);
1896  size_t pos = line.find("// global attributes");
1897  if ( pos == 0) break;
1898  }
1899 
1900  while(1) {
1901  getline( output_data_structure, line);
1902  size_t pos = line.find(" = ");
1903  if ( pos == string::npos) break;
1904 
1905  string attValue = line.substr(pos+3);
1906 
1907  // Remove any leading and trailing quotes
1908  attValue.erase(attValue.length()-2); // skip final " ;"
1909  size_t begQuote = attValue.find('"');
1910  size_t endQuote = attValue.find_last_of('"');
1911  if (begQuote == string::npos) continue; // Skip non-string global attributes
1912  attValue = attValue.substr(begQuote+1, endQuote-begQuote-1);
1913 
1914  istringstream iss(line.substr(pos+2));
1915  iss.clear();
1916  iss.str( line);
1917  iss >> skipws >> line;
1918 
1919  // Skip commented out attributes
1920  if (line.compare("//") == 0) continue;
1921 
1922  string attName;
1923  attName.assign(line.substr(1).c_str());
1924 
1925  try {
1926  geofile->putAtt(attName, attValue);
1927  }
1928  catch ( NcException& e) {
1929  e.what();
1930  cerr << "Failure creating attribute: " + attName << endl;
1931  exit(1);
1932  }
1933 
1934  } // while(1)
1935 
1936  geofile->putAtt("rta_nadir", NC_LONG, 2, rta_nadir);
1937 
1938  ngrps = 0;
1939  // Loop through groups
1940  while(1) {
1941  getline( output_data_structure, line);
1942 
1943  // Check if end of CDL file
1944  // If so then close CDL file and return
1945  if (line.substr(0,1).compare("}") == 0) {
1946  output_data_structure.close();
1947  return 0;
1948  }
1949 
1950  // Check for beginning of new group
1951  size_t pos = line.find("group:");
1952 
1953  // If found then create new group and variables
1954  if ( pos == 0) {
1955 
1956  // Parse group name
1957  istringstream iss(line.substr(6, string::npos));
1958  iss >> skipws >> line;
1959 
1960  // Create NCDF4 group
1961  ncGrps[ngrps++] = geofile->addGroup(line);
1962 
1963  int numDims=0;
1964  // int varDims[NC_MAX_DIMS];
1965  //size_t dimSize[NC_MAX_DIMS];
1966  //char dimName[NC_MAX_NAME+1];
1967  string sname;
1968  string lname;
1969  string standard_name;
1970  string units;
1971  string description;
1972  string flag_masks;
1973  string flag_meanings;
1974  double valid_min=0.0;
1975  double valid_max=0.0;
1976  double scale=1.0;
1977  double offset=0.0;
1978  double fill_value=0.0;
1979  string coordinates;
1980 
1981  vector<NcDim> varVec;
1982 
1983  int ntype=0;
1984  NcType ncType;
1985 
1986  // Loop through datasets in group
1987  getline( output_data_structure, line); // skip "variables:"
1988  while(1) {
1989  getline( output_data_structure, line);
1990  boost::trim(line);
1991 
1992  if (line.substr(0,2) == "//") continue;
1993  if (line.length() == 0) continue;
1994  if (line.substr(0,1).compare("\r") == 0) continue;
1995  if (line.substr(0,1).compare("\n") == 0) continue;
1996 
1997  if (radiance) {
1998  if (line.find("rhot_") != std::string::npos) continue;
1999  } else {
2000  if (line.find("Lt_") != std::string::npos) continue;
2001  }
2002 
2003  size_t pos = line.find(":");
2004 
2005  // No ":" found, new dataset or empty line or end-of-group
2006  if ( pos == string::npos) {
2007 
2008  if ( numDims > 0) {
2009  // Create previous dataset
2010  createField( ncGrps[ngrps-1], sname.c_str(), lname.c_str(),
2011  standard_name.c_str(), units.c_str(),
2012  description.c_str(),
2013  (void *) &fill_value,
2014  flag_masks.c_str(), flag_meanings.c_str(),
2015  valid_min, valid_max,
2016  scale, offset,
2017  ntype, varVec, coordinates);
2018 
2019  flag_masks.clear();
2020  flag_meanings.clear();
2021  units.clear();
2022  description.clear();
2023  varVec.clear();
2024  coordinates.clear();
2025  }
2026 
2027  valid_min=0.0;
2028  valid_max=0.0;
2029  scale=1.0;
2030  offset=0.0;
2031  fill_value=0.0;
2032 
2033  if (line.substr(0,10).compare("} // group") == 0) break;
2034 
2035  // Parse variable type
2036  string varType;
2037  istringstream iss(line);
2038  iss >> skipws >> varType;
2039 
2040  // Get corresponding NC variable type
2041  if ( varType.compare("char") == 0) ntype = NC_CHAR;
2042  else if ( varType.compare("byte") == 0) ntype = NC_BYTE;
2043  else if ( varType.compare("short") == 0) ntype = NC_SHORT;
2044  else if ( varType.compare("int") == 0) ntype = NC_INT;
2045  else if ( varType.compare("long") == 0) ntype = NC_INT;
2046  else if ( varType.compare("float") == 0) ntype = NC_FLOAT;
2047  else if ( varType.compare("real") == 0) ntype = NC_FLOAT;
2048  else if ( varType.compare("double") == 0) ntype = NC_DOUBLE;
2049  else if ( varType.compare("ubyte") == 0) ntype = NC_UBYTE;
2050  else if ( varType.compare("ushort") == 0) ntype = NC_USHORT;
2051  else if ( varType.compare("uint") == 0) ntype = NC_UINT;
2052  else if ( varType.compare("ulong") == 0) ntype = NC_UINT;
2053  else if ( varType.compare("int64") == 0) ntype = NC_INT64;
2054  else if ( varType.compare("uint64") == 0) ntype = NC_UINT64;
2055 
2056  // Parse short name (sname)
2057  pos = line.find("(");
2058  size_t posSname = line.substr(0, pos).rfind(" ");
2059  sname.assign(line.substr(posSname+1, pos-posSname-1));
2060 
2061  // Parse variable dimension info
2062  this->parseDims( line.substr(pos+1, string::npos), varVec);
2063  numDims = varVec.size();
2064 
2065  } else {
2066  // Parse variable attributes
2067  size_t posEql = line.find("=");
2068  size_t pos1qte = line.find("\"");
2069  size_t pos2qte = line.substr(pos1qte+1, string::npos).find("\"");
2070 
2071  string attrName = line.substr(pos+1, posEql-pos-2);
2072 
2073  // Get long_name
2074  if ( attrName.compare("long_name") == 0) {
2075  lname.assign(line.substr(pos1qte+1, pos2qte));
2076  }
2077 
2078  // Get units
2079  else if ( attrName.compare("units") == 0) {
2080  units.assign(line.substr(pos1qte+1, pos2qte));
2081  }
2082 
2083  // Get description
2084  else if ( attrName.compare("description") == 0) {
2085  description.assign(line.substr(pos1qte+1, pos2qte));
2086  }
2087 
2088  // Get _FillValue
2089  else if ( attrName.compare("_FillValue") == 0) {
2090  iss.clear();
2091  iss.str( line.substr(posEql+1, string::npos));
2092  iss >> fill_value;
2093  }
2094 
2095  // Get flag_masks
2096  else if ( attrName.compare("flag_masks") == 0) {
2097  flag_masks.assign(line.substr(pos1qte+1, pos2qte));
2098  }
2099  else if ( attrName.compare("flag_masks") == 0) {
2100  flag_masks.assign(line.substr(pos1qte+1, pos2qte));
2101  }
2102 
2103  // Get flag_meanings
2104  else if ( attrName.compare("flag_meanings") == 0) {
2105  flag_meanings.assign(line.substr(pos1qte+1, pos2qte));
2106  }
2107 
2108  // Get valid_min
2109  else if ( attrName.compare("valid_min") == 0) {
2110  iss.clear();
2111  iss.str( line.substr(posEql+1, string::npos));
2112  iss >> valid_min;
2113  }
2114 
2115  // Get valid_max
2116  else if ( attrName.compare("valid_max") == 0) {
2117  iss.clear();
2118  iss.str( line.substr(posEql+1, string::npos));
2119  iss >> valid_max;
2120  }
2121 
2122  // Get scale
2123  else if ( attrName.compare("scale_factor") == 0) {
2124  iss.clear();
2125  iss.str( line.substr(posEql+1, string::npos));
2126  iss >> scale;
2127  }
2128 
2129  // Get offset
2130  else if ( attrName.compare("add_offset") == 0) {
2131  iss.clear();
2132  iss.str( line.substr(posEql+1, string::npos));
2133  iss >> offset;
2134  }
2135 
2136  // Get coordinates
2137  else if ( attrName.compare("coordinates") == 0) {
2138  coordinates.assign(line.substr(pos1qte+1, pos2qte));
2139  }
2140 
2141  } // if ( pos == string::npos)
2142  } // datasets in group loop
2143  } // New Group loop
2144  } // Main Group loop
2145 
2146  return 0;
2147 }
2148 
2149 
2150 int geoFile::parseDims( string dimString, vector<NcDim>& varDims) {
2151 
2152  size_t curPos=0;
2153  // char dimName[NC_MAX_NAME+1];
2154  string dimName;
2155 
2156  while(1) {
2157  size_t pos = dimString.find(",", curPos);
2158  if ( pos == string::npos)
2159  pos = dimString.find(")");
2160 
2161  string varDimName;
2162  istringstream iss(dimString.substr(curPos, pos-curPos));
2163  iss >> skipws >> varDimName;
2164 
2165  for (int i=0; i<ndims; i++) {
2166 
2167  try {
2168  dimName = ncDims[i].getName();
2169  }
2170  catch ( NcException& e) {
2171  e.what();
2172  cerr << "Failure accessing dimension: " + dimName << endl;
2173  exit(1);
2174  }
2175 
2176  if ( varDimName.compare(dimName) == 0) {
2177  varDims.push_back(ncDims[i]);
2178  break;
2179  }
2180  }
2181  if ( dimString.substr(pos, 1).compare(")") == 0) break;
2182 
2183  curPos = pos + 1;
2184  }
2185 
2186  return 0;
2187 }
2188 
2189 int geoFile::close() {
2190 
2191  try {
2192  geofile->close();
2193  }
2194  catch ( NcException& e) {
2195  cout << e.what() << endl;
2196  cerr << "Failure closing: " + fileName << endl;
2197  exit(1);
2198  }
2199 
2200  return 0;
2201 }
2202 
q
data_t q
Definition: decode_rs.h:74
value
int32 value
Definition: Granule.c:1235
r
int r
Definition: decode_rs.h:73
t
data_t t[NROOTS+1]
Definition: decode_rs.h:77
geoFile::parseDims
int parseDims(std::string dimString, std::vector< netCDF::NcDim > &varDims)
Definition: geolocate_oci.cpp:2150
j
int j
Definition: decode_rs.h:73
geo_struct
Definition: common.h:25
vector
float * vector(long nl, long nh)
Definition: nrutil.c:15
day
int32_t day
Definition: atrem_corl1v2.h:308
geo_struct::rta_axis
double rta_axis[3]
Definition: common.h:35
expandEnvVar
int expandEnvVar(std::string *sValue)
Definition: ancgen.cpp:222
orb_array
float orb_array[3]
Definition: geolocate_hawkeye.h:5
quat_array
float quat_array[4]
Definition: geolocate_hawkeye.h:4
scan_ell
int scan_ell(float p[3], double sm[3][3], double coef[10])
Definition: geolocate_oci.cpp:1449
nutate
int nutate(double t, double xls, double gs, double xlm, double omega, double &dpsi, double &eps, double &epsm)
Definition: geolocate_oci.cpp:968
modis.modis_utils.geofile
geofile
Definition: modis_utils.py:368
geolocate_oci
int geolocate_oci(string l1a_filename, string geo_lut_filename, geo_struct &geoLUT, string l1b_filename, string dem_file, bool radiance, string doi, const string ephFile, string pversion)
Definition: geolocate_oci.cpp:117
l_sun
int l_sun(size_t sdim, int32_t iyr, int32_t iday, double *sec, orb_array *sunr, double *rs)
Definition: geolocate_oci.cpp:1337
C
const double C
Definition: calc_par.c:102
qtom
int qtom(float q[4], double rm[3][3])
Definition: geolocate_oci.cpp:1433
spin
void spin(double st, double *pos1, double *pos2)
read_def_eph
void read_def_eph(string ephfile, vector< double > &otime, vector< vector< double >> &posj, vector< vector< double >> &velj, double l1aEpoch)
Read a predicted PACE ephemeris file and convert vectors to ECR.
Definition: geolocate_oci.cpp:51
mat2rpy
int mat2rpy(float pos[3], float vel[3], double smat[3][3], float rpy[3])
Definition: geolocate_oci.cpp:1701
geo_struct::tilt_axis
double tilt_axis[3]
Definition: common.h:30
check_scan_times
int check_scan_times(uint32_t nscan, double *sstime, short *sfl)
Definition: common.cpp:863
geo_struct::ham_axis
double ham_axis[3]
Definition: common.h:36
color_dtdb.v
v
Definition: color_dtdb.py:430
q_interp
int q_interp(size_t n_att_rec, size_t sdim, double *tq, quat_array *q, double *time, quat_array *qi)
Definition: geolocate_oci.cpp:1249
unix2ymds
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
geo_struct::ham_at_angles
double ham_at_angles[2]
Definition: common.h:37
hico.proc_hico.valid_min
valid_min
Definition: proc_hico.py:246
pos
float32 * pos
Definition: l1_czcs_hdf.c:35
ephparms
int ephparms(double t, double &xls, double &gs, double &xlm, double &omega)
Definition: geolocate_oci.cpp:939
msec
int32 * msec
Definition: l1_czcs_hdf.c:31
geo_struct::as_planarity
double as_planarity[5]
Definition: common.h:44
jday
int32_t jday(int16_t i, int16_t j, int16_t k)
Converts a calendar date to the corresponding Julian day starting at noon on the calendar date....
Definition: jday.c:14
hico.proc_hico.valid_max
valid_max
Definition: proc_hico.py:247
dpsi
double dpsi
Definition: sun2000.c:6
nscan
int32 nscan
Definition: l1_czcs_hdf.c:19
eps
double eps
Definition: gha2000.c:3
createField
int createField(NcGroup &ncGrp, const char *sname, const char *lname, const char *standard_name, const char *units, const char *description, void *fill_value, const char *flag_masks, const char *flag_meanings, double low, double high, double scale, double offset, int nt, vector< NcDim > &varVec, string coordinates)
Definition: common.cpp:312
bzero
void bzero()
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
PI
#define PI
Definition: l3_get_org.c:6
A
const float A
Definition: calc_par.c:100
f
double precision function f(R1)
Definition: tmd.lp.f:1454
tempOut
ofstream tempOut
Definition: l1agen_oci.cpp:181
geo_struct::ham_enc_scale
double ham_enc_scale
Definition: common.h:40
hico.value_locate.x
list x
Definition: value_locate.py:64
get_oci_vecs
int get_oci_vecs(uint32_t nscan, uint16_t pdim, double as_planarity[5], double at_planarity[5], int32_t *rta_nadir, double ham_ct_angles[2], double ev_toff, int32_t spin, uint8_t hside, float *clines, double *delt, double revpsec, int32_t ppr_off, int16_t board_id, uint32_t nmcescan, int32_t *mspin, uint8_t *enc_count, float **hamenc, float **rtaenc, float **pview, double *theta, int16_t &iret)
Definition: common.cpp:195
get_history
string get_history(int argc, char *argv[])
Definition: DDOptions.cpp:496
omega
data_t omega[NROOTS+1]
Definition: decode_rs.h:77
hico.proc_hico.sena
sena
Definition: proc_hico.py:241
color_dtdb.y
y
Definition: color_dtdb.py:430
set_global_attrs
void set_global_attrs(NcFile *outfile, string history, string doi, string pversion)
Definition: global_attrs.cpp:40
rd
float rd(float x, float y, float z)
geo_struct::rta_enc_scale
double rta_enc_scale
Definition: common.h:39
viirs_ler_luts::cphi
real cphi
Definition: viirs_ler_luts.f95:49
hktgen_pace.att
int att
Definition: hktgen_pace.py:28
gha2000
int gha2000(int32_t iyr, double day, double &gha)
Definition: geolocate_oci.cpp:1048
renav_hawkeye.utcpole
utcpole
Definition: renav_hawkeye.py:294
color_dtdb.scale
int scale
Definition: color_dtdb.py:459
jd
Definition: jd.py:1
geo_struct::ham_ct_angles
double ham_ct_angles[2]
Definition: common.h:38
filename
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
orb_interp
int orb_interp(size_t n_orb_rec, size_t sdim, double *torb, orb_array *p, orb_array *v, double *time, orb_array *posi, orb_array *veli)
Definition: geolocate_oci.cpp:1139
parseDims
int parseDims(int ncid, int ndims, string dimString, int *numDims, int *dimid, int *varDims)
Definition: hawkeyeUtil.cpp:339
color_dtdb.rg
rg
Definition: color_dtdb.py:449
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
fix_mac_rpath.line
line
Definition: fix_mac_rpath.py:56
geo_struct::MCE_clock
double MCE_clock
Definition: common.h:27
get_nib_nbb
int get_nib_nbb(uint32_t ntaps, size_t *ia, uint32_t ntb[16], int16_t jagg[16], uint32_t &nib, uint32_t &nbb)
Definition: common.cpp:742
isodate2ydmsec
void isodate2ydmsec(char *date, int32_t *year, int32_t *day, int32_t *msec)
Definition: date2ydmsec.c:20
run_local.outfile
outfile
Definition: run_local.py:76
ncattredit.history
string history
Definition: ncattredit.py:30
trim
string & trim(string &s, const string &delimiters)
Definition: EnvsatUtil.cpp:29
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
tilt_interp
int tilt_interp(size_t n_tilts, size_t sdim, double *ttilt, float *tiltin, double *time, float *tilt)
Definition: geolocate_oci.cpp:1306
geoFile::createFile
int createFile(const char *filename, const char *cdlfile, size_t sdim, int *ncid, int *gid, uint32_t bbb, uint32_t rbb, int16_t pcdim, int16_t psdim, size_t nswband, int32_t *rta_nadir, bool radiance)
Definition: geolocate_oci.cpp:1804
geo_struct::master_clock
double master_clock
Definition: common.h:26
j2000_to_ecr
int j2000_to_ecr(int32_t iyr, int32_t idy, double sec, double ecmat[3][3])
Definition: geolocate_oci.cpp:819
MDN.setup.description
description
Definition: setup.py:16
oci_geonav
int oci_geonav(const char *dem_file, float pos[3], float vel[3], double smat[3][3], double coef[10], float sunr[3], float **pview, size_t npix, double *delt, float *xlat, float *xlon, short *solz, short *sola, short *senz, short *sena, short *height, float &lonwest, float &loneast, float &latsouth, float &latnorth)
Definition: geolocate_oci.cpp:1486
geo_struct::tilt_home
double tilt_home
Definition: common.h:32
xlon
float xlon[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:93
sun2000
int sun2000(size_t sdim, int32_t iyr, int32_t idy, double *sec, orb_array *sun, double *rs)
Definition: geolocate_oci.cpp:1362
pm
float pm[MODELMAX]
Definition: atrem_corl1v2.h:248
hico.proc_hico.units
units
Definition: proc_hico.py:245
get_ut1
int get_ut1(int32_t iyr, int32_t idy, double &ut1utc)
Definition: geolocate_oci.cpp:999
mapgen_overlay.norm
norm
Definition: mapgen_overlay.py:227
hico.proc_hico.sola
sola
Definition: proc_hico.py:242
fix_mac_rpath.lines
lines
Definition: fix_mac_rpath.py:35
B
const float B
Definition: calc_par.c:101
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
re
#define re
Definition: l1_czcs.c:695
geo_struct::tilt_to_oci_mech
double tilt_to_oci_mech[3][3]
Definition: common.h:33
dtype
dtype
Definition: DDataset.hpp:31
j2000_to_mod
int j2000_to_mod(int32_t iyr, int32_t idy, double sec, double j2mod[3][3])
Definition: geolocate_oci.cpp:876
time
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start time
Definition: HISTORY.txt:248
fabs
#define fabs(a)
Definition: misc.h:93
geo_struct::oci_mech_to_oci_opt
double oci_mech_to_oci_opt[3][3]
Definition: common.h:34
get_nc_height
int get_nc_height(float *xlon, float *xlat, float *senz, float *sena, float *height)
Definition: get_nc_height.c:574
ERROR_EXIT
#define ERROR_EXIT(dim)
Definition: geolocate_oci.cpp:23
global_attrs.h
qprod
int qprod(double q1[4], float q2[4], double q3[4])
Definition: geolocate_oci.cpp:1117
geo_struct::at_planarity
double at_planarity[5]
Definition: common.h:45
read_mce_tlm
int read_mce_tlm(NcFile *l1afile, geo_struct &geo_lut, NcGroup egid, uint32_t nmcescan, uint32_t nenc, int32_t &ppr_off, double &revpsec, double &secpline, int16_t &board_id, int32_t *mspin, int32_t *ot_10us, uint8_t *enc_count, float **hamenc, float **rtaenc, int16_t &iret)
Definition: common.cpp:8
s
data_t s[NROOTS]
Definition: decode_rs.h:75
ymds2unix
double ymds2unix(short year, short month, short day, double secs)
get_nut
int get_nut(int32_t iyr, int32_t idy, double xnut[3][3])
Definition: geolocate_oci.cpp:911
offset
l2prod offset
Definition: HOWTO_Add_a_product.txt:70
omf2
#define omf2
Definition: l1_czcs.c:697
geo_struct::tilt_angles
double tilt_angles[2]
Definition: common.h:31
idy
int32_t idy
Definition: atrem_corl1.h:161
tilt
int16_t * tilt
Definition: l2bin.cpp:80
standard_name
These two strings are used for the product XML output If product_id is not set then prefix is used If the last char of the name_prefix is _ then it is removed If algorithm_id is not set then name_suffix is used If the first char is _ then it is removed l2prod standard_name[0]
Definition: HOWTO_Add_a_product.txt:109
abs
#define abs(a)
Definition: misc.h:90
i
int i
Definition: decode_rs.h:71
on
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned on(as it will be in Near-Real-Time processing).
get_ev
int get_ev(double secpline, int16_t *dtype, int16_t *lines, int16_t *iagg, uint16_t &pcdim, uint16_t &psdim, double &ev_toff, float *clines, float *slines, double *deltc, double *delts, bool dark, int16_t &iret)
Definition: common.cpp:142
geo_struct::rta_nadir
int32_t rta_nadir[2]
Definition: common.h:42
iyr
int32_t iyr
Definition: atrem_corl1.h:161
mtoq
int mtoq(double rm[3][3], double q[4])
Definition: geolocate_oci.cpp:1080
k
int k
Definition: decode_rs.h:73
npix
int npix
Definition: get_cmp.c:28
geolocate_oci.h
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
viirs_ler_luts::xlat
real xlat
Definition: viirs_ler_luts.f95:44
geo_struct::sc_to_tilt
double sc_to_tilt[3][3]
Definition: common.h:29
isodate2unix
double isodate2unix(const char *isodate)
Definition: unix2isodate.c:61
count
int count
Definition: decode_rs.h:79


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster