OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
mtl_grid.c
Go to the documentation of this file.
1 /*
2  * Name: mtl_grid.c
3  *
4  * Purpose: Routines that create and use the metadata grid file for per-pixel
5  * view and sun angle generation.
6  */
7 
8 #include <stdio.h>
9 #include <math.h>
10 #include <stdlib.h>
11 #include "mtl_grid.h"
12 #define TINY 1.0e-12
13 
15  WRS2 parms, /* I: WRS2 system parameters */
16  SMETA smeta, /* I: Scene metadata structure */
17  double dellon, /* I: Path longitude adjustment (in degrees) */
18  VECTOR ecfsun, /* I: ECF sun vector */
19  int band_index, /* I: Index for current band */
20  MTL_GRID_BAND *bgrid ) /* O: Band grid structure */
21 {
22  int nrow; /* Row index */
23  int ncol; /* Column index */
24  int sca; /* SCA index */
25 
26  /* Allocate the grid */
27  bgrid->band = smeta.band_smeta[band_index].band;
28  bgrid->nsca = smeta.band_smeta[band_index].nsca;
29  /* Sudipta new addition for OLI to extract SCA and Det numbers */
30  bgrid->num_detectors = smeta.band_smeta[band_index].num_detectors;
31  /* End Sudipta new addition for OLI */
32  if ( (bgrid->sgrid = (MTL_GRID_SCA *)malloc( bgrid->nsca * sizeof( MTL_GRID_SCA ) )) == NULL )
33  {
34  printf("Unable to allocate grid for band %d\n", bgrid->band );
35  return(ERROR);
36  }
37 
38  /* Loop through the SCAs */
39  for ( sca=0; sca<bgrid->nsca; sca++ )
40  {
41  bgrid->sgrid[sca].min_line = 1.0e9;
42  bgrid->sgrid[sca].max_line = -1.0e9;
43  bgrid->sgrid[sca].min_samp = 1.0e9;
44  bgrid->sgrid[sca].max_samp = -1.0e9;
45  /* Loop through the grid */
46  for ( nrow=0; nrow<NUM_GRID_ROW; nrow++ )
47  for ( ncol=0; ncol<NUM_GRID_COL; ncol++ )
48  {
49  bgrid->sgrid[sca].pgrid[nrow][ncol].ndet = GRID_COL_BASE + (double)ncol * GRID_COL_INC;
50  bgrid->sgrid[sca].pgrid[nrow][ncol].frow = GRID_ROW_BASE + (double)nrow * GRID_ROW_INC;
51  if ( calc_grid_point( parms, smeta, dellon, band_index, sca+1, ecfsun, &bgrid->sgrid[sca].pgrid[nrow][ncol] ) != SUCCESS )
52  {
53  printf("Error generating grid point: SCA: %d, Row: %d, Col: %d\n", sca+1, nrow, ncol);
54  free( bgrid->sgrid );
55  return(ERROR);
56  }
57  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line < bgrid->sgrid[sca].min_line )
58  bgrid->sgrid[sca].min_line = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line;
59  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line > bgrid->sgrid[sca].max_line )
60  bgrid->sgrid[sca].max_line = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line;
61  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp < bgrid->sgrid[sca].min_samp )
62  bgrid->sgrid[sca].min_samp = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp;
63  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp > bgrid->sgrid[sca].max_samp )
64  bgrid->sgrid[sca].max_samp = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp;
65  }
66  /* Calculate line/sample to ndet/frow mapping */
67  if ( calc_grid_fit( &bgrid->sgrid[sca] ) != SUCCESS )
68  {
69  printf("Error calculating grid mapping coefficients for SCA %d\n", sca+1);
70  free( bgrid->sgrid );
71  return(ERROR);
72  }
73  }
74 
75  return(SUCCESS);
76 }
77 
78 
79 // need a forward declaration
81  WRS2 parms, /* I: WRS2 system parameters */
82  SMETA smeta, /* I: Scene metadata structure */
83  double dellon, /* I: Path longitude adjustment (in degrees) */
84  int band_index, /* I: Index for current band */
85  int sca, /* I: SCA number (1 to nsca) */
86  VECTOR ecfsun, /* I: ECF sun vector */
87  MTL_GRID_PT *gpt ); /* I/O: Grid point with only ndet and frow on input */
88 
89 
90 /* variant for older LS */
92  WRS2 parms, /* I: WRS2 system parameters */
93  SMETA smeta, /* I: Scene metadata structure */
94  double dellon, /* I: Path longitude adjustment (in degrees) */
95  VECTOR ecfsun, /* I: ECF sun vector */
96  int band_index, /* I: Index for current band */
97  MTL_GRID_BAND_LS *bgrid ) /* O: Band grid structure */
98 {
99  int nrow; /* Row index */
100  int ncol; /* Column index */
101  int sca = 0; /* SCA index */
102 
103  /* Allocate the grid */
104  bgrid->band = smeta.band_smeta_ls[band_index].band;
105  bgrid->nsca = smeta.band_smeta_ls[band_index].nsca;
106  if ( (bgrid->sgrid = (MTL_GRID_SCA_LS *)malloc( bgrid->nsca * sizeof( MTL_GRID_SCA_LS ) )) == NULL )
107  {
108  printf("Unable to allocate grid for band %d\n", bgrid->band );
109  return(ERROR);
110  }
111 
112  /* Loop through the SCAs */
113  for ( sca=0; sca<bgrid->nsca; sca++ )
114  {
115  bgrid->sgrid[sca].min_line = 1.0e9;
116  bgrid->sgrid[sca].max_line = -1.0e9;
117  bgrid->sgrid[sca].min_samp = 1.0e9;
118  bgrid->sgrid[sca].max_samp = -1.0e9;
119  /* Loop through the grid */
120  for ( nrow=0; nrow<NUM_GRID_ROW_LS; nrow++ )
121  for ( ncol=0; ncol<NUM_GRID_COL; ncol++ )
122  {
123  bgrid->sgrid[sca].pgrid[nrow][ncol].npix = GRID_COL_BASE + (double)ncol * GRID_COL_INC;
124  bgrid->sgrid[sca].pgrid[nrow][ncol].frow = GRID_ROW_BASE_LS + (double)nrow * GRID_ROW_INC;
125  bgrid->sgrid[sca].pgrid[nrow][ncol].ndet = GRID_COL_BASE + (double)ncol * GRID_COL_INC;
126  if ( calc_grid_point_ls( parms, smeta, dellon, band_index, sca+1, ecfsun, &bgrid->sgrid[sca].pgrid[nrow][ncol] ) != SUCCESS )
127  {
128  printf("Error generating grid point: SCA: %d, Row: %d, Col: %d\n", sca+1, nrow, ncol);
129  free( bgrid->sgrid );
130  return(ERROR);
131  }
132  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line < bgrid->sgrid[sca].min_line )
133  bgrid->sgrid[sca].min_line = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line;
134  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line > bgrid->sgrid[sca].max_line )
135  bgrid->sgrid[sca].max_line = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_line;
136  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp < bgrid->sgrid[sca].min_samp )
137  bgrid->sgrid[sca].min_samp = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp;
138  if ( bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp > bgrid->sgrid[sca].max_samp )
139  bgrid->sgrid[sca].max_samp = bgrid->sgrid[sca].pgrid[nrow][ncol].l1t_samp;
140  }
141  /* Calculate line/sample to ndet/frow mapping */
142  if ( calc_grid_fit_ls( &bgrid->sgrid[sca] ) != SUCCESS )
143  {
144  printf("Error calculating grid mapping coefficients for SCA %d\n", sca+1);
145  free( bgrid->sgrid );
146  return(ERROR);
147  }
148  }
149 
150  return(SUCCESS);
151 }
152 
154  MTL_GRID_SCA *sgrid ) /* I/O: SCA grid array to fit polynomial coefficients for */
155 {
156  int nrow; /* Row index */
157  int ncol; /* Column index */
158  double gline; /* Normalized grid line coordinate */
159  double gsamp; /* Normalized grid sample coordinate */
160  double norm[4][4]; /* Normal equations matrix */
161  double rhs1[4]; /* NDet constant vector */
162  double rhs2[4]; /* FRow constant vector */
163 
164  /* Initialize the normal equations */
165  norm[0][0] = 0.0; norm[0][1] = 0.0; norm[0][2] = 0.0; norm[0][3] = 0.0;
166  norm[1][0] = 0.0; norm[1][1] = 0.0; norm[1][2] = 0.0; norm[1][3] = 0.0;
167  norm[2][0] = 0.0; norm[2][1] = 0.0; norm[2][2] = 0.0; norm[2][3] = 0.0;
168  norm[3][0] = 0.0; norm[3][1] = 0.0; norm[3][2] = 0.0; norm[3][3] = 0.0;
169  rhs1[0] = 0.0; rhs1[1] = 0.0; rhs1[2] = 0.0; rhs1[3] = 0.0;
170  rhs2[0] = 0.0; rhs2[1] = 0.0; rhs2[2] = 0.0; rhs2[3] = 0.0;
171 
172  /* Loop through the grid */
173  for ( nrow=0; nrow<NUM_GRID_ROW; nrow++ )
174  for ( ncol=0; ncol<NUM_GRID_COL; ncol++ )
175  {
176  gline = (sgrid->pgrid[nrow][ncol].l1t_line - sgrid->min_line) / (sgrid->max_line - sgrid->min_line);
177  gsamp = (sgrid->pgrid[nrow][ncol].l1t_samp - sgrid->min_samp) / (sgrid->max_samp - sgrid->min_samp);
178  norm[0][0] += 1.0;
179  norm[0][1] += gline;
180  norm[0][2] += gsamp;
181  norm[0][3] += gline*gsamp;
182  rhs1[0] += sgrid->pgrid[nrow][ncol].ndet;
183  rhs2[0] += sgrid->pgrid[nrow][ncol].frow;
184  norm[1][1] += gline*gline;
185  norm[1][2] += gline*gsamp;
186  norm[1][3] += gline*gline*gsamp;
187  rhs1[1] += gline*sgrid->pgrid[nrow][ncol].ndet;
188  rhs2[1] += gline*sgrid->pgrid[nrow][ncol].frow;
189  norm[2][2] += gsamp*gsamp;
190  norm[2][3] += gsamp*gline*gsamp;
191  rhs1[2] += gsamp*sgrid->pgrid[nrow][ncol].ndet;
192  rhs2[2] += gsamp*sgrid->pgrid[nrow][ncol].frow;
193  norm[3][3] += gline*gsamp*gline*gsamp;
194  rhs1[3] += gline*gsamp*sgrid->pgrid[nrow][ncol].ndet;
195  rhs2[3] += gline*gsamp*sgrid->pgrid[nrow][ncol].frow;
196  }
197 
198  norm[1][0] = norm[0][1];
199  norm[2][0] = norm[0][2];
200  norm[3][0] = norm[0][3];
201  norm[2][1] = norm[1][2];
202  norm[3][1] = norm[1][3];
203  norm[3][2] = norm[2][3];
204 
205  /* Invert the normal equation matrix */
206  if ( simple_inverse( 4, norm ) != SUCCESS )
207  {
208  printf("Error inverting normal equations.\n");
209  return(ERROR);
210  }
211 
212  /* Calculate line/samp to ndet and frow coefficients */
213  for ( nrow=0; nrow<4; nrow++ )
214  {
215  sgrid->ls_to_ndet[nrow] = norm[nrow][0]*rhs1[0] + norm[nrow][1]*rhs1[1] + norm[nrow][2]*rhs1[2] + norm[nrow][3]*rhs1[3];
216  sgrid->ls_to_frow[nrow] = norm[nrow][0]*rhs2[0] + norm[nrow][1]*rhs2[1] + norm[nrow][2]*rhs2[2] + norm[nrow][3]*rhs2[3];
217  }
218 
219  return(SUCCESS);
220 }
221 
222 /* variant of above for older LS */
224  MTL_GRID_SCA_LS *sgrid ) /* I/O: SCA grid array to fit polynomial coefficients for */
225 {
226  int nrow; /* Row index */
227  int ncol; /* Column index */
228  double gline; /* Normalized grid line coordinate */
229  double gsamp; /* Normalized grid sample coordinate */
230  double norm[4][4]; /* Normal equations matrix */
231  double rhs1[4]; /* NDet constant vector */
232  double rhs2[4]; /* FRow constant vector */
233 
234  /* Initialize the normal equations */
235  norm[0][0] = 0.0; norm[0][1] = 0.0; norm[0][2] = 0.0; norm[0][3] = 0.0;
236  norm[1][0] = 0.0; norm[1][1] = 0.0; norm[1][2] = 0.0; norm[1][3] = 0.0;
237  norm[2][0] = 0.0; norm[2][1] = 0.0; norm[2][2] = 0.0; norm[2][3] = 0.0;
238  norm[3][0] = 0.0; norm[3][1] = 0.0; norm[3][2] = 0.0; norm[3][3] = 0.0;
239  rhs1[0] = 0.0; rhs1[1] = 0.0; rhs1[2] = 0.0; rhs1[3] = 0.0;
240  rhs2[0] = 0.0; rhs2[1] = 0.0; rhs2[2] = 0.0; rhs2[3] = 0.0;
241 
242  /* Loop through the grid */
243  for ( nrow=0; nrow<NUM_GRID_ROW_LS; nrow++ )
244  for ( ncol=0; ncol<NUM_GRID_COL; ncol++ )
245  {
246  gline = (sgrid->pgrid[nrow][ncol].l1t_line - sgrid->min_line) / (sgrid->max_line - sgrid->min_line);
247  gsamp = (sgrid->pgrid[nrow][ncol].l1t_samp - sgrid->min_samp) / (sgrid->max_samp - sgrid->min_samp);
248  norm[0][0] += 1.0;
249  norm[0][1] += gline;
250  norm[0][2] += gsamp;
251  norm[0][3] += gline*gsamp;
252  rhs1[0] += sgrid->pgrid[nrow][ncol].ndet;
253  rhs2[0] += sgrid->pgrid[nrow][ncol].frow;
254  norm[1][1] += gline*gline;
255  norm[1][2] += gline*gsamp;
256  norm[1][3] += gline*gline*gsamp;
257  rhs1[1] += gline*sgrid->pgrid[nrow][ncol].ndet;
258  rhs2[1] += gline*sgrid->pgrid[nrow][ncol].frow;
259  norm[2][2] += gsamp*gsamp;
260  norm[2][3] += gsamp*gline*gsamp;
261  rhs1[2] += gsamp*sgrid->pgrid[nrow][ncol].ndet;
262  rhs2[2] += gsamp*sgrid->pgrid[nrow][ncol].frow;
263  norm[3][3] += gline*gsamp*gline*gsamp;
264  rhs1[3] += gline*gsamp*sgrid->pgrid[nrow][ncol].ndet;
265  rhs2[3] += gline*gsamp*sgrid->pgrid[nrow][ncol].frow;
266  }
267 
268  norm[1][0] = norm[0][1];
269  norm[2][0] = norm[0][2];
270  norm[3][0] = norm[0][3];
271  norm[2][1] = norm[1][2];
272  norm[3][1] = norm[1][3];
273  norm[3][2] = norm[2][3];
274 
275  /* Invert the normal equation matrix */
276  if ( simple_inverse( 4, norm ) != SUCCESS )
277  {
278  printf("Error inverting normal equations.\n");
279  return(ERROR);
280  }
281 
282  /* Calculate line/samp to ndet and frow coefficients */
283  for ( nrow=0; nrow<4; nrow++ )
284  {
285  sgrid->ls_to_ndet[nrow] = norm[nrow][0]*rhs1[0] + norm[nrow][1]*rhs1[1] + norm[nrow][2]*rhs1[2] + norm[nrow][3]*rhs1[3];
286  sgrid->ls_to_frow[nrow] = norm[nrow][0]*rhs2[0] + norm[nrow][1]*rhs2[1] + norm[nrow][2]*rhs2[2] + norm[nrow][3]*rhs2[3];
287  }
288 
289  return(SUCCESS);
290 }
291 
293  int n, /* I: Dimension of matrix to invert */
294  double a[4][4] ) /* I/O: Matrix to invert */
295 {
296  int i, j, k; /* Loop indices */
297 
298  for ( k=0; k<n; k++ )
299  {
300  if ( fabs( a[k][k] ) < TINY ) return(ERROR);
301  for ( j=0; j<n; j++ )
302  if ( j != k ) a[k][j] /= a[k][k];
303  a[k][k] = 1.0 / a[k][k];
304  for ( i=0; i<n; i++ )
305  if ( i != k )
306  {
307  for ( j=0; j<n; j++ )
308  if ( j != k ) a[i][j] -= a[i][k]*a[k][j];
309  a[i][k] *= -a[k][k];
310  }
311  }
312 
313  return(SUCCESS);
314 }
316  WRS2 parms, /* I: WRS2 system parameters */
317  SMETA smeta, /* I: Scene metadata structure */
318  double dellon, /* I: Path longitude adjustment (in degrees) */
319  int band_index, /* I: Index for current band */
320  int sca, /* I: SCA number (1 to nsca) */
321  VECTOR ecfsun, /* I: ECF sun vector */
322  MTL_GRID_PT *gpt ) /* I/O: Grid point with only ndet and frow on input */
323 {
324  VECTOR pos; /* Spacecraft ECEF position vector */
325  VECTOR vel; /* Spacecraft ECEF velocity vector */
326  VECTOR slos; /* Spacecraft line-of-sight vector */
327  IAS_VECTOR gpos; /* Ground point ECEF vector */
328  VECTOR ecfview; /* ECF view vector */
329  VECTOR cursun; /* ECF sun vector rotated based on frow */
330  double gp_lat; /* Ground point latitude */
331  double gp_lon; /* Ground point longitude */
332  double proj_x; /* Ground point projection X coordinate */
333  double proj_y; /* Ground point projection Y coordinate */
334  double ecf2lsr[3][3]; /* ECEF to LSR rotation matrix */
335  double sunrot; /* Solar rotation angle due to fractional WRS row */
336  double d2r; /* Conversion from degrees to radians */
337 
338  /* Compute conversion constant */
339  d2r = atan(1.0)/45.0;
340  sunrot = 360.0*d2r*parms.cycle/(double)(parms.numpath*parms.numrow)*gpt->frow;
341 
342  /* Rotate the sun vector to the current time */
343  cursun.x = ecfsun.x*cos(sunrot) + ecfsun.y*sin(sunrot);
344  cursun.y = -ecfsun.x*sin(sunrot) + ecfsun.y*cos(sunrot);
345  cursun.z = ecfsun.z;
346 
347  /* Find the nominal ephemeris state vector */
348  pathrow_to_posvel( parms, smeta.wrs_path, dellon, (double)smeta.wrs_row+gpt->frow, &pos, &vel );
349 
350  /* Calculate the instrument LOS vector */
351  if ( calc_los( smeta.band_smeta[band_index], sca, gpt->ndet, &slos ) != SUCCESS )
352  {
353  printf("Error generating LOS vector for grid point.\n");
354  return(ERROR);
355  }
356 
357  /* Project to ground latitude/longitude */
358  if ( calc_yaw_steering_gp( parms, pos, vel, slos, smeta.roll_angle, &gp_lat, &gp_lon ) < 0 )
359  {
360  printf("Error projecting nominal scene center.\n");
361  return(ERROR);
362  }
363  gp_lat *= d2r;
364  gp_lon *= d2r;
365 
366  /* Convert lat/lon to ECEF */
367  (void)smeta_geodetic_to_ecef( gp_lat, gp_lon, 0.0, &gpos );
368 
369  /* Compute view vector */
370  slos.x = pos.x - gpos.x;
371  slos.y = pos.y - gpos.y;
372  slos.z = pos.z - gpos.z;
373  unitvec( slos, &ecfview );
374 
375  /* Construct the ECEF to LSR rotation matrix */
376  (void)smeta_geodetic_to_ecf2lsr( gp_lat, gp_lon, ecf2lsr );
377 
378  /* Convert vectors to LSR */
379  rotatevec( ecf2lsr, ecfview, &gpt->V );
380  rotatevec( ecf2lsr, cursun, &gpt->S );
381 
382  /* Apply scene map projection */
383  if ( smeta_geodetic_to_proj( smeta.projection, gp_lat, gp_lon, &proj_x, &proj_y ) != SUCCESS )
384  {
385  printf("Error converting lat/lon to projection X/Y.\n");
386  return(ERROR);
387  }
388 
389  /* Reduce to L1T line/sample */
390  if ( smeta_proj_to_l1t( &smeta, smeta.band_smeta[band_index].band, proj_x, proj_y, &gpt->l1t_line, &gpt->l1t_samp ) != SUCCESS )
391  {
392  printf("Errro converting X/Y to L1T line/sample.\n");
393  return(ERROR);
394  }
395 
396  return(SUCCESS);
397 }
398 
399 /* Variant of above for older LS */
401  WRS2 parms, /* I: WRS2 system parameters */
402  SMETA smeta, /* I: Scene metadata structure */
403  double dellon, /* I: Path longitude adjustment (in degrees) */
404  int band_index, /* I: Index for current band */
405  int sca, /* I: SCA number (1 to nsca) */
406  VECTOR ecfsun, /* I: ECF sun vector */
407  MTL_GRID_PT *gpt ) /* I/O: Grid point with only ndet and frow on input */
408 {
409  VECTOR pos; /* Spacecraft ECEF position vector */
410  VECTOR vel; /* Spacecraft ECEF velocity vector */
411  VECTOR slos; /* Spacecraft line-of-sight vector */
412  IAS_VECTOR gpos; /* Ground point ECEF vector */
413  VECTOR ecfview; /* ECF view vector */
414  VECTOR cursun; /* ECF sun vector rotated based on frow */
415  double gp_lat; /* Ground point latitude */
416  double gp_lon; /* Ground point longitude */
417  double proj_x; /* Ground point projection X coordinate */
418  double proj_y; /* Ground point projection Y coordinate */
419  double ecf2lsr[3][3]; /* ECEF to LSR rotation matrix */
420  double sunrot; /* Solar rotation angle due to fractional WRS row */
421  double d2r; /* Conversion from degrees to radians */
422 
423  /* Compute conversion constant */
424  d2r = atan(1.0)/45.0;
425  sunrot = 360.0*d2r*parms.cycle/(double)(parms.numpath*parms.numrow)*gpt->frow;
426 
427  /* Rotate the sun vector to the current time */
428  cursun.x = ecfsun.x*cos(sunrot) + ecfsun.y*sin(sunrot);
429  cursun.y = -ecfsun.x*sin(sunrot) + ecfsun.y*cos(sunrot);
430  cursun.z = ecfsun.z;
431 
432  /* Find the nominal ephemeris state vector */
433  pathrow_to_posvel( parms, smeta.wrs_path, dellon, (double)smeta.wrs_row+gpt->frow, &pos, &vel );
434 
435  /* Calculate the instrument LOS vector */
436  if ( calc_los_ls( smeta.band_smeta_ls[band_index], (double)smeta.wrs_row+gpt->frow, gpt->npix,
437  gpt->ndet, &slos ) != SUCCESS )
438  {
439  printf("Error generating LOS vector for grid point.\n");
440  return(ERROR);
441  }
442 
443  /* Project to ground latitude/longitude */
444  if ( calc_yaw_steering_gp_ls( parms, pos, vel, slos, smeta.roll_angle, &gp_lat, &gp_lon ) < 0 )
445  {
446  printf("Error projecting nominal scene center.\n");
447  return(ERROR);
448  }
449  gp_lat *= d2r;
450  gp_lon *= d2r;
451 
452  /* Convert lat/lon to ECEF */
453  (void)smeta_geodetic_to_ecef( gp_lat, gp_lon, 0.0, &gpos );
454 
455  /* Compute view vector */
456  slos.x = pos.x - gpos.x;
457  slos.y = pos.y - gpos.y;
458  slos.z = pos.z - gpos.z;
459  unitvec( slos, &ecfview );
460 
461  /* Construct the ECEF to LSR rotation matrix */
462  (void)smeta_geodetic_to_ecf2lsr( gp_lat, gp_lon, ecf2lsr );
463 
464  /* Convert vectors to LSR */
465  rotatevec( ecf2lsr, ecfview, &gpt->V );
466  rotatevec( ecf2lsr, cursun, &gpt->S );
467 
468  /* Apply scene map projection */
469  if ( smeta_geodetic_to_proj( smeta.projection, gp_lat, gp_lon, &proj_x, &proj_y ) != SUCCESS )
470  {
471  printf("Error converting lat/lon to projection X/Y.\n");
472  return(ERROR);
473  }
474 
475  /* Reduce to L1T line/sample */
476  if ( smeta_proj_to_l1t( &smeta, smeta.band_smeta_ls[band_index].band, proj_x, proj_y, &gpt->l1t_line, &gpt->l1t_samp ) != SUCCESS )
477  {
478  printf("Errro converting X/Y to L1T line/sample.\n");
479  return(ERROR);
480  }
481 
482  return(SUCCESS);
483 }
484 
486  double line, /* I: L1T line coordinate */
487  double samp, /* I: L1T sample coordinate */
488  MTL_GRID_BAND bgrid, /* I: MTL grid structure for current band */
489  /* Sudipta new addition for OLI to extract SCA and Det numbers */
490  short *sca_num, /* O: SCA containing this point */
491  short *det_num, /* O: Detector number for this point */
492  /* End Sudipta new addition for OLI */
493  double *sat_zn, /* O: Viewing zenith angle scaled to 0.01 degrees */
494  double *sat_az, /* O: Viewing azimuth angle scaled to 0.01 degrees */
495  double *sun_zn, /* O: Solar zenith angle scaled to 0.01 degrees */
496  double *sun_az ) /* O: Solar azimuth angle scaled to 0.01 degrees */
497 {
498  int sca; /* SCA index */
499  int nsca; /* Number of SCAs that see this point */
500  double sline; /* Scaled line number */
501  double ssamp; /* Scaled sample number */
502  double ndet; /* Normalized detector coordinate */
503  double frow; /* Fractional row coordinate */
504  double w0, w1, w2, w3; /* Interpolation weights */
505  double sat_zen; /* Floating point view zenith angle */
506  double sat_azm; /* Floating point view azimuth angle */
507  double sun_zen; /* Floating point sun zenith angle */
508  double sun_azm; /* Floating point sun azimuth angle */
509  double hdist; /* Horizontal component of vector */
510  int cell_row; /* Index of grid cell row */
511  int cell_col; /* Index of grid cell column */
512  VECTOR sat[2]; /* View vector */
513  VECTOR sun[2]; /* Solar vector */
514  /* Sudipta new addition for OLI to extract SCA and Det numbers */
515  short curr_sca[2]; /* Current SCA number(s) */
516  short curr_det[2]; /* Current detector number(s) */
517  /* End Sudipta new addition for OLI */
518  double r2d00; /* Units conversion factor */
519 
520  /* Scaling factor */
521  r2d00 = 45.0/atan(1.0);
522 
523  /* Initialize SCA count and loop through SCAs */
524  nsca = 0;
525  for ( sca=0; sca<bgrid.nsca; sca++ )
526  {
527 // if ( nsca > 1 ) continue;
528  if ( line < bgrid.sgrid[sca].min_line ) continue;
529  if ( line > bgrid.sgrid[sca].max_line ) continue;
530  if ( samp < bgrid.sgrid[sca].min_samp ) continue;
531  if ( samp > bgrid.sgrid[sca].max_samp ) continue;
532 
533  /* Calculate normalized detector */
534  sline = (line - bgrid.sgrid[sca].min_line) / (bgrid.sgrid[sca].max_line - bgrid.sgrid[sca].min_line);
535  ssamp = (samp - bgrid.sgrid[sca].min_samp) / (bgrid.sgrid[sca].max_samp - bgrid.sgrid[sca].min_samp);
536  ndet = bgrid.sgrid[sca].ls_to_ndet[0]
537  + bgrid.sgrid[sca].ls_to_ndet[1]*sline
538  + bgrid.sgrid[sca].ls_to_ndet[2]*ssamp
539  + bgrid.sgrid[sca].ls_to_ndet[3]*sline*ssamp;
540  if ( ndet < -1.0 || ndet > 1.0 ) continue;
541 
542  /* Calculate fractional row */
543  frow = bgrid.sgrid[sca].ls_to_frow[0]
544  + bgrid.sgrid[sca].ls_to_frow[1]*sline
545  + bgrid.sgrid[sca].ls_to_frow[2]*ssamp
546  + bgrid.sgrid[sca].ls_to_frow[3]*sline*ssamp;
547  if ( frow < -0.7 || frow > 0.7 ) continue;
548 
549  /* Sudipta new addition for OLI to extract SCA and Det numbers */
550  /* Capture the current SCA and detector numbers */
551  curr_sca[nsca] = sca + 1;
552  curr_det[nsca] = (short)floor((ndet + 1.0)/2.0 * (double)bgrid.num_detectors + 0.5);
553  /* End Sudipta new addition for OLI */
554 
555  /* Figure out which grid cell we are in */
556  cell_col = (int)floor( (ndet - GRID_COL_BASE)/GRID_COL_INC );
557  if ( cell_col > (NUM_GRID_COL-2) ) continue;
558  cell_row = (int)floor( (frow - GRID_ROW_BASE)/GRID_ROW_INC );
559  if ( cell_row > (NUM_GRID_ROW-2) ) continue;
560 
561  /* Calculate offsets within the grid cell */
562  ndet = (ndet - GRID_COL_BASE)/GRID_COL_INC - (double)cell_col;
563  frow = (frow - GRID_ROW_BASE)/GRID_ROW_INC - (double)cell_row;
564 
565  /* Interpolate vectors */
566  w0 = (1.0-ndet)*(1.0-frow);
567  w1 = (1.0-ndet)*frow;
568  w2 = ndet*(1.0-frow);
569  w3 = ndet*frow;
570  sat[nsca].x = bgrid.sgrid[sca].pgrid[cell_row][cell_col].V.x *w0
571  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].V.x *w1
572  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].V.x *w2
573  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].V.x*w3;
574  sat[nsca].y = bgrid.sgrid[sca].pgrid[cell_row][cell_col].V.y *w0
575  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].V.y *w1
576  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].V.y *w2
577  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].V.y*w3;
578  sat[nsca].z = sqrt( 1.0 - sat[nsca].x*sat[nsca].x - sat[nsca].y*sat[nsca].y );
579  sun[nsca].x = bgrid.sgrid[sca].pgrid[cell_row][cell_col].S.x *w0
580  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].S.x *w1
581  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].S.x *w2
582  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].S.x*w3;
583  sun[nsca].y = bgrid.sgrid[sca].pgrid[cell_row][cell_col].S.y *w0
584  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].S.y *w1
585  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].S.y *w2
586  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].S.y*w3;
587  sun[nsca].z = sqrt( 1.0 - sun[nsca].x*sun[nsca].x - sun[nsca].y*sun[nsca].y );
588  nsca++;
589  }
590 
591  /* Set angles to zero if point not in scene */
592  if ( nsca < 1 )
593  {
594  *sat_zn = 0;
595  *sat_az = 0;
596  *sun_zn = 0;
597  *sun_az = 0;
598  /* Sudipta new addition for OLI to extract SCA and Det numbers */
599  *sca_num = 0;
600  *det_num = 0;
601  /* End Sudipta new addition for OLI */
602  return(SUCCESS);
603  }
604 
605  /* Calculate the angles */
606  sat_zen = 0.0; sat_azm = 0.0; sun_zen = 0.0; sun_azm = 0.0;
607  for ( sca=0; sca<nsca; sca++ )
608  {
609  /* Reduce vectors to angles */
610  if ( sat[sca].z > 0.0 ) sat_zen += acos( sat[sca].z );
611  hdist = sat[sca].x*sat[sca].x + sat[sca].y*sat[sca].y;
612  if ( hdist > 0.0 ) sat_azm += atan2( sat[sca].x, sat[sca].y );
613  if ( sun[sca].z > 0.0 ) sun_zen += acos( sun[sca].z );
614  hdist = sun[sca].x*sun[sca].x + sun[sca].y*sun[sca].y;
615  if ( hdist > 0.0 ) sun_azm += atan2( sun[sca].x, sun[sca].y );
616  }
617  sat_zen /= (double)nsca;
618  sat_azm /= (double)nsca;
619  sun_zen /= (double)nsca;
620  sun_azm /= (double)nsca;
621 
622  *sat_zn = r2d00*sat_zen;
623  *sat_az = r2d00*sat_azm;
624  *sun_zn = r2d00*sun_zen;
625  *sun_az = r2d00*sun_azm;
626  /* Sudipta new addition for OLI to extract SCA and Det numbers */
627  *sca_num = curr_sca[0];
628  *det_num = curr_det[0];
629 
630  if ( nsca > 1 )
631  {
632  if ( MIN( curr_det[0], bgrid.num_detectors-curr_det[0] ) < MIN( curr_det[1], bgrid.num_detectors-curr_det[1] ) )
633  {
634  *sca_num = curr_sca[1];
635  *det_num = curr_det[1];
636  }
637  }
638  /* End Sudipta new addition for OLI */
639 
640  return(SUCCESS);
641 }
642 
643 /* varinat of smeta_angles for older LS */
645  double line, /* I: L1T line coordinate */
646  double samp, /* I: L1T sample coordinate */
647  MTL_GRID_BAND_LS bgrid, /* I: MTL grid structure for current band */
648  short *sat_zn, /* O: Viewing zenith angle scaled to 0.01 degrees */
649  short *sat_az, /* O: Viewing azimuth angle scaled to 0.01 degrees */
650  short *sun_zn, /* O: Solar zenith angle scaled to 0.01 degrees */
651  short *sun_az ) /* O: Solar azimuth angle scaled to 0.01 degrees */
652 {
653  int sca; /* SCA index */
654  int nsca; /* Number of SCAs that see this point */
655  double sline; /* Scaled line number */
656  double ssamp; /* Scaled sample number */
657  double ndet; /* Normalized detector coordinate */
658  double frow; /* Fractional row coordinate */
659  double w0, w1, w2, w3; /* Interpolation weights */
660  double sat_zen; /* Floating point view zenith angle */
661  double sat_azm; /* Floating point view azimuth angle */
662  double sun_zen; /* Floating point sun zenith angle */
663  double sun_azm; /* Floating point sun azimuth angle */
664  double hdist; /* Horizontal component of vector */
665  int cell_row; /* Index of grid cell row */
666  int cell_col; /* Index of grid cell column */
667  VECTOR sat[2]; /* View vector */
668  VECTOR sun[2]; /* Solar vector */
669  double r2d00; /* Units conversion factor */
670 
671  /* Scaling factor */
672  r2d00 = 4500.0/atan(1.0);
673 
674  /* Initialize SCA count and loop through SCAs */
675  nsca = 0;
676  for ( sca=0; sca<bgrid.nsca; sca++ )
677  {
678  if ( nsca > 1 ) continue;
679  if ( line < bgrid.sgrid[sca].min_line ) continue;
680  if ( line > bgrid.sgrid[sca].max_line ) continue;
681  if ( samp < bgrid.sgrid[sca].min_samp ) continue;
682  if ( samp > bgrid.sgrid[sca].max_samp ) continue;
683 
684  /* Calculate normalized detector */
685  sline = (line - bgrid.sgrid[sca].min_line) / (bgrid.sgrid[sca].max_line - bgrid.sgrid[sca].min_line);
686  ssamp = (samp - bgrid.sgrid[sca].min_samp) / (bgrid.sgrid[sca].max_samp - bgrid.sgrid[sca].min_samp);
687  ndet = bgrid.sgrid[sca].ls_to_ndet[0]
688  + bgrid.sgrid[sca].ls_to_ndet[1]*sline
689  + bgrid.sgrid[sca].ls_to_ndet[2]*ssamp
690  + bgrid.sgrid[sca].ls_to_ndet[3]*sline*ssamp;
691  if ( ndet < -1.0 || ndet > 1.0 ) continue;
692 
693  /* Calculate fractional row */
694  frow = bgrid.sgrid[sca].ls_to_frow[0]
695  + bgrid.sgrid[sca].ls_to_frow[1]*sline
696  + bgrid.sgrid[sca].ls_to_frow[2]*ssamp
697  + bgrid.sgrid[sca].ls_to_frow[3]*sline*ssamp;
698  if ( frow < -0.6 || frow > 0.6 ) continue;
699 
700  /* Figure out which grid cell we are in */
701  cell_col = (int)floor( (ndet - GRID_COL_BASE)/GRID_COL_INC );
702  if ( cell_col > (NUM_GRID_COL-2) ) continue;
703  cell_row = (int)floor( (frow - GRID_ROW_BASE_LS)/GRID_ROW_INC );
704  if ( cell_row > (NUM_GRID_ROW_LS-2) ) continue;
705 
706  /* Calculate offsets within the grid cell */
707  ndet = (ndet - GRID_COL_BASE)/GRID_COL_INC - (double)cell_col;
708  frow = (frow - GRID_ROW_BASE_LS)/GRID_ROW_INC - (double)cell_row;
709 
710  /* Interpolate vectors */
711  w0 = (1.0-ndet)*(1.0-frow);
712  w1 = (1.0-ndet)*frow;
713  w2 = ndet*(1.0-frow);
714  w3 = ndet*frow;
715  sat[nsca].x = bgrid.sgrid[sca].pgrid[cell_row][cell_col].V.x *w0
716  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].V.x *w1
717  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].V.x *w2
718  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].V.x*w3;
719  sat[nsca].y = bgrid.sgrid[sca].pgrid[cell_row][cell_col].V.y *w0
720  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].V.y *w1
721  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].V.y *w2
722  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].V.y*w3;
723  sat[nsca].z = sqrt( 1.0 - sat[nsca].x*sat[nsca].x - sat[nsca].y*sat[nsca].y );
724  sun[nsca].x = bgrid.sgrid[sca].pgrid[cell_row][cell_col].S.x *w0
725  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].S.x *w1
726  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].S.x *w2
727  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].S.x*w3;
728  sun[nsca].y = bgrid.sgrid[sca].pgrid[cell_row][cell_col].S.y *w0
729  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col].S.y *w1
730  + bgrid.sgrid[sca].pgrid[cell_row][cell_col+1].S.y *w2
731  + bgrid.sgrid[sca].pgrid[cell_row+1][cell_col+1].S.y*w3;
732  sun[nsca].z = sqrt( 1.0 - sun[nsca].x*sun[nsca].x - sun[nsca].y*sun[nsca].y );
733  nsca++;
734  }
735 
736  /* Set angles to zero if point not in scene */
737  if ( nsca < 1 )
738  {
739  *sat_zn = 0;
740  *sat_az = 0;
741  *sun_zn = 0;
742  *sun_az = 0;
743  return(SUCCESS);
744  }
745 
746  /* Calculate the angles */
747  sat_zen = 0.0; sat_azm = 0.0; sun_zen = 0.0; sun_azm = 0.0;
748  for ( sca=0; sca<nsca; sca++ )
749  {
750  /* Reduce vectors to angles */
751  if ( sat[sca].z > 0.0 ) sat_zen += acos( sat[sca].z );
752  hdist = sat[sca].x*sat[sca].x + sat[sca].y*sat[sca].y;
753  if ( hdist > 0.0 ) sat_azm += atan2( sat[sca].x, sat[sca].y );
754  if ( sun[sca].z > 0.0 ) sun_zen += acos( sun[sca].z );
755  hdist = sun[sca].x*sun[sca].x + sun[sca].y*sun[sca].y;
756  if ( hdist > 0.0 ) sun_azm += atan2( sun[sca].x, sun[sca].y );
757  }
758  sat_zen /= (double)nsca;
759  sat_azm /= (double)nsca;
760  sun_zen /= (double)nsca;
761  sun_azm /= (double)nsca;
762 
763  *sat_zn = (short)floor( r2d00*sat_zen + 0.5 );
764  *sat_az = (short)floor( r2d00*sat_azm + 0.5 );
765  *sun_zn = (short)floor( r2d00*sun_zen + 0.5 );
766  *sun_az = (short)floor( r2d00*sun_azm + 0.5 );
767 
768  return(SUCCESS);
769 }
double max_samp
Definition: mtl_grid.h:52
#define MIN(x, y)
Definition: rice.h:169
double frow
Definition: mtl_grid.h:27
int band
Definition: smeta.h:82
SMETA_SCENE_PROJ projection
Definition: smeta.h:130
#define SUCCESS
Definition: ObpgReadGrid.h:15
double y
Definition: mtl_geometry.h:38
int j
Definition: decode_rs.h:73
void pathrow_to_posvel(WRS2 parms, int path, double dellon, double wrsrow, VECTOR *pos, VECTOR *vel)
int calc_los(SMETA_BAND band_smeta, int sca, double ndet, VECTOR *los)
Definition: mtl_geometry.c:19
MTL_GRID_SCA_LS * sgrid
Definition: mtl_grid.h:72
int numrow
Definition: mtl_geometry.h:30
int calc_band_grid(WRS2 parms, SMETA smeta, double dellon, VECTOR ecfsun, int band_index, MTL_GRID_BAND *bgrid)
Definition: mtl_grid.c:14
MTL_GRID_PT pgrid[NUM_GRID_ROW][NUM_GRID_COL]
Definition: mtl_grid.h:41
int smeta_angles_ls(double line, double samp, MTL_GRID_BAND_LS bgrid, short *sat_zn, short *sat_az, short *sun_zn, short *sun_az)
Definition: mtl_grid.c:644
#define NULL
Definition: decode_rs.h:63
MTL_GRID_PT pgrid[NUM_GRID_ROW_LS][NUM_GRID_COL]
Definition: mtl_grid.h:53
double min_samp
Definition: mtl_grid.h:39
float32 * pos
Definition: l1_czcs_hdf.c:35
void unitvec(VECTOR a, VECTOR *b)
double max_line
Definition: mtl_grid.h:38
double min_line
Definition: mtl_grid.h:49
VECTOR V
Definition: mtl_grid.h:30
VECTOR S
Definition: mtl_grid.h:31
#define GRID_COL_BASE
Definition: mtl_grid.h:21
int calc_grid_fit_ls(MTL_GRID_SCA_LS *sgrid)
Definition: mtl_grid.c:223
int smeta_proj_to_l1t(SMETA *smeta, int band, double proj_x, double proj_y, double *l1t_line, double *l1t_samp)
int simple_inverse(int n, double a[4][4])
Definition: mtl_grid.c:292
int num_detectors
Definition: mtl_grid.h:63
double ls_to_frow[4]
Definition: mtl_grid.h:55
int wrs_path
Definition: smeta.h:124
#define NUM_GRID_ROW_LS
Definition: mtl_grid.h:16
#define TINY
Definition: mtl_grid.c:12
#define NUM_GRID_COL
Definition: mtl_grid.h:14
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 band_smeta[IAS_MAX_NBANDS]
Definition: smeta.h:132
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
int numpath
Definition: mtl_geometry.h:29
int band
Definition: smeta.h:97
double roll_angle
Definition: smeta.h:126
double ls_to_ndet[4]
Definition: mtl_grid.h:42
double l1t_line
Definition: mtl_grid.h:28
double ndet
Definition: mtl_grid.h:26
int cycle
Definition: mtl_geometry.h:31
int smeta_geodetic_to_ecf2lsr(double lat, double lon, double ecf2lsr[3][3])
#define GRID_COL_INC
Definition: mtl_grid.h:17
void rotatevec(double mat[3][3], VECTOR a, VECTOR *b)
integer, parameter double
double max_samp
Definition: mtl_grid.h:40
MTL_GRID_SCA * sgrid
Definition: mtl_grid.h:65
#define GRID_ROW_INC
Definition: mtl_grid.h:18
int calc_grid_fit(MTL_GRID_SCA *sgrid)
Definition: mtl_grid.c:153
int calc_band_grid_ls(WRS2 parms, SMETA smeta, double dellon, VECTOR ecfsun, int band_index, MTL_GRID_BAND_LS *bgrid)
Definition: mtl_grid.c:91
double x
Definition: mtl_geometry.h:37
int smeta_angles(double line, double samp, MTL_GRID_BAND bgrid, short *sca_num, short *det_num, double *sat_zn, double *sat_az, double *sun_zn, double *sun_az)
Definition: mtl_grid.c:485
int num_detectors
Definition: smeta.h:87
int nsca
Definition: smeta.h:100
SMETA_BAND_LS band_smeta_ls[IAS_MAX_NBANDS_LS]
Definition: smeta.h:133
double min_line
Definition: mtl_grid.h:37
int calc_grid_point(WRS2 parms, SMETA smeta, double dellon, int band_index, int sca, VECTOR ecfsun, MTL_GRID_PT *gpt)
Definition: mtl_grid.c:315
double npix
Definition: mtl_grid.h:25
#define fabs(a)
Definition: misc.h:93
double l1t_samp
Definition: mtl_grid.h:29
double min_samp
Definition: mtl_grid.h:51
int calc_grid_point_ls(WRS2 parms, SMETA smeta, double dellon, int band_index, int sca, VECTOR ecfsun, MTL_GRID_PT *gpt)
Definition: mtl_grid.c:400
int smeta_geodetic_to_ecef(double lat, double lon, double hgt, IAS_VECTOR *ecef)
int smeta_geodetic_to_proj(SMETA_SCENE_PROJ proj, double lat, double lon, double *proj_x, double *proj_y)
int nsca
Definition: smeta.h:85
#define GRID_ROW_BASE
Definition: mtl_grid.h:19
double z
Definition: mtl_geometry.h:39
double ls_to_ndet[4]
Definition: mtl_grid.h:54
double max_line
Definition: mtl_grid.h:50
Definition: smeta.h:118
int i
Definition: decode_rs.h:71
double ls_to_frow[4]
Definition: mtl_grid.h:43
int wrs_row
Definition: smeta.h:125
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
#define GRID_ROW_BASE_LS
Definition: mtl_grid.h:20
int calc_los_ls(SMETA_BAND_LS band_smeta, double frow, double npix, double ndet, VECTOR *los)
Definition: mtl_geometry.c:52
int k
Definition: decode_rs.h:73
#define NUM_GRID_ROW
Definition: mtl_grid.h:15
#define ERROR
Definition: ancil.h:24