/*********************************************************************************************************************************** asa_im_segment.c Version 1.5 Program to extract a smaller frame of ENVISAT ASAR IMAGE MODE Level 0 data from a long strip of Levevl 0 data set compiled on a Linux with command gcc -O2 -lm asa_im_segment.c -o asa_im_segment User will be prompted to provide information about the start and end of smaller frame file in terms of any combination of start line number / start time of day in seconds / start time in seconds wrt start of long-multi-frame file / bounding latitude && end line number/ number of lines in file/ end time of day in secs / end time in secs wrt start of long-multi-frame file / bounding latitude This code was written by taking advantage of an earlier code written to decode ASAR IM Mode Level 0 data decoder. ******* Previous code written by **************** asa_im_decoder.c : v1.0 Feb/Mar 2004, Sean M. Buckley v1.05, Mar 25, 2004, Sean M. Buckley v1.1, 17 Feb 2005, Zhenhong Li: v1.1.1, 8 Nov 2005, Vikas Gudipati, now runs correctly on 64-bit compilers. **************************************** This code written by asa_im_segment.c v0.0 June 4 2008 Vikas Gudipati and Dochul Yang. (Basic version, does not include extra padding.) v1.1 June 6 2008 Vikas Gudipati and Dochul Yang. (Converted to command line input mode. Now includes latitude based selection with padding.) v1.2 June 18 2008 Vikas Gudipati. (Corrected mistake in conversion from ECEF to ECI and vice versa) v1.3 June 21 2008 Vikas Gudipati. (Adjusted padding for large errors in 2body Orbit propagation. Padding 0.1 degrees at each end of file) (Added command prompt input parameter to specify average terrain elevation wrt WGS84) v1.4 June 25 2008 Vikas Gudipati and Dochul Yang. (Orbit propagation with J2 effect via numerical integration) (Padding due to errors in orbit propagation and zero squint assumption changed to 0.05 degrees. about 5.4 KM) v1.5 Sept 01 2008 Vikas Gudipati and Dochul Yang. (Incorporated option of using interpolated precise orbit file.) (Added flag to disable half aperture padding at both ends of raw data) ***********************************************************************************************************************************/ #include #include #include #include #include #define muE 3.9860044e+14 /* Earth's gravitational parameter*/ #define SIZE6 6 /* size of state vector and orbit element */ #define SIZE3 3 /* size of 3-D vector */ #define PI 3.14159265358979 #define TWO_PI 6.283185307179586 #define SMAJORA 6378137.0 #define SMINORA 6356752.3142 #define RTD 57.29577951 #define DTR 0.017453292 #define EROTRAD 7.292115392e-5 #define J2 0.0010826269 #define Re 6378137.0 #define STEP 0.01 #define NSTATE 8 #define NlineHeader 51 #define NlineRecord 1589 #define MaxRecLength 130 /** * The EPR_DataTypeId enumeration lists all possible data * types for field elements in ENVISAT dataset records. */ // added by Z. Li on 16/02/2005 // extracted from ESA BEAM epr_api package enum EPR_DataTypeId { /** The ID for unknown types. */ e_tid_unknown = 0, /** An array of unsigned 8-bit integers, C type is uchar* */ e_tid_uchar = 1, /** An array of signed 8-bit integers, C type is char* */ e_tid_char = 2, /** An array of unsigned 16-bit integers, C type is ushort* */ e_tid_ushort = 3, /** An array of signed 16-bit integers, C type is short* */ e_tid_short = 4, /** An array of unsigned 32-bit integers, C type is ulong* */ e_tid_ulong = 5, /** An array of signed 32-bit integers, C type is long* */ e_tid_long = 6, /** An array of 32-bit floating point numbers, C type is float* */ e_tid_float = 7, /** An array of 64-bit floating point numbers, C type is double* */ e_tid_double = 8, /** A zero-terminated ASCII string, C type is char* */ e_tid_string = 11, /** An array of unsigned character, C type is uchar* */ e_tid_spare = 13, /** A time (MJD) structure, C type is EPR_Time */ e_tid_time = 21 }; /* structures */ struct mphStruct { char product[ 62 ]; char procStage[ 1 ]; char refDoc[ 23 ]; char spare1[ 40 ]; char acquisitionStation[ 20 ]; char procCenter[ 6 ]; char procTime[ 27 ]; char softwareVer[ 14 ]; char spare2[ 40 ]; char sensingStart[ 27 ]; char sensingStop[ 27 ]; char spare3[ 40 ]; char phase[ 1 ]; int cycle; int relOrbit; int absOrbit; char stateVectorTime[ 27 ]; double deltaUt1; double xPosition; double yPosition; double zPosition; double xVelocity; double yVelocity; double zVelocity; char vectorSource[ 2 ]; char spare4[ 40 ]; char utcSbtTime[ 27 ]; unsigned int satBinaryTime; unsigned int clockStep; char spare5[ 32 ]; char leapUtc[ 27 ]; int leapSign; int leapErr; char spare6[ 40]; int productErr; int totSize; int sphSize; int numDsd; int dsdSize; int numDataSets; char spare7[ 40 ]; }; struct dsdStruct { char dsName[ 28 ]; char dsType[ 1 ]; char filename[ 62 ]; int dsOffset; int dsSize; int numDsr; int dsrSize; }; struct sphStruct { char sphDescriptor[ 28 ]; double startLat; double startLon; double stopLat; double stopLon; double satTrack; char spare1[ 50 ]; int ispErrorsSignificant; int missingIspsSignificant; int ispDiscardedSignificant; int rsSignificant; char spare2[ 50 ]; int numErrorIsps; double errorIspsThresh; int numMissingIsps; double missingIspsThresh; int numDiscardedIsps; double discardedIspsThresh; int numRsIsps; double rsThresh; char spare3[ 100 ]; char txRxPolar[ 5 ]; char swath[ 3 ]; char spare4[ 41 ]; struct dsdStruct dsd[ 4 ]; }; struct sphAuxStruct { char sphDescriptor[ 28 ]; char spare1[ 51 ]; struct dsdStruct dsd[ 1 ]; }; struct dsrTimeStruct { int days; int seconds; int microseconds; }; struct calPulseStruct { float nomAmplitude[ 32 ]; float nomPhase[ 32 ]; }; struct nomPulseStruct { float pulseAmpCoeff[ 4 ]; float pulsePhsCoeff[ 4 ]; float pulseDuration; }; struct dataConfigStruct { char echoCompMethod[ 4 ]; char echoCompRatio[ 3 ]; char echoResampFlag[ 1 ]; char initCalCompMethod[ 4 ]; char initCalCompRatio[ 3 ]; char initCalResampFlag[ 1 ]; char perCalCompMethod[ 4 ]; char perCalCompRatio[ 3 ]; char perCalResampFlag[ 1 ]; char noiseCompMethod[ 4 ]; char noiseCompRatio[ 3 ]; char noiseResampFlag[ 1 ]; }; struct swathConfigStruct { unsigned short numSampWindowsEcho[ 7 ]; unsigned short numSampWindowsInitCal[ 7 ]; unsigned short numSampWindowsPerCal[ 7 ]; unsigned short numSampWindowsNoise[ 7 ]; float resampleFactor[ 7 ]; }; struct swathIdStruct { unsigned short swathNum[ 7 ]; unsigned short beamSetNum[ 7 ]; }; struct timelineStruct { unsigned short swathNums[ 7 ]; unsigned short mValues[ 7 ]; unsigned short rValues[ 7 ]; unsigned short gValues[ 7 ]; }; /* problems begin with field 132 - check the double statement */ struct testStruct { float operatingTemp; float rxGainDroopCoeffSmb[ 16 ]; /* this needs to be converted to a double array of eight elements */ //double rxGainDroopCoeffSmb[ 8 ]; /* Something wrong here, why?*/ }; struct insGadsStruct { /* see pages 455-477 for the 142 fields associated with this gads - got length of 121712 bytes */ struct dsrTimeStruct dsrTime; unsigned int dsrLength; float radarFrequency; float sampRate; float offsetFreq; struct calPulseStruct calPulseIm0TxH1; struct calPulseStruct calPulseIm0TxV1; struct calPulseStruct calPulseIm0TxH1a; struct calPulseStruct calPulseIm0TxV1a; struct calPulseStruct calPulseIm0RxH2; struct calPulseStruct calPulseIm0RxV2; struct calPulseStruct calPulseIm0H3; struct calPulseStruct calPulseIm0V3; struct calPulseStruct calPulseImTxH1[ 7 ]; struct calPulseStruct calPulseImTxV1[ 7 ]; struct calPulseStruct calPulseImTxH1a[ 7 ]; struct calPulseStruct calPulseImTxV1a[ 7 ]; struct calPulseStruct calPulseImRxH2[ 7 ]; struct calPulseStruct calPulseImRxV2[ 7 ]; struct calPulseStruct calPulseImH3[ 7 ]; struct calPulseStruct calPulseImV3[ 7 ]; struct calPulseStruct calPulseApTxH1[ 7 ]; struct calPulseStruct calPulseApTxV1[ 7 ]; struct calPulseStruct calPulseApTxH1a[ 7 ]; struct calPulseStruct calPulseApTxV1a[ 7 ]; struct calPulseStruct calPulseApRxH2[ 7 ]; struct calPulseStruct calPulseApRxV2[ 7 ]; struct calPulseStruct calPulseApH3[ 7 ]; struct calPulseStruct calPulseApV3[ 7 ]; struct calPulseStruct calPulseWvTxH1[ 7 ]; struct calPulseStruct calPulseWvTxV1[ 7 ]; struct calPulseStruct calPulseWvTxH1a[ 7 ]; struct calPulseStruct calPulseWvTxV1a[ 7 ]; struct calPulseStruct calPulseWvRxH2[ 7 ]; struct calPulseStruct calPulseWvRxV2[ 7 ]; struct calPulseStruct calPulseWvH3[ 7 ]; struct calPulseStruct calPulseWvV3[ 7 ]; struct calPulseStruct calPulseWsTxH1[ 5 ]; struct calPulseStruct calPulseWsTxV1[ 5 ]; struct calPulseStruct calPulseWsTxH1a[ 5 ]; struct calPulseStruct calPulseWsTxV1a[ 5 ]; struct calPulseStruct calPulseWsRxH2[ 5 ]; struct calPulseStruct calPulseWsRxV2[ 5 ]; struct calPulseStruct calPulseWsH3[ 5 ]; struct calPulseStruct calPulseWsV3[ 5 ]; struct calPulseStruct calPulseGmTxH1[ 5 ]; struct calPulseStruct calPulseGmTxV1[ 5 ]; struct calPulseStruct calPulseGmTxH1a[ 5 ]; struct calPulseStruct calPulseGmTxV1a[ 5 ]; struct calPulseStruct calPulseGmRxH2[ 5 ]; struct calPulseStruct calPulseGmRxV2[ 5 ]; struct calPulseStruct calPulseGmH3[ 5 ]; struct calPulseStruct calPulseGmV3[ 5 ]; struct nomPulseStruct nomPulseIm[ 7 ]; struct nomPulseStruct nomPulseAp[ 7 ]; struct nomPulseStruct nomPulseWv[ 7 ]; struct nomPulseStruct nomPulseWs[ 5 ]; struct nomPulseStruct nomPulseGm[ 5 ]; float azPatternIs1[ 101 ]; float azPatternIs2[ 101 ]; float azPatternIs3Ss2[ 101 ]; float azPatternIs4Ss3[ 101 ]; float azPatternIs5Ss4[ 101 ]; float azPatternIs6Ss5[ 101 ]; float azPatternIs7[ 101 ]; float azPatternSs1[ 101 ]; float rangeGateBias; float rangeGateBiasGm; float adcLutI[ 255 ]; float adcLutQ[ 255 ]; char spare1[ 648 ]; float full8LutI[ 256 ]; float full8LutQ[ 256 ]; float fbaq4LutI[ 4096 ]; float fbaq3LutI[ 2048 ]; float fbaq2LutI[ 1024 ]; float fbaq4LutQ[ 4096 ]; float fbaq3LutQ[ 2048 ]; float fbaq2LutQ[ 1024 ]; float fbaq4NoAdc[ 4096 ]; float fbaq3NoAdc[ 2048 ]; float fbaq2NoAdc[ 1024 ]; float smLutI[ 16 ]; float smLutQ[ 16 ]; struct dataConfigStruct dataConfigIm; struct dataConfigStruct dataConfigAp; struct dataConfigStruct dataConfigWs; struct dataConfigStruct dataConfigGm; struct dataConfigStruct dataConfigWv; struct swathConfigStruct swathConfigIm; struct swathConfigStruct swathConfigAp; struct swathConfigStruct swathConfigWs; struct swathConfigStruct swathConfigGm; struct swathConfigStruct swathConfigWv; unsigned short perCalWindowsEc; unsigned short perCalWindowsMs; struct swathIdStruct swathIdIm; struct swathIdStruct swathIdAp; struct swathIdStruct swathIdWs; struct swathIdStruct swathIdGm; struct swathIdStruct swathIdWv; unsigned short initCalBeamSetWv; unsigned short beamSetEc; unsigned short beamSetMs; unsigned short calSeq[ 32 ]; struct timelineStruct timelineIm; struct timelineStruct timelineAp; struct timelineStruct timelineWs; struct timelineStruct timelineGm; struct timelineStruct timelineWv; unsigned short mEc; char spare2[ 44 ]; float refElevAngleIs1; float refElevAngleIs2; float refElevAngleIs3Ss2; float refElevAngleIs4Ss3; float refElevAngleIs5Ss4; float refElevAngleIs6Ss5; float refElevAngleIs7; float refElevAngleSs1; char spare3[ 64 ]; float calLoopRefIs1[ 128 ]; float calLoopRefIs2[ 128 ]; float calLoopRefIs3Ss2[ 128 ]; float calLoopRefIs4Ss3[ 128 ]; float calLoopRefIs5Ss4[ 128 ]; float calLoopRefIs6Ss5[ 128 ]; float calLoopRefIs7[ 128 ]; float calLoopRefSs1[ 128 ]; char spare4[ 5120 ]; struct testStruct im; struct testStruct ap; struct testStruct ws; struct testStruct gm; struct testStruct wv; float swstCalP2; char spare5[ 72 ]; }; // added by Z. Li on 16/02/2005 typedef enum EPR_DataTypeId EPR_EDataTypeId; typedef int boolean; typedef unsigned char uchar; typedef unsigned short ushort; typedef unsigned int uint; typedef unsigned long ulong; /* function prototypes */ struct mphStruct readMph( const char *mphPtr, const int printMphIfZero ); struct sphStruct readSph( const char *sphPtr, const int printSphIfZero, const struct mphStruct mph ); struct sphAuxStruct readSphAux( const char *sphPtr, const int printSphIfZero, const struct mphStruct mph ); struct insGadsStruct readInsGads( const char *gadsPtr, const int printInsGadsIfZero ); void printInsGads( const struct insGadsStruct ); // added by Z. Li on 16/02/2005 int is_bigendian(); void byte_swap_short(short *buffer, uint number_of_swaps); void byte_swap_ushort(ushort* buffer, uint number_of_swaps); void byte_swap_long(long *buffer, uint number_of_swaps); void byte_swap_ulong(ulong* buffer, uint number_of_swaps); void byte_swap_float(float* buffer, uint number_of_swaps); /* new byte_swap_int type added below*/ void byte_swap_int(int *buffer, uint number_of_swaps); void byte_swap_uint(uint *buffer, uint number_of_swaps); /* new byte_swap_uint type added above*/ void swap_endian_order(EPR_EDataTypeId data_type_id, void* elems, uint num_elems); void byte_swap_InsGads( struct insGadsStruct* InsGads ); double norm( double a[]); int max_int(int a, int b); int min_int(int a, int b); double dot(double a[], double b[]); void cross(double a[], double b[], double c[]); void svmul(double sm, double a[], double b[]); void svdiv(double sm, double a[], double b[]); void vadd(double a[], double b[], double c[]); void vsub(double a[], double b[], double c[]); void matvmul(double a[][3], double b[], double c[]); void rot1(double ang, double c[][3]); void rot2(double ang, double c[][3]); void rot3(double ang, double c[][3]); double MJD2GST(double MJD, double dt); void ECEF2ECI(int yr, int mo, int day, int hr, int min, double sec, double dt, double Secef[], double Seci[]); void ECI2ECEF(int yr, int mo, int day, int hr, int min, double sec, double dt, double Seci[], double Secef[]); double date2MJD(int yr, int mo, int day, int hr, int min, double sec); void orbitPropogation2body_J2Cartesian( int year, int mon, int day, int hour, int min, double sec, double delT, double sVEC[], double sVECf[]); void RK4(double RVsat0[] , double delt , double RVsatf[] ); void derivJ2(double qVEC[], double qdotVEC[]); void xyz2LatLonAlt(double rvec[] , double a, double b, double *phi, double *lam, double *h); void geolocate_radarLOS_ellipsoidWGS84(double rv[], double el, double rho, int DIR_FLAG, double ref_alt, double *lat, double *lon); int locTable(double xa[], int n, double xp); double polyint(double xa[], double ya[], int n, double xp); void orbitInterpolationPoly(int inyear, int inmonth, int inday, int inhour, int inmin, double inseconds, double delT, char orbitFileName[], double sVECout[]); /**********************************************************************************************************************************/ int main( int argc, char *argv[] ) { /* variable definitions */ FILE *imFilePtr; FILE *outFilePtr; FILE *insFilePtr; FILE *lsPtr; char imFileName[ 200 ]; char outFileName[ 200 ]; char insFileName[ 200 ]; char *mphPtr2; char *mphPtr; char *sphPtr; char *sphPtr2; char *gadsPtr; unsigned char mdsrBlockId[ 200 ]; unsigned char mdsrCheck[ 63 ]; unsigned char spare; unsigned char antennaBeamSetNumber; unsigned char mdsrLineChar[ 20000 ]; int printImMphIfZero = 1; int printImSphIfZero = 1; int printImMdsrIfZero = 1; int printInsMphIfZero = 1; int printInsSphIfZero = 1; int printInsGadsIfZero = 1; int mphSize = 1247; /* fixed size */ int outSamples = 0; int bytesRead = 0; int bytesWritten = 0; int bytesToBeWritten = 0; int i; int ii; int j; int k; int m; int n; int mdsrDsrTimeDays; int mdsrDsrTimeSeconds; int mdsrDsrTimeMicroseconds; int mdsrGsrtTimeDays; int mdsrGsrtTimeSeconds; int mdsrGsrtTimeMicroseconds; double DsrStartTime; double DsrEndTime; double dsr_rho_near_firstline, dsr_rho_far_firstline; double dsr_rho_near_lastline, dsr_rho_far_lastline; int DsrStartDay, DsrEndDay; int DsrStartDayMJD, DsrEndDayMJD; int DsrStartFlag = 0; int echoFlag; int fline; int lline; unsigned short mdsrIspLength; unsigned short mdsrIspLengthSwapped; unsigned short mdsrCrcErrs; unsigned short mdsrRsErrs; unsigned short mdsrSpare1; unsigned short mdsrPacketIdentification; unsigned short mdsrPacketSequenceControl; unsigned short mdsrPacketLength; unsigned short mdsrPacketDataHeader[ 15 ]; unsigned short priCodeword; unsigned short windowStartTimeCodeword; unsigned short windowLengthCodeword; unsigned short dataFieldHeaderLength; unsigned short cyclePacketCount; float mdsrLine[ 20000 ]; double TxPulseLength; double c = 299792458.; double pri; double windowStartTime; double windowLength; struct mphStruct mph; struct mphStruct mph2; struct mphStruct mphIns; struct sphStruct sph; struct sphStruct sph2; struct sphAuxStruct sphIns; struct insGadsStruct insGads; double ref_alt = 0.0; /* reference height of land surface above ellipsoid set to 0.0 */ double ref_elevation_angle; double rho_near_firstline, rho_far_firstline; double rho_near_lastline, rho_far_lastline; int DayMJD2000TimeFirstMDSRLine; int DayMJD2000TimeLastMDSRLine; double SecondsMJD2000TimeFirstMDSRLine; double SecondsMJD2000TimeLastMDSRLine; double MicroSecondsMJD2000TimeFirstMDSRLine; double MicroSecondsMJD2000TimeLastMDSRLine; int DayMJDTimeFirstMDSRLine; int DayMJDTimeLastMDSRLine; double TimeFirstMDSRLine; double TimeLastMDSRLine; double MaxTimeDuration; double sVEC[6]; double sVECn[6]; double sVECfirst[6]; double sVEClast[6]; double svseconds; double svMJD; double delT1; double delT2; double clat[4], clon[4]; double olat, olon; double nlat, nlon; double ave_dLat_daz; double ave_dLat_dr; double vsatv[3]; double vsat; double OrbitPropErrTime; // double OrbitPropErrDistance = 0.05*3600 *30; /* Error of 0.05 degrees in orbit propagation with 1 arc sec ~ 30 meters */ double OrbitPropErrDistance; /* Error of 0.05 degrees in orbit propagation with 1 arc sec ~ 30 meters */ double deltanLat; double start_line_TimeDaySeconds; double start_line_TimeSeconds; double start_lat; double deltaLat; double delT; double end_line_TimeDaySeconds; double end_line_TimeSeconds; double timeduration; double end_lat; double heading; double sseconds; double eseconds; double boundLat1, boundLat2; int OrbitPropErrPadding; int EdgeTargetFullCompressionPadding = 550; /* Half aperture at far range padded to facilitate full compression of frame edge targets */ int padding; int estStartLine; int StartLineErr; int estLastLine; int EndLineErr; int oflag = 0; int sflag = 0; int eflag = 0; int DIR_FLAG = -1 ; /* assumed right looking radar */ int svyear, svmonth, svday, svhour, svmin; int MaxNumLines; int svMJDday; int start_line, end_line, num_lines; int temp_line; int err_flag1 = 0, err_flag2 = 0; int sday, syear, smonth, shour, smin; int eday, eyear, emonth, ehour, emin; int ssecint, esecint; int totfilesize; int duration; int curser; int orbFileFlag = 0; int exitflag = 0; int orbdate1, orbdate2; int sensestartyyyymmdd, sensestopyyyymmdd; int status; int HApadflag = 1; int is_littlendian; char svmonthstr[4]; char swathis[3]; char smonthstr[4]; char emonthstr[4]; char shh[3]; char ehh[3]; char smm[3]; char emm[3]; char sss[3]; char ess[3]; char ssec[10]; char esec[10]; char ds_size[21]; char num_dsr[11]; char shhmmss[7]; char ehhmmss[7]; char secduration[9]; char totsize[22]; char orbFileName[63]; char orbFileDir[200]; char dummyOrbFileName[63]; char dummyphrase1[13]; char dummyphrase2[18]; char dummyphrase3[9]; char dummyphrase4[8]; char listcmd[300]; /* usage note */ printf( "\n*** asa_im_segment v1.5 by Vikas Gudipati and Dochul Yang based on asar decoder by smb ***\n\n" ); if ( (argc-1) < 8 ) { printf( "Extracts requested data segment from Envisat ASAR Image Mode Level 0 multi frame data .\n\n" ); printf( "Usage: asa_im_segment (s_line) (s_TimeofDay) (s_TimeofFile) (boundLat1) [e_line] [num_lines] [e_TimeofDay] [e_TimeofFile] [duration] [boundLat2] {ref_alt} {orbDir} {HApad}\n\n" ); printf("Quantities shown within < > are all required.\n"); printf("Any one of the quantities shown within ( ) is required to indicate start of output data file. If more than one is given, the latest option is used.\n"); printf("Any one of the quantities shown within [ ] is required to indicate end of output data file. If more than one is given, the latest option is used.\n\n"); printf( " input Image mode file (multi-frame Level 0 file)\n" ); printf( " input auxiliary instrument characterization file (ASA_INS file). See notes below for download location.\n" ); printf( " output Image mode file (single frame Level 0 file)\n" ); printf( " (s_line) starting raw data line number with respect to start of multi-frame file. (Enter integer value or '-')\n"); printf( " (s_TimeofDay) starting time of the frame in seconds from start of the day. (Enter float value or '-')\n"); printf( " (s_TimeofFile) starting time of the frame in seconds from start of multi-frame file. (Enter float value or '-')\n"); printf( " (boundLat1) first bounding latitude in decimal degrees for frame file. (Enter float value or '-')\n"); printf( " [e_line] last raw data line number with respect to start of multi-frame file. (Enter integer value or '-')\n"); printf( " [num_lines] number of raw data lines to be included in the frame file. (Enter integer value or '-')\n"); printf( " [e_TimeofDay] time of the last line to be included in the frame in seconds from start of the day. (Enter float value or '-')\n"); printf( " [e_TimeofFile] time of the last line to be included in the frame in seconds from start of multi-frame file. (Enter float value or '-')\n"); printf( " [duration] duration of data in seconds. (Enter float value or '-')\n"); printf( " [boundLat2] second bounding latitude in decimal degrees for frame file.(Enter float value or '-')\n"); printf( " {ref_alt} Reference (average) altitude of terrain surface over ellipsoid in meters.(Default: 0.0 m, Enter float value or '-')\n"); printf( " {orbDir} Full path for precise orbit directory, Enter character string value /path/to/orbitfile/directory/ or '-'\n"); printf( " {HApad} Flag to disable Half Aperture Padding at data edges when latitude bounds are used. (Enter integer value: 0 to disable, 1 to enable, Default: 1) or '-'\n\n"); printf( "Notes:\n" ); printf( "out_im_sf is a ASAR IMAGE MODE Level 0 file with headers (Updated MPH, SPH, DSDs but in same format as those in multi-frame file)\n\n" ); printf( "Data padding in case of latitude bounds based frame selection:\n" ); printf( "Extra data lines are padded to the user requested latitude bounding based estimated data set, for 2 reasons:\n" ); printf( "1) If precise orbits are not used, satellite orbit is propagated using a 4th order Runga-Kutta integrator with J2 effect. Integration error is estimated to be about 5 KM.\n" ); printf( "2) To obtain full resolution (in SLC) of targets lying near the edges of latitude bounding lines, lines amounting to half-aperture (~550 lines) are padded at start and end.\n\n" ); printf( "If only one latitude bound, and either start line or end line is given, then:\n" ); printf( "1)If start line is 'ahead' in time with respect to latitude bound, then start line will be made the end line.\n" ); printf( "2)If end line is 'behind' in time with respect to latitude bound, then end line will be made the start line.\n\n" ); printf( "No line correction is perfromed to the original dataset.\n" ); printf( "Calibration/noise lines are copied as they were in the origianl data.\n" ); printf( "In the algorithm to geolocate radar look vector, starting range is computed as (rank*pri+windowStartTime)*c/2 where rank is the number of pri between transmitted pulse and return echo\n\n" ); printf( "auxiliary data files can be found at http://envisat.esa.int/services/auxiliary_data/asar/current/\n\n"); printf( "Envisat product specifications can be found at http://envisat.esa.int/support-docs/productspecs/\n\n"); printf( "Envisat ASAR Product Handbook, Issue 1.1, 1 December 2002 can be found at http://envisat.esa.int/dataproducts/asar/CNTR6-3-6.htm#eph.asar.asardf.0pASA_IM__0P\n\n" ); return 0; } /* read in command-line arguments */ sscanf( argv[2], "%s", insFileName ); sscanf( argv[3], "%s", outFileName ); if( (strcmp(argv[4], "-") == 0) && (strcmp(argv[5], "-") == 0) && (strcmp(argv[6], "-") == 0) && (strcmp(argv[7], "-") == 0) ){ printf("\n ERROR: Enter one of the parameters (s_line, s_TimeofDay, s_TimeofFile, boundLat1) to indicate start of frame. \n"); exit(-1); } if(strcmp(argv[4], "-") != 0){ sscanf(argv[4], "%d", &start_line); sflag = 1; } if(strcmp(argv[5], "-") != 0){ sscanf(argv[5], "%lf", &start_line_TimeDaySeconds); sflag = 2; } if(strcmp(argv[6], "-") != 0){ sscanf(argv[6], "%lf", &start_line_TimeSeconds); sflag = 3; } if(strcmp(argv[7], "-") != 0){ sscanf(argv[7], "%lf", &boundLat1); sflag = 4; if( (boundLat1 > 90.0) || (boundLat1 < -90.0) ){ printf("\nLatitude bounding %12.6f is incorrect. Latitude should be between -90.0 and 90.0 Degrees.\n", boundLat1); exit(-1); } } if(strcmp(argv[8], "-") != 0){ sscanf(argv[8], "%d", &end_line); eflag = 1; } if(argc > 9){ if(strcmp(argv[9], "-") != 0){ sscanf(argv[9], "%d", &num_lines); eflag = 2; } } if(argc > 10){ if(strcmp(argv[10], "-") != 0){ sscanf(argv[10], "%lf", &end_line_TimeDaySeconds); eflag = 3; } } if(argc > 11){ if(strcmp(argv[11], "-") != 0){ sscanf(argv[11], "%lf", &end_line_TimeSeconds); eflag = 4; } } if(argc > 12){ if(strcmp(argv[12], "-") != 0){ sscanf(argv[12], "%lf", &timeduration); eflag = 5; } } if(argc > 13){ if(strcmp(argv[13], "-") != 0){ sscanf(argv[13], "%lf", &boundLat2); eflag = 6; if( (boundLat2 > 90.0) || (boundLat2 < -90.0) ){ printf("\nLatitude bounding %12.6f is incorrect. Latitude should be between -90.0 and 90.0 Degrees.\n", boundLat2); exit(-1); } } } oflag=0; for(i=8;i<(argc);i++){ oflag = oflag + (strcmp(argv[i], "-"));} if(!oflag){ printf("\n ERROR: Enter one of the parameters (e_line, num_lines, e_TimeofDay, e_TimeofFile, duration, boundLat2) to indicate end of frame. \n"); exit(-1); } if(argc > 14){ if(strcmp(argv[14], "-") != 0){ sscanf(argv[14], "%lf", &ref_alt); if(fabs(ref_alt) > 10000.0) {printf("\nWARNING: Reference elevation value %14.6f is too large. Setting reference elevation to 0.0\n", ref_alt); ref_alt =0.0;} } } if(argc > 15){ if(strcmp(argv[15], "-") != 0){ sscanf( argv[15], "%s", orbFileDir ); orbFileFlag = 1; } } if(argc > 16){ if(strcmp(argv[16], "-") != 0){ sscanf( argv[16], "%d", &HApadflag ); } } if((sflag == 4) || (eflag ==6) ){ if(HApadflag == 0) { EdgeTargetFullCompressionPadding = 0; } else if(HApadflag == 1){ EdgeTargetFullCompressionPadding = 550; } else { printf("\nUnknown value for parameter HApad: %d\n", HApadflag); printf("Should be either 1 or 0\n"); exit(-1); } } /* completed reading all required parameters, Checks on bounds of start line and end line will be done later in the code. */ /* print user suppied data to screen */ printf("\n************ User specified parameters. **********************************\n"); printf("input Image mode file (multi-frame Level 0 file) : %s\n", argv[1]); printf("input auxiliary instrument characterization data file : %s\n", argv[2]); printf("output Image mode file (single frame Level 0 file) : %s\n", argv[3]); if(sflag == 1){ printf("starting raw data line number with respect to start of multi-frame file : %s\n", argv[4]);} if(sflag == 2){ printf("starting time of the frame in seconds from start of the day : %s sec.\n", argv[5]);} if(sflag == 3){ printf("starting time of the frame in seconds from start of multi-frame file : %s sec.\n", argv[6]);} if(sflag == 4){ printf("first bounding latitude in decimal degrees for frame file : %s Deg.\n", argv[7]);} if(eflag == 1){ printf("last raw data line number with respect to start of multi-frame file : %s\n", argv[8]);} if(eflag == 2){ printf("number of raw data lines to be included in the frame file : %s\n", argv[9]);} if(eflag == 3){ printf("time of the last line to be included in the frame in seconds from start of the day : %s sec.\n", argv[10]);} if(eflag == 4){ printf("time of the last line to be included in the frame in seconds from start of multi-frame file : %s sec.\n", argv[11]);} if(eflag == 5){ printf("duration of data in seconds : %s sec.\n", argv[12]);} if(eflag == 6){ printf("second bounding latitude in decimal degrees for frame file : %s Deg.\n", argv[13]);} if(argc > 14) { printf("Reference elevation of terrain above ellipsoid : %10.2f meters\n", ref_alt);} else {printf("Using default reference elevation of terrain above ellipsoid : %10.2f meters\n", ref_alt);} if(orbFileFlag == 1) {printf("Precise orbit directory : %s\n", orbFileDir );} else {printf("Using orbit prapagation using state vector in file header \n");} if((sflag == 4) || (eflag ==6) ){printf("Half aperture padding (lines) at each end of segemented output : %d \n", EdgeTargetFullCompressionPadding);} printf("************ ***********************************\n"); printf("\n"); if (is_bigendian()) { printf("It is running under Unix or Mac...\n\n"); is_littlendian = 0; } else { printf("It is running under Linux or Windows...\n\n"); is_littlendian = 1; } /* open files */ insFilePtr = fopen( insFileName, "rb" ); if ( insFilePtr == NULL ) { printf( "*** ERROR - cannot open file: %s\n", insFileName ); printf( "\n" ); exit( -1 ); } outFilePtr = fopen( outFileName, "wb" ); if ( outFilePtr == NULL ) { printf( "*** ERROR - cannot open file: %s\n", outFileName ); printf( "\n" ); exit( -1 ); } /* read MPH of ins file */ printf( "Reading MPH of ins file...\n" ); mphPtr = ( char * ) malloc( sizeof( char ) * mphSize ); if ( mphPtr == NULL ){ printf( "ERROR - mph allocation memory\n" ); exit( -1 ); } if ( (fread( mphPtr, sizeof( char ), mphSize, insFilePtr ) ) != mphSize ){ printf( "ERROR - mph read error\n\n" ); exit( -1 ); } mphIns = readMph( mphPtr, printInsMphIfZero ); /* extract information from MPH */ free ( mphPtr ); /* read SPH from ins file */ printf( "Reading SPH from ins file...\n" ); sphPtr = ( char * ) malloc( sizeof( char ) * mphIns.sphSize ); if ( sphPtr == NULL ){ printf( "ERROR - sph allocation memory\n" ); exit( -1 ); } if ( (fread( sphPtr, sizeof( char ), mphIns.sphSize, insFilePtr ) ) != mphIns.sphSize ){ printf( "ERROR - sph read error\n\n" ); exit( -1 ); } sphIns = readSphAux( sphPtr, printInsSphIfZero, mphIns ); /* extract information from SPH */ free ( sphPtr ); /* read GADS from ins file */ printf( "Reading GADS from ins file...\n" ); if ( (fread( &insGads, sizeof( insGads ), 1, insFilePtr ) ) != 1 ){ printf( "sizeof( insGads ): %d\n", sizeof( insGads ) ); printf( "ERROR - gads read error\n\n" ); printf( "%d %d %d\n", 171648, sizeof ( insGads ), 171648-sizeof( insGads ) ); exit( -1 ); } if (is_littlendian) { byte_swap_InsGads( &insGads ); } if ( printInsGadsIfZero == 0 ) printInsGads( insGads ); fclose( insFilePtr ); /* open image file */ sscanf( argv[ 1 ], "%s", imFileName ); imFilePtr = fopen( imFileName, "rb" ); if ( imFilePtr == NULL ) { printf( "*** ERROR - cannot open file: %s\n", imFileName ); printf( "\n" ); exit( -1 ); } /* read image MPH */ printf( "Reading image MPH...\n" ); mphPtr = ( char * ) malloc( sizeof( char ) * mphSize ); if ( mphPtr == NULL ){ printf( "ERROR - mph allocation memory\n" ); exit( -1 ); } if ( (fread( mphPtr, sizeof( char ), mphSize, imFilePtr ) ) != mphSize ){ printf( "ERROR - mph read error\n\n" ); exit( -1 ); } mph = readMph( mphPtr, printImMphIfZero ); /* extract information from MPH */ /* read image SPH */ printf( "Reading image SPH...\n" ); sphPtr = ( char * ) malloc( sizeof( char ) * mph.sphSize ); if ( sphPtr == NULL ){ printf( "ERROR - sph allocation memory\n" ); exit( -1 ); } if ( (fread( sphPtr, sizeof( char ), mph.sphSize, imFilePtr ) ) != mph.sphSize ){ printf( "ERROR - sph read error\n\n" ); exit( -1 ); } sph = readSph( sphPtr, printImSphIfZero, mph ); /* extract information from SPH */ /* read the sensing start date and time of long strip line file */ sscanf(mph.sensingStart, "%2d-%3c-%4d %2d:%2d:%lf", &sday, smonthstr, &syear, &shour, &smin,&sseconds); if(strcmp(smonthstr,"JAN")==0){smonth =1;} if(strcmp(smonthstr,"FEB")==0){smonth =2;} if(strcmp(smonthstr,"MAR")==0){smonth =3;} if(strcmp(smonthstr,"APR")==0){smonth =4;} if(strcmp(smonthstr,"MAY")==0){smonth =5;} if(strcmp(smonthstr,"JUN")==0){smonth =6;} if(strcmp(smonthstr,"JUL")==0){smonth =7;} if(strcmp(smonthstr,"AUG")==0){smonth =8;} if(strcmp(smonthstr,"SEP")==0){smonth =9;} if(strcmp(smonthstr,"OCT")==0){smonth =10;} if(strcmp(smonthstr,"NOV")==0){smonth =11;} if(strcmp(smonthstr,"DEC")==0){smonth =12;} sscanf(mph.sensingStop, "%2d-%3c-%4d %2d:%2d:%lf", &eday, emonthstr, &eyear, &ehour, &emin,&eseconds); if(strcmp(emonthstr,"JAN")==0){emonth =1;} if(strcmp(emonthstr,"FEB")==0){emonth =2;} if(strcmp(emonthstr,"MAR")==0){emonth =3;} if(strcmp(emonthstr,"APR")==0){emonth =4;} if(strcmp(emonthstr,"MAY")==0){emonth =5;} if(strcmp(emonthstr,"JUN")==0){emonth =6;} if(strcmp(emonthstr,"JUL")==0){emonth =7;} if(strcmp(emonthstr,"AUG")==0){emonth =8;} if(strcmp(emonthstr,"SEP")==0){emonth =9;} if(strcmp(emonthstr,"OCT")==0){emonth =10;} if(strcmp(emonthstr,"NOV")==0){emonth =11;} if(strcmp(emonthstr,"DEC")==0){emonth =12;} sensestartyyyymmdd = syear * 10000 + smonth * 100 + sday; sensestopyyyymmdd = eyear * 10000 + emonth * 100 + eday; if(orbFileFlag == 1) { if( system("touch envisat_dor.tmp.txt") == -1) { printf("\n Cannot run system command : touch envisat_dor.tmp.txt \n"); exit(-1); } if( system("rm envisat_dor.tmp.txt") == -1) { printf("\n Cannot run system command : rm envisat_dor.tmp.txt \n"); exit(-1); } sprintf(listcmd, "ls -1 %s > envisat_dor.tmp.txt", orbFileDir); // printf("\n %s\n", listcmd); if( system(listcmd) == -1) { printf("\n Cannot run system command : ls -1 orbFileDir >> envisat_dor.tmp.txt \n"); exit(-1); } if ( (lsPtr = fopen("envisat_dor.tmp.txt", "r")) == NULL ) { printf( "*** ERROR - cannot open file: %s\n", "envisat_dor.tmp.txt" ); printf( "\n" ); exit( -1 ); } do{ // printf("\n TEST 1 \n"); if( (fgets(dummyOrbFileName,62,lsPtr)) != NULL) { // printf("\n TEST 1.5 \n"); sscanf(dummyOrbFileName,"%12c-%17c%8d%8c%8d%7c", dummyphrase1,dummyphrase2, &orbdate1, dummyphrase3, &orbdate2, dummyphrase4); // printf("\n TEST 2 \n"); if( (orbdate1 < sensestartyyyymmdd) && (orbdate2 > sensestopyyyymmdd)) { strcat(orbFileName, orbFileDir); strcat(orbFileName, dummyOrbFileName); printf("\nPrecise orbit file found: %s \n", orbFileName); exitflag = 1; } } else { exitflag = 1; printf("\nWARNING: Precise orbit file not found.\n"); printf("Listed and read all files in orbit file directory: %s\n", orbFileDir); printf("Using orbit propagator instead of precise orbit interpolation.\n"); orbFileFlag = 0; } }while(exitflag == 0); fclose(lsPtr); } if(orbFileFlag == 1) { OrbitPropErrDistance = 0.0; /* No Error in orbit propagation when precise orbit is used for interpolation */ } else { OrbitPropErrDistance = 0.05*3600 *30; /* Error of 0.05 degrees in orbit propagation with 1 arc sec ~ 30 meters */ } /* READ ANTENNA REFERENCE ELEVATION ANGLE */ strncpy(swathis,sph.swath,3); printf("\nSWATH = %s\n", swathis); if (strcmp(swathis,"IS1") == 0) ref_elevation_angle = insGads.refElevAngleIs1; if (strcmp(swathis,"IS2") == 0) ref_elevation_angle = insGads.refElevAngleIs2; if (strcmp(swathis,"IS3") == 0) ref_elevation_angle = insGads.refElevAngleIs3Ss2; if (strcmp(swathis,"SS2") == 0) ref_elevation_angle = insGads.refElevAngleIs3Ss2; if (strcmp(swathis,"IS4") == 0) ref_elevation_angle = insGads.refElevAngleIs4Ss3; if (strcmp(swathis,"SS3") == 0) ref_elevation_angle = insGads.refElevAngleIs4Ss3; if (strcmp(swathis,"IS5") == 0) ref_elevation_angle = insGads.refElevAngleIs5Ss4; if (strcmp(swathis,"SS4") == 0) ref_elevation_angle = insGads.refElevAngleIs5Ss4; if (strcmp(swathis,"IS6") == 0) ref_elevation_angle = insGads.refElevAngleIs6Ss5; if (strcmp(swathis,"SS5") == 0) ref_elevation_angle = insGads.refElevAngleIs6Ss5; if (strcmp(swathis,"IS7") == 0) ref_elevation_angle = insGads.refElevAngleIs7; if (strcmp(swathis,"SS1") == 0) ref_elevation_angle = insGads.refElevAngleSs1; heading = sph.satTrack; if(heading < 0.){heading = heading + 360.;} /* If user provides 2 latitude bounds, account for ascending/descending tracks in computing start end times */ /* default in the case only one of the two latitude bounds is given*/ if(sflag == 4){start_lat = boundLat1;} if(eflag == 6){end_lat = boundLat2;} if(heading > 270.) { printf("Satellite track heading: %10.6f Deg. (ASCENDING TRACK)\n", heading); if((sflag == 4) && (eflag == 6)) { if(boundLat1 <= boundLat2) { start_lat = boundLat1; end_lat = boundLat2; } else { start_lat = boundLat2; end_lat = boundLat1; } } } else if(heading < 270.) { printf("\nSatellite track heading: %10.6f Deg. (DESCENDING TRACK)\n", heading); if((sflag == 4) && (eflag == 6)) { if(boundLat1 <= boundLat2) { start_lat = boundLat2; end_lat = boundLat1; } else { start_lat = boundLat1; end_lat = boundLat2; } } } printf("Beam reference elevation angle: %10.6f Deg.\n", ref_elevation_angle); /* In this code elevation angle is measured with respect to local horizontal at satellite position */ ref_elevation_angle = -(90.0-ref_elevation_angle); /* LOOP OVER ALL LINES TO FIND, NEAR/FAR RANGES, MJD2000 FOR FIRST AND LAST LINE OF MDSR */ /* After setting the starting times, near and far range re-set fseek counter to start of MDSR */ bytesRead = 0; fline = 0; lline = sph.dsd[ 0 ].numDsr-2; /*parameters for penultimate and last line are read and stored in case the last line is not an echo */ /* Assumption 1: There is at max only one peiodic calibration noise pulse in any two consecutive lines. */ printf("\nReading input file ...\n"); for ( i = 0; i < sph.dsd[ 0 ].numDsr; i++ ) { bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeDays, 4, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_long, &mdsrDsrTimeDays, 1); } /* header added to the ISP by the Front End Processor (FEP) */ bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeSeconds, 4, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_ulong, &mdsrDsrTimeSeconds, 1); } bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeMicroseconds, 4, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_ulong, &mdsrDsrTimeMicroseconds, 1); } bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeDays, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeSeconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeMicroseconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrIspLength, 2, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_ushort, &mdsrIspLength, 1); } bytesRead = bytesRead + 2 * fread( &mdsrCrcErrs, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrRsErrs, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrSpare1, 2, 1, imFilePtr ); /* 6-byte ISP Packet Header */ bytesRead = bytesRead + 2 * fread( &mdsrPacketIdentification, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketSequenceControl, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketLength, 2, 1, imFilePtr ); /* 30-byte Data Field Header in Packet Data Field */ bytesRead = bytesRead + 30 * fread( &mdsrPacketDataHeader, 30, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_ushort, &mdsrPacketDataHeader, 15); } dataFieldHeaderLength = mdsrPacketDataHeader[ 0 ]; antennaBeamSetNumber = (unsigned char) ( ( mdsrPacketDataHeader[ 6 ] >> 2 ) & 63); cyclePacketCount = ( mdsrPacketDataHeader[ 7 ] & 4095); echoFlag = (unsigned char) ( ( mdsrPacketDataHeader[ 7 ] >> 15 ) & 1); priCodeword = mdsrPacketDataHeader[ 8 ]; windowStartTimeCodeword = mdsrPacketDataHeader[ 9 ]; windowLengthCodeword = mdsrPacketDataHeader[ 10 ]; fseek( imFilePtr, mdsrIspLength+1-30, SEEK_CUR ); bytesRead = bytesRead + mdsrIspLength+1-30; /* Store first and last line parameters of the long stripline file */ if((i == fline) && (echoFlag == 1)){ pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; rho_near_firstline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; //rho_far_firstline = rho_near_firstline + (outSamples/ insGads.sampRate)*c/2. ; rho_far_firstline = rho_near_firstline + ((outSamples-523)/ insGads.sampRate)*c/2. ; DayMJD2000TimeFirstMDSRLine = mdsrDsrTimeDays; SecondsMJD2000TimeFirstMDSRLine = mdsrDsrTimeSeconds; MicroSecondsMJD2000TimeFirstMDSRLine = mdsrDsrTimeMicroseconds; TimeFirstMDSRLine = SecondsMJD2000TimeFirstMDSRLine + MicroSecondsMJD2000TimeFirstMDSRLine * 1e-6; DayMJDTimeFirstMDSRLine = DayMJD2000TimeFirstMDSRLine + 51544; } else if((i == fline) && (echoFlag != 1)){ /*if first line is not echo then update fline to second line to read range parameters */ fline = fline +1; /*Only range parameters from second line are used, time parameters are still for first line even if not echo */ pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; rho_near_firstline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; //rho_far_firstline = rho_near_firstline + (outSamples/ insGads.sampRate)*c/2. ; rho_far_firstline = rho_near_firstline + ((outSamples-523)/ insGads.sampRate)*c/2. ; } if((i == lline) && (echoFlag == 1)){ pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; rho_near_lastline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; //rho_far_lastline = rho_near_lastline + (outSamples/ insGads.sampRate)*c/2. ; rho_far_lastline = rho_near_lastline + ((outSamples-523)/ insGads.sampRate)*c/2. ; DayMJD2000TimeLastMDSRLine = mdsrDsrTimeDays; SecondsMJD2000TimeLastMDSRLine = mdsrDsrTimeSeconds; MicroSecondsMJD2000TimeLastMDSRLine = mdsrDsrTimeMicroseconds; TimeLastMDSRLine = SecondsMJD2000TimeLastMDSRLine + MicroSecondsMJD2000TimeLastMDSRLine * 1e-6; DayMJDTimeLastMDSRLine = DayMJD2000TimeLastMDSRLine + 51544; lline = lline +1; /* if lline is echo, read penultimate line and update lline to lline +1 i.e. last line */ } /* if last line also echo, range parameters will be updated, if not, range parameters from penultimate line will be used*/ else if((i == lline) && (echoFlag != 1)){ /* even if not an echo, time parameters for last line are updated */ lline = lline +1; DayMJD2000TimeLastMDSRLine = mdsrDsrTimeDays; SecondsMJD2000TimeLastMDSRLine = mdsrDsrTimeSeconds; MicroSecondsMJD2000TimeLastMDSRLine = mdsrDsrTimeMicroseconds; TimeLastMDSRLine = SecondsMJD2000TimeLastMDSRLine + MicroSecondsMJD2000TimeLastMDSRLine * 1e-6; DayMJDTimeLastMDSRLine = DayMJD2000TimeLastMDSRLine + 51544; } } printf("Slant range (near and far) for first line of input file: %14.6f m, %14.6f m\n",rho_near_firstline, rho_far_firstline); printf("Slant range (near and far) for last line of input file: %14.6f m, %14.6f m\n",rho_near_lastline, rho_far_lastline); /* Total length of the long strip line file in terms of seconds and line numbers. Daughter frame file created should be of lesser length than this. */ MaxTimeDuration = TimeLastMDSRLine - TimeFirstMDSRLine; MaxNumLines = sph.dsd[ 0 ].numDsr; /* Given State vector for satellite */ sVEC[0] = mph.xPosition; sVEC[1] = mph.yPosition; sVEC[2] = mph.zPosition; sVEC[3] = mph.xVelocity; sVEC[4] = mph.yVelocity; sVEC[5] = mph.zVelocity; /* State vector time */ sscanf(mph.stateVectorTime, "%2d-%3c-%4d %2d:%2d:%lf", &svday, svmonthstr, &svyear, &svhour, &svmin,&svseconds); if(strcmp(svmonthstr,"JAN")==0){svmonth =1;} if(strcmp(svmonthstr,"FEB")==0){svmonth =2;} if(strcmp(svmonthstr,"MAR")==0){svmonth =3;} if(strcmp(svmonthstr,"APR")==0){svmonth =4;} if(strcmp(svmonthstr,"MAY")==0){svmonth =5;} if(strcmp(svmonthstr,"JUN")==0){svmonth =6;} if(strcmp(svmonthstr,"JUL")==0){svmonth =7;} if(strcmp(svmonthstr,"AUG")==0){svmonth =8;} if(strcmp(svmonthstr,"SEP")==0){svmonth =9;} if(strcmp(svmonthstr,"OCT")==0){svmonth =10;} if(strcmp(svmonthstr,"NOV")==0){svmonth =11;} if(strcmp(svmonthstr,"DEC")==0){svmonth =12;} /* 2 body propagation of state vector for first and last line times */ printf("\nSatellite vector time reported in header: \n"); printf("Year: %d , Month %d, Day %d, Hour %d, Minute %d, Seconds %14.6f \n",svyear, svmonth, svday, svhour, svmin, svseconds); printf("Satellite position vector: %14.4f %14.4f %14.4f meters", sVEC[0], sVEC[1], sVEC[2]); printf("\nSatellite velocity vector: %14.4f %14.4f %14.4f meters/sec\n",sVEC[3], sVEC[4], sVEC[5]); svMJD = date2MJD(svyear, svmonth, svday, svhour, svmin, svseconds); svMJDday = floor(svMJD); delT1= (DayMJDTimeFirstMDSRLine-svMJDday)*86400 + TimeFirstMDSRLine - ((double)svhour*3600. + (double)svmin*60. + svseconds); // printf("\n Time from reported satellite vector time to time of first line = %14.10f seconds\n", delT1); if(orbFileFlag == 1){ orbitInterpolationPoly(svyear, svmonth, svday, svhour, svmin, svseconds, delT1, orbFileName, sVECfirst); } else { orbitPropogation2body_J2Cartesian(svyear, svmonth, svday, svhour, svmin, svseconds, delT1, sVEC, sVECfirst); } delT2= (DayMJDTimeLastMDSRLine-svMJDday)*86400 + TimeLastMDSRLine - ((double)svhour*3600. + (double)svmin*60. + svseconds); // printf("\n Time from reported satellite vector time to time of last line = %14.10f \n seconds ", delT2); if(orbFileFlag == 1){ orbitInterpolationPoly(svyear, svmonth, svday, svhour, svmin, svseconds, delT2, orbFileName, sVEClast); } else { orbitPropogation2body_J2Cartesian(svyear, svmonth, svday, svhour, svmin, svseconds, delT2, sVEC, sVEClast); } /*Calculate latitudes and longitudes for the 4 corners of raw data file */ geolocate_radarLOS_ellipsoidWGS84(sVECfirst, ref_elevation_angle, rho_near_firstline, DIR_FLAG, ref_alt, &olat, &olon); clat[0] = olat; clon[0] = olon; geolocate_radarLOS_ellipsoidWGS84(sVECfirst, ref_elevation_angle, rho_far_firstline, DIR_FLAG, ref_alt, &olat, &olon); clat[1] = olat; clon[1] = olon; geolocate_radarLOS_ellipsoidWGS84(sVEClast, ref_elevation_angle, rho_near_lastline, DIR_FLAG, ref_alt, &olat, &olon); clat[2] = olat; clon[2] = olon; geolocate_radarLOS_ellipsoidWGS84(sVEClast, ref_elevation_angle, rho_far_lastline, DIR_FLAG, ref_alt, &olat, &olon); clat[3] = olat; clon[3] = olon; printf("\nApproximate geodetic co-ordinates of four corners of INPUT strip line data (first line near, first line far, last line near, last line far)\n"); for(i=0;i<4;i++){ printf("inlat[%d] = %10.6f Deg. inlon[%d] = %10.6f Deg.\n", i+1,clat[i],i+1,clon[i]); } /* Calculate average rates of latitude and longitude changes in along track and slant range direction. */ ave_dLat_daz = 0.5*((clat[2]-clat[0])/(double)MaxNumLines + (clat[3]-clat[1])/(double)MaxNumLines); ave_dLat_dr = 0.5*((clat[2]-clat[0])/outSamples + (clat[3]-clat[1])/outSamples); // printf("\n Average rate of latitude change per line: %lf\n", ave_dLat_daz); switch(sflag) { case 1: start_line_TimeDaySeconds = TimeFirstMDSRLine + (start_line-1) * pri; break; case 2: start_line = floor( (start_line_TimeDaySeconds-TimeFirstMDSRLine)/pri ) + 1; break; case 3: start_line_TimeDaySeconds = start_line_TimeSeconds + TimeFirstMDSRLine; start_line = floor( (start_line_TimeDaySeconds-TimeFirstMDSRLine)/pri ) + 1; break; case 4: deltaLat = start_lat - clat[0]; estStartLine = floor(deltaLat/ave_dLat_daz); StartLineErr = 0; /* Find the first line based on required latitude at near range through iteration */ do{ estStartLine = estStartLine - StartLineErr; delT = (DayMJDTimeFirstMDSRLine - svMJDday)*86400 + TimeFirstMDSRLine + estStartLine * pri - ((double)svhour*3600. + (double)svmin*60. + svseconds); if(orbFileFlag == 1){ orbitInterpolationPoly(svyear, svmonth, svday, svhour, svmin, svseconds, delT, orbFileName, sVECn); } else { orbitPropogation2body_J2Cartesian(svyear, svmonth, svday, svhour, svmin, svseconds, delT, sVEC, sVECn); } geolocate_radarLOS_ellipsoidWGS84(sVECn, ref_elevation_angle, rho_near_firstline, DIR_FLAG, ref_alt, &nlat, &nlon); deltanLat = (nlat - start_lat); StartLineErr = floor(deltanLat/ave_dLat_daz); }while( fabs(deltanLat) > fabs(ave_dLat_daz)); vsatv[0] =sVECn[3]; vsatv[1] =sVECn[4]; vsatv[2] =sVECn[5]; vsat = norm(vsatv); OrbitPropErrTime = OrbitPropErrDistance/vsat; OrbitPropErrPadding = floor(OrbitPropErrTime/pri); padding = OrbitPropErrPadding + EdgeTargetFullCompressionPadding; if((estStartLine < 1) || (estStartLine > MaxNumLines)) { printf("\nWARNING: Given latitude bounding: %10.6f is outside of input image area. Setting to closest boundary of input image.\n",start_lat); estStartLine = 1; } start_line = estStartLine - padding; if(start_line < 1){printf("WARNING: After padding, output file starting line is before start of input file. Setting first line of output as first line of input.\n"); start_line =1;} if( (eflag == 6) || (eflag ==2) || (eflag ==5)){ printf("\nActual number of lines padded at start of file is: %d\n", min_int(padding, (start_line-1)) ); } start_line_TimeDaySeconds = TimeFirstMDSRLine + (start_line-1) * pri; break; default: printf("\nERROR: Start line of output data has not been set.\n"); exit(-1); break; } /* end of switch(sflag) */ if( (start_line < 1) || (start_line > MaxNumLines)) { printf("\nWARNING: Starting line number: %d. Starting line number for output file creation is out of bounds. It should be between 1 and %d.\n", start_line, MaxNumLines); printf("WARNING: Starting line will be set to closest bound.\n"); if(start_line < 1) {printf("WARNING: Setting starting line number to : %d\n", 1); start_line =1;} if(start_line > MaxNumLines) {printf("WARNING: Setting starting line number to : %d\n", MaxNumLines); start_line = MaxNumLines;} } switch(eflag) { case 1: if(end_line > MaxNumLines){printf("WARNING: Output end line is after end of input file, setting last line of output as last line of input.\n"); end_line =MaxNumLines;} if((sflag == 4) && (start_line > end_line)){ start_line = start_line + 2* padding; /* pad the other way because start line will now be end line boundLat1 is 'after' the end line*/ if(start_line > MaxNumLines){printf("WARNING: After padding, output end line is after end of input file. Setting last line of output as last line of input.\n"); start_line = MaxNumLines;} printf("\nWARNING: Latitude bound provided is 'ahead' in time with respect to given end line or end time. end line will be made the start line.\n"); temp_line = end_line; end_line = start_line; start_line = temp_line; printf("\nActual number of lines padded at end of file is: %d \n", min_int(padding, (MaxNumLines-end_line)) ); } else if((sflag == 4) && (start_line < end_line)) { printf("\nActual number of lines padded at start of file is: %d\n", min_int(padding, (start_line-1)) ); } num_lines = end_line - start_line +1; /* end line is provided by user */ break; case 2: end_line = start_line + num_lines - 1; /* num_lines is provided by user */ if(end_line > MaxNumLines){printf("WARNING: Output end line is after end of input file, setting last line of output as last line of input.\n"); end_line =MaxNumLines;} break; case 3: end_line = start_line + floor( (end_line_TimeDaySeconds-start_line_TimeDaySeconds)/pri ); if(end_line > MaxNumLines){printf("WARNING: Output end line is after end of input file, setting last line of output as last line of input.\n"); end_line =MaxNumLines;} if((sflag == 4) && (start_line > end_line)){ start_line = start_line + 2* padding; /* pad the other way because start line will now be end line boundLat1 is 'after' the end line*/ if(start_line > MaxNumLines){printf("WARNING: After padding, end line is after end of input file, setting last line of output as last line of input.\n"); start_line = MaxNumLines;} printf("\nWARNING: Latitude bound provided is 'ahead' in time with respect to given end line or end time. end line will be made the start line.\n"); temp_line = end_line; end_line = start_line; start_line = temp_line; printf("\nActual number of lines padded at end of file is: %d \n", min_int(padding, (MaxNumLines-end_line)) ); } else if((sflag == 4) && (start_line < end_line)) { printf("\nActual number of lines padded at start of file is: %d\n", min_int(padding, (start_line-1)) ); } num_lines = end_line - start_line +1; break; case 4: end_line = floor(end_line_TimeSeconds / pri) + 1; if(end_line > MaxNumLines){printf("WARNING: Output end line is after end of input file, setting last line of output as last line of input.\n"); end_line =MaxNumLines;} if((sflag == 4) && (start_line > end_line)){ start_line = start_line + 2* padding; /* pad the other way because start line will now be end line boundLat1 is 'after' the end line*/ if(start_line > MaxNumLines){printf("WARNING: After padding, output end line is after end of input file, setting last line of output as last line of input.\n"); start_line = MaxNumLines;} printf("\nWARNING: Latitude bound provided is 'ahead' in time with respect to given end line or end time. end line will be made the start line.\n"); temp_line = end_line; end_line = start_line; start_line = temp_line; printf("\nActual number of lines padded at end of file is: %d \n", min_int(padding, (MaxNumLines-end_line)) ); } else if((sflag == 4) && (start_line < end_line)) { printf("\nActual number of lines padded at start of file is: %d\n", min_int(padding, (start_line-1)) ); } num_lines = end_line - start_line +1; break; case 5: num_lines = floor(timeduration / pri); end_line = start_line + num_lines -1; if(end_line > MaxNumLines){printf("WARNING: Output end line is after end of input file, setting last line of output as last line of input.\n"); end_line =MaxNumLines;} break; case 6: deltaLat = end_lat - clat[1]; estLastLine = floor(deltaLat/ave_dLat_daz); EndLineErr =0; /* Find the last line based on required latitude at far range */ do{ estLastLine = estLastLine - EndLineErr; delT = (DayMJDTimeFirstMDSRLine - svMJDday)*86400 + TimeFirstMDSRLine + estLastLine * pri - ((double)svhour*3600. + (double)svmin*60. + svseconds); if(orbFileFlag == 1){ orbitInterpolationPoly(svyear, svmonth, svday, svhour, svmin, svseconds, delT, orbFileName, sVECn); } else { orbitPropogation2body_J2Cartesian(svyear, svmonth, svday, svhour, svmin, svseconds, delT, sVEC, sVECn); } geolocate_radarLOS_ellipsoidWGS84(sVECn, ref_elevation_angle, rho_far_lastline, DIR_FLAG, ref_alt, &nlat, &nlon); deltanLat = (nlat - end_lat); EndLineErr = floor(deltanLat/ave_dLat_daz); }while( fabs(deltanLat) > fabs(ave_dLat_daz)); vsatv[0] =sVECn[3]; vsatv[1] =sVECn[4]; vsatv[2] =sVECn[5]; vsat = norm(vsatv); OrbitPropErrTime = OrbitPropErrDistance/vsat; OrbitPropErrPadding = floor(OrbitPropErrTime/pri); padding = OrbitPropErrPadding + EdgeTargetFullCompressionPadding; end_line = estLastLine; if((end_line < 0) || (end_line > MaxNumLines)) { printf("\nWARNING: Given latitude bounding: %10.6f is outside of input image area.\n", end_lat); if(end_line < 1) {printf("WARNING: Setting end line number to input line number: %d\n", 1); end_line =1;} if(end_line > MaxNumLines) {printf("WARNING: Setting end line number to input file line number: %d\n", MaxNumLines); end_line = MaxNumLines;} } /* case when only 1 Latitude bounding is given (starting line of output file is given in terms of line number or time), check that num_lines >= 1 */ if((sflag != 4) && (end_line < start_line) ){ printf("\nWARNING: Latitude bound provided is 'behind' in time with respect to given start line or start time. start line will be made the end line.\n"); temp_line = end_line; end_line = start_line; start_line = temp_line; start_line = start_line - padding; if(start_line < 1){printf("WARNING: After padding, output starting line is before start of input file, setting first line of output as first line of input.\n"); start_line =1;} printf("\nActual number of lines padded at start of file is: %d\n", min_int(padding, (start_line-1)) ); } else { end_line = end_line + padding; if(end_line > MaxNumLines){printf("WARNING: After padding, output end line is after end of input file, setting last line of output as last line of input.\n"); end_line =MaxNumLines;} printf("\nActual number of lines padded at end of file is: %d \n", min_int(padding, (MaxNumLines-end_line)) ); } num_lines = end_line - start_line +1; break; default: printf("\nERROR: End line of output data has not been set.\n"); exit(-1); break; } /* end of switch(eflag) */ if( (end_line < start_line) || (end_line > MaxNumLines)) { printf("\nERROR: Ending line number: %d. It should be between %d and %d. Re-check input parameters\n", end_line, start_line, MaxNumLines); exit(-1); } printf("\nStarting line of output file in input file: %d\n", start_line); printf("Last line of output file in input file: %d\n", end_line); printf("Number of lines in output file = %d\n", num_lines); curser = mphSize + mph.sphSize; fseek( imFilePtr, curser, SEEK_SET); bytesRead = 0; bytesToBeWritten = 0; for ( i = 1; i <= sph.dsd[ 0 ].numDsr; i++ ) { bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeDays, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeSeconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeMicroseconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeDays, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeSeconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeMicroseconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrIspLength, 2, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_ushort, &mdsrIspLength, 1); } mdsrIspLengthSwapped = mdsrIspLength; bytesRead = bytesRead + 2 * fread( &mdsrCrcErrs, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrRsErrs, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrSpare1, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketIdentification, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketSequenceControl, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketLength, 2, 1, imFilePtr ); bytesRead = bytesRead + 30 * fread( &mdsrPacketDataHeader, 30, 1, imFilePtr ); bytesRead = bytesRead + (mdsrIspLengthSwapped+1-30) * fread( &mdsrLineChar, (mdsrIspLengthSwapped+1-30), 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_long, &mdsrDsrTimeDays, 1); swap_endian_order(e_tid_ulong, &mdsrDsrTimeSeconds, 1); swap_endian_order(e_tid_ulong, &mdsrDsrTimeMicroseconds, 1); swap_endian_order(e_tid_ushort, &mdsrPacketDataHeader, 15); } dataFieldHeaderLength = mdsrPacketDataHeader[ 0 ]; antennaBeamSetNumber = (unsigned char) ( ( mdsrPacketDataHeader[ 6 ] >> 2 ) & 63); cyclePacketCount = ( mdsrPacketDataHeader[ 7 ] & 4095); echoFlag = (unsigned char) ( ( mdsrPacketDataHeader[ 7 ] >> 15 ) & 1); priCodeword = mdsrPacketDataHeader[ 8 ]; windowStartTimeCodeword = mdsrPacketDataHeader[ 9 ]; windowLengthCodeword = mdsrPacketDataHeader[ 10 ]; /* For the start and end lines store the FEP time tags to be written as sensing start and stop time in MPH of output file*/ if( (i >= start_line) && (i <= end_line) ){ if(i == start_line) { pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; dsr_rho_near_firstline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; //dsr_rho_far_firstline = dsr_rho_near_firstline + (outSamples/ insGads.sampRate)*c/2. ; dsr_rho_far_firstline = dsr_rho_near_firstline + ((outSamples-523)/ insGads.sampRate)*c/2. ; DsrStartDay = mdsrDsrTimeDays; DsrStartTime = (double)mdsrDsrTimeSeconds + (double)mdsrDsrTimeMicroseconds * 1e-6; DsrStartDayMJD = DsrStartDay + 51544; if(echoFlag == 1) {DsrStartFlag = 1;} } else if((i == start_line + 1) && (DsrStartFlag == 0)){ pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; dsr_rho_near_firstline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; // dsr_rho_far_firstline = dsr_rho_near_firstline + (outSamples/ insGads.sampRate)*c/2. ; dsr_rho_far_firstline = dsr_rho_near_firstline + ((outSamples-523)/ insGads.sampRate)*c/2. ; } if ((i == end_line-1) && (echoFlag == 1)) { pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; dsr_rho_near_lastline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; //dsr_rho_far_lastline = dsr_rho_near_lastline + (outSamples/ insGads.sampRate)*c/2. ; dsr_rho_far_lastline = dsr_rho_near_lastline + ((outSamples-523)/ insGads.sampRate)*c/2. ; } else if(i == end_line){ if(echoFlag == 1) { pri = priCodeword / insGads.sampRate; windowStartTime = windowStartTimeCodeword/ insGads.sampRate; outSamples = ((mdsrIspLength+1-30)/64)*63 + ((mdsrIspLength+1-30)%64)-1; dsr_rho_near_lastline = (insGads.timelineIm.rValues[ antennaBeamSetNumber-1 ]*pri+windowStartTime)*c/2.; //dsr_rho_far_lastline = dsr_rho_near_lastline + (outSamples/ insGads.sampRate)*c/2. ; dsr_rho_far_lastline = dsr_rho_near_lastline + ((outSamples-523)/ insGads.sampRate)*c/2. ; } DsrEndDay = mdsrDsrTimeDays; DsrEndTime = (double)mdsrDsrTimeSeconds + (double)mdsrDsrTimeMicroseconds * 1e-6; DsrEndDayMJD = DsrEndDay + 51544; } bytesToBeWritten = bytesToBeWritten + 68 + (mdsrIspLengthSwapped+1-30) ; } } printf("Slant range (near and far) for first line of output file: %14.6f m, %14.6f m\n",dsr_rho_near_firstline, dsr_rho_far_firstline); printf("Slant range (near and far) for last line of output file: %14.6f m, %14.6f m\n",dsr_rho_near_lastline, dsr_rho_far_lastline); delT1= (DsrStartDayMJD - svMJDday)*86400 + DsrStartTime - ((double)svhour*3600. + (double)svmin*60. + svseconds); // printf("\n Time from reported satellite vector time to time of first line = %14.10f seconds\n", delT1); if(orbFileFlag == 1){ orbitInterpolationPoly(svyear, svmonth, svday, svhour, svmin, svseconds, delT1, orbFileName, sVECfirst); } else { orbitPropogation2body_J2Cartesian(svyear, svmonth, svday, svhour, svmin, svseconds, delT1, sVEC, sVECfirst); } delT2= (DsrEndDayMJD - svMJDday)*86400 + DsrEndTime - ((double)svhour*3600. + (double)svmin*60. + svseconds); // printf("\n Time from reported satellite vector time to time of last line = %14.10f \n seconds ", delT2); if(orbFileFlag == 1){ orbitInterpolationPoly(svyear, svmonth, svday, svhour, svmin, svseconds, delT2, orbFileName, sVEClast); } else { orbitPropogation2body_J2Cartesian(svyear, svmonth, svday, svhour, svmin, svseconds, delT2, sVEC, sVEClast); } /*Calculate latitudes and longitudes for the 4 corners of OUTPUT data file */ geolocate_radarLOS_ellipsoidWGS84(sVECfirst, ref_elevation_angle, dsr_rho_near_firstline, DIR_FLAG, ref_alt, &olat, &olon); clat[0] = olat; clon[0] = olon; geolocate_radarLOS_ellipsoidWGS84(sVECfirst, ref_elevation_angle, dsr_rho_far_firstline, DIR_FLAG, ref_alt, &olat, &olon); clat[1] = olat; clon[1] = olon; geolocate_radarLOS_ellipsoidWGS84(sVEClast, ref_elevation_angle, dsr_rho_near_lastline, DIR_FLAG, ref_alt, &olat, &olon); clat[2] = olat; clon[2] = olon; geolocate_radarLOS_ellipsoidWGS84(sVEClast, ref_elevation_angle, dsr_rho_far_lastline, DIR_FLAG, ref_alt, &olat, &olon); clat[3] = olat; clon[3] = olon; printf("\nApproximate geodetic co-ordinates of four corners of OUTPUT data (first line near, first line far, last line near, last line far)\n"); for(i=0;i<4;i++){ printf("outlat[%d] = %10.6f Deg. outlon[%d] = %10.6f Deg.\n", i+1,clat[i],i+1,clon[i]); } printf("\nSize of ACSII header to be written to file = %d\n", mphSize + mph.sphSize); printf("Binary bytes to be written to file = %d\n", bytesToBeWritten); mphPtr2 = ( char * ) malloc( sizeof( char ) * mphSize ); if ( mphPtr2 == NULL ){ printf( "ERROR - mph allocation memory\n" ); exit( -1 ); } /* Copy MPH of long strip line file to MPH of output file */ memcpy(mphPtr2, mphPtr, mphSize); shour =0; smin = 0; /* Set sensing start time of oupput frame file */ sseconds = DsrStartTime; while(sseconds >= 60.0){ sseconds = sseconds -60.0; smin = smin +1; if(smin >= 60) { smin = smin -60; shour = shour +1; if(shour > 24) {shour = shour -24;} } } ssecint = floor(sseconds); eday = sday; eyear = syear; emonth = smonth; ehour = 0; emin = 0; /* Set sensing stop time of oupput frame file */ eseconds = DsrEndTime; while(eseconds >= 60.0){ eseconds = eseconds -60.0; emin = emin +1; if(emin >= 60) { emin = emin -60; ehour = ehour +1; if(ehour > 24) {ehour = ehour -24; eday = eday + 1;} } } esecint = floor(eseconds); duration = ehour * 3600 + emin * 60 + esecint - shour *3600 - smin * 60 - ssecint; sprintf(shh,"%02d",shour); sprintf(smm,"%02d",smin); sprintf(sss,"%02d",ssecint); sprintf(ssec,"%09.6f",sseconds); sprintf(ehh,"%02d",ehour); sprintf(emm,"%02d",emin); sprintf(ess,"%02d",esecint); sprintf(esec,"%09.6f",eseconds); strcat(shhmmss,shh); strcat(shhmmss,smm); strcat(shhmmss,sss); sprintf(secduration,"%08d",duration); /* Update PRODUCT NAME in MPH of output frame file pointer*/ memcpy(mphPtr2 + 0 + 9 + 23, shhmmss,6); memcpy(mphPtr2 + 0 + 9 + 30, secduration,8); /* Update SENSING START TIME in MPH of output frame file pointer*/ memcpy(mphPtr2 + 336 + 15 + 12, shh,2); memcpy(mphPtr2 + 336 + 15 + 15, smm,2); memcpy(mphPtr2 + 336 + 15 + 18, ssec,9); /* Update SENSING STOP TIME in MPH of output frame file pointer*/ memcpy(mphPtr2 + 380 + 14 + 12, ehh,2); memcpy(mphPtr2 + 380 + 14 + 15, emm,2); memcpy(mphPtr2 + 380 + 14 + 18, esec,9); totfilesize = bytesToBeWritten + mphSize + mph.sphSize; sprintf(totsize,"%021d",totfilesize); /* Update TOT_SIZE in MPH of output frame file pointer */ memcpy(strchr(mphPtr2 + 1066 + 0, '=')+1, totsize, 21); sphPtr2 = ( char * ) malloc( sizeof( char ) * mph.sphSize ); if ( sphPtr2 == NULL ){ printf( "ERROR - sph allocation memory\n" ); exit( -1 ); } memcpy(sphPtr2, sphPtr, mph.sphSize); sprintf(ds_size,"%020d",bytesToBeWritten); sprintf(num_dsr,"%010d",num_lines); /* Update DSD_SIZE and NUM_DSR in SPH of output frame file pointer */ memcpy(strchr(sphPtr2 + 836 + mph.dsdSize* 0 + 162 +0,'=') +2, ds_size, 20); memcpy(strchr(sphPtr2 + 836 + mph.dsdSize* 0 + 199 +0,'=') +2, num_dsr, 10); /* Write MPH of output frame file */ printf("\nWriting updated MPH of output file...\n"); if ( (fwrite( mphPtr2, sizeof ( char ), mphSize, outFilePtr) ) != mphSize) { printf(" ERROR - outFile write error\n\n"); exit(-1); } /* Write SPH of output frame file */ printf("\nWriting updated SPH of output file...\n"); if ( (fwrite( sphPtr2, sizeof ( char ), mph.sphSize, outFilePtr) ) != mph.sphSize) { printf(" ERROR - outFile write error\n\n"); exit(-1); } /* Loop over all the datasets and write the relevant lines to output file */ curser = mphSize + mph.sphSize; fseek( imFilePtr, curser, SEEK_SET); bytesRead = 0; bytesWritten = 0; printf("\nWriting records in output file...\n"); for ( i = 1; i <= sph.dsd[ 0 ].numDsr; i++ ) { bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeDays, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeSeconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrDsrTimeMicroseconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeDays, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeSeconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 4 * fread( &mdsrGsrtTimeMicroseconds, 4, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrIspLength, 2, 1, imFilePtr ); if (is_littlendian) { swap_endian_order(e_tid_ushort, &mdsrIspLength, 1); } mdsrIspLengthSwapped = mdsrIspLength; if (is_littlendian) { swap_endian_order(e_tid_ushort, &mdsrIspLength, 1); } bytesRead = bytesRead + 2 * fread( &mdsrCrcErrs, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrRsErrs, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrSpare1, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketIdentification, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketSequenceControl, 2, 1, imFilePtr ); bytesRead = bytesRead + 2 * fread( &mdsrPacketLength, 2, 1, imFilePtr ); bytesRead = bytesRead + 30 * fread( &mdsrPacketDataHeader, 30, 1, imFilePtr ); bytesRead = bytesRead + (mdsrIspLengthSwapped+1-30) * fread( &mdsrLineChar, (mdsrIspLengthSwapped+1-30), 1, imFilePtr ); if( (i >= start_line) && (i <= end_line) ){ if ( (fwrite( &mdsrDsrTimeDays, 4, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrDsrTimeSeconds, 4, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrDsrTimeMicroseconds, 4, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrGsrtTimeDays, 4, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrGsrtTimeSeconds, 4, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrGsrtTimeMicroseconds, 4, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrIspLength, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrCrcErrs, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrRsErrs, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrSpare1, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrPacketIdentification, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrPacketSequenceControl, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrPacketLength, 2, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrPacketDataHeader, 30, 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } if ( (fwrite( &mdsrLineChar, (mdsrIspLengthSwapped+1-30), 1, outFilePtr ) ) != 1 ){ printf( "ERROR - outFile write error\n\n" ); exit( -1 ); } bytesWritten = bytesWritten + 68 + (mdsrIspLengthSwapped+1-30) ; } } if(bytesToBeWritten != bytesWritten){ printf("\n Error in writing outfile. BytesToBeWritten = %d, BytesWritten = %d\n", bytesToBeWritten, bytesWritten); exit(-1);} free ( sphPtr ); free ( mphPtr ); free ( sphPtr2 ); free ( mphPtr2 ); fclose( outFilePtr ); fclose( imFilePtr ); printf( "\nDone.\n" ); return 0; } /* end main */ /**********************************************************************************************************************************/ struct mphStruct readMph( const char *mphPtr, const int printMphIfZero ) { struct mphStruct mph; if ( 1 == 0 ) { printf( "check:\n%s\n", mphPtr+1247 ); } memcpy( mph.product, mphPtr+ 0+ 9, 62 ); memcpy( mph.procStage, mphPtr+ 73+11, 1 ); memcpy( mph.refDoc, mphPtr+ 86+ 9, 23 ); memcpy( mph.spare1, mphPtr+ 120+ 0, 40 ); memcpy( mph.acquisitionStation, mphPtr+ 161+21, 20 ); memcpy( mph.procCenter, mphPtr+ 204+13, 6 ); memcpy( mph.procTime, mphPtr+ 225+11, 27 ); memcpy( mph.softwareVer, mphPtr+ 265+14, 14 ); memcpy( mph.spare2, mphPtr+ 295+ 0, 40 ); memcpy( mph.sensingStart, mphPtr+ 336+15, 27 ); memcpy( mph.sensingStop, mphPtr+ 380+14, 27 ); memcpy( mph.spare3, mphPtr+ 423+ 0, 40 ); memcpy( mph.phase, mphPtr+ 464+ 6, 1 ); mph.cycle = atoi( ( char * ) strchr( mphPtr+ 472+ 0, '=' )+1 ); mph.relOrbit = atoi( ( char * ) strchr( mphPtr+ 483+ 0, '=' )+1 ); mph.absOrbit = atoi( ( char * ) strchr( mphPtr+ 500+ 0, '=' )+1 ); memcpy( mph.stateVectorTime, mphPtr+ 517+19, 27 ); mph.deltaUt1 = atof( ( char * ) strchr( mphPtr+ 565+ 0, '=' )+1 ); mph.xPosition = atof( ( char * ) strchr( mphPtr+ 587+ 0, '=' )+1 ); mph.yPosition = atof( ( char * ) strchr( mphPtr+ 614+ 0, '=' )+1 ); mph.zPosition = atof( ( char * ) strchr( mphPtr+ 641+ 0, '=' )+1 ); mph.xVelocity = atof( ( char * ) strchr( mphPtr+ 668+ 0, '=' )+1 ); mph.yVelocity = atof( ( char * ) strchr( mphPtr+ 697+ 0, '=' )+1 ); mph.zVelocity = atof( ( char * ) strchr( mphPtr+ 726+ 0, '=' )+1 ); memcpy( mph.vectorSource, mphPtr+ 755+15, 2 ); memcpy( mph.spare4, mphPtr+ 774+ 0, 40 ); memcpy( mph.utcSbtTime, mphPtr+ 815+14, 27 ); mph.satBinaryTime = atoi( ( char * ) strchr( mphPtr+ 858+ 0, '=' )+1 ); mph.clockStep = atoi( ( char * ) strchr( mphPtr+ 886+ 0, '=' )+1 ); memcpy( mph.spare5, mphPtr+ 913+ 0, 32 ); memcpy( mph.leapUtc, mphPtr+ 946+10, 27 ); mph.leapSign = atoi( ( char * ) strchr( mphPtr+ 985+ 0, '=' )+1 ); mph.leapErr = atoi( ( char * ) strchr( mphPtr+1000+ 0, '=' )+1 ); memcpy( mph.spare6, mphPtr+1011+ 0, 40 ); mph.productErr = atoi( ( char * ) strchr( mphPtr+1052+ 0, '=' )+1 ); mph.totSize = atoi( ( char * ) strchr( mphPtr+1066+ 0, '=' )+1 ); mph.sphSize = atoi( ( char * ) strchr( mphPtr+1104+ 0, '=' )+1 ); mph.numDsd = atoi( ( char * ) strchr( mphPtr+1132+ 0, '=' )+1 ); mph.dsdSize = atoi( ( char * ) strchr( mphPtr+1152+ 0, '=' )+1 ); mph.numDataSets = atoi( ( char * ) strchr( mphPtr+1180+ 0, '=' )+1 ); memcpy( mph.spare7, mphPtr+1206+ 0, 40 ); if ( printMphIfZero == 0 ) { printf( "%s%.62s\n", "product: ", mph.product ); printf( "%s%.1s\n", "procStage: ", mph.procStage ); printf( "%s%.23s\n", "refDoc: ", mph.refDoc ); printf( "%s%.40s\n", "spare1: ", mph.spare1 ); printf( "%s%.20s\n", "acquisitionStation: ", mph.acquisitionStation ); printf( "%s%.6s\n", "procCenter: ", mph.procCenter ); printf( "%s%.27s\n", "procTime: ", mph.procTime ); printf( "%s%.14s\n", "softwareVer: ", mph.softwareVer ); printf( "%s%.40s\n", "spare2: ", mph.spare2 ); printf( "%s%.27s\n", "sensingStart: ", mph.sensingStart ); printf( "%s%.27s\n", "sensingStop: ", mph.sensingStop ); printf( "%s%.40s\n", "spare3: ", mph.spare3 ); printf( "%s%.1s\n", "phase: ", mph.phase ); printf( "%s%d\n", "cycle: ", mph.cycle ); printf( "%s%d\n", "relOrbit: ", mph.relOrbit ); printf( "%s%d\n", "absOrbit: ", mph.absOrbit ); printf( "%s%.27s\n", "stateVectorTime: ", mph.stateVectorTime ); printf( "%s%f\n", "deltaUt1: ", mph.deltaUt1 ); printf( "%s%f\n", "xPosition: ", mph.xPosition ); printf( "%s%f\n", "yPosition: ", mph.yPosition ); printf( "%s%f\n", "zPosition: ", mph.zPosition ); printf( "%s%f\n", "xVelocity: ", mph.xVelocity ); printf( "%s%f\n", "yVelocity: ", mph.yVelocity ); printf( "%s%f\n", "zVelocity: ", mph.zVelocity ); printf( "%s%.2s\n", "vectorSource: ", mph.vectorSource ); printf( "%s%.40s\n", "spare4: ", mph.spare4 ); printf( "%s%.27s\n", "utcSbtTime: ", mph.utcSbtTime ); printf( "%s%u\n", "satBinaryTime: ", mph.satBinaryTime ); printf( "%s%u\n", "clockStep: ", mph.clockStep ); printf( "%s%.32s\n", "spare5: ", mph.spare5 ); printf( "%s%.27s\n", "leapUtc: ", mph.leapUtc ); printf( "%s%d\n", "leapSign: ", mph.leapSign ); printf( "%s%d\n", "leapErr: ", mph.leapErr ); printf( "%s%.40s\n", "spare6: ", mph.spare6 ); printf( "%s%d\n", "productErr: ", mph.productErr ); printf( "%s%d\n", "totSize: ", mph.totSize ); printf( "%s%d\n", "sphSize: ", mph.sphSize ); printf( "%s%d\n", "numDsd: ", mph.numDsd ); printf( "%s%d\n", "dsdSize: ", mph.dsdSize ); printf( "%s%d\n", "numDataSets: ", mph.numDataSets ); printf( "%s%.40s\n", "spare7: ", mph.spare7 ); printf( "\n" ); } return mph; } /* end readMph */ /**********************************************************************************************************************************/ /**********************************************************************************************************************************/ struct sphStruct readSph( const char *sphPtr, const int printSphIfZero, const struct mphStruct mph ) { struct sphStruct sph; int i; memcpy( sph.sphDescriptor, sphPtr+ 0+16, 28 ); sph.startLat = atof( ( char * ) strchr( sphPtr+ 46+ 0, '=' )+1 ) * 1.e-6; sph.startLon = atof( ( char * ) strchr( sphPtr+ 78+ 0, '=' )+1 ) * 1.e-6; sph.stopLat = atof( ( char * ) strchr( sphPtr+111+ 0, '=' )+1 ) * 1.e-6; sph.stopLon = atof( ( char * ) strchr( sphPtr+142+ 0, '=' )+1 ) * 1.e-6; sph.satTrack = atof( ( char * ) strchr( sphPtr+174+ 0, '=' )+1 ); memcpy( sph.spare1, sphPtr+205+ 0, 50 ); sph.ispErrorsSignificant = atoi( ( char * ) strchr( sphPtr+256+ 0, '=' )+1 ); sph.missingIspsSignificant = atoi( ( char * ) strchr( sphPtr+281+ 0, '=' )+1 ); sph.ispDiscardedSignificant = atoi( ( char * ) strchr( sphPtr+308+ 0, '=' )+1 ); sph.rsSignificant = atoi( ( char * ) strchr( sphPtr+336+ 0, '=' )+1 ); memcpy( sph.spare2, sphPtr+353+ 0, 50 ); sph.numErrorIsps = atoi( ( char * ) strchr( sphPtr+404+ 0, '=' )+1 ); sph.errorIspsThresh = atof( ( char * ) strchr( sphPtr+431+ 0, '=' )+1 ); sph.numMissingIsps = atoi( ( char * ) strchr( sphPtr+468+ 0, '=' )+1 ); sph.missingIspsThresh = atof( ( char * ) strchr( sphPtr+497+ 0, '=' )+1 ); sph.numDiscardedIsps = atoi( ( char * ) strchr( sphPtr+536+ 0, '=' )+1 ); sph.discardedIspsThresh = atof( ( char * ) strchr( sphPtr+567+ 0, '=' )+1 ); sph.numRsIsps = atoi( ( char * ) strchr( sphPtr+608+ 0, '=' )+1 ); sph.rsThresh = atof( ( char * ) strchr( sphPtr+632+ 0, '=' )+1 ); memcpy( sph.spare3, sphPtr+661+ 0, 100 ); memcpy( sph.txRxPolar, sphPtr+762+13, 5 ); memcpy( sph.swath, sphPtr+782+ 7, 3 ); memcpy( sph.spare4, sphPtr+794+ 0, 41 ); if ( 1 == 0 ) { printf( "check:\n%s\n", sphPtr+836+ 0 ); } if ( printSphIfZero == 0 ) { printf( "%s%.28s\n", "sphDescriptor: ", sph.sphDescriptor ); printf( "%s%f\n", "startLat: ", sph.startLat ); printf( "%s%f\n", "startLon: ", sph.startLon ); printf( "%s%f\n", "stopLat: ", sph.stopLat ); printf( "%s%f\n", "stopLon: ", sph.stopLon ); printf( "%s%f\n", "satTrack: ", sph.satTrack ); printf( "%s%.50s\n", "spare1: ", sph.spare1 ); printf( "%s%d\n", "ispErrorsSignificant: ", sph.ispErrorsSignificant ); printf( "%s%d\n", "missingIspsSignificant: ", sph.missingIspsSignificant ); printf( "%s%d\n", "ispDiscardedSignificant: ", sph.ispDiscardedSignificant ); printf( "%s%d\n", "rsSignificant: ", sph.rsSignificant ); printf( "%s%.50s\n", "spare2: ", sph.spare2 ); printf( "%s%d\n", "numErrorIsps: ", sph.numErrorIsps ); printf( "%s%f\n", "errorIspsThresh: ", sph.errorIspsThresh ); printf( "%s%d\n", "numMissingIsps: ", sph.numMissingIsps ); printf( "%s%f\n", "missingIspsThresh: ", sph.missingIspsThresh ); printf( "%s%d\n", "numDiscardedIsps: ", sph.numDiscardedIsps ); printf( "%s%f\n", "discardedIspsThresh: ", sph.discardedIspsThresh ); printf( "%s%d\n", "numRsIsps: ", sph.numRsIsps ); printf( "%s%f\n", "rsThresh: ", sph.rsThresh ); printf( "%s%.100s\n", "spare3: ", sph.spare3 ); printf( "%s%.5s\n", "txRxPolar: ", sph.txRxPolar ); printf( "%s%.3s\n", "swath: ", sph.swath ); printf( "%s%.41s\n", "spare4: ", sph.spare4 ); } for ( i = 0; i < mph.numDsd; i++ ){ /* extract DSDs from SPH */ if ( i != 3 ) { /* fourth is a spare DSD - see pdf page 537 */ if (1 == 0) { printf( "check:\n%s\n", sphPtr+836+mph.dsdSize*i+ 0+ 0 ); } memcpy( sph.dsd[ i ].dsName, sphPtr+836+mph.dsdSize*i+ 0+ 9, 28 ); memcpy( sph.dsd[ i ].dsType, sphPtr+836+mph.dsdSize*i+ 39+ 8, 1 ); memcpy( sph.dsd[ i ].filename, sphPtr+836+mph.dsdSize*i+ 49+10, 62 ); sph.dsd[ i ].dsOffset = atoi( ( char * ) strchr( sphPtr+836+mph.dsdSize*i+123+ 0, '=' )+1 ); sph.dsd[ i ].dsSize = atoi( ( char * ) strchr( sphPtr+836+mph.dsdSize*i+162+ 0, '=' )+1 ); sph.dsd[ i ].numDsr = atoi( ( char * ) strchr( sphPtr+836+mph.dsdSize*i+199+ 0, '=' )+1 ); sph.dsd[ i ].dsrSize = atoi( ( char * ) strchr( sphPtr+836+mph.dsdSize*i+219+ 0, '=' )+1 ); /* write out a few things */ if ( printSphIfZero == 0 ) { printf( "%s%d%s%.28s\n", "dsd[ ", i, " ].dsName: ", sph.dsd[ i ].dsName ); printf( "%s%d%s%.1s\n", "dsd[ ", i, " ].dsType: ", sph.dsd[ i ].dsType ); printf( "%s%d%s%.62s\n", "dsd[ ", i, " ].filename: ", sph.dsd[ i ].filename ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].dsOffset: ", sph.dsd[ i ].dsOffset ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].dsSize: ", sph.dsd[ i ].dsSize ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].numDsr: ", sph.dsd[ i ].numDsr ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].dsrSize: ", sph.dsd[ i ].dsrSize ); } } } if ( printSphIfZero == 0 ) { printf( "\n" ); } return sph; } /* end readSph */ /**********************************************************************************************************************************/ /**********************************************************************************************************************************/ struct sphAuxStruct readSphAux( const char *sphPtr, const int printSphIfZero, const struct mphStruct mph ) { struct sphAuxStruct sph; int i; memcpy( sph.sphDescriptor, sphPtr+ 0+16, 28 ); memcpy( sph.spare1, sphPtr+46+ 0, 51 ); if ( printSphIfZero == 0 ) { printf( "%s%.28s\n", "sphDescriptor: ", sph.sphDescriptor ); printf( "%s%.51s\n", "spare1: ", sph.spare1 ); } for ( i = 0; i < mph.numDsd; i++ ){ /* extract DSDs from SPH */ memcpy( sph.dsd[ i ].dsName, sphPtr+ 98+mph.dsdSize*i+ 0+ 9, 28 ); memcpy( sph.dsd[ i ].dsType, sphPtr+ 98+mph.dsdSize*i+ 39+ 8, 1 ); memcpy( sph.dsd[ i ].filename, sphPtr+ 98+mph.dsdSize*i+ 49+10, 62 ); sph.dsd[ i ].dsOffset = atoi( ( char * ) strchr( sphPtr+ 98+mph.dsdSize*i+123+ 0, '=' )+1 ); sph.dsd[ i ].dsSize = atoi( ( char * ) strchr( sphPtr+ 98+mph.dsdSize*i+162+ 0, '=' )+1 ); sph.dsd[ i ].numDsr = atoi( ( char * ) strchr( sphPtr+ 98+mph.dsdSize*i+199+ 0, '=' )+1 ); sph.dsd[ i ].dsrSize = atoi( ( char * ) strchr( sphPtr+ 98+mph.dsdSize*i+219+ 0, '=' )+1 ); /* write out a few things */ if ( printSphIfZero == 0 ) { printf( "%s%d%s%.28s\n", "dsd[ ", i, " ].dsName: ", sph.dsd[ i ].dsName ); printf( "%s%d%s%.1s\n", "dsd[ ", i, " ].dsType: ", sph.dsd[ i ].dsType ); printf( "%s%d%s%.62s\n", "dsd[ ", i, " ].filename: ", sph.dsd[ i ].filename ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].dsOffset: ", sph.dsd[ i ].dsOffset ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].dsSize: ", sph.dsd[ i ].dsSize ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].numDsr: ", sph.dsd[ i ].numDsr ); printf( "%s%d%s%d\n", "dsd[ ", i, " ].dsrSize: ", sph.dsd[ i ].dsrSize ); } } if ( printSphIfZero == 0 ) { printf( "\n" ); } return sph; } /* end readSphAux */ /**********************************************************************************************************************************/ /**********************************************************************************************************************************/ void printInsGads( const struct insGadsStruct insGads ) { int i; printf( "%s%d\n", "dsrTime.days: ", insGads.dsrTime.days ); printf( "%s%d\n", "dsrTime.seconds: ", insGads.dsrTime.seconds ); printf( "%s%d\n", "dsrTime.microseconds: ", insGads.dsrTime.microseconds ); printf( "%s%d\n", "dsrLength: ", insGads.dsrLength ); printf( "%s%.9g\n", "radarFrequency: ", insGads.radarFrequency ); printf( "%s%.9g\n", "sampRate: ", insGads.sampRate ); printf( "%s%.9g\n", "offsetFreq: ", insGads.offsetFreq ); printf( "%s%.9g\n", "rangeGateBias: ", insGads.rangeGateBias ); printf( "%s%.9g\n", "rangeGateBiasGm: ", insGads.rangeGateBiasGm ); printf( "%s%f\n", "refElevAngleIs1: ", insGads.refElevAngleIs1 ); printf( "%s%f\n", "refElevAngleIs2: ", insGads.refElevAngleIs2 ); printf( "%s%f\n", "refElevAngleIs3Ss2: ", insGads.refElevAngleIs3Ss2 ); printf( "%s%f\n", "refElevAngleIs4Ss3: ", insGads.refElevAngleIs4Ss3 ); printf( "%s%f\n", "refElevAngleIs5Ss4: ", insGads.refElevAngleIs5Ss4 ); printf( "%s%f\n", "refElevAngleIs6Ss5: ", insGads.refElevAngleIs6Ss5 ); printf( "%s%f\n", "refElevAngleIs7: ", insGads.refElevAngleIs7 ); printf( "%s%f\n", "refElevAngleSs1: ", insGads.refElevAngleSs1 ); printf( "%s%.9g\n", "swstCalP2: ", insGads.swstCalP2 ); printf( "%s%u\n", "perCalWindowsEc: ", insGads.perCalWindowsEc ); printf( "%s%u\n", "perCalWindowsMs: ", insGads.perCalWindowsMs ); printf( "%s%u\n", "initCalBeamSetWv: ", insGads.initCalBeamSetWv ); printf( "%s%u\n", "beamSetEc: ", insGads.beamSetEc ); printf( "%s%u\n", "beamSetMs: ", insGads.beamSetMs ); printf( "%s%u\n", "mEc: ", insGads.mEc ); printf( ".\n" ); printf( ".\n" ); printf( ".\n" ); for ( i = 0; i < 4096; i++ ) printf( "%s%4d%s%15f %15f %15f\n", "fbaq4LutI,Q,NoAdc[ ", i, " ]: ", insGads.fbaq4LutI[ i ], insGads.fbaq4LutQ[ i ], insGads.fbaq4NoAdc[ i ] ); printf( ".\n" ); printf( ".\n" ); printf( ".\n" ); printf( "\n" ); /* exit( 0 ); */ return; } /* end printInsGads */ /**********************************************************************************************************************************/ /**********************************************************************************************************************************/ /********************************************************** ** Function: byte_swap_InsGads ** ** Purpose: Convert the bytes of struct insGadsStruct for a little endian order machine ** ** Comment: struct testStruct should be redefined in the future! ** ** Author: Zhenhong Li at UCL ** ** Created: 17/02/2005 ** ** Modified: ** ;**********************************************************/ void byte_swap_InsGads( struct insGadsStruct* InsGads ) { swap_endian_order(e_tid_long, &(*InsGads).dsrTime, 3); swap_endian_order(e_tid_ulong, &(*InsGads).dsrLength, 1); swap_endian_order(e_tid_float, &(*InsGads).radarFrequency, 1); swap_endian_order(e_tid_float, &(*InsGads).sampRate, 1); swap_endian_order(e_tid_float, &(*InsGads).offsetFreq, 1); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0TxH1, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0TxV1, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0TxH1a, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0TxV1a, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0RxH2, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0RxV2, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0H3, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseIm0V3, 64); swap_endian_order(e_tid_float, &(*InsGads).calPulseImTxH1, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImTxV1, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImTxH1a, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImTxV1a, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImRxH2, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImRxV2, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImH3, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseImV3, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApTxH1, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApTxV1, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApTxH1a, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApTxV1a, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApRxH2, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApRxV2, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApH3, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseApV3, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvTxH1, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvTxV1, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvTxH1a, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvTxV1a, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvRxH2, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvRxV2, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvH3, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWvV3, 448); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsTxH1, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsTxV1, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsTxH1a, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsTxV1a, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsRxH2, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsRxV2, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsH3, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseWsV3, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmTxH1, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmTxV1, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmTxH1a, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmTxV1a, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmRxH2, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmRxV2, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmH3, 320); swap_endian_order(e_tid_float, &(*InsGads).calPulseGmV3, 320); swap_endian_order(e_tid_float, &(*InsGads).nomPulseIm, 63); swap_endian_order(e_tid_float, &(*InsGads).nomPulseAp, 63); swap_endian_order(e_tid_float, &(*InsGads).nomPulseWv, 63); swap_endian_order(e_tid_float, &(*InsGads).nomPulseWs, 45); swap_endian_order(e_tid_float, &(*InsGads).nomPulseGm, 45); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs1, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs2, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs3Ss2, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs4Ss3, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs5Ss4, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs6Ss5, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternIs7, 101); swap_endian_order(e_tid_float, &(*InsGads).azPatternSs1, 101); swap_endian_order(e_tid_float, &(*InsGads).rangeGateBias, 1); swap_endian_order(e_tid_float, &(*InsGads).rangeGateBiasGm, 1); swap_endian_order(e_tid_float, &(*InsGads).adcLutI, 255); swap_endian_order(e_tid_float, &(*InsGads).adcLutQ, 255); swap_endian_order(e_tid_float, &(*InsGads).full8LutI, 256); swap_endian_order(e_tid_float, &(*InsGads).full8LutQ, 256); swap_endian_order(e_tid_float, &(*InsGads).fbaq4LutI, 4096); swap_endian_order(e_tid_float, &(*InsGads).fbaq3LutI, 2048); swap_endian_order(e_tid_float, &(*InsGads).fbaq2LutI, 1024); swap_endian_order(e_tid_float, &(*InsGads).fbaq4LutQ, 4096); swap_endian_order(e_tid_float, &(*InsGads).fbaq3LutQ, 2048); swap_endian_order(e_tid_float, &(*InsGads).fbaq2LutQ, 1024); swap_endian_order(e_tid_float, &(*InsGads).fbaq4NoAdc, 4096); swap_endian_order(e_tid_float, &(*InsGads).fbaq3NoAdc, 2048); swap_endian_order(e_tid_float, &(*InsGads).fbaq2NoAdc, 1024); swap_endian_order(e_tid_float, &(*InsGads).smLutI, 16); swap_endian_order(e_tid_float, &(*InsGads).smLutQ, 16); swap_endian_order(e_tid_ushort, &(*InsGads).swathConfigIm, 28); swap_endian_order(e_tid_float, &(*InsGads).swathConfigIm.resampleFactor, 7); swap_endian_order(e_tid_ushort, &(*InsGads).swathConfigAp, 28); swap_endian_order(e_tid_float, &(*InsGads).swathConfigAp.resampleFactor, 7); swap_endian_order(e_tid_ushort, &(*InsGads).swathConfigWs, 28); swap_endian_order(e_tid_float, &(*InsGads).swathConfigWs.resampleFactor, 7); swap_endian_order(e_tid_ushort, &(*InsGads).swathConfigGm, 28); swap_endian_order(e_tid_float, &(*InsGads).swathConfigGm.resampleFactor, 7); swap_endian_order(e_tid_ushort, &(*InsGads).swathConfigWv, 28); swap_endian_order(e_tid_float, &(*InsGads).swathConfigWv.resampleFactor, 7); swap_endian_order(e_tid_ushort, &(*InsGads).perCalWindowsEc, 1); swap_endian_order(e_tid_ushort, &(*InsGads).perCalWindowsMs, 1); swap_endian_order(e_tid_ushort, &(*InsGads).swathIdIm, 14); swap_endian_order(e_tid_ushort, &(*InsGads).swathIdAp, 14); swap_endian_order(e_tid_ushort, &(*InsGads).swathIdWs, 14); swap_endian_order(e_tid_ushort, &(*InsGads).swathIdGm, 14); swap_endian_order(e_tid_ushort, &(*InsGads).swathIdWv, 14); swap_endian_order(e_tid_ushort, &(*InsGads).initCalBeamSetWv, 1); swap_endian_order(e_tid_ushort, &(*InsGads).beamSetEc, 1); swap_endian_order(e_tid_ushort, &(*InsGads).beamSetMs, 1); swap_endian_order(e_tid_ushort, &(*InsGads).calSeq, 32); swap_endian_order(e_tid_ushort, &(*InsGads).timelineIm, 28); swap_endian_order(e_tid_ushort, &(*InsGads).timelineAp, 28); swap_endian_order(e_tid_ushort, &(*InsGads).timelineWs, 28); swap_endian_order(e_tid_ushort, &(*InsGads).timelineGm, 28); swap_endian_order(e_tid_ushort, &(*InsGads).timelineWv, 28); swap_endian_order(e_tid_ushort, &(*InsGads).mEc, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs1, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs2, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs3Ss2, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs4Ss3, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs5Ss4, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs6Ss5, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleIs7, 1); swap_endian_order(e_tid_float, &(*InsGads).refElevAngleSs1, 1); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs1, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs2, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs3Ss2, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs4Ss3, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs5Ss4, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs6Ss5, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefIs7, 128); swap_endian_order(e_tid_float, &(*InsGads).calLoopRefSs1, 128); //struct testStruct should be redefined in the future. swap_endian_order(e_tid_float, &(*InsGads).im, 17); swap_endian_order(e_tid_float, &(*InsGads).ap, 17); swap_endian_order(e_tid_float, &(*InsGads).ws, 17); swap_endian_order(e_tid_float, &(*InsGads).gm, 17); swap_endian_order(e_tid_float, &(*InsGads).wv, 17); swap_endian_order(e_tid_float, &(*InsGads).swstCalP2, 1); } /********************************************************** ** Function: is_bigendian ** ** Purpose: Test whether it is a bigendian machine ** ** Return values: true: 1, false: 0 ** ** Comment: ** ** Author: Eric J Fielding at JPL ** ** Created: ** ** Modified: ** ;**********************************************************/ int is_bigendian() { int bigendian, littleendian, test; unsigned char t[4]; littleendian=256; bigendian=256*256; t[0]=0; t[1]=1; t[2]=0; t[3]=0; memcpy(&test, &t[0], 4); /* printf("test: %i\n",test); */ if(test==bigendian)return(1); if(test==littleendian)return(0); printf("Error in endian test, test= %i ********\n",test); } /* * Function: byte_swap_short.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS two-byte words, * starting at address BUFFER. * * @param buffer the one element typed buffer * to convert for a little endian order machine * * @param number_of_swaps number of elements to convert * */ void byte_swap_short(short *buffer, uint number_of_swaps) { short* temp = buffer; uint swap_loop; for (swap_loop = 0, temp = buffer; swap_loop < number_of_swaps; swap_loop++, temp++) { *temp = (short)(((*temp & 0x00ff) << 8) | ((*temp & 0xff00) >> 8)); } } /* Function: byte_swap_long.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS four-byte words, * starting at address BUFFER. * * */ void byte_swap_long(long *buffer, uint number_of_swaps) { long *temp = buffer; uint swap_loop; for (swap_loop = 0, temp = buffer; swap_loop < number_of_swaps; swap_loop++, temp++) { *temp = ((*temp & 0x000000ff) << 24) | ((*temp & 0x0000ff00) << 8) | ((*temp & 0x00ff0000) >> 8) | ((*temp & 0xff000000) >> 24); } } /* ADDED THESE LINES TO TEST THE 4-BYTE INT TYPE ON 64 BIT */ /* Function: byte_swap_int.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS four-byte words, * starting at address BUFFER. * * */ void byte_swap_int(int *buffer, uint number_of_swaps) { int *temp = buffer; uint swap_loop; for (swap_loop = 0, temp = buffer; swap_loop < number_of_swaps; swap_loop++, temp++) { *temp = ((*temp & 0x000000ff) << 24) | ((*temp & 0x0000ff00) << 8) | ((*temp & 0x00ff0000) >> 8) | ((*temp & 0xff000000) >> 24); } } /* Function: byte_swap_uint.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS four-byte words, * starting at address BUFFER. * * */ void byte_swap_uint(uint *buffer, uint number_of_swaps) { uint *temp = buffer; uint swap_loop; for (swap_loop = 0, temp = buffer; swap_loop < number_of_swaps; swap_loop++, temp++) { *temp = ((*temp & 0x000000ff) << 24) | ((*temp & 0x0000ff00) << 8) | ((*temp & 0x00ff0000) >> 8) | ((*temp & 0xff000000) >> 24); } } /* ADDDED NEW LINES ABOVE */ /* ************************************************************************** */ /* Function: byte_swap_short.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS two-byte words, * starting at address BUFFER. * * @param buffer the one element typed buffer * to convert for a little endian order machine * * @param number_of_swaps number of elements to convert * */ void byte_swap_ushort(ushort* buffer, uint number_of_swaps) { byte_swap_short((short*) buffer, number_of_swaps); } /* * Function: byte_swap_ulong.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS four-byte words, * starting at address BUFFER. * * @param buffer the one element typed buffer * to convert for a little endian order machine * * @param number_of_swaps number of elements to convert * */ void byte_swap_ulong(ulong* buffer, uint number_of_swaps) { byte_swap_long((long*) buffer, number_of_swaps); } /* * Function: byte_swap_long.c */ /** * * Swaps bytes within NUMBER_OF_SWAPS four-byte words, * starting at address BUFFER. * * @param buffer the one element typed buffer * to convert for a little endian order machine * * @param number_of_swaps number of elements to convert * */ void byte_swap_float(float* buffer, uint number_of_swaps) { byte_swap_int((int*) buffer, number_of_swaps); } /* Function: epr_swap_endian_order Access: public API Changelog: 2002/02/04 mp nitial version */ /** * Converts bytes for a little endian order machine * * @param field the pointer at data reading in * */ void swap_endian_order(EPR_EDataTypeId data_type_id, void* elems, uint num_elems) { switch (data_type_id) { case e_tid_uchar: case e_tid_char: case e_tid_string: /* no conversion required */ break; case e_tid_time: byte_swap_uint((uint*)elems, 3); break; case e_tid_spare: /* no conversion required */ break; case e_tid_ushort: byte_swap_ushort((ushort*) elems, num_elems); break; case e_tid_short: byte_swap_short((short*) elems, num_elems); break; case e_tid_ulong: byte_swap_uint((uint*) elems, num_elems); break; case e_tid_long: byte_swap_int((int*) elems, num_elems); break; case e_tid_float: byte_swap_float((float*) elems, num_elems); break; case e_tid_double: printf( "swap_endian_order: DOUBLE type was not yet processed\n" ); break; default: printf( "swap_endian_order: unknown data type\n" ); } } int max_int(int a, int b){ /* find maximum of two floats */ int max; if( a > b) {max = a;} else {max = b;} return(max); } int min_int(int a, int b){ /* find minimum of two floats */ int min; if( a > b) {min = b;} else {min = a;} return(min); } double norm(double a[]){ /* norm of vector */ double mag; mag = sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]); return(mag); } double dot(double a[], double b[]){ /* dot product */ double sum; sum=a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; return(sum); } void cross(double a[], double b[], double c[]){ /* cross product */ c[0]=a[1]*b[2]-a[2]*b[1]; c[1]=a[2]*b[0]-a[0]*b[2]; c[2]=a[0]*b[1]-a[1]*b[0]; } void svdiv(double sm, double a[], double b[]) { /* vector scalar divide b = a / sm */ b[0] = a[0] / sm; b[1] = a[1] / sm; b[2] = a[2] / sm; } void svmul(double sm, double a[], double b[]) { /* vector scalar multiply b = a * sm */ b[0] = a[0] * sm; b[1] = a[1] * sm; b[2] = a[2] * sm; } void vadd(double a[], double b[], double c[]) { /* vector c = a + b */ c[0] = a[0] + b[0]; c[1] = a[1] + b[1]; c[2] = a[2] + b[2]; } void vsub(double a[], double b[], double c[]) { /* vector c = a - b */ c[0] = a[0] - b[0]; c[1] = a[1] - b[1]; c[2] = a[2] - b[2]; } void matvmul(double a[][3], double b[], double c[]){ /* matrix and vector multiplication c=a*b a:3x3 b,c:3x1 */ int i; for(i=0;i<3;i++){ c[i]=a[i][0]*b[0]+a[i][1]*b[1]+a[i][2]*b[2]; } } void rot1(double ang, double c[][3]){ /* calculation of rotation matrix about x-axis */ c[0][0]=1; c[0][1]=0; c[0][2]=0; c[1][0]=0; c[1][1]=cos(ang); c[1][2]=sin(ang); c[2][0]=0; c[2][1]=-sin(ang); c[2][2]=cos(ang); } void rot2(double ang, double c[][3]){ /* calculation of rotation matrix about y-axis */ c[0][0]=cos(ang); c[0][1]=0; c[0][2]=-sin(ang); c[1][0]=0; c[1][1]=1; c[1][2]=0; c[2][0]=sin(ang); c[2][1]=0; c[2][2]=cos(ang); } void rot3(double ang, double c[][3]){ /* calculation of rotation matrix about z-axis */ c[0][0]=cos(ang); c[0][1]=sin(ang); c[0][2]=0; c[1][0]=-sin(ang); c[1][1]=cos(ang); c[1][2]=0; c[2][0]=0; c[2][1]=0; c[2][2]=1; } double MJD2GST(double MJD, double dt){ /* dt is the seconds since the epoch time GST is the grenwich mean sidereal time Returns GST is in radians */ double JD0, we, T, GST0, GST; JD0=MJD+2400000.5; we = EROTRAD * 180/PI; /* deg/sec */ T = (JD0-2451545.0)/36525.0; GST0 = 280.46061837 + 360.98564736629*(JD0 - 2451545.0) + 0.000387933*(T*T) -(T*T*T)/38710000.0; GST = GST0 + we*dt; /* degrees */ while(GST > 360){ /* JD0 is the epoch time */ GST= GST - 360.0; } while(GST < 0){ GST = GST +360.0; } GST=GST*PI/180.0; return GST; } void ECEF2ECI(int yr, int mo, int day, int hr, int min, double sec, double dt, double Secef[] , double Seci[]){ /* convert state vector (6x1) in ECEF to state vector (6x1) in ECI at date + dt (sec)*/ int i; double MJD, GST; double Recef[3]={0}; double Vecef[3]={0}; double Reci[3]={0}; double Veci[3]={0}; double Veci_int[3]={0}; double Vecef_int[3]={0}; double w_vec[3]={0,0, EROTRAD}; double wCrossReci[3]={0}; double wCrossRecef[3]={0}; double rot3GST[3][3]={0}; // printf("\nEarth rotation vector w_vec = %10.8e %10.8e %10.8e \n", w_vec[0], w_vec[1], w_vec[2]); MJD = date2MJD(yr,mo,day,hr,min,sec); GST = MJD2GST(MJD,dt); for (i=0;i<3;i++){ Recef[i]=Secef[i]; Vecef[i]=Secef[i+3]; } rot3(-GST,rot3GST); // matvmul(rot3GST,Recef,Reci); // matvmul(rot3GST,Vecef,Veci_int); // cross(w_vec,Reci,wCrossReci); // vadd(Veci_int, wCrossReci, Veci); matvmul(rot3GST,Recef,Reci); cross(w_vec,Recef,wCrossRecef); vadd(Vecef, wCrossRecef, Vecef_int); matvmul(rot3GST,Vecef_int,Veci); for (i=0;i<3;i++){ Seci[i]=Reci[i]; Seci[i+3]=Veci[i]; } } void ECI2ECEF(int yr, int mo, int day, int hr, int min, double sec, double dt, double Seci[], double Secef[]){ /* convert state vector (6x1) in ECI to state vector (6x1) in ECEF at date + dt(sec)*/ int i; double MJD, GST; double Reci[3]={0}; double Veci[3]={0}; double Recef[3]={0}; double Vecef[3]={0}; double Vecef_int[3]={0}; double w_vec[3]={0,0, EROTRAD}; double wCrossRecef[3]={0}; double rot3GST[3][3]={0}; // printf("\nEarth rotation vector w_vec = %10.8f %10.8e %10.8e \n", w_vec[0], w_vec[1], w_vec[2]); MJD = date2MJD(yr,mo,day,hr,min,sec); GST = MJD2GST(MJD,dt); for (i=0;i<3;i++){ Reci[i]=Seci[i]; Veci[i]=Seci[i+3]; } rot3(GST,rot3GST); matvmul(rot3GST,Reci,Recef); matvmul(rot3GST,Veci,Vecef_int); cross(w_vec,Recef,wCrossRecef); vsub(Vecef_int, wCrossRecef, Vecef); for (i=0;i<3;i++){ Secef[i]=Recef[i]; Secef[i+3]=Vecef[i]; } } double date2MJD(int yr, int mo, int day, int hr, int min, double sec){ /* convert to date to MJD */ double part1,part2; double JD,MJD; part1 = 367.*(yr) - floor(7*(yr+floor((mo+9)/12.0))/4.0) + floor(275*mo/9.0) + day; part2 = 1721013.5 + ((sec/60.0+min)/60.0+hr)/24.0; JD = part1 + part2; MJD = JD - 2400000.5; return MJD; } void xyz2LatLonAlt(double rvec[] , double a, double b, double *phi, double *lam, double *h) { int i; double r; double e,ep; double E; double F; double G; double c; double s; double P; double Q; double r0; double U; double V; double z0; double xu, yu,zu; xu = rvec[0]; yu=rvec[1]; zu=rvec[2]; /* printf("\nECEF coorindates (x,y,z): %14.4f %14.4f %14.4f meters", xu, yu, zu); printf("\nsemi-major axis for reference ellipsoid: %14.4f meters", a); printf("\nsemi-minor axis for reference ellipsoid: %14.4f meters\n", b); */ e = sqrt(1- pow(b,2.0)/pow(a,2.0)); ep = (a/b) * e; r = sqrt(pow(xu,2.0) + pow(yu,2.0)); E = sqrt(pow(a,2.0) - pow(b,2.0)); F = 54*pow(b,2.0)*pow(zu,2); G = pow(r,2.0) + (1 - pow(e,2.0))*pow(zu,2.0) - pow(e,2.0)*pow(E,2.0); c = pow(e,4.0)* F * pow(r,2.0) / pow(G,3.0); s = cbrt( (1 + c + sqrt(pow(c,2.0) + 2*c)) ); P = F/ (3* pow((s + 1/s + 1),2.0) * pow(G,2.0)); Q = sqrt(1 + 2 * pow(e,4.0) * P); r0 = -(P * pow(e,2) * r/ (1 + Q)) + sqrt( 0.5*pow(a,2.0)*(1+1/Q) - P* (1-pow(e,2.0))*pow(zu,2.0)/(Q*(1+Q)) - 0.5*P*pow(r,2.0) ); U = sqrt( pow((r - pow(e,2.0)*r0),2.0) + pow(zu,2.0) ); V = sqrt( pow((r - pow(e,2.0)*r0),2.0) + (1-pow(e,2.0)) * pow(zu,2.0) ); z0 = pow(b, 2.0) * zu/(a * V); *h = U* (1 - pow(b,2.0)/(a*V)); *phi = atan( (zu + pow(ep,2.0)*z0)/r ) *RTD; *lam = atan2(yu, xu) * RTD; // printf("\nlat = %14.6f , lon = %14.6f , alt = %14.6f \n", *phi, *lam, *h); } void geolocate_radarLOS_ellipsoidWGS84(double rv[], double el, double rho, int DIR_FLAG, double ref_alt, double *lat, double *lon) { int iter; double a,b,e,ep; double rs[3], vs[3], h[3], hvec[3], rhovec[3]; double hhat[3], rshat[3], vshat[3]; double rvec[3]; double deltael; double lati, loni, eli, alt, slat, slon, salt; //printf("\nSatellite position vector: %14.4f %14.4f %14.4f meters", rs[0], rs[1], rs[2]); //printf("\nSatellite velocity vector: %14.4f %14.4f %14.4f meters/sec",vs[0], vs[1], vs[2]); //printf("\nRadar Antenna Beam Nominal elevation angle: %12.6f deg.", el); //printf("\nReference elevation: %10.2f meters", ref_alt); rs[0] = rv[0]; rs[1] = rv[1]; rs[2] = rv[2]; vs[0] = rv[3]; vs[1] = rv[4]; vs[2] = rv[5]; a = SMAJORA; b = SMINORA; eli = el*DTR; // printf("\nSatellite position vector: %14.4f %14.4f %14.4f meters", rs[0], rs[1], rs[2]); // printf("\nSatellite velocity vector: %14.4f %14.4f %14.4f meters/sec",vs[0], vs[1], vs[2]); // printf("\nRadar Antenna Beam Nominal elevation angle: %12.6f deg.", el); // printf("\nSlant range to ground point: %14.4f meters", rho); // xyz2LatLonAlt(rs, a, b, &slat, &slon, &salt); // printf("Satellite Position in geodetoc cordinates: \n"); // printf("Satellite Lat: %12.6f\n",slat); // printf("Satellite Lon: %12.6f\n",slon); // printf("Satellite altitude: %12.6f\n",salt); svdiv(norm(rs), rs, rshat); // printf("Unit satellite position vector: %14.6f %14.6f %14.6f\n", rshat[0], rshat[1], rshat[2]); svdiv(norm(vs), vs, vshat); // printf("Unit satellite velocity vector: %14.6f %14.6f %14.6f\n", vshat[0], vshat[1], vshat[2]); cross(rs,vs,h); svmul(DIR_FLAG,h,hvec); svdiv(norm(hvec), hvec,hhat); // printf("Unit cross track vector: %14.6f %14.6f %14.6f\n", hhat[0], hhat[1], hhat[2]); iter=0; do{ rhovec[0] = rho* cos(eli) * hhat[0] + rho * sin(eli) * rshat[0]; rhovec[1] = rho* cos(eli) * hhat[1] + rho * sin(eli) * rshat[1]; rhovec[2] = rho* cos(eli) * hhat[2] + rho * sin(eli) * rshat[2]; // printf("\nrho vector: %14.4f %14.4f %14.4f meters", rhovec[0], rhovec[1], rhovec[2]); rvec[0] = rs[0] + rhovec[0]; rvec[1] = rs[1] + rhovec[1]; rvec[2] = rs[2] + rhovec[2]; xyz2LatLonAlt(rvec, a, b, &lati, &loni, &alt); /* printf("lat = %12.6f\n", lati); printf("lon = %12.6f\n", loni); printf("alt = %12.6f\n", alt); */ deltael = (alt-ref_alt)/rho; eli = eli - deltael; // eli = eli - deltael; Elevation measured with respect to local horizontal place, upwards position // printf("altitude = %14.6f\n", alt); // printf("look angle = %14.6f\n", eli*RTD); iter = iter + 1; }while ( (fabs((alt-ref_alt)) > 0.0001) && (iter < 200)); *lat = lati; *lon = loni; if(iter > 199){ printf("\nERROR: cannot find intersection of radar-LOS and Earth-Ellipsoid to required precision after 100 iterations\n"); exit(-1); } // printf("\n Geodetic Co-ordinates: Latitude: %10.6f Deg. Longitude: %10.6f Deg. Altitude: %10.3f meters\n", lat, lon, alt); } void derivJ2(double qVEC[], double qdotVEC[]){ /* calculate time derivatives in cartesian orbit propagation using J2*/ int i; double dist; qdotVEC[0]=qVEC[3]; qdotVEC[1]=qVEC[4]; qdotVEC[2]=qVEC[5]; dist=sqrt(qVEC[0]*qVEC[0]+qVEC[1]*qVEC[1]+qVEC[2]*qVEC[2]); qdotVEC[3]=-muE*qVEC[0]/pow(dist,3) - (3*J2*muE*Re*Re*qVEC[0]/(2*pow(dist,5)))*(1-5*qVEC[2]*qVEC[2]/(dist*dist)); qdotVEC[4]=-muE*qVEC[1]/pow(dist,3) - (3*J2*muE*Re*Re*qVEC[1]/(2*pow(dist,5)))*(1-5*qVEC[2]*qVEC[2]/(dist*dist)); qdotVEC[5]=-muE*qVEC[2]/pow(dist,3) - (3*J2*muE*Re*Re*qVEC[2]/(2*pow(dist,5)))*(3-5*qVEC[2]*qVEC[2]/(dist*dist)); } void RK4(double RVsat0[] , double delt , double RVsatf[] ) { /* Runge-Kutta 4th order integration */ int iter, i; int nstep; double h; double temp1[6] ={0}; double temp2[6] ={0}; double temp3[6] ={0}; double temp4[6] ={0}; double y1dot[6] ={0}; double y2dot[6] ={0}; double y3dot[6] ={0}; double y4dot[6] ={0}; nstep = 0; do{ nstep = nstep +1; h = (delt/(double)nstep); }while(fabs(h) > STEP); // printf("step size: %8.6f seconds, nstep= %d",h, nstep); for (i=0;i<6;i++){temp4[i] = RVsat0[i];} for (iter=1; iter<=nstep; iter++) { derivJ2(temp4, y1dot); for (i=0;i<6;i++){ temp1[i]=y1dot[i]* h/2.0 + temp4[i]; } derivJ2(temp1,y2dot); for (i=0;i<6;i++){ temp2[i]=y2dot[i]* h/2.0 + temp4[i]; } derivJ2(temp2,y3dot); for (i=0;i<6;i++){ temp3[i]=y3dot[i]* h + temp4[i]; } derivJ2(temp3,y4dot); for (i=0;i<6;i++){ temp4[i] = temp4[i] + h/6.0*(y1dot[i]+2*y2dot[i]+2*y3dot[i]+y4dot[i]); } } for (i=0;i<6;i++){RVsatf[i] = temp4[i];} } void orbitPropogation2body_J2Cartesian(int year, int mon, int day, int hour, int min, double sec, double delT, double sVEC[], double sVECf[]) { int i; double sECI[6]={0}; /* state vector in ECI coordiantes on epoch time */ double sECIf[6] = {0}; /* state vector in ECI coordiantes after delT */ double a=0.0; /* semi-major axis */ double e=0.0; /* eccentricity */ double inc=0.0; /* inclination */ double RAAN=0.0; double w=0.0; double true=0.0; double E=0.0; double M=0.0; double n; /* mean motion */ double MJD; /* modified julian dates */ double GST; /* greenwich mean sidereal time */ /* Initial state vector in ECEF */ /* printf("\n\nOriginal State Vector in ECEF:\n\n"); for(