rtklib代码支持北斗三代电文（CNAV）解算的不完全修改

概述

因为工作原因需要使用rtklib处理北斗三代(B1C/B2a)的数据，不多说废话，直接上代码

代码修改

因为rtklib抽象后的eph_t结构已经支持北斗三代，通过eph_t结构体的code字段判断星历来源，在计算卫星位置的时候，在ephemeris.c文件作如下修改

/* broadcast ephemeris to satellite position and clock bias --------------------
* compute satellite position and clock bias with broadcast ephemeris (gps,
* galileo, qzss)
* args   : gtime_t time     I   time (gpst)
*          eph_t *eph       I   broadcast ephemeris
*          double *rs       O   satellite position (ecef) {x,y,z} (m)
*          double *dts      O   satellite clock bias (s)
*          double *var      O   satellite position and clock variance (m^2)
* return : none
* notes  : see ref [1],[7],[8]
*          satellite clock includes relativity correction without code bias
*          (tgd or bgd)
*-----------------------------------------------------------------------------*/
extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts,double *var)
{double tk,M,E,Ek,sinE,cosE,u,r,i,O,sin2u,cos2u,x,y,sinO,cosO,cosi,mu,omge;double xg,yg,zg,sino,coso;int n,sys,prn;double A,A0,deltNa,Na,N0;trace(4,"eph2pos : time=%s sat=%2d\n",time_str(time,3),eph->sat);tk=timediff(time,eph->toe);switch ((sys=satsys(eph->sat,&prn))) {case SYS_GAL: mu=MU_GAL; omge=OMGE_GAL; break;case SYS_CMP: mu=MU_CMP; omge=OMGE_CMP; break;default:      mu=MU_GPS; omge=OMGE;     break;}if(sys==SYS_CMP&&(eph->code == 7||eph->code==9))//data source (0:unknown,1:B1I,2:B1Q,3:B2I,4:B2Q,5:B3I,6:B3Q,7:B1Cd,8:B1Cp,9:B2ad,10:B2ap) */{A0 = eph->A;A = A0 + eph->Adot*tk;N0        = sqrt(mu/(A0*A0*A0));deltNa    = eph->deln +0.5*eph->ndot*tk;Na        = N0 + deltNa;M         = eph->M0 + Na*tk;} else{A = eph->A;M=eph->M0+(sqrt(mu/(eph->A*eph->A*eph->A))+eph->deln)*tk;}for (n=0,E=M,Ek=0.0;fabs(E-Ek)>RTOL_KEPLER&&n<MAX_ITER_KEPLER;n++) {Ek=E; E-=(E-eph->e*sin(E)-M)/(1.0-eph->e*cos(E));}if (n>=MAX_ITER_KEPLER) {trace(2,"eph2pos: kepler iteration overflow sat=%2d\n",eph->sat);return;}sinE=sin(E); cosE=cos(E);trace(4,"kepler: sat=%2d e=%8.5f n=%2d del=%10.3e\n",eph->sat,eph->e,n,E-Ek);u=atan2(sqrt(1.0-eph->e*eph->e)*sinE,cosE-eph->e)+eph->omg;r=A*(1.0-eph->e*cosE);i=eph->i0+eph->idot*tk;sin2u=sin(2.0*u); cos2u=cos(2.0*u);u+=eph->cus*sin2u+eph->cuc*cos2u;r+=eph->crs*sin2u+eph->crc*cos2u;i+=eph->cis*sin2u+eph->cic*cos2u;x=r*cos(u); y=r*sin(u); cosi=cos(i);/* beidou geo satellite */if (sys==SYS_CMP&&(prn<=5||prn>=59)) { /* ref [9] table 4-1 */O=eph->OMG0+eph->OMGd*tk-omge*eph->toes;sinO=sin(O); cosO=cos(O);xg=x*cosO-y*cosi*sinO;yg=x*sinO+y*cosi*cosO;zg=y*sin(i);sino=sin(omge*tk); coso=cos(omge*tk);rs[0]= xg*coso+yg*sino*COS_5+zg*sino*SIN_5;rs[1]=-xg*sino+yg*coso*COS_5+zg*coso*SIN_5;rs[2]=-yg*SIN_5+zg*COS_5;}else {O=eph->OMG0+(eph->OMGd-omge)*tk-omge*eph->toes;sinO=sin(O); cosO=cos(O);rs[0]=x*cosO-y*cosi*sinO;rs[1]=x*sinO+y*cosi*cosO;rs[2]=y*sin(i);}tk=timediff(time,eph->toc);*dts=eph->f0+eph->f1*tk+eph->f2*tk*tk;/* relativity correction */*dts-=2.0*sqrt(mu*A)*eph->e*sinE/SQR(CLIGHT);/* position and clock error variance */*var=var_uraeph(sys,eph->sva);
}