diff --git a/src/algorithms/libs/rtklib/rtklib_pntpos.cc b/src/algorithms/libs/rtklib/rtklib_pntpos.cc index a8d020b69..793216702 100644 --- a/src/algorithms/libs/rtklib/rtklib_pntpos.cc +++ b/src/algorithms/libs/rtklib/rtklib_pntpos.cc @@ -182,7 +182,7 @@ int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos, if (ionoopt == IONOOPT_BRDC) { *ion = ionmodel(time, nav->ion_gps, pos, azel); - *var = SQR(*ion*ERR_BRDCI); + *var = std::pow(*ion * ERR_BRDCI, 2.0); return 1; } /* sbas ionosphere model */ @@ -234,7 +234,7 @@ int tropcorr(gtime_t time, const nav_t *nav, const double *pos, if (tropopt == TROPOPT_SAAS || tropopt == TROPOPT_EST || tropopt == TROPOPT_ESTG) { *trp = tropmodel(time, pos, azel, REL_HUMI); - *var = SQR(ERR_SAAS / (sin(azel[1]) + 0.1)); + *var = std::pow(ERR_SAAS / (sin(azel[1]) + 0.1), 2.0); return 1; } /* sbas troposphere model */ diff --git a/src/algorithms/libs/rtklib/rtklib_rtkcmn.cc b/src/algorithms/libs/rtklib/rtklib_rtkcmn.cc index 0d741443b..0119630d6 100644 --- a/src/algorithms/libs/rtklib/rtklib_rtkcmn.cc +++ b/src/algorithms/libs/rtklib/rtklib_rtkcmn.cc @@ -48,130 +48,8 @@ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * -* -* options : -DLAPACK use LAPACK/BLAS -* -DMKL use Intel MKL -* -DTRACE enable debug trace -* -DWIN32 use WIN32 API -* -DNOCALLOC no use calloc for zero matrix -* -DIERS_MODEL use GMF instead of NMF -* -DDLL built for shared library -* -DCPUTIME_IN_GPST cputime operated in gpst -* -* references : -* [1] IS-GPS-200D, Navstar GPS Space Segment/Navigation User Interfaces, -* 7 March, 2006 -* [2] RTCA/DO-229C, Minimum operational performanc standards for global -* positioning system/wide area augmentation system airborne equipment, -* RTCA inc, November 28, 2001 -* [3] M.Rothacher, R.Schmid, ANTEX: The Antenna Exchange Format Version 1.4, -* 15 September, 2010 -* [4] A.Gelb ed., Applied Optimal Estimation, The M.I.T Press, 1974 -* [5] A.E.Niell, Global mapping functions for the atmosphere delay at radio -* wavelengths, Jounal of geophysical research, 1996 -* [6] W.Gurtner and L.Estey, RINEX The Receiver Independent Exchange Format -* Version 3.00, November 28, 2007 -* [7] J.Kouba, A Guide to using International GNSS Service (IGS) products, -* May 2009 -* [8] China Satellite Navigation Office, BeiDou navigation satellite system -* signal in space interface control document, open service signal B1I -* (version 1.0), Dec 2012 -* [9] J.Boehm, A.Niell, P.Tregoning and H.Shuh, Global Mapping Function -* (GMF): A new empirical mapping function base on numerical weather -* model data, Geophysical Research Letters, 33, L07304, 2006 -* [10] GLONASS/GPS/Galileo/Compass/SBAS NV08C receiver series BINR interface -* protocol specification ver.1.3, August, 2012 -* -* version : $Revision: 1.1 $ $Date: 2008/07/17 21:48:06 $ -* history : 2007/01/12 1.0 new -* 2007/03/06 1.1 input initial rover pos of pntpos() -* update only effective states of filter() -* fix bug of atan2() domain error -* 2007/04/11 1.2 add function antmodel() -* add gdop mask for pntpos() -* change constant MAXDTOE value -* 2007/05/25 1.3 add function execcmd(),expandpath() -* 2008/06/21 1.4 add funciton sortobs(),uniqeph(),screent() -* replace geodist() by sagnac correction way -* 2008/10/29 1.5 fix bug of ionosphereic mapping function -* fix bug of seasonal variation term of tropmapf -* 2008/12/27 1.6 add function tickget(), sleepms(), tracenav(), -* xyz2enu(), satposv(), pntvel(), covecef() -* 2009/03/12 1.7 fix bug on error-stop when localtime() returns NULL -* 2009/03/13 1.8 fix bug on time adjustment for summer time -* 2009/04/10 1.9 add function adjgpsweek(),getbits(),getbitu() -* add function geph2pos() -* 2009/06/08 1.10 add function seph2pos() -* 2009/11/28 1.11 change function pntpos() -* add function tracegnav(),tracepeph() -* 2009/12/22 1.12 change default parameter of ionos std -* valid under second for timeget() -* 2010/07/28 1.13 fix bug in tropmapf() -* added api: -* obs2code(),code2obs(),cross3(),normv3(), -* gst2time(),time2gst(),time_str(),timeset(), -* deg2dms(),dms2deg(),searchpcv(),antmodel_s(), -* tracehnav(),tracepclk(),reppath(),reppaths(), -* createdir() -* changed api: -* readpcv(), -* deleted api: -* uniqeph() -* 2010/08/20 1.14 omit to include mkl header files -* fix bug on chi-sqr(n) table -* 2010/12/11 1.15 added api: -* freeobs(),freenav(),ionppp() -* 2011/05/28 1.16 fix bug on half-hour offset by time2epoch() -* added api: -* uniqnav() -* 2012/06/09 1.17 add a leap second after 2012-6-30 -* 2012/07/15 1.18 add api setbits(),setbitu(),utc2gmst() -* fix bug on interpolation of antenna pcv -* fix bug on str2num() for string with over 256 char -* add api readblq(),satexclude(),setcodepri(), -* getcodepri() -* change api obs2code(),code2obs(),antmodel() -* 2012/12/25 1.19 fix bug on satwavelen(),code2obs(),obs2code() -* add api testsnr() -* 2013/01/04 1.20 add api gpst2bdt(),bdt2gpst(),bdt2time(),time2bdt() -* readblq(),readerp(),geterp(),crc16() -* change api eci2ecef(),sunmoonpos() -* 2013/03/26 1.21 tickget() uses clock_gettime() for linux -* 2013/05/08 1.22 fix bug on nutation coefficients for ast_args() -* 2013/06/02 1.23 add #ifdef for undefined CLOCK_MONOTONIC_RAW -* 2013/09/01 1.24 fix bug on interpolation of satellite antenna pcv -* 2013/09/06 1.25 fix bug on extrapolation of erp -* 2014/04/27 1.26 add SYS_LEO for satellite system -* add BDS L1 code for RINEX 3.02 and RTCM 3.2 -* support BDS L1 in satwavelen() -* 2014/05/29 1.27 fix bug on obs2code() to search obs code table -* 2014/08/26 1.28 fix problem on output of uncompress() for tar file -* add function to swap trace file with keywords -* 2014/10/21 1.29 strtok() -> strtok_r() in expath() for thread-safe -* add bdsmodear in procopt_default -* 2015/03/19 1.30 fix bug on interpolation of erp values in geterp() -* add leap second insertion before 2015/07/01 00:00 -* add api read_leaps() -* 2015/05/31 1.31 delte api windupcorr() -* 2015/08/08 1.32 add compile option CPUTIME_IN_GPST -* add api add_fatal() -* support usno leapsec.dat for api read_leaps() -* 2016/01/23 1.33 enable septentrio -* 2016/02/05 1.34 support GLONASS for savenav(), loadnav() -* 2016/06/11 1.35 delete trace() in reppath() to avoid deadlock -* 2016/07/01 1.36 support IRNSS -* add leap second before 2017/1/1 00:00:00 -* 2016/07/29 1.37 rename api compress() -> rtk_uncompress() -* rename api crc16() -> rtk_crc16() -* rename api crc24q() -> rtk_crc24q() -* rename api crc32() -> rtk_crc32() -* 2016/08/20 1.38 fix type incompatibility in win64 environment -* change constant _POSIX_C_SOURCE 199309 -> 199506 -* 2016/08/21 1.39 fix bug on week overflow in time2gpst()/gpst2time() -* 2016/09/05 1.40 fix bug on invalid nav data read in readnav() -* 2016/09/17 1.41 suppress warnings -* 2016/09/19 1.42 modify api deg2dms() to consider numerical error -*-----------------------------------------------------------------------------*/ + *----------------------------------------------------------------------------*/ + #include #include #include @@ -181,66 +59,66 @@ #include #include "rtklib_rtkcmn.h" -const double gpst0[]={1980,1, 6,0,0,0}; /* gps time reference */ -const double gst0 []={1999,8,22,0,0,0}; /* galileo system time reference */ -const double bdt0 []={2006,1, 1,0,0,0}; /* beidou time reference */ +const double gpst0[] = {1980, 1, 6, 0, 0, 0}; /* gps time reference */ +const double gst0 [] = {1999, 8, 22, 0, 0, 0}; /* galileo system time reference */ +const double bdt0 [] = {2006, 1, 1, 0, 0, 0}; /* beidou time reference */ -double leaps[MAXLEAPS+1][7]={ /* leap seconds (y,m,d,h,m,s,utc-gpst) */ - {2017,1,1,0,0,0,-18}, - {2015,7,1,0,0,0,-17}, - {2012,7,1,0,0,0,-16}, - {2009,1,1,0,0,0,-15}, - {2006,1,1,0,0,0,-14}, - {1999,1,1,0,0,0,-13}, - {1997,7,1,0,0,0,-12}, - {1996,1,1,0,0,0,-11}, - {1994,7,1,0,0,0,-10}, - {1993,7,1,0,0,0, -9}, - {1992,7,1,0,0,0, -8}, - {1991,1,1,0,0,0, -7}, - {1990,1,1,0,0,0, -6}, - {1988,1,1,0,0,0, -5}, - {1985,7,1,0,0,0, -4}, - {1983,7,1,0,0,0, -3}, - {1982,7,1,0,0,0, -2}, - {1981,7,1,0,0,0, -1}, +double leaps[MAXLEAPS+1][7] = { /* leap seconds (y,m,d,h,m,s,utc-gpst) */ + {2017, 1, 1, 0, 0, 0, -18}, + {2015, 7, 1, 0, 0, 0, -17}, + {2012, 7, 1, 0, 0, 0, -16}, + {2009, 1, 1, 0, 0, 0, -15}, + {2006, 1, 1, 0, 0, 0, -14}, + {1999, 1, 1, 0, 0, 0, -13}, + {1997, 7, 1, 0, 0, 0, -12}, + {1996, 1, 1, 0, 0, 0, -11}, + {1994, 7, 1, 0, 0, 0, -10}, + {1993, 7, 1, 0, 0, 0, -9}, + {1992, 7, 1, 0, 0, 0, -8}, + {1991, 1, 1, 0, 0, 0, -7}, + {1990, 1, 1, 0, 0, 0, -6}, + {1988, 1, 1, 0, 0, 0, -5}, + {1985, 7, 1, 0, 0, 0, -4}, + {1983, 7, 1, 0, 0, 0, -3}, + {1982, 7, 1, 0, 0, 0, -2}, + {1981, 7, 1, 0, 0, 0, -1}, {} }; -const prcopt_t prcopt_default={ /* defaults processing options */ - PMODE_SINGLE,0,2,SYS_GPS, /* mode,soltype,nf,navsys */ - 15.0*D2R,{}, /* elmin,snrmask */ - 0,1,1,1, /* sateph,modear,glomodear,bdsmodear */ - 5,0,10,1, /* maxout,minlock,minfix,armaxiter */ - 0,0,0,0, /* estion,esttrop,dynamics,tidecorr */ - 1,0,0,0,0, /* niter,codesmooth,intpref,sbascorr,sbassatsel */ - 0,0, /* rovpos,refpos */ - {100.0,100.0}, /* eratio[] */ - {100.0,0.003,0.003,0.0,1.0}, /* err[] */ - {30.0,0.03,0.3}, /* std[] */ - {1e-4,1e-3,1e-4,1e-1,1e-2,0.0}, /* prn[] */ - 5E-12, /* sclkstab */ - {3.0,0.9999,0.25,0.1,0.05}, /* thresar */ - 0.0,0.0,0.05, /* elmaskar,almaskhold,thresslip */ - 30.0,30.0,30.0, /* maxtdif,maxinno,maxgdop */ - {},{},{}, /* baseline,ru,rb */ +const prcopt_t prcopt_default = { /* defaults processing options */ + PMODE_SINGLE, 0, 2, SYS_GPS, /* mode, soltype, nf, navsys */ + 15.0*D2R, {}, /* elmin, snrmask */ + 0, 1, 1, 1, /* sateph, modear, glomodear, bdsmodear */ + 5, 0, 10, 1, /* maxout, minlock, minfix, armaxiter */ + 0, 0, 0, 0, /* estion, esttrop, dynamics, tidecorr */ + 1, 0, 0, 0, 0, /* niter, codesmooth, intpref, sbascorr, sbassatsel */ + 0, 0, /* rovpos, refpos */ + {100.0, 100.0}, /* eratio[] */ + {100.0, 0.003, 0.003, 0.0, 1.0}, /* err[] */ + {30.0, 0.03, 0.3}, /* std[] */ + {1e-4, 1e-3, 1e-4, 1e-1, 1e-2, 0.0}, /* prn[] */ + 5E-12, /* sclkstab */ + {3.0, 0.9999, 0.25, 0.1, 0.05}, /* thresar */ + 0.0, 0.0, 0.05, /* elmaskar, almaskhold, thresslip */ + 30.0, 30.0, 30.0, /* maxtdif, maxinno, maxgdop */ + {}, {}, {}, /* baseline, ru, rb */ {"",""}, /* anttype */ - {},{},{}, /* antdel,pcv,exsats */ - {},{},{},{},{},{},{},{},{},{} + {}, {}, {}, /* antdel, pcv, exsats */ + {}, {}, {}, {}, {}, {}, {}, {}, {}, {} }; -const solopt_t solopt_default={ /* defaults solution output options */ - SOLF_LLH,TIMES_GPST,1,3, /* posf,times,timef,timeu */ - 0,1,0,0,0,0, /* degf,outhead,outopt,datum,height,geoid */ - 0,0,0, /* solstatic,sstat,trace */ - {0.0,0.0}, /* nmeaintv */ +const solopt_t solopt_default = { /* defaults solution output options */ + SOLF_LLH, TIMES_GPST, 1, 3, /* posf, times, timef, timeu */ + 0, 1, 0, 0, 0, 0, /* degf, outhead, outopt, datum, height, geoid */ + 0, 0, 0, /* solstatic, sstat, trace */ + {0.0, 0.0}, /* nmeaintv */ " ","" /* separator/program name */ }; -const char *formatstrs[32]={ /* stream format strings */ +const char *formatstrs[32] = { /* stream format strings */ "RTCM 2", /* 0 */ "RTCM 3", /* 1 */ "NovAtel OEM6", /* 2 */ @@ -266,7 +144,7 @@ const char *formatstrs[32]={ /* stream format strings */ }; -char obscodes[][3]={ /* observation code strings */ +char obscodes[][3] = { /* observation code strings */ "" ,"1C","1P","1W","1Y", "1M","1N","1S","1L","1E", /* 0- 9 */ "1A","1B","1X","1Z","2C", "2D","2S","2L","2X","2P", /* 10-19 */ "2W","2Y","2M","2N","5I", "5Q","5X","7I","7Q","7X", /* 20-29 */ @@ -276,7 +154,7 @@ char obscodes[][3]={ /* observation code strings */ }; -unsigned char obsfreqs[]={ +unsigned char obsfreqs[] = { /* 1:L1/E1, 2:L2/B1, 3:L5/E5a/L3, 4:L6/LEX/B3, 5:E5b/B2, 6:E5(a+b), 7:S */ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, /* 0- 9 */ 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, /* 10-19 */ @@ -287,91 +165,91 @@ unsigned char obsfreqs[]={ }; -char codepris[7][MAXFREQ][16]={ /* code priority table */ +char codepris[7][MAXFREQ][16] = { /* code priority table */ /* L1/E1 L2/B1 L5/E5a/L3 L6/LEX/B3 E5b/B2 E5(a+b) S */ - {"CPYWMNSL","PYWCMNDSLX","IQX" ,"" ,"" ,"" ,"" }, /* GPS */ - {"PC" ,"PC" ,"IQX" ,"" ,"" ,"" ,"" }, /* GLO */ - {"CABXZ" ,"" ,"IQX" ,"ABCXZ" ,"IQX" ,"IQX" ,"" }, /* GAL */ - {"CSLXZ" ,"SLX" ,"IQX" ,"SLX" ,"" ,"" ,"" }, /* QZS */ - {"C" ,"" ,"IQX" ,"" ,"" ,"" ,"" }, /* SBS */ - {"IQX" ,"IQX" ,"IQX" ,"IQX" ,"IQX" ,"" ,"" }, /* BDS */ - {"" ,"" ,"ABCX" ,"" ,"" ,"" ,"ABCX"} /* IRN */ + {"CPYWMNSL", "PYWCMNDSLX", "IQX" , "" , "" , "" , "" }, /* GPS */ + {"PC" , "PC" , "IQX" , "" , "" , "" , "" }, /* GLO */ + {"CABXZ" , "" , "IQX" , "ABCXZ" , "IQX" , "IQX" , "" }, /* GAL */ + {"CSLXZ" , "SLX" , "IQX" , "SLX" , "" , "" , "" }, /* QZS */ + {"C" , "" , "IQX" , "" , "" , "" , "" }, /* SBS */ + {"IQX" , "IQX" , "IQX" , "IQX" , "IQX" , "" , "" }, /* BDS */ + {"" , "" , "ABCX" , "" , "" , "" , "ABCX"} /* IRN */ }; -fatalfunc_t *fatalfunc=NULL; /* fatal callback function */ +fatalfunc_t *fatalfunc = NULL; /* fatal callback function */ /* crc tables generated by util/gencrc ---------------------------------------*/ -const unsigned short tbl_CRC16[]={ - 0x0000,0x1021,0x2042,0x3063,0x4084,0x50A5,0x60C6,0x70E7, - 0x8108,0x9129,0xA14A,0xB16B,0xC18C,0xD1AD,0xE1CE,0xF1EF, - 0x1231,0x0210,0x3273,0x2252,0x52B5,0x4294,0x72F7,0x62D6, - 0x9339,0x8318,0xB37B,0xA35A,0xD3BD,0xC39C,0xF3FF,0xE3DE, - 0x2462,0x3443,0x0420,0x1401,0x64E6,0x74C7,0x44A4,0x5485, - 0xA56A,0xB54B,0x8528,0x9509,0xE5EE,0xF5CF,0xC5AC,0xD58D, - 0x3653,0x2672,0x1611,0x0630,0x76D7,0x66F6,0x5695,0x46B4, - 0xB75B,0xA77A,0x9719,0x8738,0xF7DF,0xE7FE,0xD79D,0xC7BC, - 0x48C4,0x58E5,0x6886,0x78A7,0x0840,0x1861,0x2802,0x3823, - 0xC9CC,0xD9ED,0xE98E,0xF9AF,0x8948,0x9969,0xA90A,0xB92B, - 0x5AF5,0x4AD4,0x7AB7,0x6A96,0x1A71,0x0A50,0x3A33,0x2A12, - 0xDBFD,0xCBDC,0xFBBF,0xEB9E,0x9B79,0x8B58,0xBB3B,0xAB1A, - 0x6CA6,0x7C87,0x4CE4,0x5CC5,0x2C22,0x3C03,0x0C60,0x1C41, - 0xEDAE,0xFD8F,0xCDEC,0xDDCD,0xAD2A,0xBD0B,0x8D68,0x9D49, - 0x7E97,0x6EB6,0x5ED5,0x4EF4,0x3E13,0x2E32,0x1E51,0x0E70, - 0xFF9F,0xEFBE,0xDFDD,0xCFFC,0xBF1B,0xAF3A,0x9F59,0x8F78, - 0x9188,0x81A9,0xB1CA,0xA1EB,0xD10C,0xC12D,0xF14E,0xE16F, - 0x1080,0x00A1,0x30C2,0x20E3,0x5004,0x4025,0x7046,0x6067, - 0x83B9,0x9398,0xA3FB,0xB3DA,0xC33D,0xD31C,0xE37F,0xF35E, - 0x02B1,0x1290,0x22F3,0x32D2,0x4235,0x5214,0x6277,0x7256, - 0xB5EA,0xA5CB,0x95A8,0x8589,0xF56E,0xE54F,0xD52C,0xC50D, - 0x34E2,0x24C3,0x14A0,0x0481,0x7466,0x6447,0x5424,0x4405, - 0xA7DB,0xB7FA,0x8799,0x97B8,0xE75F,0xF77E,0xC71D,0xD73C, - 0x26D3,0x36F2,0x0691,0x16B0,0x6657,0x7676,0x4615,0x5634, - 0xD94C,0xC96D,0xF90E,0xE92F,0x99C8,0x89E9,0xB98A,0xA9AB, - 0x5844,0x4865,0x7806,0x6827,0x18C0,0x08E1,0x3882,0x28A3, - 0xCB7D,0xDB5C,0xEB3F,0xFB1E,0x8BF9,0x9BD8,0xABBB,0xBB9A, - 0x4A75,0x5A54,0x6A37,0x7A16,0x0AF1,0x1AD0,0x2AB3,0x3A92, - 0xFD2E,0xED0F,0xDD6C,0xCD4D,0xBDAA,0xAD8B,0x9DE8,0x8DC9, - 0x7C26,0x6C07,0x5C64,0x4C45,0x3CA2,0x2C83,0x1CE0,0x0CC1, - 0xEF1F,0xFF3E,0xCF5D,0xDF7C,0xAF9B,0xBFBA,0x8FD9,0x9FF8, - 0x6E17,0x7E36,0x4E55,0x5E74,0x2E93,0x3EB2,0x0ED1,0x1EF0 +const unsigned short tbl_CRC16[] = { + 0x0000, 0x1021, 0x2042, 0x3063, 0x4084, 0x50A5, 0x60C6, 0x70E7, + 0x8108, 0x9129, 0xA14A, 0xB16B, 0xC18C, 0xD1AD, 0xE1CE, 0xF1EF, + 0x1231, 0x0210, 0x3273, 0x2252, 0x52B5, 0x4294, 0x72F7, 0x62D6, + 0x9339, 0x8318, 0xB37B, 0xA35A, 0xD3BD, 0xC39C, 0xF3FF, 0xE3DE, + 0x2462, 0x3443, 0x0420, 0x1401, 0x64E6, 0x74C7, 0x44A4, 0x5485, + 0xA56A, 0xB54B, 0x8528, 0x9509, 0xE5EE, 0xF5CF, 0xC5AC, 0xD58D, + 0x3653, 0x2672, 0x1611, 0x0630, 0x76D7, 0x66F6, 0x5695, 0x46B4, + 0xB75B, 0xA77A, 0x9719, 0x8738, 0xF7DF, 0xE7FE, 0xD79D, 0xC7BC, + 0x48C4, 0x58E5, 0x6886, 0x78A7, 0x0840, 0x1861, 0x2802, 0x3823, + 0xC9CC, 0xD9ED, 0xE98E, 0xF9AF, 0x8948, 0x9969, 0xA90A, 0xB92B, + 0x5AF5, 0x4AD4, 0x7AB7, 0x6A96, 0x1A71, 0x0A50, 0x3A33, 0x2A12, + 0xDBFD, 0xCBDC, 0xFBBF, 0xEB9E, 0x9B79, 0x8B58, 0xBB3B, 0xAB1A, + 0x6CA6, 0x7C87, 0x4CE4, 0x5CC5, 0x2C22, 0x3C03, 0x0C60, 0x1C41, + 0xEDAE, 0xFD8F, 0xCDEC, 0xDDCD, 0xAD2A, 0xBD0B, 0x8D68, 0x9D49, + 0x7E97, 0x6EB6, 0x5ED5, 0x4EF4, 0x3E13, 0x2E32, 0x1E51, 0x0E70, + 0xFF9F, 0xEFBE, 0xDFDD, 0xCFFC, 0xBF1B, 0xAF3A, 0x9F59, 0x8F78, + 0x9188, 0x81A9, 0xB1CA, 0xA1EB, 0xD10C, 0xC12D, 0xF14E, 0xE16F, + 0x1080, 0x00A1, 0x30C2, 0x20E3, 0x5004, 0x4025, 0x7046, 0x6067, + 0x83B9, 0x9398, 0xA3FB, 0xB3DA, 0xC33D, 0xD31C, 0xE37F, 0xF35E, + 0x02B1, 0x1290, 0x22F3, 0x32D2, 0x4235, 0x5214, 0x6277, 0x7256, + 0xB5EA, 0xA5CB, 0x95A8, 0x8589, 0xF56E, 0xE54F, 0xD52C, 0xC50D, + 0x34E2, 0x24C3, 0x14A0, 0x0481, 0x7466, 0x6447, 0x5424, 0x4405, + 0xA7DB, 0xB7FA, 0x8799, 0x97B8, 0xE75F, 0xF77E, 0xC71D, 0xD73C, + 0x26D3, 0x36F2, 0x0691, 0x16B0, 0x6657, 0x7676, 0x4615, 0x5634, + 0xD94C, 0xC96D, 0xF90E, 0xE92F, 0x99C8, 0x89E9, 0xB98A, 0xA9AB, + 0x5844, 0x4865, 0x7806, 0x6827, 0x18C0, 0x08E1, 0x3882, 0x28A3, + 0xCB7D, 0xDB5C, 0xEB3F, 0xFB1E, 0x8BF9, 0x9BD8, 0xABBB, 0xBB9A, + 0x4A75, 0x5A54, 0x6A37, 0x7A16, 0x0AF1, 0x1AD0, 0x2AB3, 0x3A92, + 0xFD2E, 0xED0F, 0xDD6C, 0xCD4D, 0xBDAA, 0xAD8B, 0x9DE8, 0x8DC9, + 0x7C26, 0x6C07, 0x5C64, 0x4C45, 0x3CA2, 0x2C83, 0x1CE0, 0x0CC1, + 0xEF1F, 0xFF3E, 0xCF5D, 0xDF7C, 0xAF9B, 0xBFBA, 0x8FD9, 0x9FF8, + 0x6E17, 0x7E36, 0x4E55, 0x5E74, 0x2E93, 0x3EB2, 0x0ED1, 0x1EF0 }; -const unsigned int tbl_CRC24Q[]={ - 0x000000,0x864CFB,0x8AD50D,0x0C99F6,0x93E6E1,0x15AA1A,0x1933EC,0x9F7F17, - 0xA18139,0x27CDC2,0x2B5434,0xAD18CF,0x3267D8,0xB42B23,0xB8B2D5,0x3EFE2E, - 0xC54E89,0x430272,0x4F9B84,0xC9D77F,0x56A868,0xD0E493,0xDC7D65,0x5A319E, - 0x64CFB0,0xE2834B,0xEE1ABD,0x685646,0xF72951,0x7165AA,0x7DFC5C,0xFBB0A7, - 0x0CD1E9,0x8A9D12,0x8604E4,0x00481F,0x9F3708,0x197BF3,0x15E205,0x93AEFE, - 0xAD50D0,0x2B1C2B,0x2785DD,0xA1C926,0x3EB631,0xB8FACA,0xB4633C,0x322FC7, - 0xC99F60,0x4FD39B,0x434A6D,0xC50696,0x5A7981,0xDC357A,0xD0AC8C,0x56E077, - 0x681E59,0xEE52A2,0xE2CB54,0x6487AF,0xFBF8B8,0x7DB443,0x712DB5,0xF7614E, - 0x19A3D2,0x9FEF29,0x9376DF,0x153A24,0x8A4533,0x0C09C8,0x00903E,0x86DCC5, - 0xB822EB,0x3E6E10,0x32F7E6,0xB4BB1D,0x2BC40A,0xAD88F1,0xA11107,0x275DFC, - 0xDCED5B,0x5AA1A0,0x563856,0xD074AD,0x4F0BBA,0xC94741,0xC5DEB7,0x43924C, - 0x7D6C62,0xFB2099,0xF7B96F,0x71F594,0xEE8A83,0x68C678,0x645F8E,0xE21375, - 0x15723B,0x933EC0,0x9FA736,0x19EBCD,0x8694DA,0x00D821,0x0C41D7,0x8A0D2C, - 0xB4F302,0x32BFF9,0x3E260F,0xB86AF4,0x2715E3,0xA15918,0xADC0EE,0x2B8C15, - 0xD03CB2,0x567049,0x5AE9BF,0xDCA544,0x43DA53,0xC596A8,0xC90F5E,0x4F43A5, - 0x71BD8B,0xF7F170,0xFB6886,0x7D247D,0xE25B6A,0x641791,0x688E67,0xEEC29C, - 0x3347A4,0xB50B5F,0xB992A9,0x3FDE52,0xA0A145,0x26EDBE,0x2A7448,0xAC38B3, - 0x92C69D,0x148A66,0x181390,0x9E5F6B,0x01207C,0x876C87,0x8BF571,0x0DB98A, - 0xF6092D,0x7045D6,0x7CDC20,0xFA90DB,0x65EFCC,0xE3A337,0xEF3AC1,0x69763A, - 0x578814,0xD1C4EF,0xDD5D19,0x5B11E2,0xC46EF5,0x42220E,0x4EBBF8,0xC8F703, - 0x3F964D,0xB9DAB6,0xB54340,0x330FBB,0xAC70AC,0x2A3C57,0x26A5A1,0xA0E95A, - 0x9E1774,0x185B8F,0x14C279,0x928E82,0x0DF195,0x8BBD6E,0x872498,0x016863, - 0xFAD8C4,0x7C943F,0x700DC9,0xF64132,0x693E25,0xEF72DE,0xE3EB28,0x65A7D3, - 0x5B59FD,0xDD1506,0xD18CF0,0x57C00B,0xC8BF1C,0x4EF3E7,0x426A11,0xC426EA, - 0x2AE476,0xACA88D,0xA0317B,0x267D80,0xB90297,0x3F4E6C,0x33D79A,0xB59B61, - 0x8B654F,0x0D29B4,0x01B042,0x87FCB9,0x1883AE,0x9ECF55,0x9256A3,0x141A58, - 0xEFAAFF,0x69E604,0x657FF2,0xE33309,0x7C4C1E,0xFA00E5,0xF69913,0x70D5E8, - 0x4E2BC6,0xC8673D,0xC4FECB,0x42B230,0xDDCD27,0x5B81DC,0x57182A,0xD154D1, - 0x26359F,0xA07964,0xACE092,0x2AAC69,0xB5D37E,0x339F85,0x3F0673,0xB94A88, - 0x87B4A6,0x01F85D,0x0D61AB,0x8B2D50,0x145247,0x921EBC,0x9E874A,0x18CBB1, - 0xE37B16,0x6537ED,0x69AE1B,0xEFE2E0,0x709DF7,0xF6D10C,0xFA48FA,0x7C0401, - 0x42FA2F,0xC4B6D4,0xC82F22,0x4E63D9,0xD11CCE,0x575035,0x5BC9C3,0xDD8538 +const unsigned int tbl_CRC24Q[] = { + 0x000000, 0x864CFB, 0x8AD50D, 0x0C99F6, 0x93E6E1, 0x15AA1A, 0x1933EC, 0x9F7F17, + 0xA18139, 0x27CDC2, 0x2B5434, 0xAD18CF, 0x3267D8, 0xB42B23, 0xB8B2D5, 0x3EFE2E, + 0xC54E89, 0x430272, 0x4F9B84, 0xC9D77F, 0x56A868, 0xD0E493, 0xDC7D65, 0x5A319E, + 0x64CFB0, 0xE2834B, 0xEE1ABD, 0x685646, 0xF72951, 0x7165AA, 0x7DFC5C, 0xFBB0A7, + 0x0CD1E9, 0x8A9D12, 0x8604E4, 0x00481F, 0x9F3708, 0x197BF3, 0x15E205, 0x93AEFE, + 0xAD50D0, 0x2B1C2B, 0x2785DD, 0xA1C926, 0x3EB631, 0xB8FACA, 0xB4633C, 0x322FC7, + 0xC99F60, 0x4FD39B, 0x434A6D, 0xC50696, 0x5A7981, 0xDC357A, 0xD0AC8C, 0x56E077, + 0x681E59, 0xEE52A2, 0xE2CB54, 0x6487AF, 0xFBF8B8, 0x7DB443, 0x712DB5, 0xF7614E, + 0x19A3D2, 0x9FEF29, 0x9376DF, 0x153A24, 0x8A4533, 0x0C09C8, 0x00903E, 0x86DCC5, + 0xB822EB, 0x3E6E10, 0x32F7E6, 0xB4BB1D, 0x2BC40A, 0xAD88F1, 0xA11107, 0x275DFC, + 0xDCED5B, 0x5AA1A0, 0x563856, 0xD074AD, 0x4F0BBA, 0xC94741, 0xC5DEB7, 0x43924C, + 0x7D6C62, 0xFB2099, 0xF7B96F, 0x71F594, 0xEE8A83, 0x68C678, 0x645F8E, 0xE21375, + 0x15723B, 0x933EC0, 0x9FA736, 0x19EBCD, 0x8694DA, 0x00D821, 0x0C41D7, 0x8A0D2C, + 0xB4F302, 0x32BFF9, 0x3E260F, 0xB86AF4, 0x2715E3, 0xA15918, 0xADC0EE, 0x2B8C15, + 0xD03CB2, 0x567049, 0x5AE9BF, 0xDCA544, 0x43DA53, 0xC596A8, 0xC90F5E, 0x4F43A5, + 0x71BD8B, 0xF7F170, 0xFB6886, 0x7D247D, 0xE25B6A, 0x641791, 0x688E67, 0xEEC29C, + 0x3347A4, 0xB50B5F, 0xB992A9, 0x3FDE52, 0xA0A145, 0x26EDBE, 0x2A7448, 0xAC38B3, + 0x92C69D, 0x148A66, 0x181390, 0x9E5F6B, 0x01207C, 0x876C87, 0x8BF571, 0x0DB98A, + 0xF6092D, 0x7045D6, 0x7CDC20, 0xFA90DB, 0x65EFCC, 0xE3A337, 0xEF3AC1, 0x69763A, + 0x578814, 0xD1C4EF, 0xDD5D19, 0x5B11E2, 0xC46EF5, 0x42220E, 0x4EBBF8, 0xC8F703, + 0x3F964D, 0xB9DAB6, 0xB54340, 0x330FBB, 0xAC70AC, 0x2A3C57, 0x26A5A1, 0xA0E95A, + 0x9E1774, 0x185B8F, 0x14C279, 0x928E82, 0x0DF195, 0x8BBD6E, 0x872498, 0x016863, + 0xFAD8C4, 0x7C943F, 0x700DC9, 0xF64132, 0x693E25, 0xEF72DE, 0xE3EB28, 0x65A7D3, + 0x5B59FD, 0xDD1506, 0xD18CF0, 0x57C00B, 0xC8BF1C, 0x4EF3E7, 0x426A11, 0xC426EA, + 0x2AE476, 0xACA88D, 0xA0317B, 0x267D80, 0xB90297, 0x3F4E6C, 0x33D79A, 0xB59B61, + 0x8B654F, 0x0D29B4, 0x01B042, 0x87FCB9, 0x1883AE, 0x9ECF55, 0x9256A3, 0x141A58, + 0xEFAAFF, 0x69E604, 0x657FF2, 0xE33309, 0x7C4C1E, 0xFA00E5, 0xF69913, 0x70D5E8, + 0x4E2BC6, 0xC8673D, 0xC4FECB, 0x42B230, 0xDDCD27, 0x5B81DC, 0x57182A, 0xD154D1, + 0x26359F, 0xA07964, 0xACE092, 0x2AAC69, 0xB5D37E, 0x339F85, 0x3F0673, 0xB94A88, + 0x87B4A6, 0x01F85D, 0x0D61AB, 0x8B2D50, 0x145247, 0x921EBC, 0x9E874A, 0x18CBB1, + 0xE37B16, 0x6537ED, 0x69AE1B, 0xEFE2E0, 0x709DF7, 0xF6D10C, 0xFA48FA, 0x7C0401, + 0x42FA2F, 0xC4B6D4, 0xC82F22, 0x4E63D9, 0xD11CCE, 0x575035, 0x5BC9C3, 0xDD8538 }; @@ -397,8 +275,8 @@ void fatalerr(const char *format, ...) { char msg[1024]; va_list ap; - va_start(ap,format); vsprintf(msg,format,ap); va_end(ap); - fprintf(stderr,"%s",msg); + va_start(ap, format); vsprintf(msg, format, ap); va_end(ap); + fprintf(stderr, "%s", msg); exit(-9); } @@ -411,33 +289,33 @@ void fatalerr(const char *format, ...) *-----------------------------------------------------------------------------*/ int satno(int sys, int prn) { - if (prn<=0) return 0; + if (prn <= 0) return 0; switch (sys) { case SYS_GPS: - if (prnexsats[sat-1]==1) return 1; /* excluded satellite */ - if (opt->exsats[sat-1]==2) return 0; /* included satellite */ + if (opt->exsats[sat-1] == 1) return 1; /* excluded satellite */ + if (opt->exsats[sat-1] == 2) return 0; /* included satellite */ if (!(sys&opt->navsys)) return 1; /* unselected sat sys */ } - if (sys==SYS_QZS) svh&=0xFE; /* mask QZSS LEX health */ + if (sys == SYS_QZS) svh&=0xFE; /* mask QZSS LEX health */ if (svh) { - trace(3,"unhealthy satellite: sat=%3d svh=%02X\n",sat,svh); + trace(3, "unhealthy satellite: sat=%3d svh=%02X\n", sat, svh); return 1; } return 0; @@ -604,16 +482,16 @@ int satexclude(int sat, int svh, const prcopt_t *opt) int testsnr(int base, int freq, double el, double snr, const snrmask_t *mask) { - double minsnr,a; + double minsnr, a; int i; - if (!mask->ena[base]||freq<0||freq>=NFREQ) return 0; + if (!mask->ena[base] || freq<0 || freq >= NFREQ) return 0; - a=(el*R2D+5.0)/10.0; - i=(int)floor(a); a-=i; - if (i<1) minsnr=mask->mask[freq][0]; - else if (i>8) minsnr=mask->mask[freq][8]; - else minsnr=(1.0-a)*mask->mask[freq][i-1]+a*mask->mask[freq][i]; + a = (el*R2D+5.0)/10.0; + i = (int)floor(a); a -= i; + if (i<1) minsnr = mask->mask[freq][0]; + else if (i>8) minsnr = mask->mask[freq][8]; + else minsnr = (1.0-a)*mask->mask[freq][i-1]+a*mask->mask[freq][i]; return snr>(7-i%8))&1u); + for (i = pos; i>(7-i%8))&1u); return bits; } int getbits(const unsigned char *buff, int pos, int len) { - unsigned int bits=getbitu(buff,pos,len); - if (len<=0||32<=len||!(bits&(1u<<(len-1)))) return (int)bits; + unsigned int bits = getbitu(buff, pos, len); + if (len <= 0 || 32 <= len || !(bits&(1u<<(len-1)))) return (int)bits; return (int)(bits|(~0u<>=1) + if (len <= 0 || 32>=1) { if (data&mask) buff[i/8]|=1u<<(7-i%8); else buff[i/8]&=~(1u<<(7-i%8)); } @@ -766,7 +644,7 @@ void setbitu(unsigned char *buff, int pos, int len, unsigned int data) void setbits(unsigned char *buff, int pos, int len, int data) { if (data<0) data|=1<<(len-1); else data&=~(1<<(len-1)); /* set sign bit */ - setbitu(buff,pos,len,(unsigned int)data); + setbitu(buff, pos, len, (unsigned int)data); } @@ -779,17 +657,17 @@ void setbits(unsigned char *buff, int pos, int len, int data) *-----------------------------------------------------------------------------*/ unsigned int rtk_crc32(const unsigned char *buff, int len) { - unsigned int crc=0; - int i,j; + unsigned int crc = 0; + int i, j; - trace(4,"rtk_crc32: len=%d\n",len); + trace(4, "rtk_crc32: len=%d\n", len); - for (i=0;i>1)^POLYCRC32; else crc>>=1; + if (crc&1) crc = (crc>>1)^POLYCRC32; else crc>>=1; } } return crc; @@ -805,12 +683,12 @@ unsigned int rtk_crc32(const unsigned char *buff, int len) *-----------------------------------------------------------------------------*/ unsigned int rtk_crc24q(const unsigned char *buff, int len) { - unsigned int crc=0; + unsigned int crc = 0; int i; - trace(4,"rtk_crc24q: len=%d\n",len); + trace(4, "rtk_crc24q: len=%d\n", len); - for (i=0;i>16)^buff[i]]; + for (i = 0; i>16)^buff[i]]; return crc; } @@ -824,14 +702,14 @@ unsigned int rtk_crc24q(const unsigned char *buff, int len) *-----------------------------------------------------------------------------*/ unsigned short rtk_crc16(const unsigned char *buff, int len) { - unsigned short crc=0; + unsigned short crc = 0; int i; - trace(4,"rtk_crc16: len=%d\n",len); + trace(4, "rtk_crc16: len=%d\n", len); - for (i=0;i>8)^buff[i])&0xFF]; + crc = (crc<<8)^tbl_CRC16[((crc>>8)^buff[i])&0xFF]; } return crc; } @@ -848,24 +726,24 @@ unsigned short rtk_crc16(const unsigned char *buff, int len) *-----------------------------------------------------------------------------*/ int decode_word(unsigned int word, unsigned char *data) { - const unsigned int hamming[]={ - 0xBB1F3480,0x5D8F9A40,0xAEC7CD00,0x5763E680,0x6BB1F340,0x8B7A89C0 + const unsigned int hamming[] = { + 0xBB1F3480, 0x5D8F9A40, 0xAEC7CD00, 0x5763E680, 0x6BB1F340, 0x8B7A89C0 }; - unsigned int parity=0,w; + unsigned int parity = 0, w; int i; - trace(5,"decodeword: word=%08x\n",word); + trace(5, "decodeword: word=%08x\n", word); if (word&0x40000000) word^=0x3FFFFFC0; - for (i=0;i<6;i++) + for (i = 0; i<6; i++) { - parity<<=1; - for (w=(word&hamming[i])>>6;w;w>>=1) parity^=w&1; + parity<<= 1; + for (w = (word&hamming[i])>>6; w; w>>=1) parity^=w&1; } - if (parity!=(word&0x3F)) return 0; + if (parity != (word&0x3F)) return 0; - for (i=0;i<3;i++) data[i]=(unsigned char)(word>>(22-i*8)); + for (i = 0; i<3; i++) data[i] = (unsigned char)(word>>(22-i*8)); return 1; } @@ -879,10 +757,10 @@ double *mat(int n, int m) { double *p; - if (n<=0||m<=0) return NULL; - if (!(p=(double *)malloc(sizeof(double)*n*m))) + if (n <= 0 || m <= 0) return NULL; + if (!(p = (double *)malloc(sizeof(double)*n*m))) { - fatalerr("matrix memory allocation error: n=%d,m=%d\n",n,m); + fatalerr("matrix memory allocation error: n=%d,m=%d\n", n, m); } return p; } @@ -891,16 +769,16 @@ double *mat(int n, int m) /* new integer matrix ---------------------------------------------------------- * allocate memory of integer matrix * args : int n,m I number of rows and columns of matrix - * return : matrix pointer (if n<=0 or m<=0, return NULL) + * return : matrix pointer (if n <= 0 or m <= 0, return NULL) *-----------------------------------------------------------------------------*/ int *imat(int n, int m) { int *p; - if (n<=0||m<=0) return NULL; - if (!(p=(int *)malloc(sizeof(int)*n*m))) + if (n <= 0 || m <= 0) return NULL; + if (!(p = (int *)malloc(sizeof(int)*n*m))) { - fatalerr("integer matrix memory allocation error: n=%d,m=%d\n",n,m); + fatalerr("integer matrix memory allocation error: n=%d,m=%d\n", n, m); } return p; } @@ -909,18 +787,18 @@ int *imat(int n, int m) /* zero matrix ----------------------------------------------------------------- * generate new zero matrix * args : int n,m I number of rows and columns of matrix - * return : matrix pointer (if n<=0 or m<=0, return NULL) + * return : matrix pointer (if n <= 0 or m <= 0, return NULL) *-----------------------------------------------------------------------------*/ double *zeros(int n, int m) { double *p; #if NOCALLOC - if ((p=mat(n,m))) for (n=n*m-1;n>=0;n--) p[n]=0.0; + if ((p = mat(n, m))) for (n = n*m-1; n >= 0; n--) p[n] = 0.0; #else - if (n<=0||m<=0) return NULL; - if (!(p=(double *)calloc(sizeof(double),n*m))) { - fatalerr("matrix memory allocation error: n=%d,m=%d\n",n,m); + if (n <= 0 || m <= 0) return NULL; + if (!(p = (double *)calloc(sizeof(double), n*m))) { + fatalerr("matrix memory allocation error: n=%d,m=%d\n", n, m); } #endif return p; @@ -930,14 +808,14 @@ double *zeros(int n, int m) /* identity matrix ------------------------------------------------------------- * generate new identity matrix * args : int n I number of rows and columns of matrix - * return : matrix pointer (if n<=0, return NULL) + * return : matrix pointer (if n <= 0, return NULL) *-----------------------------------------------------------------------------*/ double *eye(int n) { double *p; int i; - if ((p=zeros(n,n))) for (i=0;i=0) c+=a[n]*b[n]; + while (--n >= 0) c += a[n]*b[n]; return c; } @@ -965,7 +843,7 @@ double dot(const double *a, const double *b, int n) *-----------------------------------------------------------------------------*/ double norm(const double *a, int n) { - return sqrt(dot(a,a,n)); + return sqrt(dot(a, a, n)); } @@ -977,9 +855,9 @@ double norm(const double *a, int n) *-----------------------------------------------------------------------------*/ void cross3(const double *a, const double *b, double *c) { - 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]; + 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]; } @@ -992,10 +870,10 @@ void cross3(const double *a, const double *b, double *c) int normv3(const double *a, double *b) { double r; - if ((r=norm(a,3))<=0.0) return 0; - b[0]=a[0]/r; - b[1]=a[1]/r; - b[2]=a[2]/r; + if ((r = norm(a, 3)) <= 0.0) return 0; + b[0] = a[0]/r; + b[1] = a[1]/r; + b[2] = a[2]/r; return 1; } @@ -1009,7 +887,7 @@ int normv3(const double *a, double *b) *-----------------------------------------------------------------------------*/ void matcpy(double *A, const double *B, int n, int m) { - memcpy(A,B,sizeof(double)*n*m); + memcpy(A, B, sizeof(double)*n*m); } /* matrix routines -----------------------------------------------------------*/ @@ -1028,10 +906,10 @@ void matcpy(double *A, const double *B, int n, int m) void matmul(const char *tr, int n, int k, int m, double alpha, const double *A, const double *B, double beta, double *C) { - int lda=tr[0]=='T'?m:n,ldb=tr[1]=='T'?k:m; + int lda = tr[0] == 'T' ? m : n, ldb = tr[1] == 'T' ? k : m; - dgemm_((char *)tr,(char *)tr+1,&n,&k,&m,&alpha,(double *)A,&lda,(double *)B, - &ldb,&beta,C,&n); + dgemm_((char *)tr, (char *)tr+1, &n, &k, &m, &alpha, (double *)A, &lda, (double *)B, + &ldb, &beta, C, &n); } @@ -1044,11 +922,11 @@ void matmul(const char *tr, int n, int k, int m, double alpha, int matinv(double *A, int n) { double *work; - int info,lwork=n*16,*ipiv=imat(n,1); + int info, lwork = n*16, *ipiv = imat(n, 1); - work=mat(lwork,1); - dgetrf_(&n,&n,A,&n,ipiv,&info); - if (!info) dgetri_(&n,A,&n,ipiv,work,&lwork,&info); + work = mat(lwork, 1); + dgetrf_(&n, &n, A, &n, ipiv, &info); + if (!info) dgetri_(&n, A, &n, ipiv, work, &lwork, &info); free(ipiv); free(work); return info; } @@ -1068,13 +946,13 @@ int matinv(double *A, int n) int solve(const char *tr, const double *A, const double *Y, int n, int m, double *X) { - double *B=mat(n,n); - int info,*ipiv=imat(n,1); + double *B = mat(n, n); + int info, *ipiv = imat(n, 1); - matcpy(B,A,n,n); - matcpy(X,Y,n,m); - dgetrf_(&n,&n,B,&n,ipiv,&info); - if (!info) dgetrs_((char *)tr,&n,&m,B,&n,ipiv,X,&n,&info); + matcpy(B, A, n, n); + matcpy(X, Y, n, m); + dgetrf_(&n, &n, B, &n, ipiv, &info); + if (!info) dgetrs_((char *)tr, &n, &m, B, &n, ipiv, X, &n, &info); free(ipiv); free(B); return info; } @@ -1086,7 +964,7 @@ int solve(const char *tr, const double *A, const double *Y, int n, * least square estimation by solving normal equation (x=(A*A')^-1*A*y) * args : double *A I transpose of (weighted) design matrix (n x m) * double *y I (weighted) measurements (m x 1) - * int n,m I number of parameters and measurements (n<=m) + * int n,m I number of parameters and measurements (n <= m) * double *x O estmated parameters (n x 1) * double *Q O esimated parameters covariance matrix (n x n) * return : status (0:ok,0>:error) @@ -1100,10 +978,10 @@ int lsq(const double *A, const double *y, int n, int m, double *x, int info; if (m0.0) ix[k++]=i; - x_=mat(k,1); xp_=mat(k,1); P_=mat(k,k); Pp_=mat(k,k); H_=mat(k,m); - for (i=0;i0.0) ix[k++] = i; + x_ = mat(k, 1); xp_ = mat(k, 1); P_ = mat(k, k); Pp_ = mat(k, k); H_ = mat(k, m); + for (i = 0; i=0;s++) *p++=*s=='d'||*s=='D'?'E':*s; *p='\0'; - return sscanf(str,"%lf",&value)==1?value:0.0; + if (i<0 || (int)strlen(s)= 0; s++) *p++=*s == 'd' || *s == 'D' ? 'E' : *s; *p = '\0'; + return sscanf(str, "%lf", &value) == 1 ? value : 0.0; } @@ -1266,14 +1144,14 @@ double str2num(const char *s, int i, int n) int str2time(const char *s, int i, int n, gtime_t *t) { double ep[6]; - char str[256],*p=str; + char str[256], *p = str; - if (i<0||(int)strlen(s)=0;) *p++=*s++; *p='\0'; - if (sscanf(str,"%lf %lf %lf %lf %lf %lf",ep,ep+1,ep+2,ep+3,ep+4,ep+5)<6) + if (i<0 || (int)strlen(s)= 0;) *p++=*s++; *p = '\0'; + if (sscanf(str, "%lf %lf %lf %lf %lf %lf", ep, ep+1, ep+2, ep+3, ep+4, ep+5)<6) return -1; - if (ep[0]<100.0) ep[0]+=ep[0]<80.0?2000.0:1900.0; - *t=epoch2time(ep); + if (ep[0]<100.0) ep[0] += ep[0]<80.0 ? 2000.0 : 1900.0; + *t = epoch2time(ep); return 0; } @@ -1286,17 +1164,17 @@ int str2time(const char *s, int i, int n, gtime_t *t) *-----------------------------------------------------------------------------*/ gtime_t epoch2time(const double *ep) { - const int doy[]={1,32,60,91,121,152,182,213,244,274,305,335}; - gtime_t time={}; - int days,sec,year=(int)ep[0],mon=(int)ep[1],day=(int)ep[2]; + const int doy[] = {1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; + gtime_t time = {}; + int days, sec, year = (int)ep[0], mon = (int)ep[1], day = (int)ep[2]; - if (year<1970||2099=3?1:0); - sec=(int)floor(ep[5]); - time.time=(time_t)days*86400+(int)ep[3]*3600+(int)ep[4]*60+sec; - time.sec=ep[5]-sec; + days = (year-1970)*365+(year-1969)/4+doy[mon-1]+day-2+(year%4 == 0 && mon >= 3 ? 1 : 0); + sec = (int)floor(ep[5]); + time.time = (time_t)days*86400+(int)ep[3]*3600+(int)ep[4]*60+sec; + time.sec = ep[5]-sec; return time; } @@ -1310,21 +1188,21 @@ gtime_t epoch2time(const double *ep) *-----------------------------------------------------------------------------*/ void time2epoch(gtime_t t, double *ep) { - const int mday[]={ /* # of days in a month */ - 31,28,31,30,31,30,31,31,30,31,30,31,31,28,31,30,31,30,31,31,30,31,30,31, - 31,29,31,30,31,30,31,31,30,31,30,31,31,28,31,30,31,30,31,31,30,31,30,31 + const int mday[] = { /* # of days in a month */ + 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, + 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 }; - int days,sec,mon,day; + int days, sec, mon, day; /* leap year if year%4==0 in 1901-2099 */ - days=(int)(t.time/86400); - sec=(int)(t.time-(time_t)days*86400); - for (day=days%1461,mon=0;mon<48;mon++) + days = (int)(t.time/86400); + sec = (int)(t.time-(time_t)days*86400); + for (day = days%1461, mon = 0; mon<48; mon++) { - if (day>=mday[mon]) day-=mday[mon]; else break; + if (day >= mday[mon]) day -= mday[mon]; else break; } - ep[0]=1970+days/1461*4+mon/12; ep[1]=mon%12+1; ep[2]=day+1; - ep[3]=sec/3600; ep[4]=sec%3600/60; ep[5]=sec%60+t.sec; + ep[0] = 1970+days/1461*4+mon/12; ep[1] = mon%12+1; ep[2] = day+1; + ep[3] = sec/3600; ep[4] = sec%3600/60; ep[5] = sec%60+t.sec; } @@ -1336,11 +1214,11 @@ void time2epoch(gtime_t t, double *ep) *-----------------------------------------------------------------------------*/ gtime_t gpst2time(int week, double sec) { - gtime_t t=epoch2time(gpst0); + gtime_t t = epoch2time(gpst0); - if (sec<-1e9||1e9tm_year+1900; ep[1]=tt->tm_mon+1; ep[2]=tt->tm_mday; - ep[3]=tt->tm_hour; ep[4]=tt->tm_min; ep[5]=tt->tm_sec+tv.tv_usec*1e-6; + ep[0] = tt->tm_year+1900; ep[1] = tt->tm_mon+1; ep[2] = tt->tm_mday; + ep[3] = tt->tm_hour; ep[4] = tt->tm_min; ep[5] = tt->tm_sec+tv.tv_usec*1e-6; } - time=epoch2time(ep); + time = epoch2time(ep); #ifdef CPUTIME_IN_GPST /* cputime operated in gpst */ - time=gpst2utc(time); + time = gpst2utc(time); #endif return time; } @@ -1498,18 +1376,18 @@ void timeset(gtime_t t) /* read leap seconds table by text -------------------------------------------*/ int read_leaps_text(FILE *fp) { - char buff[256],*p; - int i,n=0,ep[6],ls; + char buff[256], *p; + int i, n = 0, ep[6], ls; rewind(fp); - while (fgets(buff,sizeof(buff),fp)&&n=13) continue; - ls[n][0]=y; - ls[n][1]=m; - ls[n][2]=d; - ls[n++][6]=(char)(19.0-tai_utc); + for (m = 1; m <= 12; m++) if (!strcmp(months[m-1], month)) break; + if (m >= 13) continue; + ls[n][0] = y; + ls[n][1] = m; + ls[n][2] = d; + ls[n++][6] = (char)(19.0-tai_utc); } - for (i=0;i0;i++) { - tu=timeadd(t,leaps[i][6]); - if (timediff(tu,epoch2time(leaps[i]))>=0.0) return tu; - } + for (i = 0; leaps[i][0]>0; i++) + { + tu = timeadd(t, leaps[i][6]); + if (timediff(tu, epoch2time(leaps[i])) >= 0.0) return tu; + } return t; } @@ -1606,9 +1485,10 @@ gtime_t utc2gpst(gtime_t t) { int i; - for (i=0;leaps[i][0]>0;i++) { - if (timediff(t,epoch2time(leaps[i]))>=0.0) return timeadd(t,-leaps[i][6]); - } + for (i = 0; leaps[i][0]>0; i++) + { + if (timediff(t, epoch2time(leaps[i])) >= 0.0) return timeadd(t, -leaps[i][6]); + } return t; } @@ -1623,7 +1503,7 @@ gtime_t utc2gpst(gtime_t t) *-----------------------------------------------------------------------------*/ gtime_t gpst2bdt(gtime_t t) { - return timeadd(t,-14.0); + return timeadd(t, -14.0); } @@ -1635,18 +1515,18 @@ gtime_t gpst2bdt(gtime_t t) *-----------------------------------------------------------------------------*/ gtime_t bdt2gpst(gtime_t t) { - return timeadd(t,14.0); + return timeadd(t, 14.0); } /* time to day and sec -------------------------------------------------------*/ double time2sec(gtime_t time, gtime_t *day) { - double ep[6],sec; - time2epoch(time,ep); - sec=ep[3]*3600.0+ep[4]*60.0+ep[5]; - ep[3]=ep[4]=ep[5]=0.0; - *day=epoch2time(ep); + double ep[6], sec; + time2epoch(time, ep); + sec = ep[3]*3600.0+ep[4]*60.0+ep[5]; + ep[3] = ep[4] = ep[5] = 0.0; + *day = epoch2time(ep); return sec; } @@ -1659,18 +1539,18 @@ double time2sec(gtime_t time, gtime_t *day) *-----------------------------------------------------------------------------*/ double utc2gmst(gtime_t t, double ut1_utc) { - const double ep2000[]={2000,1,1,12,0,0}; - gtime_t tut,tut0; - double ut,t1,t2,t3,gmst0,gmst; + const double ep2000[] = {2000, 1, 1, 12, 0, 0}; + gtime_t tut, tut0; + double ut, t1, t2, t3, gmst0, gmst; - tut=timeadd(t,ut1_utc); - ut=time2sec(tut,&tut0); - t1=timediff(tut0,epoch2time(ep2000))/86400.0/36525.0; - t2=t1*t1; t3=t2*t1; - gmst0=24110.54841+8640184.812866*t1+0.093104*t2-6.2E-6*t3; - gmst=gmst0+1.002737909350795*ut; + tut = timeadd(t, ut1_utc); + ut = time2sec(tut, &tut0); + t1 = timediff(tut0, epoch2time(ep2000))/86400.0/36525.0; + t2 = t1*t1; t3 = t2*t1; + gmst0 = 24110.54841+8640184.812866*t1+0.093104*t2-6.2E-6*t3; + gmst = gmst0+1.002737909350795*ut; - return fmod(gmst,86400.0)*PI/43200.0; /* 0 <= gmst <= 2*PI */ + return fmod(gmst, 86400.0)*PI/43200.0; /* 0 <= gmst <= 2*PI */ } @@ -1685,11 +1565,11 @@ void time2str(gtime_t t, char *s, int n) { double ep[6]; - if (n<0) n=0; else if (n>12) n=12; - if (1.0-t.sec<0.5/pow(10.0,n)) {t.time++; t.sec=0.0;}; - time2epoch(t,ep); - sprintf(s,"%04.0f/%02.0f/%02.0f %02.0f:%02.0f:%0*.*f",ep[0],ep[1],ep[2], - ep[3],ep[4],n<=0?2:n+3,n<=0?0:n,ep[5]); + if (n<0) n = 0; else if (n>12) n = 12; + if (1.0-t.sec<0.5/pow(10.0, n)) {t.time++; t.sec = 0.0;}; + time2epoch(t, ep); + sprintf(s, "%04.0f/%02.0f/%02.0f %02.0f:%02.0f:%0*.*f", ep[0], ep[1], ep[2], + ep[3], ep[4], n <= 0 ? 2:n+3, n <= 0 ? 0:n, ep[5]); } @@ -1703,7 +1583,7 @@ void time2str(gtime_t t, char *s, int n) char *time_str(gtime_t t, int n) { static char buff[64]; - time2str(t,buff,n); + time2str(t, buff, n); return buff; } @@ -1717,9 +1597,9 @@ double time2doy(gtime_t t) { double ep[6]; - time2epoch(t,ep); - ep[1]=ep[2]=1.0; ep[3]=ep[4]=ep[5]=0.0; - return timediff(t,epoch2time(ep))/86400.0+1.0; + time2epoch(t, ep); + ep[1] = ep[2] = 1.0; ep[3] = ep[4] = ep[5] = 0.0; + return timediff(t, epoch2time(ep))/86400.0+1.0; } @@ -1731,8 +1611,8 @@ double time2doy(gtime_t t) int adjgpsweek(int week) { int w; - (void)time2gpst(utc2gpst(timeget()),&w); - if (w<1560) w=1560; /* use 2009/12/1 if time is earlier than 2009/12/1 */ + (void)time2gpst(utc2gpst(timeget()), &w); + if (w<1560) w = 1560; /* use 2009/12/1 if time is earlier than 2009/12/1 */ return week+(w-week+512)/1024*1024; } @@ -1744,22 +1624,22 @@ int adjgpsweek(int week) *-----------------------------------------------------------------------------*/ unsigned int tickget(void) { - struct timespec tp={}; - struct timeval tv={}; + struct timespec tp = {}; + struct timeval tv = {}; #ifdef CLOCK_MONOTONIC_RAW /* linux kernel > 2.6.28 */ - if (!clock_gettime(CLOCK_MONOTONIC_RAW,&tp)) + if (!clock_gettime(CLOCK_MONOTONIC_RAW, &tp)) { return tp.tv_sec*1000u+tp.tv_nsec/1000000u; } else { - gettimeofday(&tv,NULL); + gettimeofday(&tv, NULL); return tv.tv_sec*1000u+tv.tv_usec/1000u; } #else - gettimeofday(&tv,NULL); + gettimeofday(&tv, NULL); return tv.tv_sec*1000u+tv.tv_usec/1000u; #endif } @@ -1773,35 +1653,35 @@ unsigned int tickget(void) void sleepms(int ms) { struct timespec ts; - if (ms<=0) return; - ts.tv_sec=(time_t)(ms/1000); - ts.tv_nsec=(long)(ms%1000*1000000); - nanosleep(&ts,NULL); + if (ms <= 0) return; + ts.tv_sec = (time_t)(ms/1000); + ts.tv_nsec = (long)(ms%1000*1000000); + nanosleep(&ts, NULL); } /* convert degree to deg-min-sec ----------------------------------------------- * convert degree to degree-minute-second * args : double deg I degree - * double *dms O degree-minute-second {deg,min,sec} + * double *dms O degree-minute-second {deg, min, sec} * int ndec I number of decimals of second * return : none *-----------------------------------------------------------------------------*/ void deg2dms(double deg, double *dms, int ndec) { - double sign=deg<0.0?-1.0:1.0,a=fabs(deg); - double unit=pow(0.1,ndec); - dms[0]=floor(a); a=(a-dms[0])*60.0; - dms[1]=floor(a); a=(a-dms[1])*60.0; - dms[2]=floor(a/unit+0.5)*unit; - if (dms[2]>=60.0) + double sign = deg<0.0 ? -1.0 : 1.0, a = fabs(deg); + double unit = pow(0.1, ndec); + dms[0] = floor(a); a = (a-dms[0])*60.0; + dms[1] = floor(a); a = (a-dms[1])*60.0; + dms[2] = floor(a/unit+0.5)*unit; + if (dms[2] >= 60.0) { - dms[2]=0.0; - dms[1]+=1.0; - if (dms[1]>=60.0) + dms[2] = 0.0; + dms[1] += 1.0; + if (dms[1] >= 60.0) { - dms[1]=0.0; - dms[0]+=1.0; + dms[1] = 0.0; + dms[0] += 1.0; } } dms[0]*=sign; @@ -1815,7 +1695,7 @@ void deg2dms(double deg, double *dms, int ndec) *-----------------------------------------------------------------------------*/ double dms2deg(const double *dms) { - double sign=dms[0]<0.0?-1.0:1.0; + double sign = dms[0]<0.0 ? -1.0 : 1.0; return sign*(fabs(dms[0])+dms[1]/60.0+dms[2]/3600.0); } @@ -1829,36 +1709,36 @@ double dms2deg(const double *dms) *-----------------------------------------------------------------------------*/ void ecef2pos(const double *r, double *pos) { - double e2=FE_WGS84*(2.0-FE_WGS84),r2=dot(r,r,2),z,zk,v=RE_WGS84,sinp; + double e2 = FE_WGS84*(2.0-FE_WGS84), r2 = dot(r, r, 2), z, zk, v = RE_WGS84, sinp; - for (z=r[2],zk=0.0;fabs(z-zk)>=1e-4;) + for (z = r[2], zk = 0.0; fabs(z-zk) >= 1e-4;) { - zk=z; - sinp=z/sqrt(r2+z*z); - v=RE_WGS84/sqrt(1.0-e2*sinp*sinp); - z=r[2]+v*e2*sinp; + zk = z; + sinp = z/sqrt(r2+z*z); + v = RE_WGS84/sqrt(1.0-e2*sinp*sinp); + z = r[2]+v*e2*sinp; } - pos[0]=r2>1e-12?atan(z/sqrt(r2)):(r[2]>0.0?PI/2.0:-PI/2.0); - pos[1]=r2>1e-12?atan2(r[1],r[0]):0.0; - pos[2]=sqrt(r2+z*z)-v; + pos[0] = r2>1e-12 ? atan(z/sqrt(r2)) : (r[2]>0.0 ? PI/2.0 : -PI/2.0); + pos[1] = r2>1e-12 ? atan2(r[1], r[0]) : 0.0; + pos[2] = sqrt(r2+z*z)-v; } /* transform geodetic to ecef position ----------------------------------------- * transform geodetic position to ecef position - * args : double *pos I geodetic position {lat,lon,h} (rad,m) + * args : double *pos I geodetic position {lat, lon,h} (rad,m) * double *r O ecef position {x,y,z} (m) * return : none * notes : WGS84, ellipsoidal height *-----------------------------------------------------------------------------*/ void pos2ecef(const double *pos, double *r) { - double sinp=sin(pos[0]),cosp=cos(pos[0]),sinl=sin(pos[1]),cosl=cos(pos[1]); - double e2=FE_WGS84*(2.0-FE_WGS84),v=RE_WGS84/sqrt(1.0-e2*sinp*sinp); + double sinp = sin(pos[0]), cosp = cos(pos[0]), sinl = sin(pos[1]), cosl = cos(pos[1]); + double e2 = FE_WGS84*(2.0-FE_WGS84), v = RE_WGS84/sqrt(1.0-e2*sinp*sinp); - r[0]=(v+pos[2])*cosp*cosl; - r[1]=(v+pos[2])*cosp*sinl; - r[2]=(v*(1.0-e2)+pos[2])*sinp; + r[0] = (v+pos[2])*cosp*cosl; + r[1] = (v+pos[2])*cosp*sinl; + r[2] = (v*(1.0-e2)+pos[2])*sinp; } @@ -1871,11 +1751,11 @@ void pos2ecef(const double *pos, double *r) *-----------------------------------------------------------------------------*/ void xyz2enu(const double *pos, double *E) { - double sinp=sin(pos[0]),cosp=cos(pos[0]),sinl=sin(pos[1]),cosl=cos(pos[1]); + double sinp = sin(pos[0]), cosp = cos(pos[0]), sinl = sin(pos[1]), cosl = cos(pos[1]); - E[0]=-sinl; E[3]=cosl; E[6]=0.0; - E[1]=-sinp*cosl; E[4]=-sinp*sinl; E[7]=cosp; - E[2]=cosp*cosl; E[5]=cosp*sinl; E[8]=sinp; + E[0] = -sinl; E[3] = cosl; E[6] = 0.0; + E[1] = -sinp*cosl; E[4] = -sinp*sinl; E[7] = cosp; + E[2] = cosp*cosl; E[5] = cosp*sinl; E[8] = sinp; } @@ -1891,7 +1771,7 @@ void ecef2enu(const double *pos, const double *r, double *e) double E[9]; xyz2enu(pos,E); - matmul("NN",3,1,3,1.0,E,r,0.0,e); + matmul("NN", 3, 1, 3, 1.0, E, r, 0.0, e); } @@ -1906,25 +1786,25 @@ void enu2ecef(const double *pos, const double *e, double *r) { double E[9]; - xyz2enu(pos,E); - matmul("TN",3,1,3,1.0,E,e,0.0,r); + xyz2enu(pos, E); + matmul("TN", 3, 1, 3, 1.0, E, e, 0.0, r); } /* transform covariance to local tangental coordinate -------------------------- * transform ecef covariance to local tangental coordinate - * args : double *pos I geodetic position {lat,lon} (rad) + * args : double *pos I geodetic position {lat, lon} (rad) * double *P I covariance in ecef coordinate * double *Q O covariance in local tangental coordinate * return : none *-----------------------------------------------------------------------------*/ void covenu(const double *pos, const double *P, double *Q) { - double E[9],EP[9]; + double E[9], EP[9]; - xyz2enu(pos,E); - matmul("NN",3,3,3,1.0,E,P,0.0,EP); - matmul("NT",3,3,3,1.0,EP,E,0.0,Q); + xyz2enu(pos, E); + matmul("NN", 3, 3, 3, 1.0, E, P, 0.0, EP); + matmul("NT", 3, 3, 3, 1.0, EP, E, 0.0, Q); } @@ -1937,18 +1817,18 @@ void covenu(const double *pos, const double *P, double *Q) *-----------------------------------------------------------------------------*/ void covecef(const double *pos, const double *Q, double *P) { - double E[9],EQ[9]; + double E[9], EQ[9]; - xyz2enu(pos,E); - matmul("TN",3,3,3,1.0,E,Q,0.0,EQ); - matmul("NN",3,3,3,1.0,EQ,E,0.0,P); + xyz2enu(pos, E); + matmul("TN", 3, 3, 3, 1.0, E, Q, 0.0, EQ); + matmul("NN", 3, 3, 3, 1.0, EQ, E, 0.0, P); } /* astronomical arguments: f={l,l',F,D,OMG} (rad) ----------------------------*/ void ast_args(double t, double *f) { - static const double fc[][5]={ /* coefficients for iau 1980 nutation */ + static const double fc[][5] = { /* coefficients for iau 1980 nutation */ { 134.96340251, 1717915923.2178, 31.8792, 0.051635, -0.00024470}, { 357.52910918, 129596581.0481, -0.5532, 0.000136, -0.00001149}, { 93.27209062, 1739527262.8478, -12.7512, -0.001037, 0.00000417}, @@ -1956,14 +1836,14 @@ void ast_args(double t, double *f) { 125.04455501, -6962890.2665, 7.4722, 0.007702, -0.00005939} }; double tt[4]; - int i,j; + int i, j; - for (tt[0]=t,i=1;i<4;i++) tt[i]=tt[i-1]*t; - for (i=0;i<5;i++) + for (tt[0] = t, i = 1; i<4; i++) tt[i] = tt[i-1]*t; + for (i = 0; i<5; i++) { - f[i]=fc[i][0]*3600.0; - for (j=0;j<4;j++) f[i]+=fc[i][j+1]*tt[j]; - f[i]=fmod(f[i]*AS2R,2.0*PI); + f[i] = fc[i][0]*3600.0; + for (j = 0; j<4; j++) f[i] += fc[i][j+1]*tt[j]; + f[i] = fmod(f[i]*AS2R, 2.0*PI); } } @@ -1971,7 +1851,7 @@ void ast_args(double t, double *f) /* iau 1980 nutation ---------------------------------------------------------*/ void nut_iau1980(double t, const double *f, double *dpsi, double *deps) { - static const double nut[106][10]={ + static const double nut[106][10] = { { 0, 0, 0, 0, 1, -6798.4, -171996, -174.2, 92025, 8.9}, { 0, 0, 2, -2, 2, 182.6, -13187, -1.6, 5736, -3.1}, { 0, 0, 2, 0, 2, 13.7, -2274, -0.2, 977, -0.5}, @@ -2080,16 +1960,16 @@ void nut_iau1980(double t, const double *f, double *dpsi, double *deps) { 0, 1, 0, 1, 0, 27.3, 1, 0.0, 0, 0.0} }; double ang; - int i,j; + int i, j; - *dpsi=*deps=0.0; + *dpsi = *deps = 0.0; - for (i=0;i<106;i++) + for (i = 0; i<106; i++) { - ang=0.0; - for (j=0;j<5;j++) ang+=nut[i][j]*f[j]; - *dpsi+=(nut[i][6]+nut[i][7]*t)*sin(ang); - *deps+=(nut[i][8]+nut[i][9]*t)*cos(ang); + ang = 0.0; + for (j = 0; j<5; j++) ang += nut[i][j]*f[j]; + *dpsi += (nut[i][6]+nut[i][7]*t)*sin(ang); + *deps += (nut[i][8]+nut[i][9]*t)*cos(ang); } *dpsi*=1e-4*AS2R; /* 0.1 mas -> rad */ *deps*=1e-4*AS2R; @@ -2109,67 +1989,67 @@ void nut_iau1980(double t, const double *f, double *dpsi, double *deps) *-----------------------------------------------------------------------------*/ void eci2ecef(gtime_t tutc, const double *erpv, double *U, double *gmst) { - const double ep2000[]={2000,1,1,12,0,0}; + const double ep2000[] = {2000, 1, 1, 12, 0, 0}; static gtime_t tutc_; - static double U_[9],gmst_; + static double U_[9], gmst_; gtime_t tgps; - double eps,ze,th,z,t,t2,t3,dpsi,deps,gast,f[5]; - double R1[9],R2[9],R3[9],R[9],W[9],N[9],P[9],NP[9]; + double eps, ze, th, z, t, t2, t3, dpsi, deps, gast, f[5]; + double R1[9], R2[9], R3[9], R[9], W[9], N[9], P[9], NP[9]; int i; - trace(4,"eci2ecef: tutc=%s\n",time_str(tutc,3)); + trace(4, "eci2ecef: tutc=%s\n", time_str(tutc, 3)); - if (fabs(timediff(tutc,tutc_))<0.01) + if (fabs(timediff(tutc, tutc_))<0.01) { /* read cache */ - for (i=0;i<9;i++) U[i]=U_[i]; - if (gmst) *gmst=gmst_; + for (i = 0; i<9; i++) U[i] = U_[i]; + if (gmst) *gmst = gmst_; return; } - tutc_=tutc; + tutc_ = tutc; /* terrestrial time */ - tgps=utc2gpst(tutc_); - t=(timediff(tgps,epoch2time(ep2000))+19.0+32.184)/86400.0/36525.0; - t2=t*t; t3=t2*t; + tgps = utc2gpst(tutc_); + t = (timediff(tgps, epoch2time(ep2000))+19.0+32.184)/86400.0/36525.0; + t2 = t*t; t3 = t2*t; /* astronomical arguments */ - ast_args(t,f); + ast_args(t, f); /* iau 1976 precession */ - ze=(2306.2181*t+0.30188*t2+0.017998*t3)*AS2R; - th=(2004.3109*t-0.42665*t2-0.041833*t3)*AS2R; - z =(2306.2181*t+1.09468*t2+0.018203*t3)*AS2R; - eps=(84381.448-46.8150*t-0.00059*t2+0.001813*t3)*AS2R; - Rz(-z,R1); Ry(th,R2); Rz(-ze,R3); - matmul("NN",3,3,3,1.0,R1,R2,0.0,R); - matmul("NN",3,3,3,1.0,R, R3,0.0,P); /* P=Rz(-z)*Ry(th)*Rz(-ze) */ + ze = (2306.2181*t+0.30188*t2+0.017998*t3)*AS2R; + th = (2004.3109*t-0.42665*t2-0.041833*t3)*AS2R; + z = (2306.2181*t+1.09468*t2+0.018203*t3)*AS2R; + eps = (84381.448-46.8150*t-0.00059*t2+0.001813*t3)*AS2R; + Rz(-z, R1); Ry(th, R2); Rz(-ze, R3); + matmul("NN", 3, 3, 3, 1.0, R1, R2, 0.0, R); + matmul("NN", 3, 3, 3, 1.0, R, R3, 0.0, P); /* P=Rz(-z)*Ry(th)*Rz(-ze) */ /* iau 1980 nutation */ - nut_iau1980(t,f,&dpsi,&deps); - Rx(-eps-deps,R1); Rz(-dpsi,R2); Rx(eps,R3); - matmul("NN",3,3,3,1.0,R1,R2,0.0,R); - matmul("NN",3,3,3,1.0,R ,R3,0.0,N); /* N=Rx(-eps)*Rz(-dspi)*Rx(eps) */ + nut_iau1980(t, f, &dpsi, &deps); + Rx(-eps-deps, R1); Rz(-dpsi, R2); Rx(eps, R3); + matmul("NN", 3, 3, 3, 1.0, R1, R2, 0.0, R); + matmul("NN", 3, 3, 3, 1.0, R , R3, 0.0, N); /* N=Rx(-eps)*Rz(-dspi)*Rx(eps) */ /* greenwich aparent sidereal time (rad) */ - gmst_=utc2gmst(tutc_,erpv[2]); - gast=gmst_+dpsi*cos(eps); - gast+=(0.00264*sin(f[4])+0.000063*sin(2.0*f[4]))*AS2R; + gmst_ = utc2gmst(tutc_, erpv[2]); + gast = gmst_+dpsi*cos(eps); + gast += (0.00264*sin(f[4])+0.000063*sin(2.0*f[4]))*AS2R; /* eci to ecef transformation matrix */ - Ry(-erpv[0],R1); Rx(-erpv[1],R2); Rz(gast,R3); - matmul("NN",3,3,3,1.0,R1,R2,0.0,W ); - matmul("NN",3,3,3,1.0,W ,R3,0.0,R ); /* W=Ry(-xp)*Rx(-yp) */ - matmul("NN",3,3,3,1.0,N ,P ,0.0,NP); - matmul("NN",3,3,3,1.0,R ,NP,0.0,U_); /* U=W*Rz(gast)*N*P */ + Ry(-erpv[0], R1); Rx(-erpv[1], R2); Rz(gast, R3); + matmul("NN", 3, 3, 3, 1.0, R1, R2, 0.0, W ); + matmul("NN", 3, 3, 3, 1.0, W , R3, 0.0, R ); /* W=Ry(-xp)*Rx(-yp) */ + matmul("NN", 3, 3, 3, 1.0, N , P , 0.0, NP); + matmul("NN", 3, 3, 3, 1.0, R , NP, 0.0, U_); /* U=W*Rz(gast)*N*P */ - for (i=0;i<9;i++) U[i]=U_[i]; - if (gmst) *gmst=gmst_; + for (i = 0; i<9; i++) U[i] = U_[i]; + if (gmst) *gmst = gmst_; - trace(5,"gmst=%.12f gast=%.12f\n",gmst_,gast); - trace(5,"P=\n"); tracemat(5,P,3,3,15,12); - trace(5,"N=\n"); tracemat(5,N,3,3,15,12); - trace(5,"W=\n"); tracemat(5,W,3,3,15,12); - trace(5,"U=\n"); tracemat(5,U,3,3,15,12); + trace(5, "gmst=%.12f gast=%.12f\n", gmst_, gast); + trace(5, "P=\n"); tracemat(5, P, 3, 3, 15, 12); + trace(5, "N=\n"); tracemat(5, N, 3, 3, 15, 12); + trace(5, "W=\n"); tracemat(5, W, 3, 3, 15, 12); + trace(5, "U=\n"); tracemat(5, U, 3, 3, 15, 12); } @@ -2178,32 +2058,32 @@ int decodef(char *p, int n, double *v) { int i; - for (i=0;inmax<=pcvs->n) + if (pcvs->nmax <= pcvs->n) { - pcvs->nmax+=256; - if (!(pcvs_pcv=(pcv_t *)realloc(pcvs->pcv,sizeof(pcv_t)*pcvs->nmax))) + pcvs->nmax += 256; + if (!(pcvs_pcv = (pcv_t *)realloc(pcvs->pcv, sizeof(pcv_t)*pcvs->nmax))) { - trace(1,"addpcv: memory allocation error\n"); - free(pcvs->pcv); pcvs->pcv=NULL; pcvs->n=pcvs->nmax=0; + trace(1, "addpcv: memory allocation error\n"); + free(pcvs->pcv); pcvs->pcv = NULL; pcvs->n = pcvs->nmax = 0; return; } - pcvs->pcv=pcvs_pcv; + pcvs->pcv = pcvs_pcv; } - pcvs->pcv[pcvs->n++]=*pcv; + pcvs->pcv[pcvs->n++] = *pcv; } @@ -2211,48 +2091,48 @@ void addpcv(const pcv_t *pcv, pcvs_t *pcvs) int readngspcv(const char *file, pcvs_t *pcvs) { FILE *fp; - static const pcv_t pcv0={}; + static const pcv_t pcv0 = {}; pcv_t pcv; double neu[3]; - int n=0; + int n = 0; char buff[256]; - if (!(fp=fopen(file,"r"))) + if (!(fp = fopen(file, "r"))) { - trace(2,"ngs pcv file open error: %s\n",file); + trace(2, "ngs pcv file open error: %s\n", file); return 0; } - while (fgets(buff,sizeof(buff),fp)) + while (fgets(buff, sizeof(buff), fp)) { - if (strlen(buff)>=62&&buff[61]=='|') continue; + if (strlen(buff) >= 62 && buff[61] == '|') continue; - if (buff[0]!=' ') n=0; /* start line */ - if (++n==1) + if (buff[0] != ' ') n = 0; /* start line */ + if (++n == 1) { - pcv=pcv0; - strncpy(pcv.type,buff,61); pcv.type[61]='\0'; + pcv = pcv0; + strncpy(pcv.type, buff, 61); pcv.type[61] = '\0'; } - else if (n==2) + else if (n == 2) { - if (decodef(buff,3,neu)<3) continue; - pcv.off[0][0]=neu[1]; - pcv.off[0][1]=neu[0]; - pcv.off[0][2]=neu[2]; + if (decodef(buff, 3, neu)<3) continue; + pcv.off[0][0] = neu[1]; + pcv.off[0][1] = neu[0]; + pcv.off[0][2] = neu[2]; } - else if (n==3) decodef(buff,10,pcv.var[0]); - else if (n==4) decodef(buff,9,pcv.var[0]+10); - else if (n==5) + else if (n == 3) decodef(buff, 10, pcv.var[0]); + else if (n == 4) decodef(buff, 9, pcv.var[0]+10); + else if (n == 5) { - if (decodef(buff,3,neu)<3) continue;; - pcv.off[1][0]=neu[1]; - pcv.off[1][1]=neu[0]; - pcv.off[1][2]=neu[2]; + if (decodef(buff, 3, neu)<3) continue;; + pcv.off[1][0] = neu[1]; + pcv.off[1][1] = neu[0]; + pcv.off[1][2] = neu[2]; } - else if (n==6) decodef(buff,10,pcv.var[1]); - else if (n==7) + else if (n == 6) decodef(buff, 10, pcv.var[1]); + else if (n == 7) { - decodef(buff,9,pcv.var[1]+10); - addpcv(&pcv,pcvs); + decodef(buff, 9, pcv.var[1]+10); + addpcv(&pcv, pcvs); } } fclose(fp); @@ -2265,75 +2145,75 @@ int readngspcv(const char *file, pcvs_t *pcvs) int readantex(const char *file, pcvs_t *pcvs) { FILE *fp; - static const pcv_t pcv0={}; + static const pcv_t pcv0 = {}; pcv_t pcv; double neu[3]; - int i,f,freq=0,state=0,freqs[]={1,2,5,6,7,8,0}; + int i, f, freq = 0, state = 0, freqs[] = {1, 2, 5, 6, 7, 8, 0}; char buff[256]; - trace(3,"readantex: file=%s\n",file); + trace(3, "readantex: file=%s\n", file); - if (!(fp=fopen(file,"r"))) + if (!(fp = fopen(file, "r"))) { - trace(2,"antex pcv file open error: %s\n",file); + trace(2, "antex pcv file open error: %s\n", file); return 0; } - while (fgets(buff,sizeof(buff),fp)) + while (fgets(buff, sizeof(buff), fp)) { - if (strlen(buff)<60||strstr(buff+60,"COMMENT")) continue; + if (strlen(buff)<60 || strstr(buff+60, "COMMENT")) continue; - if (strstr(buff+60,"START OF ANTENNA")) + if (strstr(buff+60, "START OF ANTENNA")) { - pcv=pcv0; - state=1; + pcv = pcv0; + state = 1; } - if (strstr(buff+60,"END OF ANTENNA")) + if (strstr(buff+60, "END OF ANTENNA")) { - addpcv(&pcv,pcvs); - state=0; + addpcv(&pcv, pcvs); + state = 0; } if (!state) continue; - if (strstr(buff+60,"TYPE / SERIAL NO")) + if (strstr(buff+60, "TYPE / SERIAL NO")) { - strncpy(pcv.type,buff ,20); pcv.type[20]='\0'; - strncpy(pcv.code,buff+20,20); pcv.code[20]='\0'; - if (!strncmp(pcv.code+3," ",8)) + strncpy(pcv.type, buff , 20); pcv.type[20] = '\0'; + strncpy(pcv.code, buff+20, 20); pcv.code[20] = '\0'; + if (!strncmp(pcv.code+3, " ", 8)) { - pcv.sat=satid2no(pcv.code); + pcv.sat = satid2no(pcv.code); } } - else if (strstr(buff+60,"VALID FROM")) + else if (strstr(buff+60, "VALID FROM")) { - if (!str2time(buff,0,43,&pcv.ts)) continue; + if (!str2time(buff, 0, 43, &pcv.ts)) continue; } - else if (strstr(buff+60,"VALID UNTIL")) + else if (strstr(buff+60, "VALID UNTIL")) { - if (!str2time(buff,0,43,&pcv.te)) continue; + if (!str2time(buff, 0, 43, &pcv.te)) continue; } - else if (strstr(buff+60,"START OF FREQUENCY")) + else if (strstr(buff+60, "START OF FREQUENCY")) { - if (sscanf(buff+4,"%d",&f)<1) continue; - for (i=0;in;i++) + for (i = 0; in; i++) { - pcv=pcvs->pcv+i; - trace(4,"sat=%2d type=%20s code=%s off=%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n", - pcv->sat,pcv->type,pcv->code,pcv->off[0][0],pcv->off[0][1], - pcv->off[0][2],pcv->off[1][0],pcv->off[1][1],pcv->off[1][2]); + pcv = pcvs->pcv+i; + trace(4, "sat=%2d type=%20s code=%s off=%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n", + pcv->sat, pcv->type, pcv->code, pcv->off[0][0], pcv->off[0][1], + pcv->off[0][2], pcv->off[1][0], pcv->off[1][1], pcv->off[1][2]); } return stat; } @@ -2396,41 +2276,41 @@ pcv_t *searchpcv(int sat, const char *type, gtime_t time, const pcvs_t *pcvs) { pcv_t *pcv; - char buff[MAXANT],*types[2],*p; - int i,j,n=0; + char buff[MAXANT], *types[2], *p; + int i, j, n = 0; - trace(3,"searchpcv: sat=%2d type=%s\n",sat,type); + trace(3, "searchpcv: sat=%2d type=%s\n", sat, type); if (sat) { /* search satellite antenna */ - for (i=0;in;i++) + for (i = 0; in; i++) { - pcv=pcvs->pcv+i; - if (pcv->sat!=sat) continue; - if (pcv->ts.time!=0&&timediff(pcv->ts,time)>0.0) continue; - if (pcv->te.time!=0&&timediff(pcv->te,time)<0.0) continue; + pcv = pcvs->pcv+i; + if (pcv->sat != sat) continue; + if (pcv->ts.time != 0 && timediff(pcv->ts, time)>0.0) continue; + if (pcv->te.time != 0 && timediff(pcv->te, time)<0.0) continue; return pcv; } } else { - strcpy(buff,type); - for (p=strtok(buff," ");p&&n<2;p=strtok(NULL," ")) types[n++]=p; - if (n<=0) return NULL; + strcpy(buff, type); + for (p = strtok(buff, " "); p && n<2; p = strtok(NULL, " ")) types[n++] = p; + if (n <= 0) return NULL; /* search receiver antenna with radome at first */ - for (i=0;in;i++) + for (i = 0; in; i++) { - pcv=pcvs->pcv+i; - for (j=0;jtype,types[j])) break; - if (j>=n) return pcv; + pcv = pcvs->pcv+i; + for (j = 0; jtype, types[j])) break; + if (j >= n) return pcv; } /* search receiver antenna without radome */ - for (i=0;in;i++) + for (i = 0; in; i++) { - pcv=pcvs->pcv+i; - if (strstr(pcv->type,types[0])!=pcv->type) continue; + pcv = pcvs->pcv+i; + if (strstr(pcv->type, types[0]) != pcv->type) continue; - trace(2,"pcv without radome is used type=%s\n",type); + trace(2, "pcv without radome is used type=%s\n", type); return pcv; } } @@ -2452,33 +2332,33 @@ void readpos(const char *file, const char *rcv, double *pos) static double poss[2048][3]; static char stas[2048][16]; FILE *fp; - int i,j,len,np=0; - char buff[256],str[256]; + int i, j, len, np = 0; + char buff[256], str[256]; - trace(3,"readpos: file=%s\n",file); + trace(3, "readpos: file=%s\n", file); - if (!(fp=fopen(file,"r"))) + if (!(fp = fopen(file, "r"))) { - fprintf(stderr,"reference position file open error : %s\n",file); + fprintf(stderr, "reference position file open error : %s\n", file); return; } - while (np<2048&&fgets(buff,sizeof(buff),fp)) + while (np<2048 && fgets(buff, sizeof(buff), fp)) { - if (buff[0]=='%'||buff[0]=='#') continue; - if (sscanf(buff,"%lf %lf %lf %s",&poss[np][0],&poss[np][1],&poss[np][2], + if (buff[0] == '%' || buff[0] == '#') continue; + if (sscanf(buff, "%lf %lf %lf %s", &poss[np][0], &poss[np][1], &poss[np][2], str)<4) continue; - strncpy(stas[np],str,15); stas[np++][15]='\0'; + strncpy(stas[np], str, 15); stas[np++][15] = '\0'; } fclose(fp); - len=(int)strlen(rcv); - for (i=0;in>=erp->nmax) + if (erp->n >= erp->nmax) { - erp->nmax=erp->nmax<=0?128:erp->nmax*2; - erp_data=(erpd_t *)realloc(erp->data,sizeof(erpd_t)*erp->nmax); + erp->nmax = erp->nmax <= 0 ? 128 : erp->nmax*2; + erp_data = (erpd_t *)realloc(erp->data, sizeof(erpd_t)*erp->nmax); if (!erp_data) { - free(erp->data); erp->data=NULL; erp->n=erp->nmax=0; + free(erp->data); erp->data = NULL; erp->n = erp->nmax = 0; fclose(fp); return 0; } - erp->data=erp_data; + erp->data = erp_data; } - erp->data[erp->n].mjd=v[0]; - erp->data[erp->n].xp=v[1]*1e-6*AS2R; - erp->data[erp->n].yp=v[2]*1e-6*AS2R; - erp->data[erp->n].ut1_utc=v[3]*1e-7; - erp->data[erp->n].lod=v[4]*1e-7; - erp->data[erp->n].xpr=v[12]*1e-6*AS2R; - erp->data[erp->n++].ypr=v[13]*1e-6*AS2R; + erp->data[erp->n].mjd = v[0]; + erp->data[erp->n].xp = v[1]*1e-6*AS2R; + erp->data[erp->n].yp = v[2]*1e-6*AS2R; + erp->data[erp->n].ut1_utc = v[3]*1e-7; + erp->data[erp->n].lod = v[4]*1e-7; + erp->data[erp->n].xpr = v[12]*1e-6*AS2R; + erp->data[erp->n++].ypr = v[13]*1e-6*AS2R; } fclose(fp); return 1; @@ -2604,60 +2484,60 @@ int readerp(const char *file, erp_t *erp) *-----------------------------------------------------------------------------*/ int geterp(const erp_t *erp, gtime_t time, double *erpv) { - const double ep[]={2000,1,1,12,0,0}; - double mjd,day,a; - int i,j,k; + const double ep[] = {2000, 1, 1, 12, 0, 0}; + double mjd, day, a; + int i, j, k; - trace(4,"geterp:\n"); + trace(4, "geterp:\n"); - if (erp->n<=0) return 0; + if (erp->n <= 0) return 0; - mjd=51544.5+(timediff(gpst2utc(time),epoch2time(ep)))/86400.0; + mjd = 51544.5+(timediff(gpst2utc(time), epoch2time(ep)))/86400.0; - if (mjd<=erp->data[0].mjd) + if (mjd <= erp->data[0].mjd) { - day=mjd-erp->data[0].mjd; - erpv[0]=erp->data[0].xp +erp->data[0].xpr*day; - erpv[1]=erp->data[0].yp +erp->data[0].ypr*day; - erpv[2]=erp->data[0].ut1_utc-erp->data[0].lod*day; - erpv[3]=erp->data[0].lod; + day = mjd-erp->data[0].mjd; + erpv[0] = erp->data[0].xp +erp->data[0].xpr*day; + erpv[1] = erp->data[0].yp +erp->data[0].ypr*day; + erpv[2] = erp->data[0].ut1_utc-erp->data[0].lod*day; + erpv[3] = erp->data[0].lod; return 1; } - if (mjd>=erp->data[erp->n-1].mjd) + if (mjd >= erp->data[erp->n-1].mjd) { - day=mjd-erp->data[erp->n-1].mjd; - erpv[0]=erp->data[erp->n-1].xp +erp->data[erp->n-1].xpr*day; - erpv[1]=erp->data[erp->n-1].yp +erp->data[erp->n-1].ypr*day; - erpv[2]=erp->data[erp->n-1].ut1_utc-erp->data[erp->n-1].lod*day; - erpv[3]=erp->data[erp->n-1].lod; + day = mjd-erp->data[erp->n-1].mjd; + erpv[0] = erp->data[erp->n-1].xp +erp->data[erp->n-1].xpr*day; + erpv[1] = erp->data[erp->n-1].yp +erp->data[erp->n-1].ypr*day; + erpv[2] = erp->data[erp->n-1].ut1_utc-erp->data[erp->n-1].lod*day; + erpv[3] = erp->data[erp->n-1].lod; return 1; } - for (j=0,k=erp->n-1;jn-1; jdata[i].mjd) k=i; else j=i; + i = (j+k)/2; + if (mjddata[i].mjd) k = i; else j = i; } - if (erp->data[j].mjd==erp->data[j+1].mjd) + if (erp->data[j].mjd == erp->data[j+1].mjd) { - a=0.5; + a = 0.5; } else { - a=(mjd-erp->data[j].mjd)/(erp->data[j+1].mjd-erp->data[j].mjd); + a = (mjd-erp->data[j].mjd)/(erp->data[j+1].mjd-erp->data[j].mjd); } - erpv[0]=(1.0-a)*erp->data[j].xp +a*erp->data[j+1].xp; - erpv[1]=(1.0-a)*erp->data[j].yp +a*erp->data[j+1].yp; - erpv[2]=(1.0-a)*erp->data[j].ut1_utc+a*erp->data[j+1].ut1_utc; - erpv[3]=(1.0-a)*erp->data[j].lod +a*erp->data[j+1].lod; + erpv[0] = (1.0-a)*erp->data[j].xp +a*erp->data[j+1].xp; + erpv[1] = (1.0-a)*erp->data[j].yp +a*erp->data[j+1].yp; + erpv[2] = (1.0-a)*erp->data[j].ut1_utc+a*erp->data[j+1].ut1_utc; + erpv[3] = (1.0-a)*erp->data[j].lod +a*erp->data[j+1].lod; return 1; } /* compare ephemeris ---------------------------------------------------------*/ int cmpeph(const void *p1, const void *p2) { - eph_t *q1=(eph_t *)p1,*q2=(eph_t *)p2; - return q1->ttr.time!=q2->ttr.time?(int)(q1->ttr.time-q2->ttr.time): - (q1->toe.time!=q2->toe.time?(int)(q1->toe.time-q2->toe.time): + eph_t *q1 = (eph_t *)p1, *q2 = (eph_t *)p2; + return q1->ttr.time != q2->ttr.time ? (int)(q1->ttr.time-q2->ttr.time) : + (q1->toe.time != q2->toe.time ? (int)(q1->toe.time-q2->toe.time) : q1->sat-q2->sat); } @@ -2668,41 +2548,41 @@ void uniqeph(nav_t *nav) eph_t *nav_eph; int i,j; - trace(3,"uniqeph: n=%d\n",nav->n); + trace(3, "uniqeph: n=%d\n", nav->n); - if (nav->n<=0) return; + if (nav->n <= 0) return; - qsort(nav->eph,nav->n,sizeof(eph_t),cmpeph); + qsort(nav->eph, nav->n, sizeof(eph_t), cmpeph); - for (i=1,j=0;in;i++) + for (i = 1, j = 0; in; i++) { - if (nav->eph[i].sat!=nav->eph[j].sat|| - nav->eph[i].iode!=nav->eph[j].iode) + if (nav->eph[i].sat != nav->eph[j].sat || + nav->eph[i].iode != nav->eph[j].iode) { - nav->eph[++j]=nav->eph[i]; + nav->eph[++j] = nav->eph[i]; } } - nav->n=j+1; + nav->n = j+1; - if (!(nav_eph=(eph_t *)realloc(nav->eph,sizeof(eph_t)*nav->n))) + if (!(nav_eph = (eph_t *)realloc(nav->eph, sizeof(eph_t)*nav->n))) { - trace(1,"uniqeph malloc error n=%d\n",nav->n); - free(nav->eph); nav->eph=NULL; nav->n=nav->nmax=0; + trace(1, "uniqeph malloc error n=%d\n", nav->n); + free(nav->eph); nav->eph = NULL; nav->n = nav->nmax = 0; return; } - nav->eph=nav_eph; - nav->nmax=nav->n; + nav->eph = nav_eph; + nav->nmax = nav->n; - trace(4,"uniqeph: n=%d\n",nav->n); + trace(4, "uniqeph: n=%d\n", nav->n); } /* compare glonass ephemeris -------------------------------------------------*/ int cmpgeph(const void *p1, const void *p2) { - geph_t *q1=(geph_t *)p1,*q2=(geph_t *)p2; - return q1->tof.time!=q2->tof.time?(int)(q1->tof.time-q2->tof.time): - (q1->toe.time!=q2->toe.time?(int)(q1->toe.time-q2->toe.time): + geph_t *q1 = (geph_t *)p1, *q2 = (geph_t *)p2; + return q1->tof.time != q2->tof.time ? (int)(q1->tof.time-q2->tof.time) : + (q1->toe.time != q2->toe.time ? (int)(q1->toe.time-q2->toe.time) : q1->sat-q2->sat); } @@ -2711,44 +2591,44 @@ int cmpgeph(const void *p1, const void *p2) void uniqgeph(nav_t *nav) { geph_t *nav_geph; - int i,j; + int i, j; - trace(3,"uniqgeph: ng=%d\n",nav->ng); + trace(3, "uniqgeph: ng=%d\n", nav->ng); - if (nav->ng<=0) return; + if (nav->ng <= 0) return; - qsort(nav->geph,nav->ng,sizeof(geph_t),cmpgeph); + qsort(nav->geph, nav->ng, sizeof(geph_t), cmpgeph); - for (i=j=0;ing;i++) + for (i = j = 0; ing; i++) { - if (nav->geph[i].sat!=nav->geph[j].sat|| - nav->geph[i].toe.time!=nav->geph[j].toe.time|| - nav->geph[i].svh!=nav->geph[j].svh) + if (nav->geph[i].sat != nav->geph[j].sat || + nav->geph[i].toe.time != nav->geph[j].toe.time || + nav->geph[i].svh != nav->geph[j].svh) { - nav->geph[++j]=nav->geph[i]; + nav->geph[++j] = nav->geph[i]; } } - nav->ng=j+1; + nav->ng = j+1; - if (!(nav_geph=(geph_t *)realloc(nav->geph,sizeof(geph_t)*nav->ng))) + if (!(nav_geph = (geph_t *)realloc(nav->geph, sizeof(geph_t)*nav->ng))) { - trace(1,"uniqgeph malloc error ng=%d\n",nav->ng); - free(nav->geph); nav->geph=NULL; nav->ng=nav->ngmax=0; + trace(1, "uniqgeph malloc error ng=%d\n", nav->ng); + free(nav->geph); nav->geph = NULL; nav->ng = nav->ngmax = 0; return; } - nav->geph=nav_geph; - nav->ngmax=nav->ng; + nav->geph = nav_geph; + nav->ngmax = nav->ng; - trace(4,"uniqgeph: ng=%d\n",nav->ng); + trace(4, "uniqgeph: ng=%d\n", nav->ng); } /* compare sbas ephemeris ----------------------------------------------------*/ int cmpseph(const void *p1, const void *p2) { - seph_t *q1=(seph_t *)p1,*q2=(seph_t *)p2; - return q1->tof.time!=q2->tof.time?(int)(q1->tof.time-q2->tof.time): - (q1->t0.time!=q2->t0.time?(int)(q1->t0.time-q2->t0.time): + seph_t *q1 = (seph_t *)p1, *q2 = (seph_t *)p2; + return q1->tof.time != q2->tof.time ? (int)(q1->tof.time-q2->tof.time) : + (q1->t0.time != q2->t0.time ? (int)(q1->t0.time-q2->t0.time) : q1->sat-q2->sat); } @@ -2757,34 +2637,34 @@ int cmpseph(const void *p1, const void *p2) void uniqseph(nav_t *nav) { seph_t *nav_seph; - int i,j; + int i, j; - trace(3,"uniqseph: ns=%d\n",nav->ns); + trace(3, "uniqseph: ns=%d\n", nav->ns); - if (nav->ns<=0) return; + if (nav->ns <= 0) return; - qsort(nav->seph,nav->ns,sizeof(seph_t),cmpseph); + qsort(nav->seph, nav->ns, sizeof(seph_t), cmpseph); - for (i=j=0;ins;i++) + for (i = j = 0; ins; i++) { - if (nav->seph[i].sat!=nav->seph[j].sat|| - nav->seph[i].t0.time!=nav->seph[j].t0.time) + if (nav->seph[i].sat != nav->seph[j].sat || + nav->seph[i].t0.time != nav->seph[j].t0.time) { - nav->seph[++j]=nav->seph[i]; + nav->seph[++j] = nav->seph[i]; } } - nav->ns=j+1; + nav->ns = j+1; - if (!(nav_seph=(seph_t *)realloc(nav->seph,sizeof(seph_t)*nav->ns))) + if (!(nav_seph = (seph_t *)realloc(nav->seph, sizeof(seph_t)*nav->ns))) { - trace(1,"uniqseph malloc error ns=%d\n",nav->ns); - free(nav->seph); nav->seph=NULL; nav->ns=nav->nsmax=0; + trace(1, "uniqseph malloc error ns=%d\n", nav->ns); + free(nav->seph); nav->seph = NULL; nav->ns = nav->nsmax = 0; return; } - nav->seph=nav_seph; - nav->nsmax=nav->ns; + nav->seph = nav_seph; + nav->nsmax = nav->ns; - trace(4,"uniqseph: ns=%d\n",nav->ns); + trace(4, "uniqseph: ns=%d\n", nav->ns); } @@ -2795,9 +2675,9 @@ void uniqseph(nav_t *nav) *-----------------------------------------------------------------------------*/ void uniqnav(nav_t *nav) { - int i,j; + int i, j; - trace(3,"uniqnav: neph=%d ngeph=%d nseph=%d\n",nav->n,nav->ng,nav->ns); + trace(3, "uniqnav: neph=%d ngeph=%d nseph=%d\n", nav->n, nav->ng, nav->ns); /* unique ephemeris */ uniqeph (nav); @@ -2805,9 +2685,9 @@ void uniqnav(nav_t *nav) uniqseph(nav); /* update carrier wave length */ - for (i=0;ilam[i][j]=satwavelen(i+1,j,nav); + nav->lam[i][j] = satwavelen(i+1, j, nav); } } @@ -2815,10 +2695,10 @@ void uniqnav(nav_t *nav) /* compare observation data -------------------------------------------------*/ int cmpobs(const void *p1, const void *p2) { - obsd_t *q1=(obsd_t *)p1,*q2=(obsd_t *)p2; - double tt=timediff(q1->time,q2->time); - if (fabs(tt)>DTTOL) return tt<0?-1:1; - if (q1->rcv!=q2->rcv) return (int)q1->rcv-(int)q2->rcv; + obsd_t *q1 = (obsd_t *)p1, *q2 = (obsd_t *)p2; + double tt = timediff(q1->time, q2->time); + if (fabs(tt)>DTTOL) return tt<0 ? -1 : 1; + if (q1->rcv != q2->rcv) return (int)q1->rcv-(int)q2->rcv; return (int)q1->sat-(int)q2->sat; } @@ -2830,31 +2710,31 @@ int cmpobs(const void *p1, const void *p2) *-----------------------------------------------------------------------------*/ int sortobs(obs_t *obs) { - int i,j,n; + int i, j, n; - trace(3,"sortobs: nobs=%d\n",obs->n); + trace(3, "sortobs: nobs=%d\n", obs->n); - if (obs->n<=0) return 0; + if (obs->n <= 0) return 0; - qsort(obs->data,obs->n,sizeof(obsd_t),cmpobs); + qsort(obs->data, obs->n, sizeof(obsd_t), cmpobs); /* delete duplicated data */ - for (i=j=0;in;i++) + for (i = j = 0; in; i++) { - if (obs->data[i].sat!=obs->data[j].sat|| - obs->data[i].rcv!=obs->data[j].rcv|| - timediff(obs->data[i].time,obs->data[j].time)!=0.0) + if (obs->data[i].sat != obs->data[j].sat || + obs->data[i].rcv != obs->data[j].rcv || + timediff(obs->data[i].time, obs->data[j].time) != 0.0) { - obs->data[++j]=obs->data[i]; + obs->data[++j] = obs->data[i]; } } - obs->n=j+1; + obs->n = j+1; - for (i=n=0;in;i=j,n++) + for (i = n = 0; in; i = j, n++) { - for (j=i+1;jn;j++) + for (j = i+1; jn; j++) { - if (timediff(obs->data[j].time,obs->data[i].time)>DTTOL) break; + if (timediff(obs->data[j].time, obs->data[i].time)>DTTOL) break; } } return n; @@ -2871,9 +2751,9 @@ int sortobs(obs_t *obs) *-----------------------------------------------------------------------------*/ int screent(gtime_t time, gtime_t ts, gtime_t te, double tint) { - return (tint<=0.0||fmod(time2gpst(time,NULL)+DTTOL,tint)<=DTTOL*2.0)&& - (ts.time==0||timediff(time,ts)>=-DTTOL)&& - (te.time==0||timediff(time,te)< DTTOL); + return (tint <= 0.0 || fmod(time2gpst(time, NULL)+DTTOL, tint) <= DTTOL*2.0) && + (ts.time == 0 || timediff(time, ts) >= -DTTOL) && + (te.time == 0 || timediff(time, te)< DTTOL); } @@ -2886,70 +2766,70 @@ int screent(gtime_t time, gtime_t ts, gtime_t te, double tint) int readnav(const char *file, nav_t *nav) { FILE *fp; - eph_t eph0={}; - geph_t geph0={}; - char buff[4096],*p; - long toe_time,tof_time,toc_time,ttr_time; - int i,sat,prn; + eph_t eph0 = {}; + geph_t geph0 = {}; + char buff[4096], *p; + long toe_time, tof_time, toc_time, ttr_time; + int i, sat, prn; - trace(3,"loadnav: file=%s\n",file); + trace(3, "loadnav: file=%s\n", file); - if (!(fp=fopen(file,"r"))) return 0; + if (!(fp = fopen(file, "r"))) return 0; - while (fgets(buff,sizeof(buff),fp)) + while (fgets(buff, sizeof(buff), fp)) { - if (!strncmp(buff,"IONUTC",6)) + if (!strncmp(buff, "IONUTC", 6)) { - for (i=0;i<8;i++) nav->ion_gps[i]=0.0; - for (i=0;i<4;i++) nav->utc_gps[i]=0.0; - nav->leaps=0; - sscanf(buff,"IONUTC,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%d", - &nav->ion_gps[0],&nav->ion_gps[1],&nav->ion_gps[2],&nav->ion_gps[3], - &nav->ion_gps[4],&nav->ion_gps[5],&nav->ion_gps[6],&nav->ion_gps[7], - &nav->utc_gps[0],&nav->utc_gps[1],&nav->utc_gps[2],&nav->utc_gps[3], + for (i = 0; i<8; i++) nav->ion_gps[i] = 0.0; + for (i = 0; i<4; i++) nav->utc_gps[i] = 0.0; + nav->leaps = 0; + sscanf(buff, "IONUTC,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%d", + &nav->ion_gps[0], &nav->ion_gps[1], &nav->ion_gps[2], &nav->ion_gps[3], + &nav->ion_gps[4], &nav->ion_gps[5], &nav->ion_gps[6], &nav->ion_gps[7], + &nav->utc_gps[0], &nav->utc_gps[1], &nav->utc_gps[2], &nav->utc_gps[3], &nav->leaps); continue; } - if ((p=strchr(buff,','))) *p='\0'; else continue; - if (!(sat=satid2no(buff))) continue; - if (satsys(sat,&prn)==SYS_GLO) + if ((p = strchr(buff, ','))) *p = '\0'; else continue; + if (!(sat = satid2no(buff))) continue; + if (satsys(sat, &prn) == SYS_GLO) { - nav->geph[prn-1]=geph0; - nav->geph[prn-1].sat=sat; - toe_time=tof_time=0; - sscanf(p+1,"%d,%d,%d,%d,%d,%ld,%ld,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf," - "%lf,%lf,%lf,%lf", - &nav->geph[prn-1].iode,&nav->geph[prn-1].frq,&nav->geph[prn-1].svh, - &nav->geph[prn-1].sva,&nav->geph[prn-1].age, - &toe_time,&tof_time, - &nav->geph[prn-1].pos[0],&nav->geph[prn-1].pos[1],&nav->geph[prn-1].pos[2], - &nav->geph[prn-1].vel[0],&nav->geph[prn-1].vel[1],&nav->geph[prn-1].vel[2], - &nav->geph[prn-1].acc[0],&nav->geph[prn-1].acc[1],&nav->geph[prn-1].acc[2], - &nav->geph[prn-1].taun ,&nav->geph[prn-1].gamn ,&nav->geph[prn-1].dtaun); - nav->geph[prn-1].toe.time=toe_time; - nav->geph[prn-1].tof.time=tof_time; + nav->geph[prn-1] = geph0; + nav->geph[prn-1].sat = sat; + toe_time = tof_time = 0; + sscanf(p+1, "%d,%d,%d,%d,%d,%ld,%ld,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf," + "%lf,%lf,%lf,%lf", + &nav->geph[prn-1].iode, &nav->geph[prn-1].frq, &nav->geph[prn-1].svh, + &nav->geph[prn-1].sva, &nav->geph[prn-1].age, + &toe_time, &tof_time, + &nav->geph[prn-1].pos[0], &nav->geph[prn-1].pos[1], &nav->geph[prn-1].pos[2], + &nav->geph[prn-1].vel[0], &nav->geph[prn-1].vel[1], &nav->geph[prn-1].vel[2], + &nav->geph[prn-1].acc[0], &nav->geph[prn-1].acc[1], &nav->geph[prn-1].acc[2], + &nav->geph[prn-1].taun , &nav->geph[prn-1].gamn , &nav->geph[prn-1].dtaun); + nav->geph[prn-1].toe.time = toe_time; + nav->geph[prn-1].tof.time = tof_time; } else { - nav->eph[sat-1]=eph0; - nav->eph[sat-1].sat=sat; - toe_time=toc_time=ttr_time=0; - sscanf(p+1,"%d,%d,%d,%d,%ld,%ld,%ld,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf," - "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%d,%d", - &nav->eph[sat-1].iode,&nav->eph[sat-1].iodc,&nav->eph[sat-1].sva , + nav->eph[sat-1] = eph0; + nav->eph[sat-1].sat = sat; + toe_time = toc_time = ttr_time = 0; + sscanf(p+1, "%d,%d,%d,%d,%ld,%ld,%ld,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf," + "%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%d,%d", + &nav->eph[sat-1].iode, &nav->eph[sat-1].iodc, &nav->eph[sat-1].sva , &nav->eph[sat-1].svh , - &toe_time,&toc_time,&ttr_time, - &nav->eph[sat-1].A ,&nav->eph[sat-1].e ,&nav->eph[sat-1].i0 , - &nav->eph[sat-1].OMG0,&nav->eph[sat-1].omg ,&nav->eph[sat-1].M0 , - &nav->eph[sat-1].deln,&nav->eph[sat-1].OMGd,&nav->eph[sat-1].idot, - &nav->eph[sat-1].crc ,&nav->eph[sat-1].crs ,&nav->eph[sat-1].cuc , - &nav->eph[sat-1].cus ,&nav->eph[sat-1].cic ,&nav->eph[sat-1].cis , - &nav->eph[sat-1].toes,&nav->eph[sat-1].fit ,&nav->eph[sat-1].f0 , - &nav->eph[sat-1].f1 ,&nav->eph[sat-1].f2 ,&nav->eph[sat-1].tgd[0], - &nav->eph[sat-1].code, &nav->eph[sat-1].flag); - nav->eph[sat-1].toe.time=toe_time; - nav->eph[sat-1].toc.time=toc_time; - nav->eph[sat-1].ttr.time=ttr_time; + &toe_time, &toc_time, &ttr_time, + &nav->eph[sat-1].A , &nav->eph[sat-1].e , &nav->eph[sat-1].i0 , + &nav->eph[sat-1].OMG0, &nav->eph[sat-1].omg , &nav->eph[sat-1].M0 , + &nav->eph[sat-1].deln, &nav->eph[sat-1].OMGd, &nav->eph[sat-1].idot, + &nav->eph[sat-1].crc , &nav->eph[sat-1].crs , &nav->eph[sat-1].cuc , + &nav->eph[sat-1].cus , &nav->eph[sat-1].cic , &nav->eph[sat-1].cis , + &nav->eph[sat-1].toes, &nav->eph[sat-1].fit , &nav->eph[sat-1].f0 , + &nav->eph[sat-1].f1 , &nav->eph[sat-1].f2 , &nav->eph[sat-1].tgd[0], + &nav->eph[sat-1].code, &nav->eph[sat-1].flag); + nav->eph[sat-1].toe.time = toe_time; + nav->eph[sat-1].toc.time = toc_time; + nav->eph[sat-1].ttr.time = ttr_time; } } fclose(fp); @@ -2957,52 +2837,52 @@ int readnav(const char *file, nav_t *nav) } -int savenav(const char *file, const nav_t *nav) +int savenav(const char *file, const nav_t *nav) { FILE *fp; int i; char id[32]; - trace(3,"savenav: file=%s\n",file); + trace(3, "savenav: file=%s\n", file); - if (!(fp=fopen(file,"w"))) return 0; + if (!(fp = fopen(file, "w"))) return 0; - for (i=0;ieph[i].ttr.time==0) continue; - satno2id(nav->eph[i].sat,id); - fprintf(fp,"%s,%d,%d,%d,%d,%d,%d,%d,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E," + if (nav->eph[i].ttr.time == 0) continue; + satno2id(nav->eph[i].sat, id); + fprintf(fp, "%s,%d,%d,%d,%d,%d,%d,%d,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E," "%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E," "%.14E,%.14E,%.14E,%.14E,%.14E,%d,%d\n", - id,nav->eph[i].iode,nav->eph[i].iodc,nav->eph[i].sva , - nav->eph[i].svh ,(int)nav->eph[i].toe.time, - (int)nav->eph[i].toc.time,(int)nav->eph[i].ttr.time, - nav->eph[i].A ,nav->eph[i].e ,nav->eph[i].i0 ,nav->eph[i].OMG0, - nav->eph[i].omg ,nav->eph[i].M0 ,nav->eph[i].deln,nav->eph[i].OMGd, - nav->eph[i].idot,nav->eph[i].crc,nav->eph[i].crs ,nav->eph[i].cuc , - nav->eph[i].cus ,nav->eph[i].cic,nav->eph[i].cis ,nav->eph[i].toes, - nav->eph[i].fit ,nav->eph[i].f0 ,nav->eph[i].f1 ,nav->eph[i].f2 , - nav->eph[i].tgd[0],nav->eph[i].code,nav->eph[i].flag); + id, nav->eph[i].iode, nav->eph[i].iodc, nav->eph[i].sva , + nav->eph[i].svh , (int)nav->eph[i].toe.time, + (int)nav->eph[i].toc.time, (int)nav->eph[i].ttr.time, + nav->eph[i].A , nav->eph[i].e , nav->eph[i].i0 , nav->eph[i].OMG0, + nav->eph[i].omg , nav->eph[i].M0 , nav->eph[i].deln, nav->eph[i].OMGd, + nav->eph[i].idot, nav->eph[i].crc, nav->eph[i].crs , nav->eph[i].cuc , + nav->eph[i].cus , nav->eph[i].cic, nav->eph[i].cis , nav->eph[i].toes, + nav->eph[i].fit , nav->eph[i].f0 , nav->eph[i].f1 , nav->eph[i].f2 , + nav->eph[i].tgd[0], nav->eph[i].code, nav->eph[i].flag); } - for (i=0;igeph[i].tof.time==0) continue; - satno2id(nav->geph[i].sat,id); - fprintf(fp,"%s,%d,%d,%d,%d,%d,%d,%d,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E," + if (nav->geph[i].tof.time == 0) continue; + satno2id(nav->geph[i].sat, id); + fprintf(fp, "%s,%d,%d,%d,%d,%d,%d,%d,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E," "%.14E,%.14E,%.14E,%.14E,%.14E,%.14E\n", - id,nav->geph[i].iode,nav->geph[i].frq,nav->geph[i].svh, - nav->geph[i].sva,nav->geph[i].age,(int)nav->geph[i].toe.time, + id, nav->geph[i].iode, nav->geph[i].frq, nav->geph[i].svh, + nav->geph[i].sva, nav->geph[i].age, (int)nav->geph[i].toe.time, (int)nav->geph[i].tof.time, - nav->geph[i].pos[0],nav->geph[i].pos[1],nav->geph[i].pos[2], - nav->geph[i].vel[0],nav->geph[i].vel[1],nav->geph[i].vel[2], - nav->geph[i].acc[0],nav->geph[i].acc[1],nav->geph[i].acc[2], - nav->geph[i].taun,nav->geph[i].gamn,nav->geph[i].dtaun); + nav->geph[i].pos[0], nav->geph[i].pos[1], nav->geph[i].pos[2], + nav->geph[i].vel[0], nav->geph[i].vel[1], nav->geph[i].vel[2], + nav->geph[i].acc[0], nav->geph[i].acc[1], nav->geph[i].acc[2], + nav->geph[i].taun, nav->geph[i].gamn, nav->geph[i].dtaun); } fprintf(fp,"IONUTC,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E,%.14E," "%.14E,%.14E,%.14E,%d", - nav->ion_gps[0],nav->ion_gps[1],nav->ion_gps[2],nav->ion_gps[3], - nav->ion_gps[4],nav->ion_gps[5],nav->ion_gps[6],nav->ion_gps[7], - nav->utc_gps[0],nav->utc_gps[1],nav->utc_gps[2],nav->utc_gps[3], + nav->ion_gps[0], nav->ion_gps[1], nav->ion_gps[2], nav->ion_gps[3], + nav->ion_gps[4], nav->ion_gps[5], nav->ion_gps[6], nav->ion_gps[7], + nav->utc_gps[0], nav->utc_gps[1], nav->utc_gps[2], nav->utc_gps[3], nav->leaps); fclose(fp); @@ -3017,7 +2897,7 @@ int savenav(const char *file, const nav_t *nav) *-----------------------------------------------------------------------------*/ void freeobs(obs_t *obs) { - free(obs->data); obs->data=NULL; obs->n=obs->nmax=0; + free(obs->data); obs->data = NULL; obs->n = obs->nmax = 0; } @@ -3033,14 +2913,14 @@ void freeobs(obs_t *obs) *-----------------------------------------------------------------------------*/ void freenav(nav_t *nav, int opt) { - if (opt&0x01) {free(nav->eph ); nav->eph =NULL; nav->n =nav->nmax =0;} - if (opt&0x02) {free(nav->geph); nav->geph=NULL; nav->ng=nav->ngmax=0;} - if (opt&0x04) {free(nav->seph); nav->seph=NULL; nav->ns=nav->nsmax=0;} - if (opt&0x08) {free(nav->peph); nav->peph=NULL; nav->ne=nav->nemax=0;} - if (opt&0x10) {free(nav->pclk); nav->pclk=NULL; nav->nc=nav->ncmax=0;} - if (opt&0x20) {free(nav->alm ); nav->alm =NULL; nav->na=nav->namax=0;} - if (opt&0x40) {free(nav->tec ); nav->tec =NULL; nav->nt=nav->ntmax=0;} - if (opt&0x80) {free(nav->fcb ); nav->fcb =NULL; nav->nf=nav->nfmax=0;} + if (opt&0x01) {free(nav->eph ); nav->eph = NULL; nav->n = nav->nmax = 0;} + if (opt&0x02) {free(nav->geph); nav->geph = NULL; nav->ng = nav->ngmax = 0;} + if (opt&0x04) {free(nav->seph); nav->seph = NULL; nav->ns = nav->nsmax = 0;} + if (opt&0x08) {free(nav->peph); nav->peph = NULL; nav->ne = nav->nemax = 0;} + if (opt&0x10) {free(nav->pclk); nav->pclk = NULL; nav->nc = nav->ncmax = 0;} + if (opt&0x20) {free(nav->alm ); nav->alm = NULL; nav->na = nav->namax = 0;} + if (opt&0x40) {free(nav->tec ); nav->tec = NULL; nav->nt = nav->ntmax = 0;} + if (opt&0x80) {free(nav->fcb ); nav->fcb = NULL; nav->nf = nav->nfmax = 0;} } @@ -3106,7 +2986,7 @@ void freenav(nav_t *nav, int opt) // va_list ap; // // /* print error message to stderr */ -// if (level<=1) { +// if (level <= 1) { // va_start(ap,format); vfprintf(stderr,format,ap); va_end(ap); // } // if (!fp_trace||level>level_trace) return; @@ -3271,7 +3151,7 @@ void traceb (int level, const unsigned char *p, int n) {} *-----------------------------------------------------------------------------*/ int execcmd(const char *cmd) { - trace(3,"execcmd: cmd=%s\n",cmd); + trace(3, "execcmd: cmd=%s\n", cmd); return system(cmd); } @@ -3284,34 +3164,34 @@ int execcmd(const char *cmd) *-----------------------------------------------------------------------------*/ void createdir(const char *path) { - char buff[1024],*p; + char buff[1024], *p; - tracet(3,"createdir: path=%s\n",path); + tracet(3, "createdir: path=%s\n", path); - strcpy(buff,path); - if (!(p=strrchr(buff,FILEPATHSEP))) return; - *p='\0'; + strcpy(buff, path); + if (!(p = strrchr(buff, FILEPATHSEP))) return; + *p = '\0'; - mkdir(buff,0777); + mkdir(buff, 0777); } /* replace string ------------------------------------------------------------*/ int repstr(char *str, const char *pat, const char *rep) { - int len=(int)strlen(pat); - char buff[1024],*p,*q,*r; + int len = (int)strlen(pat); + char buff[1024], *p, *q, *r; - for (p=str,r=buff;*p;p=q+len) + for (p = str, r = buff; *p; p = q+len) { - if (!(q=strstr(p,pat))) break; - strncpy(r,p,q-p); - r+=q-p; - r+=sprintf(r,"%s",rep); + if (!(q = strstr(p, pat))) break; + strncpy(r, p, q-p); + r += q-p; + r += sprintf(r, "%s", rep); } - if (p<=str) return 0; - strcpy(r,p); - strcpy(str,buff); + if (p <= str) return 0; + strcpy(r, p); + strcpy(str, buff); return 1; } @@ -3347,42 +3227,42 @@ int repstr(char *str, const char *pat, const char *rep) int reppath(const char *path, char *rpath, gtime_t time, const char *rov, const char *base) { - double ep[6],ep0[6]={2000,1,1,0,0,0}; - int week,dow,doy,stat=0; + double ep[6], ep0[6] = {2000, 1, 1, 0, 0, 0}; + int week, dow, doy, stat = 0; char rep[64]; - strcpy(rpath,path); + strcpy(rpath, path); - if (!strstr(rpath,"%")) return 0; - if (*rov ) stat|=repstr(rpath,"%r",rov ); - if (*base) stat|=repstr(rpath,"%b",base); - if (time.time!=0) + if (!strstr(rpath, "%")) return 0; + if (*rov ) stat|=repstr(rpath, "%r", rov ); + if (*base) stat|=repstr(rpath, "%b", base); + if (time.time != 0) { - time2epoch(time,ep); - ep0[0]=ep[0]; - dow=(int)floor(time2gpst(time,&week)/86400.0); - doy=(int)floor(timediff(time,epoch2time(ep0))/86400.0)+1; - sprintf(rep,"%02d", ((int)ep[3]/3)*3); stat|=repstr(rpath,"%ha",rep); - sprintf(rep,"%02d", ((int)ep[3]/6)*6); stat|=repstr(rpath,"%hb",rep); - sprintf(rep,"%02d", ((int)ep[3]/12)*12); stat|=repstr(rpath,"%hc",rep); - sprintf(rep,"%04.0f",ep[0]); stat|=repstr(rpath,"%Y",rep); - sprintf(rep,"%02.0f",fmod(ep[0],100.0)); stat|=repstr(rpath,"%y",rep); - sprintf(rep,"%02.0f",ep[1]); stat|=repstr(rpath,"%m",rep); - sprintf(rep,"%02.0f",ep[2]); stat|=repstr(rpath,"%d",rep); - sprintf(rep,"%02.0f",ep[3]); stat|=repstr(rpath,"%h",rep); - sprintf(rep,"%02.0f",ep[4]); stat|=repstr(rpath,"%M",rep); - sprintf(rep,"%02.0f",floor(ep[5])); stat|=repstr(rpath,"%S",rep); - sprintf(rep,"%03d", doy); stat|=repstr(rpath,"%n",rep); - sprintf(rep,"%04d", week); stat|=repstr(rpath,"%W",rep); - sprintf(rep,"%d", dow); stat|=repstr(rpath,"%D",rep); - sprintf(rep,"%c", 'a'+(int)ep[3]); stat|=repstr(rpath,"%H",rep); - sprintf(rep,"%02d", ((int)ep[4]/15)*15); stat|=repstr(rpath,"%t",rep); + time2epoch(time, ep); + ep0[0] = ep[0]; + dow = (int)floor(time2gpst(time, &week)/86400.0); + doy = (int)floor(timediff(time, epoch2time(ep0))/86400.0)+1; + sprintf(rep, "%02d", ((int)ep[3]/3)*3); stat|=repstr(rpath, "%ha", rep); + sprintf(rep, "%02d", ((int)ep[3]/6)*6); stat|=repstr(rpath, "%hb", rep); + sprintf(rep, "%02d", ((int)ep[3]/12)*12); stat|=repstr(rpath, "%hc", rep); + sprintf(rep, "%04.0f", ep[0]); stat|=repstr(rpath, "%Y", rep); + sprintf(rep, "%02.0f", fmod(ep[0], 100.0)); stat|=repstr(rpath, "%y", rep); + sprintf(rep, "%02.0f", ep[1]); stat|=repstr(rpath, "%m", rep); + sprintf(rep, "%02.0f", ep[2]); stat|=repstr(rpath, "%d", rep); + sprintf(rep, "%02.0f", ep[3]); stat|=repstr(rpath, "%h", rep); + sprintf(rep, "%02.0f", ep[4]); stat|=repstr(rpath, "%M", rep); + sprintf(rep, "%02.0f", floor(ep[5])); stat|=repstr(rpath, "%S", rep); + sprintf(rep, "%03d", doy); stat|=repstr(rpath, "%n", rep); + sprintf(rep, "%04d", week); stat|=repstr(rpath, "%W", rep); + sprintf(rep, "%d", dow); stat|=repstr(rpath, "%D", rep); + sprintf(rep, "%c", 'a'+(int)ep[3]); stat|=repstr(rpath, "%H", rep); + sprintf(rep, "%02d", ((int)ep[4]/15)*15); stat|=repstr(rpath, "%t", rep); } - else if (strstr(rpath,"%ha")||strstr(rpath,"%hb")||strstr(rpath,"%hc")|| - strstr(rpath,"%Y" )||strstr(rpath,"%y" )||strstr(rpath,"%m" )|| - strstr(rpath,"%d" )||strstr(rpath,"%h" )||strstr(rpath,"%M" )|| - strstr(rpath,"%S" )||strstr(rpath,"%n" )||strstr(rpath,"%W" )|| - strstr(rpath,"%D" )||strstr(rpath,"%H" )||strstr(rpath,"%t" )) + else if (strstr(rpath, "%ha") || strstr(rpath, "%hb") || strstr(rpath, "%hc") || + strstr(rpath, "%Y" ) || strstr(rpath, "%y" ) || strstr(rpath, "%m" ) || + strstr(rpath, "%d" ) || strstr(rpath, "%h" ) || strstr(rpath, "%M" ) || + strstr(rpath, "%S" ) || strstr(rpath, "%n" ) || strstr(rpath, "%W" ) || + strstr(rpath, "%D" ) || strstr(rpath, "%H" ) || strstr(rpath, "%t" )) { return -1; /* no valid time */ } @@ -3408,26 +3288,26 @@ int reppaths(const char *path, char *rpath[], int nmax, gtime_t ts, gtime_t te, const char *rov, const char *base) { gtime_t time; - double tow,tint=86400.0; - int i,n=0,week; + double tow, tint = 86400.0; + int i, n = 0, week; - trace(3,"reppaths: path =%s nmax=%d rov=%s base=%s\n",path,nmax,rov,base); + trace(3, "reppaths: path =%s nmax=%d rov=%s base=%s\n", path, nmax, rov, base); - if (ts.time==0||te.time==0||timediff(ts,te)>0.0) return 0; + if (ts.time == 0 || te.time == 0 || timediff(ts, te)>0.0) return 0; - if (strstr(path,"%S")||strstr(path,"%M")||strstr(path,"%t")) tint=900.0; - else if (strstr(path,"%h")||strstr(path,"%H")) tint=3600.0; + if (strstr(path, "%S") || strstr(path, "%M") || strstr(path, "%t")) tint = 900.0; + else if (strstr(path, "%h") || strstr(path, "%H")) tint = 3600.0; - tow=time2gpst(ts,&week); - time=gpst2time(week,floor(tow/tint)*tint); + tow = time2gpst(ts, &week); + time = gpst2time(week, floor(tow/tint)*tint); - while (timediff(time,te)<=0.0&&nng;i++) + for (i = 0; ing; i++) { - if (nav->geph[i].sat!=sat) continue; + if (nav->geph[i].sat != sat) continue; return SPEED_OF_LIGHT/(freq_glo[frq]+dfrq_glo[frq]*nav->geph[i].frq); } } - else if (frq==2) + else if (frq == 2) { /* L3 */ return SPEED_OF_LIGHT/FREQ3_GLO; } } - else if (sys==SYS_BDS) + else if (sys == SYS_BDS) { - if (frq==0) return SPEED_OF_LIGHT/FREQ1_BDS; /* B1 */ - else if (frq==1) return SPEED_OF_LIGHT/FREQ2_BDS; /* B2 */ - else if (frq==2) return SPEED_OF_LIGHT/FREQ3_BDS; /* B3 */ + if (frq == 0) return SPEED_OF_LIGHT/FREQ1_BDS; /* B1 */ + else if (frq == 1) return SPEED_OF_LIGHT/FREQ2_BDS; /* B2 */ + else if (frq == 2) return SPEED_OF_LIGHT/FREQ3_BDS; /* B3 */ } else { - if (frq==0) return SPEED_OF_LIGHT/FREQ1; /* L1/E1 */ - else if (frq==1) return SPEED_OF_LIGHT/FREQ2; /* L2 */ - else if (frq==2) return SPEED_OF_LIGHT/FREQ5; /* L5/E5a */ - else if (frq==3) return SPEED_OF_LIGHT/FREQ6; /* L6/LEX */ - else if (frq==4) return SPEED_OF_LIGHT/FREQ7; /* E5b */ - else if (frq==5) return SPEED_OF_LIGHT/FREQ8; /* E5a+b */ - else if (frq==6) return SPEED_OF_LIGHT/FREQ9; /* S */ + if (frq == 0) return SPEED_OF_LIGHT/FREQ1; /* L1/E1 */ + else if (frq == 1) return SPEED_OF_LIGHT/FREQ2; /* L2 */ + else if (frq == 2) return SPEED_OF_LIGHT/FREQ5; /* L5/E5a */ + else if (frq == 3) return SPEED_OF_LIGHT/FREQ6; /* L6/LEX */ + else if (frq == 4) return SPEED_OF_LIGHT/FREQ7; /* E5b */ + else if (frq == 5) return SPEED_OF_LIGHT/FREQ8; /* E5a+b */ + else if (frq == 6) return SPEED_OF_LIGHT/FREQ9; /* S */ } return 0.0; } @@ -3493,10 +3373,10 @@ double geodist(const double *rs, const double *rr, double *e) double r; int i; - if (norm(rs,3)-RE_WGS84) { - ecef2enu(pos,e,enu); - az=dot(enu,enu,2)<1e-12?0.0:atan2(enu[0],enu[1]); - if (az<0.0) az+=2*PI; - el=asin(enu[2]); + ecef2enu(pos, e, enu); + az = dot(enu, enu, 2)<1e-12 ? 0.0 : atan2(enu[0], enu[1]); + if (az<0.0) az += 2*PI; + el = asin(enu[2]); } - if (azel) {azel[0]=az; azel[1]=el;} + if (azel) {azel[0] = az; azel[1] = el;} return el; } @@ -3536,29 +3416,29 @@ double satazel(const double *pos, const double *e, double *azel) *-----------------------------------------------------------------------------*/ void dops(int ns, const double *azel, double elmin, double *dop) { - double H[4*MAXSAT],Q[16],cosel,sinel; - int i,n; + double H[4*MAXSAT], Q[16], cosel, sinel; + int i, n; - for (i=0;i<4;i++) dop[i]=0.0; - for (i=n=0;i 0.416) phi= 0.416; - else if (phi<-0.416) phi=-0.416; - lam=pos[1]/PI+psi*sin(azel[0])/cos(phi*PI); + phi = pos[0]/PI+psi*cos(azel[0]); + if (phi> 0.416) phi = 0.416; + else if (phi<-0.416) phi = -0.416; + lam = pos[1]/PI+psi*sin(azel[0])/cos(phi*PI); /* geomagnetic latitude (semi-circle) */ - phi+=0.064*cos((lam-1.617)*PI); + phi += 0.064*cos((lam-1.617)*PI); /* local time (s) */ - tt=43200.0*lam+time2gpst(t,&week); - tt-=floor(tt/86400.0)*86400.0; /* 0<=tt<86400 */ + tt = 43200.0*lam+time2gpst(t, &week); + tt -= floor(tt/86400.0)*86400.0; /* 0 <= tt<86400 */ /* slant factor */ - f=1.0+16.0*pow(0.53-azel[1]/PI,3.0); + f = 1.0+16.0*pow(0.53-azel[1]/PI, 3.0); /* ionospheric delay */ - amp=ion[0]+phi*(ion[1]+phi*(ion[2]+phi*ion[3])); - per=ion[4]+phi*(ion[5]+phi*(ion[6]+phi*ion[7])); - amp=amp< 0.0? 0.0:amp; - per=per<72000.0?72000.0:per; - x=2.0*PI*(tt-50400.0)/per; + amp = ion[0]+phi*(ion[1]+phi*(ion[2]+phi*ion[3])); + per = ion[4]+phi*(ion[5]+phi*(ion[6]+phi*ion[7])); + amp = amp< 0.0 ? 0.0 : amp; + per = per<72000.0 ? 72000.0 : per; + x = 2.0*PI*(tt-50400.0)/per; - return SPEED_OF_LIGHT*f*(fabs(x)<1.57?5E-9+amp*(1.0+x*x*(-0.5+x*x/24.0)):5E-9); + return SPEED_OF_LIGHT*f*(fabs(x)<1.57 ? 5E-9+amp*(1.0+x*x*(-0.5+x*x/24.0)) : 5E-9); } @@ -3622,7 +3502,7 @@ double ionmodel(gtime_t t, const double *ion, const double *pos, *-----------------------------------------------------------------------------*/ double ionmapf(const double *pos, const double *azel) { - if (pos[2]>=HION) return 1.0; + if (pos[2] >= HION) return 1.0; return 1.0/cos(asin((RE_WGS84+pos[2])/(RE_WGS84+HION)*sin(PI/2.0-azel[1]))); } @@ -3641,23 +3521,23 @@ double ionmapf(const double *pos, const double *azel) double ionppp(const double *pos, const double *azel, double re, double hion, double *posp) { - double cosaz,rp,ap,sinap,tanap; + double cosaz, rp, ap, sinap, tanap; - rp=re/(re+hion)*cos(azel[1]); - ap=PI/2.0-azel[1]-asin(rp); - sinap=sin(ap); - tanap=tan(ap); - cosaz=cos(azel[0]); - posp[0]=asin(sin(pos[0])*cos(ap)+cos(pos[0])*sinap*cosaz); + rp = re/(re+hion)*cos(azel[1]); + ap = PI/2.0-azel[1]-asin(rp); + sinap = sin(ap); + tanap = tan(ap); + cosaz = cos(azel[0]); + posp[0] = asin(sin(pos[0])*cos(ap)+cos(pos[0])*sinap*cosaz); - if ((pos[0]> 70.0*D2R&& tanap*cosaz>tan(PI/2.0-pos[0]))|| - (pos[0]<-70.0*D2R&&-tanap*cosaz>tan(PI/2.0+pos[0]))) + if ((pos[0]> 70.0*D2R && tanap*cosaz>tan(PI/2.0-pos[0])) || + (pos[0]<-70.0*D2R && -tanap*cosaz>tan(PI/2.0+pos[0]))) { - posp[1]=pos[1]+PI-asin(sinap*sin(azel[0])/cos(posp[0])); + posp[1] = pos[1]+PI-asin(sinap*sin(azel[0])/cos(posp[0])); } else { - posp[1]=pos[1]+asin(sinap*sin(azel[0])/cos(posp[0])); + posp[1] = pos[1]+asin(sinap*sin(azel[0])/cos(posp[0])); } return 1.0/sqrt(1.0-rp*rp); } @@ -3674,22 +3554,22 @@ double ionppp(const double *pos, const double *azel, double re, double tropmodel(gtime_t time, const double *pos, const double *azel, double humi) { - const double temp0=15.0; /* temparature at sea level */ - double hgt,pres,temp,e,z,trph,trpw; + const double temp0 = 15.0; /* temparature at sea level */ + double hgt, pres, temp, e, z, trph, trpw; - if (pos[2]<-100.0||1e44) return coef[4]; return coef[i-1]*(1.0-lat/15.0+i)+coef[i]*(lat/15.0-i); } @@ -3705,7 +3585,7 @@ double interpc(const double coef[], double lat) double mapf(double el, double a, double b, double c) { - double sinel=sin(el); + double sinel = sin(el); return (1.0+a/(1.0+b/(1.0+c)))/(sinel+(a/(sinel+b/(sinel+c)))); } @@ -3715,7 +3595,7 @@ double nmf(gtime_t time, const double pos[], const double azel[], { /* ref [5] table 3 */ /* hydro-ave-a,b,c, hydro-amp-a,b,c, wet-a,b,c at latitude 15,30,45,60,75 */ - const double coef[][5]={ + const double coef[][5] = { { 1.2769934E-3, 1.2683230E-3, 1.2465397E-3, 1.2196049E-3, 1.2045996E-3}, { 2.9153695E-3, 2.9152299E-3, 2.9288445E-3, 2.9022565E-3, 2.9024912E-3}, { 62.610505E-3, 62.837393E-3, 63.721774E-3, 63.824265E-3, 64.258455E-3}, @@ -3728,33 +3608,33 @@ double nmf(gtime_t time, const double pos[], const double azel[], { 1.4275268E-3, 1.5138625E-3, 1.4572752E-3, 1.5007428E-3, 1.7599082E-3}, { 4.3472961e-2, 4.6729510E-2, 4.3908931e-2, 4.4626982E-2, 5.4736038E-2} }; - const double aht[]={ 2.53E-5, 5.49E-3, 1.14E-3}; /* height correction */ + const double aht[] = { 2.53E-5, 5.49E-3, 1.14E-3}; /* height correction */ - double y,cosy,ah[3],aw[3],dm,el=azel[1],lat=pos[0]*R2D,hgt=pos[2]; + double y, cosy, ah[3], aw[3], dm, el = azel[1], lat = pos[0]*R2D, hgt = pos[2]; int i; - if (el<=0.0) + if (el <= 0.0) { - if (mapfw) *mapfw=0.0; + if (mapfw) *mapfw = 0.0; return 0.0; } /* year from doy 28, added half a year for southern latitudes */ - y=(time2doy(time)-28.0)/365.25+(lat<0.0?0.5:0.0); + y = (time2doy(time)-28.0)/365.25+(lat<0.0 ? 0.5 : 0.0); - cosy=cos(2.0*PI*y); - lat=fabs(lat); + cosy = cos(2.0*PI*y); + lat = fabs(lat); - for (i=0;i<3;i++) + for (i = 0; i<3; i++) { - ah[i]=interpc(coef[i ],lat)-interpc(coef[i+3],lat)*cosy; - aw[i]=interpc(coef[i+6],lat); + ah[i] = interpc(coef[i ], lat)-interpc(coef[i+3], lat)*cosy; + aw[i] = interpc(coef[i+6], lat); } /* ellipsoidal height is used instead of height above sea level */ - dm=(1.0/sin(el)-mapf(el,aht[0],aht[1],aht[2]))*hgt/1e3; + dm = (1.0/sin(el)-mapf(el, aht[0], aht[1], aht[2]))*hgt/1e3; - if (mapfw) *mapfw=mapf(el,aw[0],aw[1],aw[2]); + if (mapfw) *mapfw = mapf(el, aw[0], aw[1], aw[2]); - return mapf(el,ah[0],ah[1],ah[2])+dm; + return mapf(el, ah[0], ah[1], ah[2])+dm; } #endif /* !IERS_MODEL */ @@ -3775,31 +3655,31 @@ double tropmapf(gtime_t time, const double pos[], const double azel[], double *mapfw) { #ifdef IERS_MODEL - const double ep[]={2000,1,1,12,0,0}; - double mjd,lat,lon,hgt,zd,gmfh,gmfw; + const double ep[] = {2000, 1, 1, 12, 0, 0}; + double mjd, lat, lon, hgt, zd, gmfh, gmfw; #endif - trace(4,"tropmapf: pos=%10.6f %11.6f %6.1f azel=%5.1f %4.1f\n", - pos[0]*R2D,pos[1]*R2D,pos[2],azel[0]*R2D,azel[1]*R2D); + trace(4, "tropmapf: pos=%10.6f %11.6f %6.1f azel=%5.1f %4.1f\n", + pos[0]*R2D, pos[1]*R2D, pos[2], azel[0]*R2D, azel[1]*R2D); - if (pos[2]<-1000.0||pos[2]>20000.0) + if (pos[2]<-1000.0 || pos[2]>20000.0) { - if (mapfw) *mapfw=0.0; + if (mapfw) *mapfw = 0.0; return 0.0; } #ifdef IERS_MODEL - mjd=51544.5+(timediff(time,epoch2time(ep)))/86400.0; - lat=pos[0]; - lon=pos[1]; - hgt=pos[2]-geoidh(pos); /* height in m (mean sea level) */ - zd =PI/2.0-azel[1]; + mjd = 51544.5+(timediff(time, epoch2time(ep)))/86400.0; + lat = pos[0]; + lon = pos[1]; + hgt = pos[2]-geoidh(pos); /* height in m (mean sea level) */ + zd = PI/2.0-azel[1]; /* call GMF */ - gmf_(&mjd,&lat,&lon,&hgt,&zd,&gmfh,&gmfw); + gmf_(&mjd, &lat, &lon, &hgt, &zd, &gmfh, &gmfw); - if (mapfw) *mapfw=gmfw; + if (mapfw) *mapfw = gmfw; return gmfh; #else - return nmf(time,pos,azel,mapfw); /* NMF */ + return nmf(time, pos, azel, mapfw); /* NMF */ #endif } @@ -3807,9 +3687,9 @@ double tropmapf(gtime_t time, const double pos[], const double azel[], /* interpolate antenna phase center variation --------------------------------*/ double interpvar(double ang, const double *var) { - double a=ang/5.0; /* ang=0-90 */ - int i=(int)a; - if (i<0) return var[0]; else if (i>=18) return var[18]; + double a = ang/5.0; /* ang=0-90 */ + int i = (int)a; + if (i<0) return var[0]; else if (i >= 18) return var[18]; return var[i]*(1.0-a+i)+var[i+1]*(a-i); } @@ -3826,22 +3706,22 @@ double interpvar(double ang, const double *var) void antmodel(const pcv_t *pcv, const double *del, const double *azel, int opt, double *dant) { - double e[3],off[3],cosel=cos(azel[1]); - int i,j; + double e[3], off[3], cosel = cos(azel[1]); + int i, j; - trace(4,"antmodel: azel=%6.1f %4.1f opt=%d\n",azel[0]*R2D,azel[1]*R2D,opt); + trace(4, "antmodel: azel=%6.1f %4.1f opt=%d\n", azel[0]*R2D, azel[1]*R2D, opt); - e[0]=sin(azel[0])*cosel; - e[1]=cos(azel[0])*cosel; - e[2]=sin(azel[1]); + e[0] = sin(azel[0])*cosel; + e[1] = cos(azel[0])*cosel; + e[2] = sin(azel[1]); - for (i=0;ioff[i][j]+del[j]; + for (j = 0; j<3; j++) off[j] = pcv->off[i][j]+del[j]; - dant[i]=-dot(off,e,3)+(opt?interpvar(90.0-azel[1]*R2D,pcv->var[i]):0.0); + dant[i] = -dot(off, e, 3)+(opt ? interpvar(90.0-azel[1]*R2D, pcv->var[i]) : 0.0); } - trace(5,"antmodel: dant=%6.3f %6.3f\n",dant[0],dant[1]); + trace(5, "antmodel: dant=%6.3f %6.3f\n", dant[0], dant[1]); } @@ -3856,62 +3736,62 @@ void antmodel_s(const pcv_t *pcv, double nadir, double *dant) { int i; - trace(4,"antmodel_s: nadir=%6.1f\n",nadir*R2D); + trace(4, "antmodel_s: nadir=%6.1f\n", nadir*R2D); - for (i=0;ivar[i]); + dant[i] = interpvar(nadir*R2D*5.0, pcv->var[i]); } - trace(5,"antmodel_s: dant=%6.3f %6.3f\n",dant[0],dant[1]); + trace(5, "antmodel_s: dant=%6.3f %6.3f\n", dant[0], dant[1]); } /* sun and moon position in eci (ref [4] 5.1.1, 5.2.1) -----------------------*/ void sunmoonpos_eci(gtime_t tut, double *rsun, double *rmoon) { - const double ep2000[]={2000,1,1,12,0,0}; - double t,f[5],eps,Ms,ls,rs,lm,pm,rm,sine,cose,sinp,cosp,sinl,cosl; + const double ep2000[] = {2000, 1, 1, 12, 0, 0}; + double t, f[5], eps, Ms, ls, rs, lm, pm, rm, sine, cose, sinp, cosp, sinl, cosl; - trace(4,"sunmoonpos_eci: tut=%s\n",time_str(tut,3)); + trace(4, "sunmoonpos_eci: tut=%s\n", time_str(tut, 3)); - t=timediff(tut,epoch2time(ep2000))/86400.0/36525.0; + t = timediff(tut, epoch2time(ep2000))/86400.0/36525.0; /* astronomical arguments */ - ast_args(t,f); + ast_args(t, f); /* obliquity of the ecliptic */ - eps=23.439291-0.0130042*t; - sine=sin(eps*D2R); cose=cos(eps*D2R); + eps = 23.439291-0.0130042*t; + sine = sin(eps*D2R); cose = cos(eps*D2R); /* sun position in eci */ if (rsun) { - Ms=357.5277233+35999.05034*t; - ls=280.460+36000.770*t+1.914666471*sin(Ms*D2R)+0.019994643*sin(2.0*Ms*D2R); - rs=AU*(1.000140612-0.016708617*cos(Ms*D2R)-0.000139589*cos(2.0*Ms*D2R)); - sinl=sin(ls*D2R); cosl=cos(ls*D2R); - rsun[0]=rs*cosl; - rsun[1]=rs*cose*sinl; - rsun[2]=rs*sine*sinl; + Ms = 357.5277233+35999.05034*t; + ls = 280.460+36000.770*t+1.914666471*sin(Ms*D2R)+0.019994643*sin(2.0*Ms*D2R); + rs = AU*(1.000140612-0.016708617*cos(Ms*D2R)-0.000139589*cos(2.0*Ms*D2R)); + sinl = sin(ls*D2R); cosl = cos(ls*D2R); + rsun[0] = rs*cosl; + rsun[1] = rs*cose*sinl; + rsun[2] = rs*sine*sinl; - trace(5,"rsun =%.3f %.3f %.3f\n",rsun[0],rsun[1],rsun[2]); + trace(5, "rsun =%.3f %.3f %.3f\n", rsun[0], rsun[1], rsun[2]); } /* moon position in eci */ if (rmoon) { - lm=218.32+481267.883*t+6.29*sin(f[0])-1.27*sin(f[0]-2.0*f[3])+ + lm = 218.32+481267.883*t+6.29*sin(f[0])-1.27*sin(f[0]-2.0*f[3])+ 0.66*sin(2.0*f[3])+0.21*sin(2.0*f[0])-0.19*sin(f[1])-0.11*sin(2.0*f[2]); - pm=5.13*sin(f[2])+0.28*sin(f[0]+f[2])-0.28*sin(f[2]-f[0])- + pm = 5.13*sin(f[2])+0.28*sin(f[0]+f[2])-0.28*sin(f[2]-f[0])- 0.17*sin(f[2]-2.0*f[3]); - rm=RE_WGS84/sin((0.9508+0.0518*cos(f[0])+0.0095*cos(f[0]-2.0*f[3])+ + rm = RE_WGS84/sin((0.9508+0.0518*cos(f[0])+0.0095*cos(f[0]-2.0*f[3])+ 0.0078*cos(2.0*f[3])+0.0028*cos(2.0*f[0]))*D2R); - sinl=sin(lm*D2R); cosl=cos(lm*D2R); - sinp=sin(pm*D2R); cosp=cos(pm*D2R); - rmoon[0]=rm*cosp*cosl; - rmoon[1]=rm*(cose*cosp*sinl-sine*sinp); - rmoon[2]=rm*(sine*cosp*sinl+cose*sinp); + sinl = sin(lm*D2R); cosl = cos(lm*D2R); + sinp = sin(pm*D2R); cosp = cos(pm*D2R); + rmoon[0] = rm*cosp*cosl; + rmoon[1] = rm*(cose*cosp*sinl-sine*sinp); + rmoon[2] = rm*(sine*cosp*sinl+cose*sinp); - trace(5,"rmoon=%.3f %.3f %.3f\n",rmoon[0],rmoon[1],rmoon[2]); + trace(5, "rmoon=%.3f %.3f %.3f\n", rmoon[0], rmoon[1], rmoon[2]); } } @@ -3929,22 +3809,22 @@ void sunmoonpos(gtime_t tutc, const double *erpv, double *rsun, double *rmoon, double *gmst) { gtime_t tut; - double rs[3],rm[3],U[9],gmst_; + double rs[3], rm[3], U[9], gmst_; - trace(4,"sunmoonpos: tutc=%s\n",time_str(tutc,3)); + trace(4, "sunmoonpos: tutc=%s\n", time_str(tutc, 3)); - tut=timeadd(tutc,erpv[2]); /* utc -> ut1 */ + tut = timeadd(tutc, erpv[2]); /* utc -> ut1 */ /* sun and moon position in eci */ - sunmoonpos_eci(tut,rsun?rs:NULL,rmoon?rm:NULL); + sunmoonpos_eci(tut, rsun ? rs : NULL, rmoon ? rm : NULL); /* eci to ecef transformation matrix */ - eci2ecef(tutc,erpv,U,&gmst_); + eci2ecef(tutc, erpv, U, &gmst_); /* sun and moon postion in ecef */ - if (rsun ) matmul("NN",3,1,3,1.0,U,rs,0.0,rsun ); - if (rmoon) matmul("NN",3,1,3,1.0,U,rm,0.0,rmoon); - if (gmst ) *gmst=gmst_; + if (rsun ) matmul("NN", 3, 1, 3, 1.0, U, rs, 0.0, rsun ); + if (rmoon) matmul("NN", 3, 1, 3, 1.0, U, rm, 0.0, rmoon); + if (gmst ) *gmst = gmst_; } @@ -3956,28 +3836,28 @@ void sunmoonpos(gtime_t tutc, const double *erpv, double *rsun, *-----------------------------------------------------------------------------*/ void csmooth(obs_t *obs, int ns) { - double Ps[2][MAXSAT][NFREQ]={},Lp[2][MAXSAT][NFREQ]={},dcp; - int i,j,s,r,n[2][MAXSAT][NFREQ]={}; + double Ps[2][MAXSAT][NFREQ] = {}, Lp[2][MAXSAT][NFREQ] = {}, dcp; + int i, j, s, r, n[2][MAXSAT][NFREQ] = {}; obsd_t *p; - trace(3,"csmooth: nobs=%d,ns=%d\n",obs->n,ns); + trace(3, "csmooth: nobs=%d,ns=%d\n", obs->n, ns); - for (i=0;in;i++) + for (i = 0; in; i++) { - p=&obs->data[i]; s=p->sat; r=p->rcv; - for (j=0;jdata[i]; s = p->sat; r = p->rcv; + for (j = 0; jP[j]==0.0||p->L[j]==0.0) continue; - if (p->LLI[j]) n[r-1][s-1][j]=0; - if (n[r-1][s-1][j]==0) Ps[r-1][s-1][j]=p->P[j]; + if (s <= 0 || MAXSATP[j] == 0.0 || p->L[j] == 0.0) continue; + if (p->LLI[j]) n[r-1][s-1][j] = 0; + if (n[r-1][s-1][j] == 0) Ps[r-1][s-1][j] = p->P[j]; else { - dcp=lam_carr[j]*(p->L[j]-Lp[r-1][s-1][j]); - Ps[r-1][s-1][j]=p->P[j]/ns+(Ps[r-1][s-1][j]+dcp)*(ns-1)/ns; + dcp = lam_carr[j]*(p->L[j]-Lp[r-1][s-1][j]); + Ps[r-1][s-1][j] = p->P[j]/ns+(Ps[r-1][s-1][j]+dcp)*(ns-1)/ns; } - if (++n[r-1][s-1][j]P[j]=0.0; else p->P[j]=Ps[r-1][s-1][j]; - Lp[r-1][s-1][j]=p->L[j]; + if (++n[r-1][s-1][j]P[j] = 0.0; else p->P[j] = Ps[r-1][s-1][j]; + Lp[r-1][s-1][j] = p->L[j]; } } } @@ -3993,56 +3873,56 @@ void csmooth(obs_t *obs, int ns) *-----------------------------------------------------------------------------*/ int rtk_uncompress(const char *file, char *uncfile) { - int stat=0; - char *p,cmd[2048]="",tmpfile[1024]="",buff[1024],*fname,*dir=(char*)""; + int stat = 0; + char *p, cmd[2048] = "", tmpfile[1024] = "", buff[1024], *fname, *dir = (char*)""; - trace(3,"rtk_uncompress: file=%s\n",file); + trace(3, "rtk_uncompress: file=%s\n", file); - strcpy(tmpfile,file); - if (!(p=strrchr(tmpfile,'.'))) return 0; + strcpy(tmpfile, file); + if (!(p = strrchr(tmpfile, '.'))) return 0; /* uncompress by gzip */ - if (!strcmp(p,".z" )||!strcmp(p,".Z" )|| - !strcmp(p,".gz" )||!strcmp(p,".GZ" )|| - !strcmp(p,".zip")||!strcmp(p,".ZIP")) + if (!strcmp(p, ".z" ) || !strcmp(p, ".Z" ) || + !strcmp(p, ".gz" ) || !strcmp(p, ".GZ" ) || + !strcmp(p, ".zip") || !strcmp(p, ".ZIP")) { - strcpy(uncfile,tmpfile); uncfile[p-tmpfile]='\0'; - sprintf(cmd,"gzip -f -d -c \"%s\" > \"%s\"",tmpfile,uncfile); + strcpy(uncfile, tmpfile); uncfile[p-tmpfile] = '\0'; + sprintf(cmd, "gzip -f -d -c \"%s\" > \"%s\"", tmpfile, uncfile); if (execcmd(cmd)) { remove(uncfile); return -1; } - strcpy(tmpfile,uncfile); - stat=1; + strcpy(tmpfile, uncfile); + stat = 1; } /* extract tar file */ - if ((p=strrchr(tmpfile,'.'))&&!strcmp(p,".tar")) + if ((p = strrchr(tmpfile, '.')) && !strcmp(p, ".tar")) { - strcpy(uncfile,tmpfile); uncfile[p-tmpfile]='\0'; - strcpy(buff,tmpfile); - fname=buff; - if ((p=strrchr(buff,'/'))) + strcpy(uncfile, tmpfile); uncfile[p-tmpfile] = '\0'; + strcpy(buff, tmpfile); + fname = buff; + if ((p = strrchr(buff, '/'))) { - *p='\0'; dir=fname; fname=p+1; + *p = '\0'; dir = fname; fname = p+1; } - sprintf(cmd,"tar -C \"%s\" -xf \"%s\"",dir,tmpfile); + sprintf(cmd, "tar -C \"%s\" -xf \"%s\"", dir, tmpfile); if (execcmd(cmd)) { if (stat) remove(tmpfile); return -1; } if (stat) remove(tmpfile); - stat=1; + stat = 1; } /* extract hatanaka-compressed file by cnx2rnx */ - else if ((p=strrchr(tmpfile,'.'))&&strlen(p)>3&&(*(p+3)=='d'||*(p+3)=='D')) + else if ((p = strrchr(tmpfile, '.')) && strlen(p)>3 && (*(p+3) == 'd' || *(p+3) == 'D')) { - strcpy(uncfile,tmpfile); - uncfile[p-tmpfile+3]=*(p+3)=='D'?'O':'o'; - sprintf(cmd,"crx2rnx < \"%s\" > \"%s\"",tmpfile,uncfile); + strcpy(uncfile, tmpfile); + uncfile[p-tmpfile+3] = *(p+3) == 'D' ? 'O' : 'o'; + sprintf(cmd, "crx2rnx < \"%s\" > \"%s\"", tmpfile, uncfile); if (execcmd(cmd)) { @@ -4051,9 +3931,9 @@ int rtk_uncompress(const char *file, char *uncfile) return -1; } if (stat) remove(tmpfile); - stat=1; + stat = 1; } - trace(3,"rtk_uncompress: stat=%d\n",stat); + trace(3, "rtk_uncompress: stat=%d\n", stat); return stat; } @@ -4068,50 +3948,49 @@ int rtk_uncompress(const char *file, char *uncfile) *-----------------------------------------------------------------------------*/ int expath(const char *path, char *paths[], int nmax) { - int i,j,n=0; + int i, j, n = 0; char tmp[1024]; struct dirent *d; DIR *dp; - const char *file=path; - char dir[1024]="",s1[1024],s2[1024],*p,*q,*r; + const char *file = path; + char dir[1024] = "", s1[1024], s2[1024], *p, *q, *r; - trace(3,"expath : path=%s nmax=%d\n",path,nmax); + trace(3, "expath : path=%s nmax=%d\n", path, nmax); //TODO: Fix invalid conversion from ‘const char*’ to ‘char*’ - //if ((p=strrchr(path,'/'))||(p=strrchr(path,'\\'))) { + //if ((p=strrchr(path,'/')) || (p=strrchr(path,'\\'))) { // file=p+1; strncpy(dir,path,p-path+1); dir[p-path+1]='\0'; //} - if (!(dp=opendir(*dir?dir:"."))) return 0; - while ((d=readdir(dp))) + if (!(dp = opendir(*dir ? dir : "."))) return 0; + while ((d = readdir(dp))) { - if (*(d->d_name)=='.') continue; - sprintf(s1,"^%s$",d->d_name); - sprintf(s2,"^%s$",file); - for (p=s1;*p;p++) *p=(char)tolower((int)*p); - for (p=s2;*p;p++) *p=(char)tolower((int)*p); + if (*(d->d_name) == '.') continue; + sprintf(s1, "^%s$", d->d_name); + sprintf(s2, "^%s$", file); + for (p = s1; *p;p++) *p = (char)tolower((int)*p); + for (p = s2; *p;p++) *p = (char)tolower((int)*p); - for (p=s1,q=strtok_r(s2,"*",&r);q;q=strtok_r(NULL,"*",&r)) + for (p = s1, q = strtok_r(s2, "*", &r); q; q = strtok_r(NULL, "*", &r)) { - if ((p=strstr(p,q))) p+=strlen(q); else break; + if ((p = strstr(p, q))) p += strlen(q); else break; } - if (p&&nd_name); + if (p && nd_name); } closedir(dp); /* sort paths in alphabetical order */ - for (i=0;i0) + if (strcmp(paths[i], paths[j])>0) { - strcpy(tmp,paths[i]); - strcpy(paths[i],paths[j]); - strcpy(paths[j],tmp); + strcpy(tmp, paths[i]); + strcpy(paths[i], paths[j]); + strcpy(paths[j], tmp); } } } - for (i=0;i