OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
hawkeye_methods.cpp
Go to the documentation of this file.
1 #include <fstream>
2 #include <iostream>
3 #include <sstream>
4 
5 #include <gsl/gsl_blas.h>
6 #include <gsl/gsl_linalg.h>
7 #include <gsl/gsl_matrix_double.h>
8 #include <string.h>
9 #include <string>
10 #include "hawkeye_methods.h"
11 //#include "hawkeyeUtil.h"
12 #include "netcdf.h"
13 
14 #include <cstdio>
15 #include <cstdlib>
16 #include <vector>
17 #include "spline.h"
18 #include <math.h>
19 #include <timeutils.h>
20 
21 #define RADEG 57.29577951
22 
23 using namespace std;
24 
25 
26 
27 
28 //this function ONLY WORKS WITH MONOTONIC INCREASING X VARIABLE!!!
29 int latlon_interp_dist(size_t num_gridlines, double *dist_ai, float *lati, float *loni,double *dist, float *lati2, float *loni2){
30  // std::vector<double> X(5),Y(5);
31 // double f = INFINITY;
32  std::vector<double> dscan;
33  std::vector<double> latpix,lonpix;
34 // X[0]=0.1; X[1]=0.4; X[2]=1.2; X[3]=1.8; X[4]=2.0;
35  // Y[0]=0.1; Y[1]=0.7; Y[2]=0.6; Y[3]=1.1; Y[4]=0.9;
36 
37 //big loop for each pixel across-track we create a time series---
38  // for(int j=0;j<num_pixels;j++){
39 
40 
41  for(size_t i=0;i<num_gridlines-1;i++){
42  dscan.push_back(dist_ai[i]);
43  lati[i]+=90;
44  latpix.push_back(lati[i]);
45  if (loni[i]<0.0) loni[i]+=360;
46  lonpix.push_back(loni[i]);
47  // for (auto it = latpix.begin(); it != latpix.end(); ++it)
48  // cout << *it << " ";
49  }
50  //interpolate pix time series--
51  for(size_t k=0;k<num_gridlines-1;k++){
52  tk::spline s1,s2;
53  s1.set_boundary(tk::spline::second_deriv, 0.0,
54  tk::spline::second_deriv, 0.0);
55  s2.set_boundary(tk::spline::second_deriv, 0.0,
56  tk::spline::second_deriv, 0.0);
57 
58  s1.set_points(dscan,latpix);
59  s2.set_points(dscan,lonpix);
60  lati2[k]=s1(dist[k]);
61  loni2[k]=s2(dist[k]);
62  // cout<<"lat.."<<lati2[k]<<"lon.."<<loni2[k]<<endl;
63 
64  if (loni2[k]>180.) loni2[k]-=360;
65  lati2[k]-=90;
66  }
67 
68 
69  dscan.clear();
70  latpix.clear();
71  lonpix.clear();
72  return(0);
73 }
74 
75 
76 int latlon_interp2(size_t num_gridlines,size_t num_pixels, double *gtime, float *lati, float *loni,double *gtime2, float *lati2, float *loni2){
77  // std::vector<double> X(5),Y(5);
78 // double f = INFINITY;
79  std::vector<double> tscan;
80  std::vector<double> latpix,lonpix;
81 // X[0]=0.1; X[1]=0.4; X[2]=1.2; X[3]=1.8; X[4]=2.0;
82  // Y[0]=0.1; Y[1]=0.7; Y[2]=0.6; Y[3]=1.1; Y[4]=0.9;
83 
84 //big loop for each pixel across-track we create a time series---
85  // for(int j=0;j<num_pixels;j++){
86  // int j=(num_pixels/2+1/2)-1;//nad pixel
87  for(size_t i=0;i<num_gridlines;i++){
88  tscan.push_back(gtime[i]);
89  lati[i]+=90;
90  latpix.push_back(lati[i]);
91  if (loni[i]<0.0) loni[i]+=360;
92  lonpix.push_back(loni[i]);
93  // for (auto it = latpix.begin(); it != latpix.end(); ++it)
94  // cout << *it << " ";
95  }
96  //interpolate pix time series--
97  for(size_t k=0;k<num_gridlines;k++){
98  tk::spline s1,s2;
99  s1.set_points(tscan,latpix);
100  s2.set_points(tscan,lonpix);
101  lati2[k]=s1(gtime2[k]);
102  loni2[k]=s2(gtime2[k]);
103  if (loni2[k]>180.) loni2[k]-=360;
104  lati2[k]-=90;
105  }
106 
107 
108  tscan.clear();
109  latpix.clear();
110  lonpix.clear();
111  return(0);
112 }
113 
114 
115 
116 int latlon_interp1pix(size_t num_gridlines, size_t gd_index,double *tgrid, float *lat_nad, float *lon_nad, float *latipix, float *lonipix){
117 
118 
119  std::vector<double> tscan;
120  std::vector<double> latpix,lonpix;
121 
122  cout<<"interpolating 1 nadir position at the time..latlon_interp1pix."<<endl;
123 
124  for(size_t i=0;i<num_gridlines;i++){
125  tscan.push_back(tgrid[i]);
126  lat_nad[i]+=90;
127  latpix.push_back(lat_nad[i]);
128  if (lon_nad[i]<0.0) lon_nad[i]+=360;
129  lonpix.push_back(lon_nad[i]);
130  }
131  //interpolate pix time series--
132 
133  tk::spline s1,s2;
134  s1.set_boundary(tk::spline::second_deriv, 0.0,
135  tk::spline::second_deriv, 0.0);
136  s2.set_boundary(tk::spline::second_deriv, 0.0,
137  tk::spline::second_deriv, 0.0);
138 
139  s1.set_points(tscan,latpix);
140  s2.set_points(tscan,lonpix);
141  cout<<"ok filling latpix/lonpix vectors..."<<endl;
142  cout<<"gd.."<<gd_index<<"time.."<<tgrid[gd_index]<<"latnad.."<<lat_nad[gd_index]<<"lonnad.."<<lon_nad[gd_index]<<"lat.."<<s1(tgrid[gd_index])<<"..."<<"lon.."<<s2(tgrid[gd_index])<<endl;
143 
144 
145  *latipix=s1(tgrid[gd_index]);
146  *lonipix=s2(tgrid[gd_index]);
147  if (*lonipix>180.) *lonipix-=360;
148  *latipix-=90;
149 
150  tscan.clear();
151  latpix.clear();
152  lonpix.clear();
153  return 0;
154  }
155 
156 
157 
158 int latlon_interp_vec(size_t n_orb_rec, size_t num_gridlines,double *torb, double *latorb, double *lonorb,double *tmgv, float *lati, float *loni){
159 
160  double elon1,elon2,elat1,elat2,thres=5;
161  std::vector<double> tscan,tlat,tlon;
162  std::vector<double> latpix,lonpix,latclean,lonclean;
163  std::vector<int> latidx,lonidx;
164  double latemp=0.,lontemp=0.;
165 
166 //big loop for each pixel across-track we create a time series---
167 
168  for(size_t i=0;i<n_orb_rec;i++){
169  tscan.push_back(torb[i]);
170  latemp=latorb[i]+90;
171  latpix.push_back(latemp);
172 
173  if (lonorb[i]<0.0){
174  lontemp=lonorb[i]+360;}
175  else lontemp=lonorb[i];
176 
177  lonpix.push_back(lontemp);
178  }
179 
180  //interpolate pix time series--
181  for(size_t k=0;k<num_gridlines;k++){
182  tk::spline s1,s2;
183  s1.set_boundary(tk::spline::second_deriv, 0.0,
184  tk::spline::second_deriv, 0.0);
185  s2.set_boundary(tk::spline::second_deriv, 0.0,
186  tk::spline::second_deriv, 0.0);
187 
188  s1.set_points(tscan,latpix);
189  s2.set_points(tscan,lonpix);
190  lati[k]=s1(tmgv[k]);
191  loni[k]=s2(tmgv[k]);
192  if (loni[k]>180.) loni[k]-=360;
193  lati[k]-=90;
194 
195 
196  //check for interpolation artifacts and remove values------
197  //longitude
198 
199  if(k>=1 && k+1<num_gridlines){
200  elon1=abs(loni[k]-loni[k-1]);
201  elon2=abs(loni[k+1]-loni[k]);
202  elat1=abs(lati[k]-lati[k-1]);
203  elat2=abs(lati[k+1]-lati[k]);
204 
205  if(abs(100*elon2/elon1)<=100+thres && abs(100*elon2/elon1)>=100-thres){
206  tlon.push_back(tmgv[k]);
207  lonidx.push_back(k);
208  lonclean.push_back(loni[k]);}
209  if(abs(100*elat2/elat1)<=100+thres && abs(100*elat2/elat1)>=100-thres){
210  tlat.push_back(tmgv[k]);
211  latidx.push_back(k);
212  latclean.push_back(lati[k]);}
213  }//check lat lon difference
214 
215  }//end gridline
216 
217  tscan.clear();
218  tlat.clear();
219  tlon.clear();
220  latpix.clear();
221  lonpix.clear();
222  latidx.clear();
223  lonidx.clear();
224  latclean.clear();
225  lonclean.clear();
226  return 0;
227  }
228 
229 
230 
231 
232 
233 
234 
235 
236 int latlon_interp(size_t n_orb_rec, size_t num_gridlines,size_t num_pixels, double *torb, float **lat, float **lon,double *time, float *lati, float *loni){
237  // std::vector<double> X(5),Y(5);
238  double elon1,elon2,elat1,elat2,thres=5;
239 // double f = INFINITY;
240  std::vector<double> tscan,tlat,tlon;
241  std::vector<double> latpix,lonpix,latclean,lonclean;
242  // std::vector<double> erelon,erelat;
243  std::vector<int> latidx,lonidx;
244 // int k0lat,k0lon,c=0;
245 // X[0]=0.1; X[1]=0.4; X[2]=1.2; X[3]=1.8; X[4]=2.0;
246  // Y[0]=0.1; Y[1]=0.7; Y[2]=0.6; Y[3]=1.1; Y[4]=0.9;
247 
248 //big loop for each pixel across-track we create a time series---
249  // for(int j=0;j<num_pixels;j++){
250  int j=(num_pixels-1)/2;//nad pixel
251  cout<<"interpolating based on nad pixels....index #.."<<j<<endl;
252  for(size_t i=0;i<n_orb_rec;i++){
253  tscan.push_back(torb[i]);
254  lat[i][j]+=90;
255  latpix.push_back(lat[i][j]);
256 
257  if (lon[i][j]<0.0) lon[i][j]+=360;
258  lonpix.push_back(lon[i][j]);
259  }
260 
261  //interpolate pix time series--
262  for(size_t k=0;k<num_gridlines;k++){
263  tk::spline s1,s2;
264  s1.set_boundary(tk::spline::second_deriv, 0.0,
265  tk::spline::second_deriv, 0.0);
266  s2.set_boundary(tk::spline::second_deriv, 0.0,
267  tk::spline::second_deriv, 0.0);
268 
269  s1.set_points(tscan,latpix);
270  s2.set_points(tscan,lonpix);
271  lati[k]=s1(time[k]);
272  loni[k]=s2(time[k]);
273  if (loni[k]>180.) loni[k]-=360;
274  lati[k]-=90;
275 
276 
277  //check for interpolation artifacts and remove values------
278  //longitude
279 
280  if(k>=1 && k+1<num_gridlines){
281  elon1=abs(loni[k]-loni[k-1]);
282  elon2=abs(loni[k+1]-loni[k]);
283  elat1=abs(lati[k]-lati[k-1]);
284  elat2=abs(lati[k+1]-lati[k]);
285 
286  if(abs(100*elon2/elon1)<=100+thres && abs(100*elon2/elon1)>=100-thres){
287  tlon.push_back(time[k]);
288  lonidx.push_back(k);
289  lonclean.push_back(loni[k]);}
290  if(abs(100*elat2/elat1)<=100+thres && abs(100*elat2/elat1)>=100-thres){
291  tlat.push_back(time[k]);
292  latidx.push_back(k);
293  latclean.push_back(lati[k]);}
294 
295  }//check lat lon difference
296 
297  }//end gridline
298 
299 
300 /* cout<<"tlat.size."<<tlat.size()<<"latclean size.."<<latclean.size()<<"tlon.size."<<tlon.size()<<"lonclean size.."<<lonclean.size()<<endl;
301  //second interpolation after removing weird values
302  for(int k=0;k<num_gridlines;k++){
303  tk::spline s3,s4;
304  s3.set_points(tlat,latclean);
305  s4.set_points(tlon,lonclean);
306  lati[k]=s3(time[k]);
307  loni[k]=s4(time[k]);
308  }
309 
310 
311 
312  for(int j=0;j<lonidx.size();j++) cout<<"lonidx.."<<lonidx.operator[](j)<<"lonclean.."<<lonclean.operator[](j)<<endl;
313  for(int j=0;j<latidx.size();j++) cout<<"latidx.."<<latidx.operator[](j)<<"latclean.."<<latclean.operator[](j)<<endl;
314 
315  exit(1);
316 */
317 
318  tscan.clear();
319  tlat.clear();
320  tlon.clear();
321  latpix.clear();
322  lonpix.clear();
323  latidx.clear();
324  lonidx.clear();
325  latclean.clear();
326  lonclean.clear();
327  // erelon.clear();
328  // erelat.clear();
329  // cout<<"pix.."<<j+1<<endl;
330 
331  //************* backup vectors ****************
332  /* ofstream myfile ("/accounts/mamontes/images/OCIS/lati2.txt");
333  if (myfile.is_open())
334  {
335  myfile << "lati interp mean v.\n";
336  // myfile << "along-track ground positions.\n";
337  for(int count = 0; count <num_gridlines; count ++){
338  myfile << lati[count][j] << "," ;
339  }
340  myfile.close();
341  }
342  else cout << "Unable to open file";
343 
344  ofstream myfile2 ("/accounts/mamontes/images/OCIS/loni2.txt");
345  if (myfile2.is_open())
346  {
347  myfile2 << "loni interp mean v.\n";
348  // myfile << "along-track ground positions.\n";
349  for(int count = 0; count <num_gridlines; count ++){
350  myfile2 << loni[count][j] << "," ;
351  }
352  myfile2.close();
353  }
354  else cout << "Unable to open file";
355 
356 
357  ofstream myfile3 ("/accounts/mamontes/images/OCIS/gtime.txt");
358  if (myfile3.is_open())
359  {
360  myfile3 << "gtime....\n";
361  // myfile << "along-track ground positions.\n";
362  for(int count = 0; count <num_gridlines; count ++){
363  myfile3 << time[count] << "," ;
364  }
365  myfile3.close();
366  }
367  else cout << "Unable to open file";
368 
369 
370  ofstream myfile4 ("/accounts/mamontes/images/OCIS/torb.txt");
371  if (myfile4.is_open())
372  {
373  myfile4 << "torb.\n";
374  // myfile << "along-track ground positions.\n";
375  for(int count = 0; count <n_orb_rec; count ++){
376  myfile4 << torb[count] << "," ;
377  }
378  myfile4.close();
379  }
380  else cout << "Unable to open file";
381 
382  ofstream myfile5 ("/accounts/mamontes/images/OCIS/latorb.txt");
383  if (myfile5.is_open())
384  {
385  myfile5 << "latorb.\n";
386  // myfile << "along-track ground positions.\n";
387  for(int count = 0; count <n_orb_rec; count ++){
388  myfile5 << lat[count][j] << "," ;
389  }
390  myfile5.close();
391  }
392  else cout << "Unable to open file";
393 
394  ofstream myfile6 ("/accounts/mamontes/images/OCIS/lonorb.txt");
395  if (myfile6.is_open())
396  {
397  myfile6 << "lonorb.\n";
398  // myfile << "along-track ground positions.\n";
399  for(int count = 0; count <n_orb_rec; count ++){
400  myfile6 << lon[count][j] << "," ;
401  }
402  myfile6.close();
403  }
404  else cout << "Unable to open file";
405 
406  */
407 
408 
409 
410  // exit(1);
411  // } //end pixel
412 
413  // tk::spline s;
414  // s.set_points(X,Y);
415  // s.set_points(tscan,latpix); // currently it is required that X is already sorted
416  // double x=1.5;
417  // printf("spline at %f is %f\n", x, s(x));
418 // for(int k=0;k<num_gridlines;k++){
419 // cout<<"gd.."<<k+1<<"pix..."<<j+1<<"lati.."<<lati[k][j]<<"loni..."<<loni[k][j]<<endl;
420 // }
421 
422 // exit(1);
423  return 0;
424  }
425 
426 
427 
428 int orb_interp2(size_t n_orb_rec, size_t sdim, double *torb, orb_array2 *p,
429  orb_array2 *v, double *time, orb_array2 *posi, orb_array2 *veli) {
430  // Purpose: Interpolate orbit position and velocity vectors to scan line
431  // times
432  //
433  //
434  // Calling Arguments:
435  //
436  // Name Type I/O Description
437  // -------- ---- --- -----------
438  // torb(*) double I Array of orbit vector times in seconds of day
439  // p(3,*) float I Array of orbit position vectors for
440  // each time torb
441  // v(3,*) float I Array of orbit velocity vectors for
442  // each time torb
443  // time(*) double I Array of time in seconds of day
444  // for every scan line
445  // pi(3,*) float O Array of interpolated positions
446  // vi(3,*) float O Array of interpolated velocities
447  //
448  //
449  // By: Frederick S. Patt, GSC, August 10, 1993
450  //
451  // Notes: Method uses cubic polynomial to match positions and velocities
452  // at input data points.
453  //
454  // Modification History:
455  //
456  // Created IDL version from FORTRAN code. F.S. Patt, SAIC, November 29, 2006
457  //
458 
459  double a0[3], a1[3], a2[3], a3[3];
460 
461  /*
462  // Make sure that first orbit vector precedes first scan line
463  k = where (time lt torb(0))
464  if (k(0) ne -1) then begin
465  posi(*,k) = 0.0
466  veli(*,k) = 0.0
467  orbfl(k) = 1
468  print, 'Scan line times before available orbit data'
469  i1 = max(k) + 1
470  endif
471  */
472 
473  size_t ind = 0;
474  for (size_t i = 0; i < sdim; i++) {
475  // Find input orbit vectors bracketing scan
476  for (size_t j = ind; j < n_orb_rec; j++) {
477  if (time[i] > torb[j] && time[i] <= torb[j + 1]) {
478  ind = j;
479  break;
480  }
481  }
482 
483  // Set up cubic interpolation
484  double dt = torb[ind + 1] - torb[ind];
485  for (size_t j = 0; j < 3; j++) {
486  a0[j] = p[ind][j];
487  a1[j] = v[ind][j] * dt;
488  a2[j] = 3 * p[ind + 1][j] - 3 * p[ind][j] - 2 * v[ind][j] * dt -
489  v[ind + 1][j] * dt;
490  a3[j] = 2 * p[ind][j] - 2 * p[ind + 1][j] + v[ind][j] * dt +
491  v[ind + 1][j] * dt;
492  }
493 
494  // Interpolate orbit position and velocity components to scan line time
495  double x = (time[i] - torb[ind]) / dt;
496  double x2 = x * x;
497  double x3 = x2 * x;
498  for (size_t j = 0; j < 3; j++) {
499  posi[i][j] = a0[j] + a1[j] * x + a2[j] * x2 + a3[j] * x3;
500  veli[i][j] = (a1[j] + 2 * a2[j] * x + 3 * a3[j] * x2) / dt;
501  }
502  } // i-loop
503 
504  return 0;
505 }
506 
507 
508 
509 int orb_interp(size_t n_orb_rec, size_t sdim, double *torb, orb_array *p,
510  orb_array *v, double *time, orb_array *posi, orb_array *veli) {
511  // Purpose: Interpolate orbit position and velocity vectors to scan line
512  // times
513  //
514  //
515  // Calling Arguments:
516  //
517  // Name Type I/O Description
518  // -------- ---- --- -----------
519  // torb(*) double I Array of orbit vector times in seconds of day
520  // p(3,*) float I Array of orbit position vectors for
521  // each time torb
522  // v(3,*) float I Array of orbit velocity vectors for
523  // each time torb
524  // time(*) double I Array of time in seconds of day
525  // for every scan line
526  // pi(3,*) float O Array of interpolated positions
527  // vi(3,*) float O Array of interpolated velocities
528  //
529  //
530  // By: Frederick S. Patt, GSC, August 10, 1993
531  //
532  // Notes: Method uses cubic polynomial to match positions and velocities
533  // at input data points.
534  //
535  // Modification History:
536  //
537  // Created IDL version from FORTRAN code. F.S. Patt, SAIC, November 29, 2006
538  //
539 
540  double a0[3], a1[3], a2[3], a3[3];
541 
542  /*
543  // Make sure that first orbit vector precedes first scan line
544  k = where (time lt torb(0))
545  if (k(0) ne -1) then begin
546  posi(*,k) = 0.0
547  veli(*,k) = 0.0
548  orbfl(k) = 1
549  print, 'Scan line times before available orbit data'
550  i1 = max(k) + 1
551  endif
552  */
553 
554  size_t ind = 0;
555  for (size_t i = 0; i < sdim; i++) {
556  // Find input orbit vectors bracketing scan
557  for (size_t j = ind; j < n_orb_rec; j++) {
558  if (time[i] > torb[j] && time[i] <= torb[j + 1]) {
559  ind = j;
560  break;
561  }
562  }
563 
564  // Set up cubic interpolation
565  double dt = torb[ind + 1] - torb[ind];
566  for (size_t j = 0; j < 3; j++) {
567  a0[j] = p[ind][j];
568  a1[j] = v[ind][j] * dt;
569  a2[j] = 3 * p[ind + 1][j] - 3 * p[ind][j] - 2 * v[ind][j] * dt -
570  v[ind + 1][j] * dt;
571  a3[j] = 2 * p[ind][j] - 2 * p[ind + 1][j] + v[ind][j] * dt +
572  v[ind + 1][j] * dt;
573  }
574 
575  // Interpolate orbit position and velocity components to scan line time
576  double x = (time[i] - torb[ind]) / dt;
577  double x2 = x * x;
578  double x3 = x2 * x;
579  for (size_t j = 0; j < 3; j++) {
580  posi[i][j] = a0[j] + a1[j] * x + a2[j] * x2 + a3[j] * x3;
581  veli[i][j] = (a1[j] + 2 * a2[j] * x + 3 * a3[j] * x2) / dt;
582  }
583  } // i-loop
584 
585  return 0;
586 }
587 
588 
589 int j2000_to_ecr(int32_t iyr, int32_t idy, double sec, double ecmat[3][3]) {
590  // Get J2000 to ECEF transformation matrix
591 
592  // Arguments:
593 
594  // Name Type I/O Description
595  // --------------------------------------------------------
596  // iyr I I Year, four digits
597  // idy I I Day of year
598  // sec R*8 I Seconds of day
599  // ecmat(3,3) R O J2000 to ECEF matrix
600 
601  // Get transformation from J2000 to mean-of-date inertial
602  double j2mod[3][3];
603  j2000_to_mod(iyr, idy, sec, j2mod);
604 
605  // Get nutation and UT1-UTC (once per run)
606  double xnut[3][3], ut1utc;
607  get_nut(iyr, idy, xnut);
608  get_ut1(iyr, idy, ut1utc);
609 
610  // Compute Greenwich hour angle for time of day
611  double day = idy + (sec + ut1utc) / 86400;
612  double gha, gham[3][3];
613  gha2000(iyr, day, gha);
614 
615  gham[0][0] = cos(gha / RADEG);
616  gham[1][1] = cos(gha / RADEG);
617  gham[2][2] = 1;
618  gham[0][1] = sin(gha / RADEG);
619  gham[1][0] = -sin(gha / RADEG);
620 
621  gham[0][2] = 0;
622  gham[2][0] = 0;
623  gham[1][2] = 0;
624  gham[2][1] = 0;
625 
626  // Combine all transformations
627  gsl_matrix_view A = gsl_matrix_view_array(&gham[0][0], 3, 3); // gham
628  gsl_matrix_view B = gsl_matrix_view_array(&xnut[0][0], 3, 3); // xnut
629  gsl_matrix *C = gsl_matrix_alloc(3, 3);
630 
631  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, &A.matrix, &B.matrix, 0.0, C);
632 
633  gsl_matrix_view D = gsl_matrix_view_array(&j2mod[0][0], 3, 3); // j2mod
634  gsl_matrix *E = gsl_matrix_alloc(3, 3);
635  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, C, &D.matrix, 0.0, E);
636  double *ptr_E = gsl_matrix_ptr(E, 0, 0);
637 
638  memcpy(ecmat, ptr_E, 9 * sizeof(double));
639 
640  gsl_matrix_free(C);
641  gsl_matrix_free(E);
642 
643  return 0;
644 }
645 
646 int j2000_to_mod(int32_t iyr, int32_t idy, double sec, double j2mod[3][3]) {
647  // Get J2000 to MOD (precession) transformation
648 
649  // Arguments:
650 
651  // Name Type I/O Description
652  // --------------------------------------------------------
653  // iyr I I Year, four digits
654  // idy I I Day of year
655  // sec R*8 I Seconds of day
656  // j2mod(3,3) R O J2000 to MOD matrix
657 
658  int16_t iyear = iyr;
659  int16_t iday = idy;
660 
661  double t = jday(iyear, 1, iday) - (double)2451545.5 + sec / 86400;
662  t /= 36525;
663 
664  double zeta0 = t * (2306.2181 + 0.302 * t + 0.018 * t * t) / RADEG / 3600;
665  double thetap = t * (2004.3109 - 0.4266 * t - 0.04160 * t * t) / RADEG / 3600;
666  double xip = t * (2306.2181 + 1.095 * t + 0.018 * t * t) / RADEG / 3600;
667 
668  j2mod[0][0] = -sin(zeta0) * sin(xip) + cos(zeta0) * cos(xip) * cos(thetap);
669  j2mod[0][1] = -cos(zeta0) * sin(xip) - sin(zeta0) * cos(xip) * cos(thetap);
670  j2mod[0][2] = -cos(xip) * sin(thetap);
671  j2mod[1][0] = sin(zeta0) * cos(xip) + cos(zeta0) * sin(xip) * cos(thetap);
672  j2mod[1][1] = cos(zeta0) * cos(xip) - sin(zeta0) * sin(xip) * cos(thetap);
673  j2mod[1][2] = -sin(xip) * sin(thetap);
674  j2mod[2][0] = cos(zeta0) * sin(thetap);
675  j2mod[2][1] = -sin(zeta0) * sin(thetap);
676  j2mod[2][2] = cos(thetap);
677 
678  return 0;
679 }
680 
681 int get_nut(int32_t iyr, int32_t idy, double xnut[3][3]) {
682  int16_t iyear = iyr;
683  int16_t iday = idy;
684 
685  double t = jday(iyear, 1, iday) - (double)2451545.5;
686 
687  double xls, gs, xlm, omega;
688  double dpsi, eps, epsm;
689  ephparms(t, xls, gs, xlm, omega);
690  nutate(t, xls, gs, xlm, omega, dpsi, eps, epsm);
691 
692  xnut[0][0] = cos(dpsi / RADEG);
693  xnut[1][0] = -sin(dpsi / RADEG) * cos(epsm / RADEG);
694  xnut[2][0] = -sin(dpsi / RADEG) * sin(epsm / RADEG);
695  xnut[0][1] = sin(dpsi / RADEG) * cos(eps / RADEG);
696  xnut[1][1] = cos(dpsi / RADEG) * cos(eps / RADEG) * cos(epsm / RADEG) +
697  sin(eps / RADEG) * sin(epsm / RADEG);
698  xnut[2][1] = cos(dpsi / RADEG) * cos(eps / RADEG) * sin(epsm / RADEG) -
699  sin(eps / RADEG) * cos(epsm / RADEG);
700  xnut[0][2] = sin(dpsi / RADEG) * sin(eps / RADEG);
701  xnut[1][2] = cos(dpsi / RADEG) * sin(eps / RADEG) * cos(epsm / RADEG) -
702  cos(eps / RADEG) * sin(epsm / RADEG);
703  xnut[2][2] = cos(dpsi / RADEG) * sin(eps / RADEG) * sin(epsm / RADEG) +
704  cos(eps / RADEG) * cos(epsm / RADEG);
705 
706  return 0;
707 }
708 
709  int ephparms(double t, double &xls, double &gs, double &xlm, double &omega) {
710  // This subroutine computes ephemeris parameters used by other Mission
711  // Operations routines: the solar mean longitude and mean anomaly, and
712  // the lunar mean longitude and mean ascending node. It uses the model
713  // referenced in The Astronomical Almanac for 1984, Section S
714  // (Supplement) and documented for the SeaWiFS Project in "Constants
715  // and Parameters for SeaWiFS Mission Operations", in TBD. These
716  // parameters are used to compute the solar longitude and the nutation
717  // in longitude and obliquity.
718 
719  // Sun Mean Longitude
720  xls = (double)280.46592 + ((double)0.9856473516) * t;
721  xls = fmod(xls, (double)360);
722 
723  // Sun Mean Anomaly
724  gs = (double)357.52772 + ((double)0.9856002831) * t;
725  gs = fmod(gs, (double)360);
726 
727  // Moon Mean Longitude
728  xlm = (double)218.31643 + ((double)13.17639648) * t;
729  xlm = fmod(xlm, (double)360);
730 
731  // Ascending Node of Moon's Mean Orbit
732  omega = (double)125.04452 - ((double)0.0529537648) * t;
733  omega = fmod(omega, (double)360);
734 
735  return 0;
736 }
737 
738 int nutate(double t, double xls, double gs, double xlm, double omega,
739  double &dpsi, double &eps, double &epsm) {
740  // This subroutine computes the nutation in longitude and the obliquity
741  // of the ecliptic corrected for nutation. It uses the model referenced
742  // in The Astronomical Almanac for 1984, Section S (Supplement) and
743  // documented for the SeaWiFS Project in "Constants and Parameters for
744  // SeaWiFS Mission Operations", in TBD. These parameters are used to
745  // compute the apparent time correction to the Greenwich Hour Angle and
746  // for the calculation of the geocentric Sun vector. The input
747  // ephemeris parameters are computed using subroutine ephparms. Terms
748  // are included to 0.1 arcsecond.
749 
750  // Nutation in Longitude
751  dpsi = -17.1996 * sin(omega / RADEG) + 0.2062 * sin(2. * omega / RADEG) -
752  1.3187 * sin(2. * xls / RADEG) + 0.1426 * sin(gs / RADEG) -
753  0.2274 * sin(2. * xlm / RADEG);
754 
755  // Mean Obliquity of the Ecliptic
756  epsm = (double)23.439291 - ((double)3.560e-7) * t;
757 
758  // Nutation in Obliquity
759  double deps = 9.2025 * cos(omega / RADEG) + 0.5736 * cos(2. * xls / RADEG);
760 
761  // True Obliquity of the Ecliptic
762  eps = epsm + deps / 3600;
763 
764  dpsi = dpsi / 3600;
765 
766  return 0;
767 }
768 
769 int get_ut1(int32_t iyr, int32_t idy, double &ut1utc) {
770  int16_t iyear = iyr;
771  int16_t iday = idy;
772 
773  static int32_t ijd[25000];
774  static double ut1[25000];
775  string utcpole = "$OCVARROOT/var/modis/utcpole.dat";
776  static bool first = true;
777  int k = 0;
778  if (first) {
779  string line;
781  istringstream istr;
782 
783  ifstream utcpfile(utcpole.c_str());
784  if (utcpfile.is_open()) {
785  getline(utcpfile, line);
786  getline(utcpfile, line);
787  while (getline(utcpfile, line)) {
788  // cout << line << '\n';
789  istr.clear();
790  istr.str(line.substr(0, 5));
791  istr >> ijd[k];
792  istr.clear();
793  istr.str(line.substr(44, 9));
794  istr >> ut1[k];
795  k++;
796  }
797  ijd[k] = 0;
798  utcpfile.close();
799  first = false;
800  } else {
801  cout << utcpole.c_str() << " not found" << endl;
802  exit(1);
803  }
804  } // if (first)
805 
806  k = 0;
807  int mjd = jday(iyear, 1, iday) - 2400000;
808  while (ijd[k] > 0) {
809  if (mjd == ijd[k]) {
810  ut1utc = ut1[k];
811  break;
812  }
813  k++;
814  }
815 
816  return 0;
817 }
818 
819 int gha2000(int32_t iyr, double day, double &gha) {
820  // This subroutine computes the Greenwich hour angle in degrees for the
821  // input time. It uses the model referenced in The Astronomical Almanac
822  // for 1984, Section S (Supplement) and documented for the SeaWiFS
823  // Project in "Constants and Parameters for SeaWiFS Mission Operations",
824  // in TBD. It includes the correction to mean sidereal time for nutation
825  // as well as precession.
826  //
827 
828  // Compute days since J2000
829  int16_t iday = day;
830  double fday = day - iday;
831  int jd = jday(iyr, 1, iday);
832  double t = jd - (double)2451545.5 + fday;
833 
834  // Compute Greenwich Mean Sidereal Time (degrees)
835  double gmst = (double)100.4606184 + (double)0.9856473663 * t +
836  (double)2.908e-13 * t * t;
837 
838  // Check if need to compute nutation correction for this day
839  double xls, gs, xlm, omega;
840  double dpsi, eps, epsm;
841  ephparms(t, xls, gs, xlm, omega);
842  nutate(t, xls, gs, xlm, omega, dpsi, eps, epsm);
843 
844  // Include apparent time correction and time-of-day
845  gha = gmst + dpsi * cos(eps / RADEG) + fday * 360;
846  gha = fmod(gha, (double)360);
847 
848  return 0;
849 }
850 
851 
852 int expandEnvVar(std::string *sValue) {
853  if ( (*sValue).find_first_of( "$" ) == string::npos) return 0;
854  string::size_type posEndIdx = (*sValue).find_first_of( "/" );
855  if ( posEndIdx == string::npos) return 0;
856  char *envVar_str = getenv((*sValue).substr( 1, posEndIdx-1 ).c_str());
857  if (envVar_str == 0x0) {
858  printf("Environment variable: %s not defined.\n", envVar_str);
859  exit(1);
860  }
861  *sValue = envVar_str + (*sValue).substr( posEndIdx);
862 
863  return 0;
864 }
865 
866 
867 
868 
869 
870 
871 
872 
873 
874 
875 
876 
877 
878 
879 
880 
881 
882 
883 
884 
885 
886 
887 
888 
889 
890 
891 
892 
893 
894 
895 
896 
897 
898 
899 
900 
901 
902 
903 
904 
905 
906 
907 
908 
909 
910 
911 
912 
913 
914 
915 
916 
917 
918 
919 
920 
921 
922 
data_t t[NROOTS+1]
Definition: decode_rs.h:77
int j
Definition: decode_rs.h:73
int32_t day
int j2000_to_mod(int32_t iyr, int32_t idy, double sec, double j2mod[3][3])
int latlon_interp_dist(size_t num_gridlines, double *dist_ai, float *lati, float *loni, double *dist, float *lati2, float *loni2)
int orb_interp2(size_t n_orb_rec, size_t sdim, double *torb, orb_array2 *p, orb_array2 *v, double *time, orb_array2 *posi, orb_array2 *veli)
int latlon_interp_vec(size_t n_orb_rec, size_t num_gridlines, double *torb, double *latorb, double *lonorb, double *tmgv, float *lati, float *loni)
float orb_array[3]
const double C
int ephparms(double t, double &xls, double &gs, double &xlm, double &omega)
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
float * lat
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
@ string
int nutate(double t, double xls, double gs, double xlm, double omega, double &dpsi, double &eps, double &epsm)
int expandEnvVar(std::string *sValue)
int latlon_interp2(size_t num_gridlines, size_t num_pixels, double *gtime, float *lati, float *loni, double *gtime2, float *lati2, float *loni2)
int j2000_to_ecr(int32_t iyr, int32_t idy, double sec, double ecmat[3][3])
int get_nut(int32_t iyr, int32_t idy, double xnut[3][3])
data_t omega[NROOTS+1]
Definition: decode_rs.h:77
int get_ut1(int32_t iyr, int32_t idy, double &ut1utc)
Definition: jd.py:1
int gha2000(int32_t iyr, double day, double &gha)
integer, parameter double
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)
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
int latlon_interp(size_t n_orb_rec, size_t num_gridlines, size_t num_pixels, double *torb, float **lat, float **lon, double *time, float *lati, float *loni)
float * lon
int32_t idy
Definition: atrem_corl1.h:161
#define abs(a)
Definition: misc.h:90
int i
Definition: decode_rs.h:71
int32_t iyr
Definition: atrem_corl1.h:161
int latlon_interp1pix(size_t num_gridlines, size_t gd_index, double *tgrid, float *lat_nad, float *lon_nad, float *latipix, float *lonipix)
#define RADEG
int k
Definition: decode_rs.h:73
float p[MODELMAX]
Definition: atrem_corl1.h:131
double orb_array2[3]