OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
GEO_get_GRing_points.c
Go to the documentation of this file.
1 /* GEO_cumulate_GRing() and GEO_get_GRing_points() */
2 
3 #include <float.h>
4 #include "GEO_input.h"
5 #include "GEO_earth.h"
6 #include "imsl_wrap.h"
7 #include "PGS_SMF.h"
8 #include "PGS_MET.h"
9 #include "PGS_CSC.h"
10 #include "GEO_geo.h"
11 #include "PGS_MODIS_35251.h"
12 #include "smfio.h"
13 
14 static double tmin=DBL_MAX;
15 static double smin=DBL_MAX; /* Upper and lower bounds on*/
16 static double tmax=-DBL_MAX;
17 static double smax=-DBL_MAX; /* projected points. */
18 static PGSt_double project[3][3]; /* Projection axes. */
19 
20 PGSt_SMF_status GEO_cumulate_GRing(
21  GEO_param_struct const * const params,
22  int32 const num_frames,
23  frame_state_struct const sc_ev_frame_state[],
24  unsigned char frame_flags[DETECTORS_1KM][MAX_FRAMES],
25  double ecr_frame_position[DETECTORS_1KM][MAX_FRAMES][3]
26 )
27 
28 /**************************************************************************
29 !C
30 
31 !Description: Routine in the output group of the Level-1A geolocation
32  software to cumulate the information needed to determine
33  the GRing for a granule.
34 
35 !Input Parameters:
36  params Geolocation parameters
37  num_frames number of frames in this scan of data
38  sc_ev_frame_state An array containing, among other things, the
39  spacecraft position at the time of each frame.
40  frame_flags An array indicating any problems that may have
41  occurred while geolocating a pixel.
42  ecr_frame_position The ECR coordinates for each pixel.
43 
44 !Output Parameters:
45  None
46 
47 Return Parameters:
48  MODIS_E_BAD_INPUT_ARG If there's any problem with the inputs
49  MODIS_E_GEO If the matrix inversion fails
50  PGS_S_SUCCESS Otherwise
51 
52 Externally Defined:
53  DBL_EPSILON <float.h>
54  DBL_MAX <float.h>
55  DETECTORS_1KM "GEO_geo.h"
56  IMSL_INVERSE_USER "imsl.h"
57  IMSL_INVERSE_ONLY "imsl.h"
58  MAX_FRAMES "GEO_geo.h"
59  MODIS_E_BAD_INPUT_ARG "PGS_MODIS_35251.h"
60  MODIS_E_GEO "PGS_MODIS_35251.h"
61  NO_ELLIPSE_INTERSECT "GEO_geo.h"
62  PGS_S_SUCCESS "PGS_SMF.h"
63 
64 Called by:
65  GEO_locate_one_scan
66 
67 Routines Called:
68  imsl_d_lin_sol_gen
69  imsl_error_code
70  modsmf
71  PGS_CSC_crossProduct
72  PGS_CSC_dotProduct
73  PGS_CSC_Norm
74  PGS_CSC_quatRotate
75 
76 !Revision History:
77  * $Log: GEO_get_GRing_points.c,v $
78  * Revision 6.2 2010/06/18 21:10:37 kuyper
79  * Removed leading 'const' qualifiers from parameters that are pointers to
80  * arrays.
81  * Corrected error message for invalid input arguments.
82  *
83  * Revision 6.1 2009/05/31 20:12:27 kuyper
84  * Changes to GEO_cumulate_GRing():
85  * Use DETECTORS_1KM and MAX_FRAMES where needed.
86  * Changed names with "sample" to use "frame" instead.
87  * Removed validation of params->num_detectors; no longer needed.
88  * Expanded message buffer.
89  *
90  * Revision 4.3 2003/10/27 01:07:12 vlin
91  * expanded message buffer size.
92  *
93  * Revision 4.2 2003/08/01 18:22:42 kuyper
94  * Corrected error message format.
95  *
96  * Revision 4.1 2003/07/28 15:19:34 kuyper
97  * Changed to use imsl_wrap.h.
98  * Relaxed radius test for division by zero.
99  * Corrected test for an error return from imsl_d_lin_sol_gen().
100  * Corrected to cast pointer arguments for "%p" to (void*).
101  *
102  * Revision 3.8 2001/12/06 16:11:35 kuyper
103  * Relaxed division by zero guards.
104  *
105 
106 Requirements:
107  PR03-F-4.1-1
108  PR03-F-4.4-1
109  PR03-I-1
110  PR03-I-2
111  PR03-S-1
112 
113 !Team-unique Header:
114 
115  This software is developed by the MODIS Science Data Support
116  Team for the National Aeronautics and Space Administration,
117  Goddard Space Flight Center, under contract NAS5-32373.
118 
119 References and Credits
120 
121 !END**************************************************************************/
122 
123 #define T_HALF 150.0 /* Half the length of a normal granule (seconds) */
124 {
125  static char filefunc[] = __FILE__ ", GEO_cumulate_GRing";
126  static double inverse[3][3]; /* Inverse of project[][]. */
127 
128  PGSt_SMF_status retval=PGS_S_SUCCESS; /* Value to be returned.*/
129  int det, frame; /* detector number, frame number */
130  int endframe; /* Frame number of last good pixel. */
131  PGSt_double positionECR[3]; /* Spacecraft position. */
132  PGSt_double velocityECR[3]; /* Spacecraft velocity. */
133  PGSt_double quat[4]; /* Rotation quaternion. */
134  PGSt_double temp[3]; /* Temporary storage for a vector. */
135  double radius; /* magnitude of the positionECR vector */
136  double pcrossv; /* magnitude of the cross-product of the*/
137  /* positionECR and velocityECR vectors */
138  double theta; /* Rotation angel for T_HALF. */
139  double length=0.0; /* Length of a vector. */
140  double matrix[3][3]; /* Projection axes. */
141  double rts[3]; /* projected components of ecr_frame_position*/
142  int row, col;
143  char msgbuf[PGS_SMF_MAX_MSGBUF_SIZE]="";
144 
145  if (params == NULL || sc_ev_frame_state == NULL || frame_flags == NULL
146  || ecr_frame_position == NULL || num_frames > MAX_FRAMES)
147  {
148  sprintf(msgbuf, "params:%p, num_frames:%ld, sc_ev_frame_state:%p, "
149  "frame_flags:%p, ecr_frame_position:%p",
150  (void*)params, (long)num_frames, (void*)sc_ev_frame_state,
151  (void*)frame_flags, (void*)ecr_frame_position);
152  modsmf(MODIS_E_BAD_INPUT_ARG, msgbuf, filefunc);
153 
154  return MODIS_E_BAD_INPUT_ARG;
155  }
156 
157  for (frame = 0; frame < num_frames; frame++)
158  {
159  for (det = 0; det < DETECTORS_1KM; det++)
160  {
161  if (frame_flags[det][frame] < (unsigned char) NO_ELLIPSE_INTERSECT)
162  /* the pixel was successfully geolocated */
163  {
164  if (tmin >= DBL_MAX)
165  { /* this is first pixel processed for this granule. */
166  /* Determine the projection matrix */
167 
168  for (col=0; col<3; col++) {
169  positionECR[col] =
170  (PGSt_double)sc_ev_frame_state[frame].positionECR[col];
171  velocityECR[col] =
172  (PGSt_double)sc_ev_frame_state[frame].velocityECR[col];
173  }
174  PGS_CSC_crossProduct(positionECR, velocityECR, quat+1);
175  radius = PGS_CSC_Norm(positionECR);
176  pcrossv = PGS_CSC_Norm(quat+1);
177 
178  if (radius < 0.5 * params->orbit_valid_params.position_mag_limit[0] ||
179  pcrossv < params->orbit_valid_params.ang_mom_limit[0])
180  {
181  sprintf(msgbuf, "sc_ev_frame_state[%d].positionECR too small",
182  frame);
183  modsmf(MODIS_E_BAD_INPUT_ARG, msgbuf, filefunc);
184 
185  retval = MODIS_E_BAD_INPUT_ARG;
186  break; /* break out of detector loop. This problem will affect
187  * every pixel in same frame. */
188  }
189 
190  memcpy(project[1], velocityECR, sizeof(project[1]));
191 
192  /* Search for last good frame in same scanline */
193  for(endframe=num_frames-1;
194  endframe>frame && frame_flags[det][endframe]>=NO_ELLIPSE_INTERSECT;
195  endframe--);
196  if(endframe > frame)
197  { /* Normal case: found. */
198  for (col=0; col<3; col++)
199  project[2][col] = ecr_frame_position[det][endframe][col]
200  - ecr_frame_position[det][frame][col];
201  /* Calculate normal to projection plane from scan and track */
202  PGS_CSC_crossProduct(project[1], project[2], project[0]);
203  length = PGS_CSC_Norm(project[0]);
204  }
205  if(length < DBL_EPSILON)
206  { /* length was either not calculated, or is too short. */
207  /* Either there was exactly one good pixel in the scanline, or the
208  * yaw was exactly 90 degrees. Either case should be pretty rare.
209  * Substitute the rotation axis for the scan direction. */
210  memcpy(project[0], positionECR, sizeof(project[0]));
211  memcpy(project[2], quat+1, sizeof(project[2]));
212  }
213  else if(PGS_CSC_dotProduct(project[0], positionECR, 3) < 0.0)
214  { /* spacecraft flying backwards; unlikely but possible. */
215  for (col=0; col<3; col++)
216  { /* Reverse two axes, to maintain intended orientation. */
217  project[0][col] *= -1.0;
218  project[2][col] *= -1.0;
219  }
220  }
221 
222  /* Rotate the vectors forward corresponding to T_HALF seconds of
223  * orbital rotation. */
224  theta = T_HALF*pcrossv/(radius*radius);
225  quat[0] = cos(0.5*theta);
226  for (col=1; col<4; col++)
227  quat[col] *= sin(0.5*theta)/pcrossv;
228  for(row=0; row<3; row++)
229  {
230  if(PGS_CSC_quatRotate(quat, project[row], temp)!=PGS_S_SUCCESS)
231  {
232  modsmf(MODIS_E_GEO, "PGS_CSC_quatRotate()", filefunc);
233  break;
234  }
235  length = PGS_CSC_Norm(temp);
236  if(length < DBL_EPSILON)
237  {
238  modsmf(MODIS_E_GEO, "PGS_CSC_Norm()", filefunc);
239  break;
240  }
241  for(col=0; col<3; col++)
242  project[row][col] = temp[col]/length;
243  }
244  if(row<3)
245  continue;
246 
247  for (row=0; row<3; row++)
248  for (col=0; col<3; col++)
249  matrix[row][col] = (double)project[row][col];
250 
251  /* Invert matrix. */
253  IMSL_INVERSE_USER, inverse[0],
254  IMSL_INVERSE_ONLY,
255  0);
256  if(imsl_error_code())
257  { /* Matrix singular or not positive definite */
258  modsmf(MODIS_E_GEO, "imsl_d_lin_sol_gen()", filefunc);
259  continue;
260  }
261  } /* if tmin >= DBL_MAX ended */
262 
263  for (row=0; row<3; row++)
264  {
265  rts[row] = 0.0;
266  for (col=0; col<3; col++)
267  rts[row] += inverse[col][row] *
268  ecr_frame_position[det][frame][col];
269  }
270 
271  if (rts[0] > (0.01 * params->orbit_valid_params.position_mag_limit[0]))
272  {
273  rts[1] /= rts[0];
274  rts[2] /= rts[0];
275  tmin = (tmin < rts[1]) ? tmin : rts[1];
276  tmax = (tmax > rts[1]) ? tmax : rts[1];
277  smin = (smin < rts[2]) ? smin : rts[2];
278  smax = (smax > rts[2]) ? smax : rts[2];
279  }
280  else
281  {
282  sprintf(msgbuf,"ecr_frame_position[%d][%d]: rts[0]=%g",
283  det, frame, rts[0]);
284  modsmf(MODIS_E_BAD_INPUT_ARG, msgbuf, filefunc);
285  retval = MODIS_E_BAD_INPUT_ARG;
286  }
287  } /* end of if the pixel was successfully geolocated */
288  } /* for det loop ended */
289  } /* for frame loop ended */
290  return retval;
291 }
292 
293 /******************************************************************************/
294 
295 PGSt_SMF_status GEO_get_GRing_points(
296  GEO_GRing_struct * const GRing_points)
297 
298 /************************************************************************
299 !C
300 
301 !Description:
302  Routine in Output group of the Level-1A geolocation software to
303  calculate the G-Ring polygon vertices for the ECS Core Metadata.
304  At this point, GEO_cumulate_GRing() has already projected each ground
305  location from the center of the earth onto a plane that's roughly
306  oriented with the granule. It has determined the bounding box in that
307  plane containing all of those projected points. This routine projects
308  the corners of that bounding box toward the center of the earth,
309  stopping at the surface of the earth. As result, the great circle arcs
310  that connect those points correspond precisely to the edges of the
311  bounding box, and will therefore enclose all of the ground locations.
312 
313 !Input Parameters:
314  None
315 
316 !Output Parameters:
317  GRing_points G-Ring polygon latitudes and longitudes
318 
319 Return Parameters:
320  MODIS_E_BAD_INPUT_ARG If GRing_points is NULL
321  MODIS_E_GEO If projection failed for any vertex
322  MODIS_W_NO_GEO If no geolocatable pixels were found
323  PGS_S_SUCCESS If projection was successful
324 
325 Global variables: None
326 
327 Externally Defined:
328  DBL_MAX <float.h>
329  MODIS_E_GEO "PGS_MODIS_35251.h"
330  MODIS_E_BAD_INPUT__ARG "PGS_MODIS_35251.h"
331  MODIS_W_NO_GEO "PGS_MODIS_35251.h"
332  EARTH_MODEL "GEO_earth.h"
333  PGS_S_SUCCESS "PGS_SMF.h"
334  PGSCSC_W_HIT_EARTH "PGS_CSC.h"
335  RAD2DEG "GEO_geo.h"
336 
337 Called by:
338  GEO_write_granule_metadata
339 
340 Routines Called:
341  modsmf Logs a status message
342  PGS_CSC_GrazingRay Projects from a point along a line of
343  sight to the surface of the Earth.
344  PGS_CSC_ECRtoGEO Converts from ECR to Geodetic coordinates.
345 
346 !Revision History:
347  Please see prologue in "GEO_cumulate_GRing.c"
348 
349 Requirements:
350  PR03-F-4.1-1
351  PR03-F-4.4-1
352  PR03-I-1
353  PR03-I-2
354  PR03-S-1
355 
356 !Team-unique Header:
357  This software is developed by the MODIS Science Data Support
358  Team for the National Aeronautics and Space Administration,
359  Goddard Space Flight Center, under contract NAS5-32373.
360 
361 !END
362 ***************************************************************************/
363 
364 {
365 #define CORNERS 4
366 #define BORDER 1e-5 /* roughly 1/2 the angular width in radians
367  of a pixel at nadir */
368  PGSt_double corner_point[CORNERS][3], posECR[3];
369  PGSt_double length, ray[3];
370  PGSt_double missAlt, slantRg;
371  PGSt_double posNEAR[3], posSURF[3];
372  int cor, j;
373  GEO_GRing_struct new_points;
374  PGSt_SMF_status rtnfun;
375 
376  if (GRing_points == NULL)
377  {
378  modsmf(MODIS_E_BAD_INPUT_ARG,".",
379  "GEO_get_GRing_points.c, GEO_get_GRing_points");
380  return MODIS_E_BAD_INPUT_ARG;
381  }
382 
383  new_points = *GRing_points;
384 
385  if (tmin >= DBL_MAX)
386  {
387  /* There were no geolocatable pixels in the granule. This is a normal
388  * occurrance in the face of bad data or an empty granule; not a defect in
389  * this program.*/
390  modsmf(MODIS_W_NO_GEO,"", "GEO_get_GRing_points.c, GEO_get_GRing_points");
391  return MODIS_W_NO_GEO;
392  }
393 
394 /* The following adjustments guarantee that even if there's only one
395  geolocatable pixel in the granule, the GRing will cover a finite area -
396  GRings with zero area are not permitted by the archiving system. */
397 
398  tmax += BORDER;
399  smax += BORDER;
400  tmin -= BORDER;
401  smin -= BORDER;
402 
403  for (j=0; j<3; j++) {
404  corner_point[0][j] = project[0][j] + tmin*project[1][j] + smin*project[2][j];
405  corner_point[1][j] = project[0][j] + tmin*project[1][j] + smax*project[2][j];
406  corner_point[2][j] = project[0][j] + tmax*project[1][j] + smax*project[2][j];
407  corner_point[3][j] = project[0][j] + tmax*project[1][j] + smin*project[2][j];
408  }
409 
410  for (cor = 0; cor < CORNERS; cor++)
411  {
412  length = 1.0;
413  for (j=0; j<3; j++) {
414  ray[j] = -corner_point[cor][j];
415  posECR[j] = 7.0e6 * corner_point[cor][j];
416  if (fabs(ray[j]) >= length)
417  length = fabs(ray[j]);
418  }
419 
420  if (length > 1.0)
421  for (j=0; j<3; j++)
422  ray[j] /=length;
423 
424  rtnfun = PGS_CSC_GrazingRay(EARTH_MODEL, posECR, ray,
425  &new_points.latitude[cor], &new_points.longitude[cor],
426  &missAlt, &slantRg, posNEAR, posSURF);
427  /* project along a ray parrallel to the corner's
428  position vector, up to the earth's surface,
429  starting at the corner's position */
430 
431  if (rtnfun != PGSCSC_W_HIT_EARTH)
432  {
433  modsmf(MODIS_E_GEO,"PGS_CSC_GrazingRay()",
434  "GEO_get_GRing_points.c, GEO_get_GRing_points");
435  return MODIS_E_GEO;
436  }
437 
438  rtnfun = PGS_CSC_ECRtoGEO(posSURF, EARTH_MODEL,
439  &new_points.longitude[cor], &new_points.latitude[cor], &missAlt);
440 
441  if (rtnfun != PGS_S_SUCCESS)
442  {
443  modsmf(MODIS_E_GEO,"PGS_CSC_ECRtoGEO()",
444  "GEO_get_GRing_points.c, GEO_get_GRing_points");
445  return MODIS_E_GEO;
446  }
447  new_points.latitude[cor] *= RAD2DEG;
448  new_points.longitude[cor] *= RAD2DEG;
449 
450  } /* for cor loop ended */
451 
452  tmin=DBL_MAX;
453  smin=DBL_MAX;
454  tmax=-DBL_MAX;
455  smax=-DBL_MAX;
456 
457  *GRing_points = new_points;
458 
459  return PGS_S_SUCCESS;
460 }
#define T_HALF
int j
Definition: decode_rs.h:73
#define NULL
Definition: decode_rs.h:63
#define MODIS_E_BAD_INPUT_ARG
#define MAX_FRAMES
Definition: GEO_geo.h:79
PGSt_SMF_status GEO_cumulate_GRing(GEO_param_struct const *const params, int32 const num_frames, frame_state_struct const sc_ev_frame_state[], unsigned char frame_flags[DETECTORS_1KM][MAX_FRAMES], double ecr_frame_position[DETECTORS_1KM][MAX_FRAMES][3])
double *IMSL_DECL imsl_d_lin_sol_gen(int const dim, double matrix[], double not_used[],...)
const unsigned char NO_ELLIPSE_INTERSECT
float ** matrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:60
#define CORNERS
#define MODIS_E_GEO
#define BORDER
#define RAD2DEG
Definition: GEO_geo.h:173
long IMSL_DECL imsl_error_code(void)
Definition: imsl_error.c:27
#define DETECTORS_1KM
Definition: GEO_geo.h:85
#define EARTH_MODEL
Definition: GEO_earth.h:118
#define MODIS_W_NO_GEO
PGSt_SMF_status GEO_get_GRing_points(GEO_GRing_struct *const GRing_points)
#define fabs(a)
Definition: misc.h:93
#define DBL_MAX
Definition: make_L3_v1.1.c:60
void radius(double A)
Definition: proj_report.c:132