OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
mtl_geometry.c
Go to the documentation of this file.
1 /*
2  * Name: mtl_geometry.c
3  *
4  * Purpose:
5  * Source file containing the functions used to perform geometric operations
6  * on the data extracted from the MTL file.
7  */
8 
9 #include "mtl_geometry.h"
10 
11 /*
12  * Name: calc_los
13  *
14  * Purpose:
15  * Source file containing the function used to compute OLI and TIRS line-of-sight
16  * vectors and transform them to the spacecraft coordinate system.
17  */
18 
19 int calc_los(
20  SMETA_BAND band_smeta, /* I: Metadata band structure */
21  int sca, /* I: SCA number */
22  double ndet, /* I: Normalized detector coordinate */
23  VECTOR *los) /* O: Line of sight unit vector */ {
24  VECTOR ilos; /* Instrument LOS */
25  VECTOR slos; /* Spacecraft LOS */
26 
27  /* Check SCA number */
28  if (sca < 1 || sca > band_smeta.nsca) {
29  printf("Invalid SCA %d requested.\n", sca);
30  return (ERROR);
31  }
32 
33  /* Calculate the instrument LOS vector */
34  ilos.x = band_smeta.legendre[sca - 1][0][0]
35  + band_smeta.legendre[sca - 1][0][1] * ndet
36  + band_smeta.legendre[sca - 1][0][2] * (1.5 * ndet * ndet - 0.5)
37  + band_smeta.legendre[sca - 1][0][3] * ndet * (2.5 * ndet * ndet - 1.5);
38  ilos.y = band_smeta.legendre[sca - 1][1][0]
39  + band_smeta.legendre[sca - 1][1][1] * ndet
40  + band_smeta.legendre[sca - 1][1][2] * (1.5 * ndet * ndet - 0.5)
41  + band_smeta.legendre[sca - 1][1][3] * ndet * (2.5 * ndet * ndet - 1.5);
42  ilos.z = 1.0;
43  unitvec(ilos, &slos);
44 
45  /* Rotate to spacecraft coordinates */
46  rotatevec(band_smeta.align, slos, los);
47 
48  return (SUCCESS);
49 }
50 
51 /* Different varant for older LS sensors */
52 int calc_los_ls(
53  SMETA_BAND_LS band_smeta, /* I: Metadata band structure */
54  double frow, /* I: Fractional row */
55  double npix, /* I: Sample coordinate */
56  double ndet, /* I: Normalized detector coordinate */
57  VECTOR *los) /* O: Line of sight unit vector */ {
58  VECTOR ilos; /* Instrument LOS */
59  VECTOR slos; /* Spacecraft LOS */
60  double along_angle; /* Along scan angle */
61  double across_angle; /* Across scan angle */
62  //double scan_number; /* Current scan number */
63 
64  /* Calculate the instrument LOS vector */
65  //scan_number = band_smeta.nscans * (1.0 + ndet)/2.0;
66 
67  across_angle = 0.0;
68  along_angle = -band_smeta.fascan + (1.0 + npix) / 2.0 * (band_smeta.fascan + band_smeta.rascan) - band_smeta.aoffset;
69 
70  ilos.x = sin(across_angle);
71  ilos.y = -cos(across_angle) * sin(along_angle);
72  ilos.z = cos(across_angle) * cos(along_angle);
73  unitvec(ilos, &slos);
74 
75  /* Rotate to spacecraft coordinates */
76  rotatevec(band_smeta.align, slos, los);
77 
78  return (SUCCESS);
79 }
80 
81 /******************************************************************************
82 NAME: calc_yaw_steering_gp
83 
84 PURPOSE:
85 Calculates the ground point intersection latitude and longitude for a specified
86 input instrument LOS vector, spacecraft ephemeris, and roll angle. Uses the
87 yaw steering pointing logic to do the calculation.
88 
89 RETURN VALUE:
90 None.
91 
92 NOTES:
93 
94 ALGORITHM REFERENCES:
95  ******************************************************************************
96  Property of the U.S. Government
97  USGS EROS Data Center
98  ******************************************************************************/
99 
100 int calc_yaw_steering_gp(
101  WRS2 parms, /* I: WRS2 system parameters */
102  VECTOR pos, /* I: Position vector */
103  VECTOR vel, /* I: Velocity vector */
104  VECTOR i_los, /* I: Instrument line-of-sight unit vector */
105  double roll, /* I: Off-nadir roll angle (in degrees) */
106  double *gp_lat, /* O: Ground point geodetic latitude (in degrees) */
107  double *gp_lon) /* O: Ground point geodetic longitude (in degrees) */ {
108  double a[3][3];
109  VECTOR omega_los, /* LOS angular rate vector */
110  u_los, /* LOS vector */
111  o_x_u, /* angular rate cross LOS */
112  rp, /* Earth point vector */
113  vp, /* Earth point velocity */
114  er, /* Earth rotation vector */
115  er_x_rp, /* Earth rotation cross position */
116  vgt, /* Target velocity */
117  u_x_vgt, /* LOS cross target velocity */
118  x, /* YSF X axis */
119  y; /* YSF Y axis */
120  double e_a, /* Earth semimajor axis */
121  e_b, /* Earth semiminor axis */
122  e_rate, /* Earth rotation rate */
123  d2r, /* Degrees to radians conversion */
124  scale, /* Orbit to Earth scale factor */
125  rsc, /* Orbital radius */
126  rho, /* Orbital altitude */
127  rho_dot, /* Radial velocity */
128  angmo, /* Scalar angular momentum */
129  sang, /* Rescaled (for ellipsoid) roll angle */
130  cosr, /* Cosine of roll angle */
131  sinr; /* Sine of roll angle */
132 
133  /* Initialize Earth constants */
134  e_a = parms.semimajor;
135  e_b = parms.semiminor;
136  e_rate = parms.inertialvel;
137  d2r = atan(1.0) / 45.0;
138 
139  /* Compute roll trig functions */
140  cosr = cos(roll * d2r);
141  sinr = sin(roll * d2r);
142 
143  /* Compute orbital radius and angular momentum */
144  rsc = vecnorm(pos);
145  crossprod(pos, vel, &omega_los);
146  angmo = vecnorm(omega_los);
147 
148  /* Construct the LOS */
149  u_los.x = -cosr * pos.x / rsc + sinr * omega_los.x / angmo;
150  u_los.y = -cosr * pos.y / rsc + sinr * omega_los.y / angmo;
151  u_los.z = -cosr * pos.z / rsc + sinr * omega_los.z / angmo;
152 
153  /* Compute LOS angular rate */
154  omega_los.x = cosr * angmo / (rsc * rsc)*(cosr * omega_los.x / angmo + sinr * pos.x / rsc);
155  omega_los.y = cosr * angmo / (rsc * rsc)*(cosr * omega_los.y / angmo + sinr * pos.y / rsc);
156  omega_los.z = cosr * angmo / (rsc * rsc)*(cosr * omega_los.z / angmo + sinr * pos.z / rsc);
157 
158  /* Compute orbit to Earth scale factor */
159  scale = pos.x / e_a * pos.x / e_a + pos.y / e_a * pos.y / e_a + pos.z / e_b * pos.z / e_b;
160  sang = u_los.x / e_a * u_los.x / e_a + u_los.y / e_a * u_los.y / e_a + u_los.z / e_b * u_los.z / e_b;
161  sang = (-pos.x / e_a * u_los.x / e_a - pos.y / e_a * u_los.y / e_a - pos.z / e_b * u_los.z / e_b) / sqrt(sang * scale);
162  sang = MIN(MAX(-1.0, sang), 1.0);
163  sang = acos(sang);
164  cosr = cos(sang);
165  sinr = sin(sang);
166  scale = 1.0 / scale - sinr*sinr;
167  if (scale > 0.0) {
168  scale = sqrt(scale);
169  } else {
170  printf("Roll angle too large, LOS misses Earth\n");
171  return (-1);
172  }
173  rho = rsc * (cosr - scale);
174  rho_dot = dotprod(pos, vel) / rsc * (cosr + sinr * sinr / scale);
175 
176  /* LOS rate cross LOS */
177  crossprod(omega_los, u_los, &o_x_u);
178 
179  /* Construct ground target point */
180  rp.x = pos.x + rho * u_los.x;
181  rp.y = pos.y + rho * u_los.y;
182  rp.z = pos.z + rho * u_los.z;
183  vp.x = vel.x + rho_dot * u_los.x + rho * o_x_u.x;
184  vp.y = vel.y + rho_dot * u_los.y + rho * o_x_u.y;
185  vp.z = vel.z + rho_dot * u_los.z + rho * o_x_u.z;
186 
187  /* Construct the Earth rotation vector */
188  er.x = 0.0;
189  er.y = 0.0;
190  er.z = e_rate;
191 
192  /* Compute er cross rp and ground target velocity */
193  crossprod(er, rp, &er_x_rp);
194  vgt.x = vp.x - er_x_rp.x;
195  vgt.y = vp.y - er_x_rp.y;
196  vgt.z = vp.z - er_x_rp.z;
197  crossprod(u_los, vgt, &u_x_vgt);
198  unitvec(u_x_vgt, &y);
199  crossprod(y, u_los, &x);
200 
201  /* Construct ECI2YSF orientation matrix */
202  a[0][0] = x.x;
203  a[0][1] = x.y;
204  a[0][2] = x.z;
205  a[1][0] = y.x;
206  a[1][1] = y.y;
207  a[1][2] = y.z;
208  a[2][0] = u_los.x;
209  a[2][1] = u_los.y;
210  a[2][2] = u_los.z;
211 
212  /* Rotate input LOS vector using the orientation matrix */
213  u_los.x = a[0][0] * i_los.x + a[1][0] * i_los.y + a[2][0] * i_los.z;
214  u_los.y = a[0][1] * i_los.x + a[1][1] * i_los.y + a[2][1] * i_los.z;
215  u_los.z = a[0][2] * i_los.x + a[1][2] * i_los.y + a[2][2] * i_los.z;
216 
217  /* Compute orbit to Earth scale factor */
218  scale = pos.x / e_a * pos.x / e_a + pos.y / e_a * pos.y / e_a + pos.z / e_b * pos.z / e_b;
219  sang = u_los.x / e_a * u_los.x / e_a + u_los.y / e_a * u_los.y / e_a + u_los.z / e_b * u_los.z / e_b;
220  sang = (-pos.x / e_a * u_los.x / e_a - pos.y / e_a * u_los.y / e_a - pos.z / e_b * u_los.z / e_b) / sqrt(sang * scale);
221  sang = MIN(MAX(-1.0, sang), 1.0);
222  sang = acos(sang);
223  cosr = cos(sang);
224  sinr = sin(sang);
225  scale = 1.0 / scale - sinr*sinr;
226  if (scale > 0.0) {
227  scale = sqrt(scale);
228  } else {
229  printf("Roll angle too large, LOS misses Earth");
230  return (-1);
231  }
232  rho = rsc * (cosr - scale);
233 
234  /* Construct ground target point */
235  rp.x = pos.x + rho * u_los.x;
236  rp.y = pos.y + rho * u_los.y;
237  rp.z = pos.z + rho * u_los.z;
238 
239  /* Convert ground point to lat/lon */
240  *gp_lat = atan(rp.z / sqrt(rp.x * rp.x + rp.y * rp.y) * e_a / e_b * e_a / e_b) / d2r;
241  *gp_lon = atan2(rp.y, rp.x) / d2r;
242 
243  return (0);
244 }
245 
246 /******************************************************************************
247 NAME: calc_yaw_steering_gp_ls
248 
249 PURPOSE:
250 Calculates the ground point intersection latitude and longitude for a specified
251 input instrument LOS vector, spacecraft ephemeris, and roll angle. Uses the
252 yaw steering pointing logic to do the calculation.
253 
254 RETURN VALUE:
255 None.
256 
257 NOTES:
258 
259 ALGORITHM REFERENCES:
260  ******************************************************************************
261  Property of the U.S. Government
262  USGS EROS Data Center
263  ******************************************************************************/
264 
265 int calc_yaw_steering_gp_ls(
266  WRS2 parms, /* I: WRS2 system parameters */
267  VECTOR pos, /* I: Position vector */
268  VECTOR vel, /* I: Velocity vector */
269  VECTOR i_los, /* I: Instrument line-of-sight unit vector */
270  double roll, /* I: Off-nadir roll angle (in degrees) */
271  double *gp_lat, /* O: Ground point geodetic latitude (in degrees) */
272  double *gp_lon) /* O: Ground point geodetic longitude (in degrees) */ {
273  double a[3][3];
274  VECTOR omega_los, /* LOS angular rate vector */
275  u_los, /* LOS vector */
276  o_x_u, /* angular rate cross LOS */
277  rp, /* Earth point vector */
278  vp, /* Earth point velocity */
279  er, /* Earth rotation vector */
280  er_x_rp, /* Earth rotation cross position */
281  vgt, /* Target velocity */
282  u_x_vgt, /* LOS cross target velocity */
283  x, /* YSF X axis */
284  y; /* YSF Y axis */
285  double e_a, /* Earth semimajor axis */
286  e_b, /* Earth semiminor axis */
287  e_rate, /* Earth rotation rate */
288  d2r, /* Degrees to radians conversion */
289  scale, /* Orbit to Earth scale factor */
290  rsc, /* Orbital radius */
291  rho, /* Orbital altitude */
292  rho_dot, /* Radial velocity */
293  angmo, /* Scalar angular momentum */
294  sang, /* Rescaled (for ellipsoid) roll angle */
295  cosr, /* Cosine of roll angle */
296  sinr; /* Sine of roll angle */
297 
298  /* Initialize Earth constants */
299  e_a = parms.semimajor;
300  e_b = parms.semiminor;
301  e_rate = parms.inertialvel;
302  d2r = atan(1.0) / 45.0;
303 
304  /* Compute roll trig functions */
305  cosr = cos(roll * d2r);
306  sinr = sin(roll * d2r);
307 
308  /* Compute orbital radius and angular momentum */
309  rsc = vecnorm(pos);
310  crossprod(pos, vel, &omega_los);
311  angmo = vecnorm(omega_los);
312 
313  /* Construct the LOS */
314  u_los.x = -cosr * pos.x / rsc + sinr * omega_los.x / angmo;
315  u_los.y = -cosr * pos.y / rsc + sinr * omega_los.y / angmo;
316  u_los.z = -cosr * pos.z / rsc + sinr * omega_los.z / angmo;
317 
318  /* Compute LOS angular rate */
319  omega_los.x = cosr * angmo / (rsc * rsc)*(cosr * omega_los.x / angmo + sinr * pos.x / rsc);
320  omega_los.y = cosr * angmo / (rsc * rsc)*(cosr * omega_los.y / angmo + sinr * pos.y / rsc);
321  omega_los.z = cosr * angmo / (rsc * rsc)*(cosr * omega_los.z / angmo + sinr * pos.z / rsc);
322 
323  /* Compute orbit to Earth scale factor */
324  scale = pos.x / e_a * pos.x / e_a + pos.y / e_a * pos.y / e_a + pos.z / e_b * pos.z / e_b;
325  sang = u_los.x / e_a * u_los.x / e_a + u_los.y / e_a * u_los.y / e_a + u_los.z / e_b * u_los.z / e_b;
326  sang = (-pos.x / e_a * u_los.x / e_a - pos.y / e_a * u_los.y / e_a - pos.z / e_b * u_los.z / e_b) / sqrt(sang * scale);
327  sang = MIN(MAX(-1.0, sang), 1.0);
328  sang = acos(sang);
329  cosr = cos(sang);
330  sinr = sin(sang);
331  scale = 1.0 / scale - sinr*sinr;
332  if (scale > 0.0) {
333  scale = sqrt(scale);
334  } else {
335  printf("Roll angle too large, LOS misses Earth\n");
336  return (-1);
337  }
338  rho = rsc * (cosr - scale);
339  rho_dot = dotprod(pos, vel) / rsc * (cosr + sinr * sinr / scale);
340 
341  /* LOS rate cross LOS */
342  crossprod(omega_los, u_los, &o_x_u);
343 
344  /* Construct ground target point */
345  rp.x = pos.x + rho * u_los.x;
346  rp.y = pos.y + rho * u_los.y;
347  rp.z = pos.z + rho * u_los.z;
348  vp.x = vel.x + rho_dot * u_los.x + rho * o_x_u.x;
349  vp.y = vel.y + rho_dot * u_los.y + rho * o_x_u.y;
350  vp.z = vel.z + rho_dot * u_los.z + rho * o_x_u.z;
351 
352  /* Construct the Earth rotation vector */
353  er.x = 0.0;
354  er.y = 0.0;
355  er.z = e_rate;
356 
357  /* Compute er cross rp and ground target velocity */
358  crossprod(er, rp, &er_x_rp);
359  vgt.x = vp.x; // - er_x_rp.x;
360  vgt.y = vp.y; // - er_x_rp.y;
361  vgt.z = vp.z; // - er_x_rp.z;
362  crossprod(u_los, vgt, &u_x_vgt);
363  unitvec(u_x_vgt, &y);
364  crossprod(y, u_los, &x);
365 
366  /* Construct ECI2YSF orientation matrix */
367  a[0][0] = x.x;
368  a[0][1] = x.y;
369  a[0][2] = x.z;
370  a[1][0] = y.x;
371  a[1][1] = y.y;
372  a[1][2] = y.z;
373  a[2][0] = u_los.x;
374  a[2][1] = u_los.y;
375  a[2][2] = u_los.z;
376 
377  /* Rotate input LOS vector using the orientation matrix */
378  u_los.x = a[0][0] * i_los.x + a[1][0] * i_los.y + a[2][0] * i_los.z;
379  u_los.y = a[0][1] * i_los.x + a[1][1] * i_los.y + a[2][1] * i_los.z;
380  u_los.z = a[0][2] * i_los.x + a[1][2] * i_los.y + a[2][2] * i_los.z;
381 
382  /* Compute orbit to Earth scale factor */
383  scale = pos.x / e_a * pos.x / e_a + pos.y / e_a * pos.y / e_a + pos.z / e_b * pos.z / e_b;
384  sang = u_los.x / e_a * u_los.x / e_a + u_los.y / e_a * u_los.y / e_a + u_los.z / e_b * u_los.z / e_b;
385  sang = (-pos.x / e_a * u_los.x / e_a - pos.y / e_a * u_los.y / e_a - pos.z / e_b * u_los.z / e_b) / sqrt(sang * scale);
386  sang = MIN(MAX(-1.0, sang), 1.0);
387  sang = acos(sang);
388  cosr = cos(sang);
389  sinr = sin(sang);
390  scale = 1.0 / scale - sinr*sinr;
391  if (scale > 0.0) {
392  scale = sqrt(scale);
393  } else {
394  printf("Roll angle too large, LOS misses Earth");
395  return (-1);
396  }
397  rho = rsc * (cosr - scale);
398 
399  /* Construct ground target point */
400  rp.x = pos.x + rho * u_los.x;
401  rp.y = pos.y + rho * u_los.y;
402  rp.z = pos.z + rho * u_los.z;
403 
404  /* Convert ground point to lat/lon */
405  *gp_lat = atan(rp.z / sqrt(rp.x * rp.x + rp.y * rp.y) * e_a / e_b * e_a / e_b) / d2r;
406  *gp_lon = atan2(rp.y, rp.x) / d2r;
407 
408  return (0);
409 }
410 
411 /*
412  * Name: frame_scene.c
413  *
414  * Purpose: Calculate the nominal scene frame and adjust the longitude
415  * to fit the observed metadata frame.
416  *
417  */
418 
419 int frame_scene(
420  WRS2 wrsparms, /* I: Structure of WRS-2 parameters */
421  SMETA smeta, /* I: Input scene metadata structure */
422  double *dellon, /* O: WRS longitude offset */
423  double *delrow) /* O: Scene length in fractional WRS rows */ {
424  int band; /* Current band number */
425  int band_index; /* Index of current band in metadata structure */
426  int sca; /* Current SCA number */
427  double ndet; /* Normalized detector coordinate */
428  double start_row; /* Fractional row at scene start */
429  double stop_row; /* Fractional row at scene end */
430  double max_line; /* Maximum product line number */
431  double max_samp; /* Maximum product sample number */
432  double l1t_line[2]; /* Ground point L1T line coordinate */
433  double l1t_samp[2]; /* Ground point L1t sample coordinate */
434  double ul_ls[2]; /* Upper left corner line/sample */
435  double ur_ls[2]; /* Upper right corner line/sample */
436  double lr_ls[2]; /* Lower right corner line/sample */
437  double ll_ls[2]; /* Lower left corner line/sample */
438  double dr_right; /* Scene size on right side */
439  double dr_left; /* Scene size on left side */
440  double sc_lat; /* Scene center latitude */
441  double sc_lon; /* Scene center longitude */
442  double dlon; /* Longitude offset */
443  double radius; /* Earth radius in parallel of latitude */
444  double e2; /* Earth eccentricity squared */
445  double d2r; /* Conversion from degrees to radians */
446 
447  /* Compute conversion constant */
448  d2r = atan(1.0) / 45.0;
449 
450  /* Construct the nominal scene frame */
451  /* Find the upper left corner */
452  band = 9;
453  sca = 1;
454  ndet = -1.0;
455  /* Sudipta new addition for OLI to extract SCA and Det numbers */
456  if ((band_index = smeta_band_number_to_index(&smeta, band)) < 0) {
457  band = 10;
458  sca = 3;
459  ndet = 1.0; /* TIRS only */
460  }
461  if ((band_index = smeta_band_number_to_index(&smeta, band)) < 0) {
462  printf("Scene framing bands 9 and 10 are not available.\n");
463  return (ERROR);
464  }
465  /* End Sudipta new addition for OLI */
466  start_row = (double) smeta.wrs_row - 0.7;
467  if (project_point(wrsparms, smeta, *dellon, band, sca, start_row, ndet, &l1t_line[0], &l1t_samp[0]) != SUCCESS) {
468  printf("Error projecting scene UL corner to L1T line/sample.\n");
469  smeta_release_projection();
470  return (ERROR);
471  }
472  start_row = (double) smeta.wrs_row - 0.6;
473  if (project_point(wrsparms, smeta, *dellon, band, sca, start_row, ndet, &l1t_line[1], &l1t_samp[1]) != SUCCESS) {
474  printf("Error projecting scene UL corner to L1T line/sample.\n");
475  smeta_release_projection();
476  return (ERROR);
477  }
478  start_row = (double) smeta.wrs_row - (0.7 * l1t_line[1] - 0.6 * l1t_line[0]) / (l1t_line[1] - l1t_line[0]);
479  if (project_point(wrsparms, smeta, *dellon, band, sca, start_row, ndet, &ul_ls[0], &ul_ls[1]) != SUCCESS) {
480  printf("Error projecting scene UL corner to L1T line/sample.\n");
481  smeta_release_projection();
482  return (ERROR);
483  }
484 
485  /* Find the lower right corner */
486  band = 9;
487  sca = 14;
488  ndet = 1.0;
489  band_index = smeta_band_number_to_index(&smeta, band);
490  max_line = (double) (smeta.band_smeta[band_index].l1t_lines - 1);
491  max_samp = (double) (smeta.band_smeta[band_index].l1t_samps - 1);
492 
493  stop_row = (double) smeta.wrs_row + 0.6;
494  if (project_point(wrsparms, smeta, *dellon, band, sca, stop_row, ndet, &l1t_line[0], &l1t_samp[0]) != SUCCESS) {
495  printf("Error projecting scene LR corner to L1T line/sample.\n");
496  smeta_release_projection();
497  return (ERROR);
498  }
499  stop_row = (double) smeta.wrs_row + 0.7;
500  if (project_point(wrsparms, smeta, *dellon, band, sca, stop_row, ndet, &l1t_line[1], &l1t_samp[1]) != SUCCESS) {
501  printf("Error projecting scene LR corner to L1T line/sample.\n");
502  smeta_release_projection();
503  return (ERROR);
504  }
505  stop_row = (double) smeta.wrs_row + (0.6 * (l1t_line[1] - max_line) + 0.7 * (max_line - l1t_line[0])) / (l1t_line[1] - l1t_line[0]);
506  if (project_point(wrsparms, smeta, *dellon, band, sca, stop_row, ndet, &lr_ls[0], &lr_ls[1]) != SUCCESS) {
507  printf("Error projecting scene LR corner to L1T line/sample.\n");
508  smeta_release_projection();
509  return (ERROR);
510  }
511 
512  /* Find provisional upper right corner */
513  if (project_point(wrsparms, smeta, *dellon, band, sca, start_row, ndet, &ur_ls[0], &ur_ls[1]) != SUCCESS) {
514  printf("Error projecting scene UR corner to L1T line/sample.\n");
515  smeta_release_projection();
516  return (ERROR);
517  }
518 
519  /* Find the corner of the last odd SCA */
520  sca = 13;
521  if (project_point(wrsparms, smeta, *dellon, band, sca, start_row, ndet, &l1t_line[0], &l1t_samp[0]) != SUCCESS) {
522  printf("Error projecting scene UR corner to L1T line/sample.\n");
523  smeta_release_projection();
524  return (ERROR);
525  }
526 
527  /* Calculate the size of the scene, in rows */
528  dr_right = (start_row - stop_row)
529  * ((ur_ls[0] - l1t_line[0])*(ul_ls[1] - l1t_samp[0]) - (ur_ls[1] - l1t_samp[0])*(ul_ls[0] - l1t_line[0]))
530  / ((ur_ls[1] - lr_ls[1])*(ul_ls[0] - l1t_line[0]) - (ur_ls[0] - lr_ls[0])*(ul_ls[1] - l1t_samp[0]));
531  dr_right = stop_row - start_row - dr_right;
532 
533  /* Locate the provisional LL corner */
534  sca = 1;
535  ndet = -1.0;
536  if (project_point(wrsparms, smeta, *dellon, band, sca, stop_row, ndet, &ll_ls[0], &ll_ls[1]) != SUCCESS) {
537  printf("Error projecting scene LL corner to L1T line/sample.\n");
538  smeta_release_projection();
539  return (ERROR);
540  }
541 
542  /* Find the corner of the first even SCA */
543  sca = 2;
544  if (project_point(wrsparms, smeta, *dellon, band, sca, stop_row, ndet, &l1t_line[0], &l1t_samp[0]) != SUCCESS) {
545  printf("Error projecting scene LL corner to L1T line/sample.\n");
546  smeta_release_projection();
547  return (ERROR);
548  }
549 
550  /* Calculate the size of the scene, in rows */
551  dr_left = (start_row - stop_row)
552  * ((ll_ls[0] - l1t_line[0])*(lr_ls[1] - l1t_samp[0]) - (ll_ls[1] - l1t_samp[0])*(lr_ls[0] - l1t_line[0]))
553  / ((ll_ls[1] - ul_ls[1])*(lr_ls[0] - l1t_line[0]) - (ll_ls[0] - ul_ls[0])*(lr_ls[1] - l1t_samp[0]));
554  dr_left = stop_row - start_row - dr_left;
555 
556  /* Combine the results */
557  *delrow = (dr_left + dr_right) / 2.0;
558 
559  /* Locate the UR corner from the LR corner */
560  sca = 14;
561  ndet = 1.0;
562  if (project_point(wrsparms, smeta, *dellon, band, sca, stop_row - (*delrow), ndet, &ur_ls[0], &ur_ls[1]) != SUCCESS) {
563  printf("Error projecting scene UR corner to L1T line/sample.\n");
564  smeta_release_projection();
565  return (ERROR);
566  }
567 
568  /* Locate the LL corner from the UL corner */
569  sca = 1;
570  ndet = -1.0;
571  if (project_point(wrsparms, smeta, *dellon, band, sca, start_row + (*delrow), ndet, &ll_ls[0], &ll_ls[1]) != SUCCESS) {
572  printf("Error projecting scene LL corner to L1T line/sample.\n");
573  smeta_release_projection();
574  return (ERROR);
575  }
576 
577  /* Calculate longitude offset in meters */
578  dlon = smeta.band_smeta[band_index].pixsize * (max_samp - MAX(MAX(ul_ls[1], ur_ls[1]), MAX(ll_ls[1], lr_ls[1]))
579  - MIN(MIN(ul_ls[1], ur_ls[1]), MIN(ll_ls[1], lr_ls[1]))) / 2.0;
580 
581  /* Convert to longitude offset in degrees */
582  pathrow_to_latlon(wrsparms, smeta.wrs_path, (double) smeta.wrs_row, &sc_lat, &sc_lon);
583  e2 = 1.0 - wrsparms.semiminor / wrsparms.semimajor * wrsparms.semiminor / wrsparms.semimajor;
584  radius = cos(d2r * sc_lat) * wrsparms.semimajor / sqrt(1.0 - e2 * sin(d2r * sc_lat) * sin(d2r * sc_lat));
585  *dellon = *dellon + dlon / radius / d2r;
586 
587  return (SUCCESS);
588 }
589 
590 /* another variant for frame_scene for older LS */
591 int frame_scene_ls(
592  WRS2 wrsparms, /* I: Structure of WRS-2 parameters */
593  SMETA smeta, /* I: Input scene metadata structure */
594  double *dellon, /* O: WRS longitude offset */
595  double *delrow) /* O: Scene length in fractional WRS rows */ {
596  int band; /* Current band number */
597  int band_index; /* Index of current band in metadata structure */
598  double nsamp; /* Sample coordinate */
599  double nline; /* Line coordinate */
600  double start_row; /* Fractional row at scene start */
601  double stop_row; /* Fractional row at scene end */
602  double max_line; /* Maximum product line number */
603  double max_samp; /* Maximum product sample number */
604  double l1t_line[2]; /* Ground point L1T line coordinate */
605  double l1t_samp[2]; /* Ground point L1t sample coordinate */
606  double ul_ls[2]; /* Upper left corner line/sample */
607  double ur_ls[2]; /* Upper right corner line/sample */
608  double lr_ls[2]; /* Lower right corner line/sample */
609  double ll_ls[2]; /* Lower left corner line/sample */
610  double sc_lat; /* Scene center latitude */
611  double sc_lon; /* Scene center longitude */
612  double dlon; /* Longitude offset */
613  double radius; /* Earth radius in parallel of latitude */
614  double e2; /* Earth eccentricity squared */
615  double d2r; /* Conversion from degrees to radians */
616 
617  /* Compute conversion constant */
618  d2r = atan(1.0) / 45.0;
619 
620  /* Construct the nominal scene frame */
621  /* Find the upper left corner */
622  /* if( smeta.sensor_id == IAS_MSS &&
623  (smeta.spacecraft_id == IAS_L1 || smeta.spacecraft_id == IAS_L2 || smeta.spacecraft_id == IAS_L3) )
624  band = 4;
625  else if( smeta.sensor_id == IAS_MSS &&
626  (smeta.spacecraft_id == IAS_L4 || smeta.spacecraft_id == IAS_L5) )
627  band = 1;
628  else if( (smeta.spacecraft_id == IAS_L7) || (smeta.spacecraft_id == IAS_L5) ) */
629  band = 4;
630  /* else if( smeta.spacecraft_id == IAS_L4 )
631  band = 0;
632  else
633  {
634  printf("Error in sensor id.\n");
635  smeta_release_projection();
636  return(ERROR);
637  } */
638 
639  nsamp = -1.0;
640  nline = -1.0;
641  start_row = (double) smeta.wrs_row - 0.6;
642  if (project_point_ls(wrsparms, smeta, *dellon, band, start_row, nsamp, nline, &l1t_line[0], &l1t_samp[0]) != SUCCESS) {
643  printf("Error projecting scene UL corner to L1T line/sample.\n");
644  smeta_release_projection();
645  return (ERROR);
646  }
647  start_row = (double) smeta.wrs_row - 0.5;
648  if (project_point_ls(wrsparms, smeta, *dellon, band, start_row, nsamp, nline, &l1t_line[1], &l1t_samp[1]) != SUCCESS) {
649  printf("Error projecting scene UL corner to L1T line/sample.\n");
650  smeta_release_projection();
651  return (ERROR);
652  }
653  start_row = (double) smeta.wrs_row - (0.6 * l1t_line[1] - 0.5 * l1t_line[0]) / (l1t_line[1] - l1t_line[0]);
654  // start_row = (double)smeta.wrs_row - (0.8*l1t_line[1] - 0.5*l1t_line[0])/(l1t_line[1] - l1t_line[0]);
655  if (project_point_ls(wrsparms, smeta, *dellon, band, start_row, nsamp, nline, &ul_ls[0], &ul_ls[1]) != SUCCESS) {
656  printf("Error projecting scene UL corner to L1T line/sample.\n");
657  smeta_release_projection();
658  return (ERROR);
659  }
660 
661  /* Find the lower right corner */
662  nsamp = +1.0;
663  nline = +1.0;
664  band_index = smeta_band_number_to_index(&smeta, band);
665  max_line = (double) (smeta.band_smeta_ls[band_index].l1t_lines - 1);
666  max_samp = (double) (smeta.band_smeta_ls[band_index].l1t_samps - 1);
667 
668  stop_row = (double) smeta.wrs_row + 0.5;
669  if (project_point_ls(wrsparms, smeta, *dellon, band, stop_row, nsamp, nline, &l1t_line[0], &l1t_samp[0]) != SUCCESS) {
670  printf("Error projecting scene LR corner to L1T line/sample.\n");
671  smeta_release_projection();
672  return (ERROR);
673  }
674  stop_row = (double) smeta.wrs_row + 0.6;
675  if (project_point_ls(wrsparms, smeta, *dellon, band, stop_row, nsamp, nline, &l1t_line[1], &l1t_samp[1]) != SUCCESS) {
676  printf("Error projecting scene LR corner to L1T line/sample.\n");
677  smeta_release_projection();
678  return (ERROR);
679  }
680  stop_row = (double) smeta.wrs_row + (0.5 * (l1t_line[1] - max_line) + 0.6 * (max_line - l1t_line[0])) / (l1t_line[1] - l1t_line[0]);
681  // stop_row = (double)smeta.wrs_row + (0.5*(l1t_line[1] - max_line) + 0.8*(max_line - l1t_line[0]))/(l1t_line[1] - l1t_line[0]);
682  if (project_point_ls(wrsparms, smeta, *dellon, band, stop_row, nsamp, nline, &lr_ls[0], &lr_ls[1]) != SUCCESS) {
683  printf("Error projecting scene LR corner to L1T line/sample.\n");
684  smeta_release_projection();
685  return (ERROR);
686  }
687 
688  /* Find provisional upper right corner */
689  nsamp = +1.0;
690  nline = -1.0;
691  if (project_point_ls(wrsparms, smeta, *dellon, band, start_row, nsamp, nline, &ur_ls[0], &ur_ls[1]) != SUCCESS) {
692  printf("Error projecting scene UR corner to L1T line/sample.\n");
693  smeta_release_projection();
694  return (ERROR);
695  }
696 
697  /* Locate the provisional LL corner */
698  nsamp = -1.0;
699  nline = +1.0;
700  if (project_point_ls(wrsparms, smeta, *dellon, band, stop_row, nsamp, nline, &ll_ls[0], &ll_ls[1]) != SUCCESS) {
701  printf("Error projecting scene LL corner to L1T line/sample.\n");
702  smeta_release_projection();
703  return (ERROR);
704  }
705 
706  /* Combine the results */
707  *delrow = stop_row - start_row;
708 
709  /* Calculate longitude offset in meters */
710  dlon = smeta.band_smeta_ls[band_index].pixsize * (max_samp - MAX(MAX(ul_ls[1], ur_ls[1]), MAX(ll_ls[1], lr_ls[1]))
711  - MIN(MIN(ul_ls[1], ur_ls[1]), MIN(ll_ls[1], lr_ls[1]))) / 2.0;
712 
713  /* Convert to longitude offset in degrees */
714  pathrow_to_latlon(wrsparms, smeta.wrs_path, (double) smeta.wrs_row, &sc_lat, &sc_lon);
715  e2 = 1.0 - wrsparms.semiminor / wrsparms.semimajor * wrsparms.semiminor / wrsparms.semimajor;
716  radius = cos(d2r * sc_lat) * wrsparms.semimajor / sqrt(1.0 - e2 * sin(d2r * sc_lat) * sin(d2r * sc_lat));
717  *dellon = *dellon + dlon / radius / d2r;
718 
719  return (SUCCESS);
720 }
721 
722 int project_point(
723  WRS2 wrsparms, /* I: Structure of WRS-2 parameters */
724  SMETA smeta, /* I: Input scene metadata structure */
725  double dellon, /* I: WRS longitude offset */
726  int band, /* I: Band number (user band) */
727  int sca, /* I: SCA number (1 to nsca) */
728  double frow, /* I: Fractional WRS row coordinate */
729  double ndet, /* I: Normalized detector coordinate */
730  double *l1t_line, /* O: Projected line coordinate */
731  double *l1t_samp) /* O: Projected sample coordinate */ {
732  VECTOR pos; /* Spacecraft ECEF position vector */
733  VECTOR vel; /* Spacecraft ECEF velocity vector */
734  VECTOR slos; /* Spacecraft line-of-sight vector */
735  int band_index; /* Index of current band in metadata structure */
736  double gp_lat; /* Ground point latitude */
737  double gp_lon; /* Ground point longitude */
738  double proj_x; /* Ground point projection X coordinate */
739  double proj_y; /* Ground point projection Y coordinate */
740  double d2r; /* Conversion from degrees to radians */
741 
742  /* Compute conversion constant */
743  d2r = atan(1.0) / 45.0;
744 
745  /* Find the nominal ephemeris state vector */
746  pathrow_to_posvel(wrsparms, smeta.wrs_path, dellon, frow, &pos, &vel);
747 
748  /* Find the current band in the metadata */
749  band_index = smeta_band_number_to_index(&smeta, band);
750 
751  /* Calculate the instrument LOS vector */
752  if (calc_los(smeta.band_smeta[band_index], sca, ndet, &slos) != SUCCESS) {
753  printf("Error generating LOS vector.\n");
754  return (ERROR);
755  }
756 
757  /* Project to ground latitude/longitude */
758  if (calc_yaw_steering_gp(wrsparms, pos, vel, slos, smeta.roll_angle, &gp_lat, &gp_lon) < 0) {
759  printf("Error projecting nominal scene center.\n");
760  return (ERROR);
761  }
762  gp_lat *= d2r;
763  gp_lon *= d2r;
764 
765  /* Apply scene map projection */
766  if (smeta_geodetic_to_proj(smeta.projection, gp_lat, gp_lon, &proj_x, &proj_y) != SUCCESS) {
767  printf("Error converting scene center lat/lon to projection X/Y.\n");
768  return (ERROR);
769  }
770 
771  /* Reduce to L1T line/sample */
772  if (smeta_proj_to_l1t(&smeta, band, proj_x, proj_y, l1t_line, l1t_samp) != SUCCESS) {
773  printf("Errro converting scene center X/Y to L1T line/sample.\n");
774  return (ERROR);
775  }
776 
777  return (SUCCESS);
778 }
779 
780 /* Another variant for project_point for older LS */
781 int project_point_ls(
782  WRS2 wrsparms, /* I: Structure of WRS-2 parameters */
783  SMETA smeta, /* I: Input scene metadata structure */
784  double dellon, /* I: WRS longitude offset */
785  int band, /* I: Band number (user band) */
786  double frow, /* I: Fractional WRS row coordinate */
787  double npix, /* I: Normalized detector coordinate */
788  double ndet, /* I: Line detector coordinate */
789  double *l1t_line, /* O: Projected line coordinate */
790  double *l1t_samp) /* O: Projected sample coordinate */ {
791  VECTOR pos; /* Spacecraft ECEF position vector */
792  VECTOR vel; /* Spacecraft ECEF velocity vector */
793  VECTOR slos; /* Spacecraft line-of-sight vector */
794  int band_index; /* Index of current band in metadata structure */
795  double gp_lat; /* Ground point latitude */
796  double gp_lon; /* Ground point longitude */
797  double proj_x; /* Ground point projection X coordinate */
798  double proj_y; /* Ground point projection Y coordinate */
799  double d2r; /* Conversion from degrees to radians */
800 
801  /* Compute conversion constant */
802  d2r = atan(1.0) / 45.0;
803 
804  /* Find the nominal ephemeris state vector */
805  pathrow_to_posvel(wrsparms, smeta.wrs_path, dellon, frow, &pos, &vel);
806 
807  /* Find the current band in the metadata */
808  band_index = smeta_band_number_to_index(&smeta, band);
809 
810  /* Calculate the instrument LOS vector */
811  if (calc_los_ls(smeta.band_smeta_ls[band_index], frow, npix, ndet, &slos) != SUCCESS) {
812  printf("Error generating LOS vector.\n");
813  return (ERROR);
814  }
815 
816  /* Project to ground latitude/longitude */
817  if (calc_yaw_steering_gp_ls(wrsparms, pos, vel, slos, smeta.roll_angle, &gp_lat, &gp_lon) < 0) {
818  printf("Error projecting nominal scene center.\n");
819  return (ERROR);
820  }
821  gp_lat *= d2r;
822  gp_lon *= d2r;
823 
824  /* Apply scene map projection */
825  if (smeta_geodetic_to_proj(smeta.projection, gp_lat, gp_lon, &proj_x, &proj_y) != SUCCESS) {
826  printf("Error converting scene center lat/lon to projection X/Y.\n");
827  return (ERROR);
828  }
829 
830  /* Reduce to L1T line/sample */
831  if (smeta_proj_to_l1t(&smeta, band, proj_x, proj_y, l1t_line, l1t_samp) != SUCCESS) {
832  printf("Errro converting scene center X/Y to L1T line/sample.\n");
833  return (ERROR);
834  }
835 
836  return (SUCCESS);
837 }
838 
839 /*
840  * Name: get_ecfsun.c
841  *
842  * Purpose: Convert the metadata sun angles into an ECEF solar vector.
843  * As a first approximation the sun is assumed not to change
844  * position during the scene as the actual change should be less
845  * than 10 arc-minutes.
846  *
847  */
848 
849 int get_ecfsun(
850  WRS2 wrsparms, /* Structure of WRS-2 parameters */
851  SMETA smeta, /* Input scene metadata structure */
852  double dellon, /* WRS longitude offset */
853  VECTOR *ecfsun) /* Sun vector in ECEF coordinates */ {
854  VECTOR pos; /* Spacecraft ECEF position vector */
855  VECTOR vel; /* Spacecraft ECEF velocity vector */
856  VECTOR los; /* Spacecraft line-of-sight vector */
857  double drow; /* Fractional WRS row */
858  double gp_lat; /* Ground point latitude */
859  double gp_lon; /* Ground point longitude */
860  double d2r; /* Conversion from degrees to radians */
861  double e2; /* WGS84 eccentricity squared */
862  double ecf2lsr[3][3]; /* ECEF to Local Space Rectangular transformation matrix */
863  VECTOR lsrsun; /* Sun vector in LSR coordinates */
864 
865  /* Compute conversion constant */
866  d2r = atan(1.0) / 45.0;
867  e2 = 1.0 - wrsparms.semiminor / wrsparms.semimajor * wrsparms.semiminor / wrsparms.semimajor;
868 
869  /* Calculate the scene center latitude/longitude */
870  drow = (double) smeta.wrs_row;
871  pathrow_to_posvel(wrsparms, smeta.wrs_path, dellon, drow, &pos, &vel);
872  los.x = 0.0;
873  los.y = 0.0;
874  los.z = 1.0; /* boresight */
875  if (calc_yaw_steering_gp(wrsparms, pos, vel, los, smeta.roll_angle, &gp_lat, &gp_lon) < 0) {
876  printf("Error projecting nominal scene center.\n");
877  return (ERROR);
878  }
879  gp_lat *= d2r;
880  gp_lon *= d2r;
881 
882  /* Convert latitude to geocentric */
883  /* This looks wrong but this is how the metadata angles
884  are calculated, so we need to back out the geocentric
885  - geodetic effect */
886  gp_lat = atan(tan(gp_lat) * (1.0 - e2));
887 
888  /* Compute ECEF solar vector */
889  lsrsun.x = cos(d2r * smeta.sun_elevation) * sin(d2r * smeta.sun_azimuth);
890  lsrsun.y = cos(d2r * smeta.sun_elevation) * cos(d2r * smeta.sun_azimuth);
891  lsrsun.z = sin(d2r * smeta.sun_elevation);
892  (void) smeta_geodetic_to_ecf2lsr(gp_lat, gp_lon, ecf2lsr);
893  ecfsun->x = ecf2lsr[0][0] * lsrsun.x + ecf2lsr[1][0] * lsrsun.y + ecf2lsr[2][0] * lsrsun.z;
894  ecfsun->y = ecf2lsr[0][1] * lsrsun.x + ecf2lsr[1][1] * lsrsun.y + ecf2lsr[2][1] * lsrsun.z;
895  ecfsun->z = ecf2lsr[0][2] * lsrsun.x + ecf2lsr[1][2] * lsrsun.y + ecf2lsr[2][2] * lsrsun.z;
896 
897  return (SUCCESS);
898 }
899 
900 int get_ecfsun_ls(
901  WRS2 wrsparms, /* Structure of WRS-2 parameters */
902  SMETA smeta, /* Input scene metadata structure */
903  double dellon, /* WRS longitude offset */
904  VECTOR *ecfsun) /* Sun vector in ECEF coordinates */ {
905  VECTOR pos; /* Spacecraft ECEF position vector */
906  VECTOR vel; /* Spacecraft ECEF velocity vector */
907  VECTOR los; /* Spacecraft line-of-sight vector */
908  double drow; /* Fractional WRS row */
909  double gp_lat; /* Ground point latitude */
910  double gp_lon; /* Ground point longitude */
911  double d2r; /* Conversion from degrees to radians */
912  double e2; /* WGS84 eccentricity squared */
913  double ecf2lsr[3][3]; /* ECEF to Local Space Rectangular transformation matrix */
914  VECTOR lsrsun; /* Sun vector in LSR coordinates */
915 
916  /* Compute conversion constant */
917  d2r = atan(1.0) / 45.0;
918  e2 = 1.0 - wrsparms.semiminor / wrsparms.semimajor * wrsparms.semiminor / wrsparms.semimajor;
919 
920  /* Calculate the scene center latitude/longitude */
921  drow = (double) smeta.wrs_row;
922  pathrow_to_posvel(wrsparms, smeta.wrs_path, dellon, drow, &pos, &vel);
923  los.x = 0.0;
924  los.y = 0.0;
925  los.z = 1.0; /* boresight */
926  if (calc_yaw_steering_gp_ls(wrsparms, pos, vel, los, smeta.roll_angle, &gp_lat, &gp_lon) < 0) {
927  printf("Error projecting nominal scene center.\n");
928  return (ERROR);
929  }
930  gp_lat *= d2r;
931  gp_lon *= d2r;
932 
933  /* Convert latitude to geocentric */
934  /* This looks wrong but this is how the metadata angles
935  are calculated, so we need to back out the geocentric
936  - geodetic effect */
937  gp_lat = atan(tan(gp_lat) * (1.0 - e2));
938 
939  /* Compute ECEF solar vector */
940  lsrsun.x = cos(d2r * smeta.sun_elevation) * sin(d2r * smeta.sun_azimuth);
941  lsrsun.y = cos(d2r * smeta.sun_elevation) * cos(d2r * smeta.sun_azimuth);
942  lsrsun.z = sin(d2r * smeta.sun_elevation);
943  (void) smeta_geodetic_to_ecf2lsr(gp_lat, gp_lon, ecf2lsr);
944  ecfsun->x = ecf2lsr[0][0] * lsrsun.x + ecf2lsr[1][0] * lsrsun.y + ecf2lsr[2][0] * lsrsun.z;
945  ecfsun->y = ecf2lsr[0][1] * lsrsun.x + ecf2lsr[1][1] * lsrsun.y + ecf2lsr[2][1] * lsrsun.z;
946  ecfsun->z = ecf2lsr[0][2] * lsrsun.x + ecf2lsr[1][2] * lsrsun.y + ecf2lsr[2][2] * lsrsun.z;
947 
948  return (SUCCESS);
949 }
950 
951 /*
952  * Name: pathrow_to_latlon.c
953  *
954  * Purpose: Computes the latitude and longitude corresponding to an input WRS-2
955  * path/row including fractional row.
956  *
957  */
958 
959 void pathrow_to_latlon(
960  WRS2 parms, /* I: WRS2 system parameters */
961  int path, /* I: WRS2 path */
962  double wrsrow, /* I: WRS2 fractional row */
963  double *lat, /* O: Nominal latitude */
964  double *lon) /* O: Nominal longitude */ {
965  double cta; /* Central travel angle */
966  double dnodelon; /* Longitude of the descending node */
967  double dlon; /* Longitude offset from descending node */
968  double gc_lat; /* Geocentric latitude */
969  double d2r; /* Degrees to radians conversion */
970 
971  /* Initialize constants */
972  d2r = atan(1.0) / 45.0;
973 
974  /* Calculate CTA from row */
975  cta = 360.0 * (wrsrow - (double) parms.row0lat) / (double) parms.numrow;
976 
977  /* Calculate latitude from CTA */
978  gc_lat = asin(-sin(d2r * cta) * sin(d2r * parms.orbitincl));
979  *lat = atan(tan(gc_lat) * parms.semimajor / parms.semiminor * parms.semimajor / parms.semiminor) / d2r;
980 
981  /* Calculate longitude of descending node from path and CTA */
982  dnodelon = 360.0 + 180.0 + parms.path1long - (double) (path - 1)*360.0 / (double) parms.numpath;
983  dnodelon = fmod(dnodelon, 360.0) - 180.0 - cta * (double) parms.cycle / (double) parms.numpath;
984 
985  /* Calculate the longitude offset */
986  dlon = atan2(tan(d2r * gc_lat) / tan(parms.orbitincl), cos(d2r * cta) / cos(d2r * gc_lat)) / d2r;
987  *lon = dnodelon - dlon;
988  if (*lon > 180.0) *lon -= 360.0;
989  if (*lon < -180.0) *lon += 360.0;
990 
991  return;
992 }
993 
994 /*
995  * Name: pathrow_to_posvel.c
996  *
997  * Purpose: Computes the nominal satellite position and velocity (ECEF) vectors
998  * that correspond to the specified WRS path/row including fractional row
999  * and path longitude adjustments.
1000  *
1001  */
1002 
1003 void pathrow_to_posvel(
1004  WRS2 parms, /* I: WRS2 system parameters */
1005  int path, /* I: WRS2 path */
1006  double dellon, /* I: Path longitude adjustment (in degrees) */
1007  double wrsrow, /* I: WRS2 fractional row */
1008  VECTOR *pos, /* O: Nominal position vector */
1009  VECTOR *vel) /* O: Nominal velocity vector */ {
1010  VECTOR on; /* Orbit normal vector */
1011  VECTOR dnode; /* Descending node vector */
1012  VECTOR y; /* In-plane normal vector */
1013  double cta; /* Central travel angle */
1014  double dnodelon; /* Longitude of the descending node */
1015  double orate; /* Orbital angular rate */
1016  double d2r; /* Degrees to radians conversion */
1017 
1018  /* Initialize constants */
1019  d2r = atan(1.0) / 45.0;
1020 
1021  /* Calculate CTA from row */
1022  cta = 360.0 * (wrsrow - (double) parms.row0lat) / (double) parms.numrow;
1023 
1024  /* Calculate longitude of descending node from path and CTA */
1025  dnodelon = 360.0 + 180.0 + parms.path1long - (double) (path - 1)*360.0 / (double) parms.numpath + dellon;
1026  dnodelon = fmod(dnodelon, 360.0) - 180.0 - cta * (double) parms.cycle / (double) parms.numpath;
1027 
1028  /* Compute the orbit normal vector */
1029  on.x = -sin(d2r * parms.orbitincl) * sin(d2r * dnodelon);
1030  on.y = sin(d2r * parms.orbitincl) * cos(d2r * dnodelon);
1031  on.z = cos(d2r * parms.orbitincl);
1032 
1033  /* Compute the descending node vector */
1034  dnode.x = cos(d2r * dnodelon);
1035  dnode.y = sin(d2r * dnodelon);
1036  dnode.z = 0.0;
1037  crossprod(on, dnode, &y);
1038 
1039  /* Construct the position vector */
1040  pos->x = (dnode.x * cos(d2r * cta) + y.x * sin(d2r * cta)) * parms.orbitrad;
1041  pos->y = (dnode.y * cos(d2r * cta) + y.y * sin(d2r * cta)) * parms.orbitrad;
1042  pos->z = (dnode.z * cos(d2r * cta) + y.z * sin(d2r * cta)) * parms.orbitrad;
1043 
1044  /* Construct the velocity vector */
1045  orate = 360.0 * (double) parms.numpath / (double) parms.cycle / 86400.0 * d2r;
1046  on.x *= orate;
1047  on.y *= orate;
1048  on.z *= orate;
1049  crossprod(on, *pos, vel);
1050 
1051  return;
1052 }
1053 
1054 /*
1055  * Name: vector_math.c
1056  *
1057  * Purpose: Code units to calculate vector functions.
1058  *
1059  */
1060 
1061 double dotprod(/* Returns dot product of two VECTOR structures */
1062  VECTOR a, /* I: First input vector */
1063  VECTOR b) /* I: Second input vector */ {
1064  return ( a.x * b.x + a.y * b.y + a.z * b.z);
1065 }
1066 
1067 void crossprod(/* Computes vector cross product */
1068  VECTOR a, /* I: First input vector */
1069  VECTOR b, /* I: Second input vector */
1070  VECTOR *c) /* O: Output vector */ {
1071  c->x = a.y * b.z - a.z * b.y;
1072  c->y = a.z * b.x - a.x * b.z;
1073  c->z = a.x * b.y - a.y * b.x;
1074  return;
1075 }
1076 
1077 double vecnorm(/* Returns vector magnitude */
1078  VECTOR a) /* I: Input vector */ {
1079  return ( sqrt(dotprod(a, a)));
1080 }
1081 
1082 void unitvec(/* Normalizes input vector */
1083  VECTOR a, /* I: Input vector */
1084  VECTOR *b) /* O: Normalized output vector */ {
1085  double mag; /* Vector magnitude */
1086 
1087  mag = vecnorm(a);
1088  if (mag > 0.0) {
1089  b->x = a.x / mag;
1090  b->y = a.y / mag;
1091  b->z = a.z / mag;
1092  } else {
1093  b->x = a.x;
1094  b->y = a.y;
1095  b->z = a.z;
1096  }
1097  return;
1098 }
1099 
1100 void rotatevec(
1101  double mat[3][3], /* 3-by-3 rotation matrix */
1102  VECTOR a, /* Original vector */
1103  VECTOR *b) /* Rotated vector */ {
1104  b->x = mat[0][0] * a.x + mat[0][1] * a.y + mat[0][2] * a.z;
1105  b->y = mat[1][0] * a.x + mat[1][1] * a.y + mat[1][2] * a.z;
1106  b->z = mat[2][0] * a.x + mat[2][1] * a.y + mat[2][2] * a.z;
1107 
1108  return;
1109 }
dotprod
double dotprod(VECTOR a, VECTOR b)
Definition: mtl_geometry.c:1061
frame_scene_ls
int frame_scene_ls(WRS2 wrsparms, SMETA smeta, double *dellon, double *delrow)
Definition: mtl_geometry.c:591
MAX
#define MAX(A, B)
Definition: swl0_utils.h:26
MIN
#define MIN(x, y)
Definition: rice.h:169
SMETA::projection
SMETA_SCENE_PROJ projection
Definition: smeta.h:130
SMETA_BAND_LS
Definition: smeta.h:95
SUCCESS
#define SUCCESS
Definition: ObpgReadGrid.h:15
smeta_band_number_to_index
int smeta_band_number_to_index(SMETA *smeta, int band)
Definition: smeta_geometry.c:205
crossprod
void crossprod(VECTOR a, VECTOR b, VECTOR *c)
Definition: mtl_geometry.c:1067
VECTOR::y
double y
Definition: mtl_geometry.h:38
pathrow_to_posvel
void pathrow_to_posvel(WRS2 parms, int path, double dellon, double wrsrow, VECTOR *pos, VECTOR *vel)
Definition: mtl_geometry.c:1003
SMETA_BAND_LS::rascan
double rascan
Definition: smeta.h:107
calc_los
int calc_los(SMETA_BAND band_smeta, int sca, double ndet, VECTOR *los)
Definition: mtl_geometry.c:19
WRS2::inertialvel
double inertialvel
Definition: mtl_geometry.h:25
WRS2::numrow
int numrow
Definition: mtl_geometry.h:30
SMETA_BAND_LS::l1t_samps
int l1t_samps
Definition: smeta.h:99
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
SMETA::sun_elevation
double sun_elevation
Definition: smeta.h:129
WRS2::semimajor
double semimajor
Definition: mtl_geometry.h:23
pos
float32 * pos
Definition: l1_czcs_hdf.c:35
unitvec
void unitvec(VECTOR a, VECTOR *b)
Definition: mtl_geometry.c:1082
lat
float * lat
Definition: l1_hmodis_hdf.c:756
SMETA_BAND_LS::fascan
double fascan
Definition: smeta.h:106
vecnorm
double vecnorm(VECTOR a)
Definition: mtl_geometry.c:1077
smeta_proj_to_l1t
int smeta_proj_to_l1t(SMETA *smeta, int band, double proj_x, double proj_y, double *l1t_line, double *l1t_samp)
Definition: smeta_geometry.c:73
names::if
character(len=1000) if
Definition: names.f90:13
SMETA::sun_azimuth
double sun_azimuth
Definition: smeta.h:128
SMETA_BAND_LS::aoffset
double aoffset
Definition: smeta.h:111
SMETA_BAND::l1t_samps
int l1t_samps
Definition: smeta.h:84
smeta_release_projection
void smeta_release_projection()
Definition: smeta_geometry.c:196
WRS2::orbitincl
double orbitincl
Definition: mtl_geometry.h:26
pathrow_to_latlon
void pathrow_to_latlon(WRS2 parms, int path, double wrsrow, double *lat, double *lon)
Definition: mtl_geometry.c:959
SMETA_BAND::legendre
double legendre[IAS_MAX_NSCAS][2][4]
Definition: smeta.h:91
SMETA_BAND::align
double align[3][3]
Definition: smeta.h:90
SMETA_BAND::l1t_lines
int l1t_lines
Definition: smeta.h:83
mtl_geometry.h
SMETA::wrs_path
int wrs_path
Definition: smeta.h:124
color_dtdb.path
string path
Definition: color_dtdb.py:221
hico.value_locate.x
list x
Definition: value_locate.py:64
color_dtdb.y
y
Definition: color_dtdb.py:430
calc_yaw_steering_gp
int calc_yaw_steering_gp(WRS2 parms, VECTOR pos, VECTOR vel, VECTOR i_los, double roll, double *gp_lat, double *gp_lon)
Definition: mtl_geometry.c:100
SMETA::band_smeta
SMETA_BAND band_smeta[IAS_MAX_NBANDS]
Definition: smeta.h:132
WRS2::semiminor
double semiminor
Definition: mtl_geometry.h:24
calc_yaw_steering_gp_ls
int calc_yaw_steering_gp_ls(WRS2 parms, VECTOR pos, VECTOR vel, VECTOR i_los, double roll, double *gp_lat, double *gp_lon)
Definition: mtl_geometry.c:265
WRS2::numpath
int numpath
Definition: mtl_geometry.h:29
color_dtdb.scale
int scale
Definition: color_dtdb.py:459
SMETA_BAND_LS::pixsize
double pixsize
Definition: smeta.h:104
SMETA::roll_angle
double roll_angle
Definition: smeta.h:126
SMETA_BAND
Definition: smeta.h:80
WRS2::cycle
int cycle
Definition: mtl_geometry.h:31
frame_scene
int frame_scene(WRS2 wrsparms, SMETA smeta, double *dellon, double *delrow)
Definition: mtl_geometry.c:419
smeta_geodetic_to_ecf2lsr
int smeta_geodetic_to_ecf2lsr(double lat, double lon, double ecf2lsr[3][3])
Definition: smeta_geometry.c:48
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
rotatevec
void rotatevec(double mat[3][3], VECTOR a, VECTOR *b)
Definition: mtl_geometry.c:1100
WRS2::row0lat
int row0lat
Definition: mtl_geometry.h:32
WRS2::path1long
double path1long
Definition: mtl_geometry.h:28
VECTOR::x
double x
Definition: mtl_geometry.h:37
science_parameters::d2r
real, parameter d2r
Definition: science_parameters.f90:60
project_point_ls
int project_point_ls(WRS2 wrsparms, SMETA smeta, double dellon, int band, double frow, double npix, double ndet, double *l1t_line, double *l1t_samp)
Definition: mtl_geometry.c:781
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
get_ecfsun
int get_ecfsun(WRS2 wrsparms, SMETA smeta, double dellon, VECTOR *ecfsun)
Definition: mtl_geometry.c:849
SMETA::band_smeta_ls
SMETA_BAND_LS band_smeta_ls[IAS_MAX_NBANDS_LS]
Definition: smeta.h:133
get_ecfsun_ls
int get_ecfsun_ls(WRS2 wrsparms, SMETA smeta, double dellon, VECTOR *ecfsun)
Definition: mtl_geometry.c:900
lon
float * lon
Definition: l1_hmodis_hdf.c:755
smeta_geodetic_to_proj
int smeta_geodetic_to_proj(SMETA_SCENE_PROJ proj, double lat, double lon, double *proj_x, double *proj_y)
Definition: smeta_geometry.c:158
WRS2::orbitrad
double orbitrad
Definition: mtl_geometry.h:27
c
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
SMETA_BAND::nsca
int nsca
Definition: smeta.h:85
radius
void radius(double A)
Definition: proj_report.c:132
VECTOR::z
double z
Definition: mtl_geometry.h:39
SMETA_BAND_LS::align
double align[3][3]
Definition: smeta.h:105
SMETA
Definition: smeta.h:118
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).
SMETA::wrs_row
int wrs_row
Definition: smeta.h:125
a
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
calc_los_ls
int calc_los_ls(SMETA_BAND_LS band_smeta, double frow, double npix, double ndet, VECTOR *los)
Definition: mtl_geometry.c:52
npix
int npix
Definition: get_cmp.c:27
SMETA_BAND_LS::l1t_lines
int l1t_lines
Definition: smeta.h:98
project_point
int project_point(WRS2 wrsparms, SMETA smeta, double dellon, int band, int sca, double frow, double ndet, double *l1t_line, double *l1t_samp)
Definition: mtl_geometry.c:722
ERROR
#define ERROR
Definition: ancil.h:24
SMETA_BAND::pixsize
double pixsize
Definition: smeta.h:89


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

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 Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster