rinex.c

一个很不错关于GPS中RINEX文件的编码！

rinex.c

void gps_dat(fp,ifp,xhidat,odat_t,odat_dt,odat_n,In0,rt,dt,odat,
             id_odat,iodat,odat1,iodat1,
             sdat,id_sdat,isdat,sdat_n,ddat,id_ddat,iddat,ddat_n,
             tdat,id_tdat,itdat,tdat_n,ref_sv,
             slips_t,id_slips,ion_r,iion_r
             ,icon,x0,Height,borb,itborb,Tr1,singlex,singlet,singlx,singlt)
FILE *fp[Ista];
long ifp,*In0,ref_sv;
long xhidat[2*Isat*Ista];
double rt,dt;
double odat[Isato*I10*Ista*Iepoches],odat1[Isato*I10*Ista*Iepoches];
double odat_t[Iepoches],odat_dt[Ista*Iepoches];
long odat_n[Iepoches],id_odat[I10*Ista*Iepoches];
unsigned int iodat[Isato*I10*Ista*Iepoches],iodat1[Isato*I10*Ista*Iepoches];
double sdat[Isato*I10*Ista*Iepoches],ddat[Isato*I10*Ista*Iepoches],
       tdat[Isato*I10*Ista*Iepoches];
long id_sdat[I10*Ista*Iepoches],id_ddat[I10*Ista*Iepoches],
     id_tdat[I10*Ista*Iepoches],
     sdat_n[Iepoches],ddat_n[Iepoches],tdat_n[Iepoches];
unsigned int isdat[Isato*I10*Ista*Iepoches],iddat[Isato*I10*Ista*Iepoches],
             itdat[Isato*I10*Ista*Iepoches];
double slips_t[Islips],ion_r[Isato*10*Ista*Iepoches];
long id_slips[Islips];
unsigned int iion_r[Isato*10*Ista*Iepoches];
int icon[I10*Isat];
double x0[3*Ista],Height[Ista];
long itborb[8]; /* orbit */
double borb[Isat*I96*8];
double Tr1[I3*Ista*I96]; /* tide */
double singlex[3*Ista*Iepoches],singlet[Ista*Iepoches];
double singlx[3*Ista],singlt[Ista];
{    int hour,minute;
     double sec;
     long IN0,k,j,j1;
     if(O_gpsdat==7) printf("gps_dat_b ");
     IN0=*In0;
/* 此处，决定是否需要继续读取数据 */
     if(IN0==Iswitch || IN0==0) 
      { if(O_gpsdat==7) printf("%d ",IN0);       /* 得到原始数据 */
        getdat_o(fp,ifp,xhidat,odat_t,odat_dt,odat_n,&IN0,rt,dt,odat,
                 id_odat,iodat);
        if(O_gpsdat==7) printf("singlep ");
/* 单点定位 */
        singlep(icon,x0,Height,xhidat,odat_t,odat_dt,odat_n,id_odat,iodat,odat,
               borb,itborb,Tr1,singlex,singlet);
        if(O_gpsdat==7) printf("getdat_o ");
/* 处理毫秒问题 */
        odat_tc(ifp,xhidat,odat_t,odat1,odat_dt,odat_n,odat,
                id_odat,iodat,iodat1,singlet);
        if(O_gpsdat==7) printf("odat_tc ");
/* 组成单差、双差和三差 */
        sdtdat(ifp,xhidat,odat_t,odat_dt,odat_n,odat1,id_odat,iodat1,
               sdat,id_sdat,isdat,sdat_n,ddat,id_ddat,iddat,ddat_n,
               tdat,id_tdat,itdat,tdat_n,ref_sv);
        if(O_gpsdat==7) printf("sdtdat ");
/* 跳周检验 */
        cycleslips_test(ifp,xhidat,odat_n,odat1,odat_t,id_odat,iodat1,
                         slips_t,id_slips,ion_r,iion_r);
        if(O_gpsdat==7) printf("cycle_slips "); }
     if(IN0==0) IN0+=1; 
     *In0=IN0;
/* 为检查输出该历元的数据 */
     k=IN0; 					/* 设置 singlx and singlt */
     for(j1=0;j1<icon[0];j1++) 
      { for(j=0;j<3;j++) singlx[j+j1*3]=singlex[j+j1*3+k*icon[0]*3];
        singlt[j1]=singlet[j1+k*Ista]*1.e-6; }
     if(O_gpsdat!=0&&O_gpsdat!=-1) goto gps1d_dat;
     if(odat_n[k]!=0) 
      { t_hms(odat_t[k],&hour,&minute,&sec);
        printf("odat at time %2d %2d %4.1f ",hour,minute,sec);
        printf(", odat_n = %3d\n",odat_n[k]);}
     for (j1=0;j1<odat_n[k];j1++) 
      { printf("%8d ",id_odat[j1+I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf("%15.3f",odat1[j+j1*Isato+Isato*I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf(" %d",iodat1[j+j1*Isato+Isato*I10*Ista*k]);
        printf("\n"); }
     gps1d_dat:; if(O_gpsdat!=1&&O_gpsdat!=-1) goto gps2d_dat;
     if(sdat_n[k]!=0) 
      { t_hms(odat_t[k],&hour,&minute,&sec);
        printf("sdat at time %2d %2d %4.1f ",hour,minute,sec);
        printf(", sdat_n = %3d\n",sdat_n[k]); }
     for (j1=0;j1<sdat_n[k];j1++) 
      { printf("%8d ",id_sdat[j1+I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf("%15.3lf",sdat[j+j1*Isato+Isato*I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf(" %d",isdat[j+j1*Isato+Isato*I10*Ista*k]);
        printf("\n"); }
     gps2d_dat:; if(O_gpsdat!=2&&O_gpsdat!=-1) goto gps3d_dat;
     if(ddat_n[k]!=0) 
      { t_hms(odat_t[k],&hour,&minute,&sec);
        printf("ddat at time %2d %2d %4.1f ",hour,minute,sec);
        printf(", ddat_n = %3d\n",ddat_n[k]); }
     for (j1=0;j1<ddat_n[k];j1++) 
      { printf("%8d ",id_ddat[j1+I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf("%15.3lf",ddat[j+j1*Isato+Isato*I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf(" %d",iddat[j+j1*Isato+Isato*I10*Ista*k]);
        printf("\n"); }
     gps3d_dat:; if(O_gpsdat!=3&&O_gpsdat!=-1) goto cycle_slips;
     if(tdat_n[k]!=0) 
      { t_hms(odat_t[k],&hour,&minute,&sec);
        printf("tdat at time %2d %2d %4.1f ",hour,minute,sec);
        printf(", tdat_n = %3d",tdat_n[k]);
        printf(", dt = %f\n",odat_t[k]-odat_t[k-1]); }
     for (j1=0;j1<tdat_n[k];j1++) 
      { printf("%8d ",id_tdat[j1+I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf("%15.3lf",tdat[j+j1*Isato+Isato*I10*Ista*k]);
        for(j=0;j<Ob_num;j++)
        printf(" %d",itdat[j+j1*Isato+Isato*I10*Ista*k]);
        printf("\n"); }
     cycle_slips:; if(O_gpsdat!=-1 && O_gpsdat!=4) goto single1;
     for (j1=0,j=0;j<Islip;j++) 
      { if(fabs(slips_t[j]-odat_t[k])<0.1) 
         { if(j1==0) 
           printf("cycle slips id,t vectors : j,id_slips[],slips_t[]\n");
           j1+=1;
           printf("%4d %8d %8.1f\n",j+1,id_slips[j],slips_t[j]);} }
     single1:; if(O_gpsdat<-1 || O_gpsdat>5) goto outend;
     printf("sta   singlx and singlt\n"); 
     for(j1=0;j1<icon[0];j1++) 
      { printf(" %d  ",j1+1);
        for(j=0;j<3;j++) printf("%14.3f ",singlx[j+j1*3]);
        printf("%14.8f \n",singlt[j1]); }
     outend:;
     if(O_gpsdat==7) printf("gps_dat_end\n");
}
void cycleslips_test(ista,xhidat,odat_n,odat1,odat_t,id_odat,iodat1,
                     slips_t,id_slips,ion_r,iion_r) 
long ista,xhidat[Ista*Isat*2],odat_n[Iepoches];
double odat1[Isato*I10*Ista*Iepoches],odat_t[Iepoches];
long id_odat[I10*Ista*Iepoches];
unsigned int iodat1[Isato*I10*Ista*Iepoches];
double slips_t[Islips],ion_r[Isato*10*Ista*Iepoches];
long id_slips[Islips];
unsigned int iion_r[Isato*10*Ista*Iepoches];
/* 函数 cycle_slips_test ：如果存在跳周时，应用非差数据检查 */
{
    int i,i1,i2,j,j1,j2,k,k0,k1,k2,LL1,LL2,Co,D,Islip0;
    double LambdaL1,LambdaL2,t0,t1,t2,y0,y1,y2,dint,dt;
/* 计算参差 */ 
    Islip0=Islip; dt=odat_t[2]-odat_t[1];
for (k=1;k<Iepoches-3;k++) 			/* 每一个历元 */
 { if(O_cycle==0) printf("%8.1f %5.1f\n",odat_t[k],Cycleslip);
       for(k1=0;k1<odat_n[k];k1++) 		/* 每一个卫星 */
        { for(j=0;j<5;j++) iion_r[j+k1*Isato+Isato*10*Ista*k] = 0;
          /* 求出最后和下一个历元的有关数据，即 k0,k1 */
          k0=-1;k2=-1;
          for(k0=0;k0<odat_n[k-1];k0++)
          if(id_odat[k1+I10*Ista*k]==id_odat[k0+I10*Ista*(k-1)]) goto foundk0;
          foundk0:;
          if(k>=2) for(k2=0;k2<odat_n[k-2];k2++)
          if(id_odat[k1+I10*Ista*k]==id_odat[k2+I10*Ista*(k-2)]) goto foundk2;
          foundk2:;
          if(k0==-1) goto nextsat;
/* 双差情况 */
          if(L1+L2>=100) goto phasecode;
          if(iodat1[L1+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[L2+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[L1+k0*Isato+Isato*I10*Ista*(k-1)]==1 &&
          iodat1[L2+k0*Isato+Isato*I10*Ista*(k-1)]==1) {
             ion_r[0+k1*Isato+Isato*10*Ista*k] = 
                odat1[L1+k1*Isato+Isato*I10*Ista*k]
                - odat1[L2+k1*Isato+Isato*I10*Ista*k]*(77./60)
                - (odat1[L1+k0*Isato+Isato*I10*Ista*(k-1)]
                - odat1[L2+k0*Isato+Isato*I10*Ista*(k-1)]*(77./60));
             iion_r[0+k1*Isato+Isato*10*Ista*k] = 1; }
          phasecode:;
          if(L1<100){LL1=L1;LambdaL1=Lambda1;}else{LL1=L2;LambdaL1=Lambda2;}  
          if(L2<100){LL2=L2;LambdaL2=Lambda2;}else{LL2=L1;LambdaL2=Lambda1;} 
          if(Code==0) 
           { if(P1<100) Co=P1; else Co=P2;
             if(Co>=100) Co=C1; if(Co>=100) Co=C2;}
          else 
           { if(Code==1) Co=P1; if(Code==2) Co=P2;
             if(Code==3) Co=C1; if(Code==4) Co=C2; }
          if(LL1+Co>=100) goto phasedoppler;
          if(iodat1[LL1+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[Co+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[LL1+k0*Isato+Isato*I10*Ista*(k-1)]==1 &&
          iodat1[Co+k0*Isato+Isato*I10*Ista*(k-1)]==1) 
          { ion_r[1+k1*Isato+Isato*10*Ista*k] = 
                (odat1[LL1+k1*Isato+Isato*I10*Ista*k]
                - odat1[LL1+k0*Isato+Isato*I10*Ista*(k-1)])
                - (odat1[Co+k1*Isato+Isato*I10*Ista*k]
                - odat1[Co+k0*Isato+Isato*I10*Ista*(k-1)])/LambdaL1;
             iion_r[1+k1*Isato+Isato*10*Ista*k] = 1; }
          if(iodat1[LL2+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[Co+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[LL2+k0*Isato+Isato*I10*Ista*(k-1)]==1 &&
          iodat1[Co+k0*Isato+Isato*I10*Ista*(k-1)]==1) 
           { ion_r[2+k1*Isato+Isato*10*Ista*k] = 
                (odat1[LL2+k1*Isato+Isato*I10*Ista*k]
                - odat1[LL2+k0*Isato+Isato*I10*Ista*(k-1)])
                - (odat1[Co+k1*Isato+Isato*I10*Ista*k]
                - odat1[Co+k0*Isato+Isato*I10*Ista*(k-1)])/LambdaL2;
             iion_r[2+k1*Isato+Isato*10*Ista*k] = 1; }
          phasedoppler:;
          if(D1<100) D=D1; else D=D2;
          if(D>=100) goto nextsat;
          if(k<1) goto nextsat;
          if(iodat1[LL1+k1*Isato+Isato*I10*Ista*k]==1 &&
          iodat1[D+k1*Isato+Isato*I10*Ista*k]==1) 
          {  t2=odat_t[k]; t1=odat_t[k-1];
             if(k2!=-1) t0=odat_t[k-2]; else t0 = t1;
             y2=odat1[D+k1*Isato+Isato*I10*Ista*k]; 
             if(iodat1[D+k0*Isato+Isato*I10*Ista*(k-1)]==1)
             y1=odat1[D+k0*Isato+Isato*I10*Ista*(k-1)]; 
             y0=y1;
             if(k2!=-1&&iodat1[D+k2*Isato+Isato*I10*Ista*(k-2)]==1)
             y0=odat1[D+k2*Isato+Isato*I10*Ista*(k-2)]; 
/* 固定D 并凑整 */
             dint=line2_integrate(t0,t1,t2,y0,y1,y2);
             if(Doppler_sign==-1) dint=-dint; /* 由于多普勒反号 */
             ion_r[3+k1*Isato+Isato*10*Ista*k] = 
                (odat1[LL1+k1*Isato+Isato*I10*Ista*k] - 
                odat1[LL1+k0*Isato+Isato*I10*Ista*(k-1)]) - dint;
             iion_r[3+k1*Isato+Isato*10*Ista*k] = 1;
             ion_r[4+k1*Isato+Isato*10*Ista*k] = 
                (odat1[LL2+k1*Isato+Isato*I10*Ista*k] - 
                odat1[LL2+k0*Isato+Isato*I10*Ista*(k-1)])
                - dint*(LambdaL1/LambdaL2);
             iion_r[4+k1*Isato+Isato*10*Ista*k] = 1; }
          nextsat:;  } }
/* 检查新观测值，跳周，消去观测值 */
if(Islip==0) 
{ k=1;
       Islip = odat_n[k];
       for(k1=0;k1<odat_n[k];k1++) 
       {  id_slips[k1]=id_odat[k1+I10*Ista*k]*100;
          slips_t[k1]=odat_t[k]; } }
for(k=1;k<Iepoches-3;k++) 		/* 检查新观测值 */
{ for(k1=0;k1<odat_n[k];k1++) 
  { for(j=0,k0=0;k0<odat_n[k-1];k0++) 
    { if(id_odat[k1+I10*Ista*k]==id_odat[k0+I10*Ista*(k-1)]) j=1; }
          if(j==0) 
        {    id_slips[Islip]=id_odat[k1+I10*Ista*k]*100;
             slips_t[Islip]=odat_t[k];
             Islip +=1; } }
/* 经检验，如(D<100)应用相位多普勒法，否则应用双频法决定是否存在跳周 */
       if(D<90&&M_cyctest==1) 
       { for(k1=0;k1<odat_n[k];k1++) 
         { j=0;
           if(iion_r[3+k1*Isato+Isato*10*Ista*k]==1) 
            { if(fabs(ion_r[3+k1*Isato+Isato*10*Ista*k])>Cycleslip)
              { id_slips[Islip]=id_odat[k1+I10*Ista*k]*100+1;
                slips_t[Islip]=odat_t[k];
                Islip +=1; j=1; } }
             if(iion_r[4+k1*Isato+Isato*10*Ista*k]==1) 
            { if(fabs(ion_r[4+k1*Isato+Isato*10*Ista*k])>Cycleslip)
              { if(j==0) 
                { id_slips[Islip]=id_odat[k1+I10*Ista*k]*100+2;
                  slips_t[Islip]=odat_t[k];
                  Islip +=1; }
             else
              { id_slips[Islip-1]+=2; } } } } } 
   else
   { for(k1=0;k1<odat_n[k];k1++) 
     { if(iion_r[0+k1*Isato+Isato*10*Ista*k]==1) 
       { if(fabs(ion_r[0+k1*Isato+Isato*10*Ista*k])>Cycleslip/**0.07*/)
         { id_slips[Islip]=id_odat[k1+I10*Ista*k]*100+3;
           slips_t[Islip]=odat_t[k];
           Islip +=1; } } } }
/* 检查消失的观测值 */
       for(k0=0;k0<odat_n[k-1];k0++) 
        { for(j=0,k1=0;k1<odat_n[k];k1++) 
          { if(id_odat[k1+I10*Ista*k]==id_odat[k0+I10*Ista*(k-1)]) j=1;}
          if(j==0) 
           { id_slips[Islip]=id_odat[k0+I10*Ista*(k-1)]*100+10;
             slips_t[Islip]=odat_t[k];
             Islip +=1; }} }
/* 自检：跳周允许，但在消失前新观测值不允许两次 */
    repeat:;
for(k=0;k<Islip-1;k++) 
{  j=id_slips[k]/100; j1=id_slips[k]-j*100; if(j1!=4) goto nextk;
       k2=k; 
       for(k1=k+1;k1<Islip;k1++) 
        { i=id_slips[k1]/100; i1=id_slips[k1]-i*100;
          if(i1==10) goto nextk;
          if(i1==4 && i==j) 			/* 同一id 的两个观测值 */
           { for(i2=k2+1;i2<Islip;i2++) 
             {  id_slips[i2-1]=id_slips[i2];
                slips_t[i2-1]=slips_t[i2]; } 
             Islip -=1; 
             goto repeat; } } 
       nextk:;
    }
/* 自检：新观测值，跳周，在相同时间不允许 */
    repeat1:;
for(k=Islip-1;k>0;k--) 
{  j=id_slips[k]/100; j1=id_slips[k]-j*100; 
       k2=k; 
       for(k1=k-1;k1>=0;k1--) 
        { i=id_slips[k1]/100; i1=id_slips[k1]-i*100;
          if(fabs(slips_t[k]-slips_t[k1])>0.4) goto nextkk;
          /* 相同时间和相同id有两个观测值 */
          if(i==j) 
           { if(j1<=i1) k2=k1;
             for(i2=k2+1;i2<Islip;i2++) 
              { id_slips[i2-1]=id_slips[i2];
                slips_t[i2-1]=slips_t[i2]; }