|
ocssw
1.0
|
00001 #include "PGS_CSC.h" 00002 #include "PGS_MODIS_35251.h" 00003 #include "GEO_earth.h" 00004 #include "GEO_global_arrays.h" 00005 #include "GEO_util.h" 00006 #include "smfio.h" 00007 00008 int GEO_ellip_position( 00009 int const scan_number, 00010 int const sample_number, 00011 int const num_detectors, 00012 double u_inst[][3], 00013 double ecr_sc_sample_position[][3], 00014 double ecr_sc_velocity[][3], 00015 double T_inst2ecr[][3][3], 00016 double ecr_sample_position[][MAX_PADDED][3], 00017 double ellip_sample_position[][MAX_PADDED][3], 00018 unsigned char sample_flags[][MAX_PADDED], 00019 double sample_view_vec[][MAX_PADDED][3] 00020 ) 00021 /* 00022 !C***************************************************************************** 00023 !Description: 00024 Routine in Earth Location group of the Level-1A geolocation software to 00025 calculate the ellipsoidal position of a pixel. It transforms each view 00026 vector to the ECR frame and solves for the magnitude of the vector 00027 (See Geolocation ATBD section 3.1.3.2). Finally, it computes the 00028 geodetic latitude and longitude for the pixel and stores them in the 00029 global array. If the view vector for a detector does not intersect the 00030 earth ellipsoid, no ecr_sample_position nor ellip_sample_position will 00031 be determined, and the sample_view_vec will not be rescaled for that 00032 detector. 00033 00034 GEO_ellip_position initializes the terrain sample position validity 00035 sample_flags, setting NO_ELLIPSE_INTERSECT for detectors with view 00036 vectors that do not intersect the earth ellipsoid. The 00037 INVALID_INPUT_DATA flag is not set if a routine failure occurs, it 00038 being set for each detector in the frame in GEO_locate_one_scan. 00039 00040 !Input Parameters: 00041 scan_number The scan number 00042 sample_number The sample number within the scan 00043 num_detectors Number of detectors for the sample 00044 u_inst View vectors for each detector in the 00045 instrument coordinate system (meters). 00046 ecr_sc_sample_position The spacecraft position. 00047 ecr_sc_velocity The spacecraft velocity. 00048 T_inst2ecr Matrix for conversion from instrument to ECR 00049 coordinates. 00050 00051 !Output Parameters: 00052 ecr_sample_position Ground positions in ECR coordinates. 00053 ellip_sample_position Ground positions without terrain correction. 00054 sample_flags Error flags for each sample. 00055 sample_view_vec View vector from satellite to ground position. 00056 00057 Return Values: 00058 FAIL If arguments are invalid or PGS_CSC_GetEarthFigure() 00059 fails. 00060 SUCCESS Otherwise 00061 00062 Externally Defined: 00063 DETECTORS_QKM "GEO_geo.h" 00064 FAIL "GEO_basic.h" 00065 INVALID_INPUT_DATA "GEO_geo.h" 00066 MAX_SCAN_NUMBER "GEO_geo.h" 00067 MODIS_E_BAD_INPUT_ARG "PGS_MODIS_35251.h" 00068 MODIS_E_BAD_VIEW_VEC "PGS_MODIS_35251.h" 00069 MODIS_E_GEO "PGS_MODIS_35251.h" 00070 MODIS_W_SCAN_OFF "PGS_MODIS_35251.h" 00071 NO_ELLIPSE_INTERSECT "GEO_geo.h" 00072 MAX_PADDED "GEO_geo.h" 00073 PGSd_GEO_ERROR_VALUE "PGS_TD.h" 00074 PGS_S_SUCCESS "PGS_SMF.h" 00075 SUCCESS "GEO_basic.h" 00076 00077 Called by: GEO_earth_location 00078 00079 Routines called: 00080 PGS_CSC_ECRtoGEO "PGS_CSC.h" 00081 PGS_CSC_GetEarthFigure "PGS_CSC.h" 00082 modsmf "smfio.h" 00083 00084 !Revision History: 00085 $Log: GEO_ellip_position.c,v $ 00086 Revision 6.4 2010/06/29 18:45:31 kuyper 00087 Replaced floating point equality test with an inequality. 00088 00089 Revision 6.3 2009/05/31 19:02:17 kuyper 00090 Corrected to initialize p_prime_sq, up_prod, u_prime_sq where needed. 00091 Improved error messages. 00092 00093 Revision 6.2 2009/05/15 17:00:18 xgeng 00094 Minor correction on error messages. 00095 00096 Xu Geng (xu.geng@saic.com) 00097 Revision 6.1 2009/05/06 19:34:43 xgeng 00098 Code updated to match the PDL change made in revision 6.1, 6.2 and 6.3. 00099 00100 Revision 4.4 2003/09/03 17:17:04 vlin 00101 corrected a typo. 00102 00103 Revision 4.3 2003/08/28 15:56:51 kuyper 00104 Corrected prolog. 00105 00106 Revision 4.2 2003/08/22 14:24:21 kuyper 00107 Removed T_inst2ecr and sc_attitude globals. 00108 00109 Revision 4.1 2003/08/08 19:36:28 vlin 00110 Calls to GEO_interp_ephemeris_attitude() and GEO_get_T_inst2ecr() are moved 00111 up to a higher level. As a result, ecr_sc_position, ecr_sc_velocity, 00112 scan_start_time, and T_inst2ecr are now inputs rather than outputs. 00113 Removed intitialization of sc_attitude. 00114 Shifted messages to MODIS_E_GEO, where applicable. 00115 00116 Revision 2.4 1999/01/29 20:14:01 kuyper 00117 Removed max_extrap from GEO_interp_ephemeris() argument list. 00118 00119 4/25/95 00120 Ruiming Chen 00121 Finished coding. 00122 00123 !Team-unique Header: 00124 This software is developed by the MODIS Science Data Support 00125 Team for the National Aeronautics and Space Administration, 00126 Goddard Space Flight Center, under contract NAS5-32373. 00127 00128 Requirements: 00129 PR03-F-3.2.2-1 00130 PR03-F-3.2.3-1 00131 PR03-F-3.2.3-2 00132 PR03-F-3.2.4-1 00133 PR03-I-1 00134 PR03-I-2 00135 00136 !END************************************************************************* 00137 */ 00138 00139 { 00140 static PGSt_double equat=-1.0; /* equatorial radius of the earth */ 00141 static PGSt_double polar=-1.0; /* polar radius of the earth. */ 00142 PGSt_double p_prime[3]; /* rescaled sc position */ 00143 double p_prime_sq = 0.0; /* square of p_prime_length */ 00144 PGSt_SMF_status PGS_error_code; /* SDP toolkit function return value */ 00145 int det; /* detector index */ 00146 int ecr; /* iteration parameter */ 00147 char msgbuf[128] = ""; 00148 static char filefunc[] = __FILE__ ", GEO_ellip_position"; 00149 00150 /**************************************************************************** 00151 calculate from here 00152 ****************************************************************************/ 00153 00154 if( scan_number < 0 || scan_number >= MAX_SCAN_NUMBER || 00155 num_detectors < 0 || num_detectors > DETECTORS_QKM || 00156 sample_number < 0 || sample_number >= MAX_PADDED || 00157 u_inst == NULL || ecr_sc_sample_position == NULL || 00158 ecr_sc_velocity == NULL || T_inst2ecr == NULL || 00159 ecr_sample_position == NULL || ellip_sample_position == NULL || 00160 sample_flags == NULL || sample_view_vec == NULL ) 00161 { 00162 sprintf(msgbuf, "scan:%d detectors:%d sample:%d u_inst:%p " 00163 "ecr_sc_sample_position:%p ecr_sc_velocity:%p " 00164 "T_inst2ecr:%p ecr_sample_position:%p " 00165 "ellip_sample_position:%p sample_flags:%p " 00166 "sample_view_vec:%p", 00167 scan_number, num_detectors, sample_number, (void *)u_inst, 00168 (void *)ecr_sc_sample_position, (void *)ecr_sc_velocity, 00169 (void *)T_inst2ecr, (void *)ecr_sample_position, 00170 (void *)ellip_sample_position, (void *)sample_flags, 00171 (void *)sample_view_vec 00172 ); 00173 00174 modsmf(MODIS_E_BAD_INPUT_ARG, msgbuf, filefunc); 00175 return FAIL; 00176 } 00177 00178 if (ecr_sc_sample_position[sample_number][0] >= PGSd_GEO_ERROR_VALUE) { 00179 for (det = 0; det < num_detectors; det++) 00180 sample_flags[det][sample_number] = INVALID_INPUT_DATA; 00181 return SUCCESS; 00182 } 00183 00184 for (det=0; det<num_detectors ; det++) 00185 sample_flags[det][sample_number] = 0; 00186 00187 if (equat < 0.0) { 00188 /* the equatorial and polar earth radii have not been set */ 00189 PGS_error_code = PGS_CSC_GetEarthFigure(EARTH_MODEL, &equat, &polar); 00190 if (PGS_error_code!=PGS_S_SUCCESS) 00191 { 00192 modsmf(MODIS_E_GEO, "PGS_CSC_GetEarthFigure(\"" EARTH_MODEL "\")", 00193 filefunc); 00194 return FAIL; 00195 } 00196 } 00197 00198 /* Rescale orbit position vector */ 00199 p_prime[0] = ecr_sc_sample_position[sample_number][0] / equat; 00200 p_prime[1] = ecr_sc_sample_position[sample_number][1] / equat; 00201 p_prime[2] = ecr_sc_sample_position[sample_number][2] / polar; 00202 00203 for (ecr = 0; ecr < 3; ecr ++) 00204 p_prime_sq += p_prime[ecr]*p_prime[ecr]; 00205 00206 /* Loop through detectors in sample */ 00207 for (det = 0; det < num_detectors; det++) 00208 { 00209 PGSt_double u_prime[3]; /* rescaled view vec */ 00210 double up_prod_sq; /* square of up_prod */ 00211 double sqrt_arg; /* argument of square root in the calculation of d */ 00212 double d; /* scaling factor of the viewing vector */ 00213 double up_prod = 0.0; /* product between u_prime and p_prime */ 00214 double u_prime_sq = 0.0; /* square of u_prime_length */ 00215 00216 /* Transform the viewing vector to ECR */ 00217 for (ecr = 0; ecr < 3; ecr ++) { 00218 int inst; /* iteration parameter */ 00219 sample_view_vec[det][sample_number][ecr] = 0.0; 00220 00221 for (inst = 0; inst < 3; inst ++) { 00222 sample_view_vec[det][sample_number][ecr] += 00223 T_inst2ecr[sample_number][ecr][inst]*u_inst[det][inst]; 00224 } 00225 } 00226 00227 /* Rescale viewing vector */ 00228 u_prime[0] = sample_view_vec[det][sample_number][0] / equat; 00229 u_prime[1] = sample_view_vec[det][sample_number][1] / equat; 00230 u_prime[2] = sample_view_vec[det][sample_number][2] / polar; 00231 00232 /* Compute dot product of u_prime and p_prime 00233 and squared magnitudes of u_prime and p_prime */ 00234 for (ecr = 0; ecr < 3; ecr ++) { 00235 up_prod += u_prime[ecr]*p_prime[ecr]; 00236 u_prime_sq += u_prime[ecr]*u_prime[ecr]; 00237 } 00238 00239 up_prod_sq = up_prod * up_prod; 00240 00241 /* Check for negative under square root (scan off of Earth) */ 00242 sqrt_arg = up_prod_sq - u_prime_sq*(p_prime_sq - 1.0); 00243 if (sqrt_arg < 0.0 || up_prod>=0.0) { /* No real ellipsoid intersection*/ 00244 /* call SDP function to report event */ 00245 modsmf(MODIS_W_SCAN_OFF, "", filefunc); 00246 sample_flags[det][sample_number] |= NO_ELLIPSE_INTERSECT; 00247 } 00248 else if ( (-up_prod - sqrt(sqrt_arg))/DBL_MAX >= u_prime_sq) 00249 { /* Calculation of d will produce overflow. */ 00250 modsmf(MODIS_E_BAD_VIEW_VEC, "", filefunc); 00251 sample_flags[det][sample_number] |= INVALID_INPUT_DATA; 00252 } 00253 else 00254 { 00255 PGSt_double latitude; /* latitude of the ellipsoid position */ 00256 PGSt_double longitude; /* longitude of the ellipsoid position */ 00257 PGSt_double altitude; /* altitude of the ellipsoid position */ 00258 00259 d = (-up_prod - sqrt(sqrt_arg)) / u_prime_sq; 00260 00261 /* calculate the geocentric position vector */ 00262 for (ecr = 0; ecr < 3; ++ecr) { 00263 sample_view_vec[det][sample_number][ecr] *= d; 00264 ecr_sample_position[det][sample_number][ecr] = 00265 ecr_sc_sample_position[sample_number][ecr] + 00266 sample_view_vec[det][sample_number][ecr]; 00267 } 00268 00269 /* calculate the geodetic lat and lon */ 00270 00271 if (PGS_CSC_ECRtoGEO( ecr_sample_position[det][sample_number], 00272 EARTH_MODEL, &longitude, &latitude, &altitude) != PGS_S_SUCCESS) 00273 { 00274 modsmf(MODIS_E_GEO, "PGS_CSC_ECRtoGEO(\"" EARTH_MODEL "\")", filefunc); 00275 sample_flags[det][sample_number] |= INVALID_INPUT_DATA; 00276 continue; 00277 } 00278 00279 /* store latitude and longitude in global variable */ 00280 ellip_sample_position[det][sample_number][LAT_IDX] = latitude; 00281 ellip_sample_position[det][sample_number][LON_IDX] = longitude; 00282 ellip_sample_position[det][sample_number][HT_IDX] = 0.0; 00283 } 00284 } 00285 00286 return SUCCESS; 00287 } 00288
1.7.6.1