NASA Ocean Color

ocssw V2022

 /* 
  *Revision History:
  * $Log: GEO_hires.c,v $
  * Revision 6.5  2010/12/14 21:12:18  kuyper
  * Split samp into osamp and psamp, thereby simplifying the logic, enabling
  *   identification and correction of some off-by-one errors.
  *
  * Revision 6.4  2010/06/18 20:38:31  kuyper
  * Corrected to initialize 'triplets'.
  *
  * Revision 6.3  2009/06/12 16:01:46  kuyper
  * Corrected sign error in formula for the interpolation factor in the scan
  *   direction.
  * Corrected handling of the case where the first call to GEO_calculate_offsets
  *   fails.
  *
  * Revision 6.2  2009/05/22 18:15:20  kuyper
  * Removed design notes; the ones in the PDL are sufficient.
  * Added sample_flags to invalid input arguments message.
  * Corrected several misuses of N_samp.
  *
  * Revision 6.1  2009/04/28 18:05:37  kuyper
  * Initial revision.
  *
  * James Kuyper James.R.Kuyper@nasa.gov
  */
  
 #include "PGS_CSC.h"
 #include "smfio.h"
 #include "GEO_earth.h"
 #include "GEO_geo.h"
 #include "PGS_MODIS_35251.h"
  
 #define KILO 1000.0
  
 static double GEO_triple_product(
     double x[3],
     double y[3],
     double z[3]
 )
 {
     return x[0] * (y[1] * z[2] - y[2] * z[1]) +
        x[1] * (y[2] * z[0] - y[0] * z[2]) +
        x[2] * (y[0] * z[1] - y[1] * z[0]);
 }
  
 static int GEO_calculate_offsets(
     double      hires_scale,
     double      num,
     double      den,
     double      nst,
     double      nds,
     double      ndt,
     double      nqs,
     double      nqt,
     const double    delta[],
     const double    quad[],
     const double    normal[],
     const double    scan[],
     const double    track[],
     int8        sth[]
 )
  
  
 {
     int ind;        /* Loop counter */
     double doff[3];     /* temporary offsets */
     double den1, den2;  /* possible denominators */
     const double max_offset = (-0.5 - HIRES_FVALUE) * hires_scale;
  
     if(fabs(den) <= fabs(num)/DBL_MAX)
     return 0;
  
     doff[0] = num/den;
     if(fabs(doff[0]) >= max_offset)
     return 0;
  
     den1 = doff[0]*nqs - nst;
     den2 = doff[0]*nqt;
     if(fabs(den1) > fabs(den2))
     {
     num = nds;
     den = den1;
     }
     else
     {
     num = ndt - doff[0]*nst;
     den = den2;
     }
  
     if(fabs(den) <= fabs(num)/DBL_MAX)
     return 0;
  
     doff[1] = num/den;
  
     doff[2] = 0.0;
     for(ind=0; ind<3; ind++)
     doff[2] += normal[ind] * (delta[ind] - doff[0]*scan[ind] -
         doff[1]*track[ind] - doff[0]*doff[1]*quad[ind]) / (KILO * KILO);
  
     if(fabs(doff[2]) >= max_offset)
     return 0;
  
     for(ind=0; ind<3; ind++)
     sth[ind] = (int)floor(doff[ind]/hires_scale + 0.5);
  
     return 1;
 }
  
  
 PGSt_SMF_status GEO_hires(
         uint16  N_samp,
     int padded_samples,
     double  hires_scale,
         double  terrain_sample_position[][MAX_PADDED][3],
         double  ecr_sample_position[][MAX_PADDED][3],
         double  ecr_frame_position[][MAX_FRAMES][3],
     uint8   sample_flags[][MAX_PADDED],
     uint8   frame_flags[][MAX_FRAMES],
     int8    hires_offsets[][DETECTORS_QKM][SAMPLES_QKM]
 )
 {
     /* loop counters to identify which high resolution pixel is being
      * processed.
      */
     int hdet;
     int samples = padded_samples - N_samp + 1;
     int frames = (samples+1)/N_samp;
  
     if(!terrain_sample_position || !ecr_sample_position || !ecr_frame_position
     || !sample_flags || !frame_flags || !hires_offsets)
     {
     char filefunc[] = __FILE__ ", GEO_hires()";
     char msgbuf[PGS_SMF_MAX_MSGBUF_SIZE];
  
     sprintf(msgbuf, "terrain_sample_position=%p, "
         "ecr_sample_position=%p, ecr_frame_position=%p, sample_flags=%p, "
         "frame_flags=%p, hires_offsets=%p", (void*)terrain_sample_position,
         (void*)ecr_sample_position, (void*) ecr_frame_position,
         (void*)sample_flags, (void*)frame_flags, (void*)hires_offsets);
     modsmf(MODIS_E_BAD_INPUT_ARG, msgbuf, filefunc);
  
     return MODIS_E_BAD_INPUT_ARG;
     }
  
     for(hdet = 0; hdet < DETECTORS_1KM * N_samp; hdet++)
     {
     int ind;  /* loop counter to index across coordinates. */
     /* The addition of N_samp cancels final -1, mathematically , but in C
      * it is needed to ensure the same behavior. */
     int det = (hdet + N_samp - N_samp/2)/N_samp - 1;
     double dt;  /* Fractional offset in track direction. */
     int osamp;  /* output sample number */
  
     if(det < 0)
         det = 0;
     else if(det > DETECTORS_1KM - 2)
         det = DETECTORS_1KM - 2;
  
     dt = (hdet + 0.5)/N_samp - (det + 0.5);
  
     /* The bracketing detectors are det and det+1 */
     for(osamp =0; osamp < samples; osamp++)
     {   /* Determine which 1km frames to use for the    */
         /* estimate, and the corresponding values of ds.*/
         double  estimated[3], delta[3], scan[3], track[3];
         double normal[3], quad[3];
         double nds, nqs, nst, ndt, nqt;  /* triple products */
         double a, b, c, q;  /* quadratic coefficients. */
         int8 sth[2][3]; /* tentative s/t/h triplets. */
         int triplets=0; /* number of triplets selected. */
         double ds;      /* Fractional offset in scan direction */
         /* sample number in padded hi-resolution arrays. */
         int psamp = osamp + N_samp -1;
         int frame = osamp/N_samp;
  
         if(frame > frames - 2)
         frame = frames - 2;
         ds = osamp/(double)N_samp - frame;
         /* The bracketing frames are frame and frame+1 */
  
         if(sample_flags[hdet ][psamp  ] >= INVALID_INPUT_DATA || 
         frame_flags[det  ][frame  ] >= INVALID_INPUT_DATA ||
         frame_flags[det  ][frame+1] >= INVALID_INPUT_DATA ||
         frame_flags[det+1][frame  ] >= INVALID_INPUT_DATA ||
         frame_flags[det+1][frame+1] >= INVALID_INPUT_DATA)
         {
         for(ind = 0; ind<3; ind++) /* hires_offsets is NOT padded */
             hires_offsets[ind][hdet][osamp] = HIRES_FVALUE;
         }
  
         /*Perform bilinear interpolation/extrapolation */
  
         for(ind =0; ind<3; ind++)
         {
         estimated[ind] =
             (1 - dt)*((1 - ds)*ecr_frame_position[det  ][frame  ][ind] +
                    ds *ecr_frame_position[det  ][frame+1][ind])+
              dt *((1 - ds)*ecr_frame_position[det+1][frame  ][ind] +
                    ds *ecr_frame_position[det+1][frame+1][ind]);
         delta[ind] =
             ecr_sample_position[hdet][psamp][ind] - estimated[ind];
         scan[ind] = 
             (1 - dt)*(ecr_frame_position[det  ][frame+1][ind]  -
                   ecr_frame_position[det  ][frame  ][ind]) +
              dt *(ecr_frame_position[det+1][frame+1][ind]  -
                   ecr_frame_position[det+1][frame  ][ind]);
         track[ind] = 
             (1 - ds)*(ecr_frame_position[det+1][frame  ][ind]  -
                   ecr_frame_position[det  ][frame  ][ind]) +
              ds *(ecr_frame_position[det+1][frame+1][ind]  -
                   ecr_frame_position[det  ][frame+1][ind]);
         quad[ind] = ecr_frame_position[det  ][frame  ][ind]-
                 ecr_frame_position[det  ][frame+1][ind]-
                 ecr_frame_position[det+1][frame  ][ind]+
                 ecr_frame_position[det+1][frame+1][ind];
         }
  
         normal[0] = KILO *
         cos(terrain_sample_position[hdet][psamp][LAT_IDX]) *
         cos(terrain_sample_position[hdet][psamp][LON_IDX]);
         normal[1] = KILO *
         cos(terrain_sample_position[hdet][psamp][LAT_IDX]) *
         sin(terrain_sample_position[hdet][psamp][LON_IDX]);
         normal[2] = KILO *
         sin(terrain_sample_position[hdet][psamp][LAT_IDX]);
  
         nds = GEO_triple_product(normal, delta, scan);
         ndt = GEO_triple_product(normal, delta, track);
         nst = GEO_triple_product(normal, scan, track);
         nqs = GEO_triple_product(normal, quad, scan);
         nqt = GEO_triple_product(normal, quad, track);
  
         /* Calculate quadratic coefficients */
         a = nst * nqs;
         b = nds * nqt - nst * nst - ndt * nqs;
         c = ndt * nst;
  
         if(b*b >= 4.0*a*c)
         {
         q = -0.5*(b + (b<0 ? -1 : 1)*sqrt(b*b - 4.0*a*c));
  
         triplets = GEO_calculate_offsets(hires_scale, q, a, nst, nds,
             ndt, nqs, nqt, delta, quad, normal, scan, track, sth[0]);
         triplets += GEO_calculate_offsets(hires_scale, c, q, nst, nds,
             ndt, nqs, nqt, delta, quad, normal, scan, track,
             sth[triplets]);
         }
  
         if(triplets > 1 && 
         sth[1][0]*sth[1][0]+sth[1][1]*sth[1][1]+sth[1][2]*sth[1][2] <
         sth[0][0]*sth[0][0]+sth[0][1]*sth[0][1]+sth[0][2]*sth[0][2]) 
             memcpy(sth[0], sth[1], sizeof sth[0]);
  
         if(triplets > 0)
         {
         for(ind=0; ind < 3; ind++) /* hires_offsets is NOT padded */
             hires_offsets[ind][hdet][osamp] = sth[0][ind];
         }
         else
         {
         for(ind=0; ind < 3; ind++) /* hires_offsets is NOT padded */
             hires_offsets[ind][hdet][osamp] = HIRES_FVALUE;
         }
     }
     }
  
     return PGS_S_SUCCESS;
 }
  

data_t q

Definition: decode_rs.h:74

GEO_earth.h

MODIS_E_BAD_INPUT_ARG

#define MODIS_E_BAD_INPUT_ARG

Definition: PGS_MODIS_35251.h:4

MAX_FRAMES

#define MAX_FRAMES

Definition: GEO_geo.h:79

scan

MOD_PR01 Production producing one five minute granule of output data in each run It can be configured to produce as many as three five minute granules per run Each execution with one construction record and one date file for each dataset In normal these are created by which splits them out of the hour datasets For LANCE they are created by which merges all session MODIS L0 datasets overlapping the requested time and extracts from the merged data those packets which fall within that time period Each scan of data is stored in the L1A granule that covers the start time of that scan

Definition: MOD_PR01_pr.txt:19

MAX_PADDED

#define MAX_PADDED

Definition: GEO_geo.h:84

hico.value_locate.x

list x

Definition: value_locate.py:64

color_dtdb.y

Definition: color_dtdb.py:430

delta

const double delta

Definition: VcstCmnConsts.h:523

generalauxtype::double

integer, parameter double

Definition: GeneralAuxType.f90:27

KILO

#define KILO

Definition: GEO_hires.c:34

DETECTORS_1KM

#define DETECTORS_1KM

Definition: GEO_geo.h:85

LON_IDX

@ LON_IDX

Definition: VcstGeoParameters.h:161

PGS_MODIS_35251.h

data_t b[NROOTS+1]

Definition: decode_rs.h:77

LAT_IDX

@ LAT_IDX

Definition: VcstGeoParameters.h:160

HIRES_FVALUE

#define HIRES_FVALUE

Definition: GEO_geo.h:157

fabs

#define fabs(a)

Definition: misc.h:93

den

data_t den

Definition: decode_rs.h:74

GEO_hires

PGSt_SMF_status GEO_hires(uint16 N_samp, int padded_samples, double hires_scale, double terrain_sample_position[][MAX_PADDED][3], double ecr_sample_position[][MAX_PADDED][3], double ecr_frame_position[][MAX_FRAMES][3], uint8 sample_flags[][MAX_PADDED], uint8 frame_flags[][MAX_FRAMES], int8 hires_offsets[][DETECTORS_QKM][SAMPLES_QKM])

Definition: GEO_hires.c:205

GEO_geo.h

DBL_MAX

#define DBL_MAX

Definition: make_L3_v1.1.c:60

INVALID_INPUT_DATA

const unsigned char INVALID_INPUT_DATA

Definition: VcstGeoParameters.h:148

SAMPLES_QKM

#define SAMPLES_QKM

Definition: GEO_geo.h:83

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

DETECTORS_QKM

#define DETECTORS_QKM

Definition: GEO_geo.h:87

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