📄 rinex.c
字号:
Islip -=1;
goto repeat1; }}
nextkk:;
}
/* 自检：不允许同样两个 id_slips[] 和 slips_t[] */
repeat2:;
for(k=Islip-1;k>0;k--)
{ for(k1=k-1;k1>=0;k1--)
{ if(fabs(slips_t[k]-slips_t[k1])<0.4)
{ if(id_slips[k]==id_slips[k1]) /* 找到两个相同的，取消第 k个元素 */
{ for(i2=k+1;i2<Islip;i2++)
{ id_slips[i2-1]=id_slips[i2];
slips_t[i2-1]=slips_t[i2]; }
Islip -=1;
goto repeat2; }}}}
/* 自检：在相邻时间上不允许两个'相同'的id_slips[]，即 id 30700 - 29065 and id 30703 - 29066 */
if(Islip>=Islips) {printf("dimension of Islip error\n"); exit(0);} /* Islip 检查 */
if(O_cycle <= 1 && Islip>Islip0) /* 为控制输出数据 */
{ printf("cycle slips id and time vectors : i,id_slips[],slips_t[]\n");
for (i=0;i<Islip;i++)
printf("%4d %8d %8.1f\n",i+1,id_slips[i],slips_t[i]); }
if(O_cycle <= 0)
{ printf("cycle slips id and time vectors : i,id_slips[],slips_t[]\n");
for (i=0;i<Islip;i++)
printf("%4d %8d %8.1f\n",i+1,id_slips[i],slips_t[i]); }
endof_f:;}
double line2_integrate(t0,t1,t2,y0,y1,y2)
/* 此函数应用t0 t1 t2 的多普勒数据预测相位从t1到t2的增加，比较其相位增加以检查是否存在跳周。为此f(t)从t1到t2凑整，t1到t2 的平均值为凑整值(t1,t2) f(t)*dt/(t2-t1) */
double t0,t1,t2,y0,y1,y2;
{ double b,c=0.,var;
if(fabs(t2-t1)<1.e-2 && fabs(t1-t0)<1.e-2) return(y2);
if(fabs(t1-t0)<1.e-2) {b=(y2-y1)/(t2-t1);c=0.;}
else
{ c = ((y2-y1)/(t2-t1) - (y1-y0)/(t1-t0))/(t2-t0);
b = (y2-y1)/(t2-t1) - c*(t2-t1); }
var = y1*(t2-t1)+(b/2.)*(t2-t1)*(t2-t1)
+(c/3.)*(t2-t1)*(t2-t1)*(t2-t1);
return(var);}
void id_obs(xhidat,iC1,iC2,iL1,iL2,iP1,iP2,iD1,iD2)
/* 利用xhidat[i] 设置 C1, C2, P1, P2, L1, L2, D1, D2 的值，1=L1，2=L2，3=C1 */
long xhidat[Ista*Isat*2];
int *iC1,*iC2,*iL1,*iL2,*iP1,*iP2,*iD1,*iD2;
{ int i,j,i1;
*iC1=100;*iC2=100;*iL1=100;*iL2=100;*iP1=100;*iP2=100;*iD1=100;*iD2=100;
j = (int)xhidat[0];
for(i=1;i<=j;i++)
{ i1 = xhidat[i];
switch(i1)
{ case 1: { *iL1 = i-1;break;}
case 2: { *iL2 = i-1;break;}
case 3: { *iC1 = i-1;break;}
case 4: { *iC2 = i-1;break;}
case 5: { *iP1 = i-1;break;}
case 6: { *iP2 = i-1;break;}
case 7: { *iD1 = i-1;break;}
case 8: { *iD2 = i-1;break;}}}}
void sdtdat(ista,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)
/* 此函数利用 GPS RINEX 数据或接收机经时间改正的数据，组成单差sd 、双差dd和三差td */
long ista,ref_sv;
long xhidat[Ista*Isat*2];
long odat_n[Iepoches];
double odat1[Isato*I10*Ista*Iepoches],odat_t[Iepoches];
double odat_dt[Ista*Iepoches];
long id_odat[I10*Ista*Iepoches];
unsigned int 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];
{ long i,i1,i2,j,j1,j2,j3,k,k1,kp1,k2,ob_num;
double rit,sec;
int hour,minute,I1,I2,J1,J2;
ob_num = (int)xhidat[0]; /*printf("%d\n",ob_num);*/
/* 组成单差数据 ( 径向基线情况) */
for(k=0;k<Iepoches;k++) /* 每一历元 */
{ for(i=0,k1=0;k1<odat_n[k];k1++) /* 每一卫星 */
{ id_i12j12(&id_odat[k1+I10*Ista*k],&I1,&I2,&J1,&J2);
i1 = I1-1; j1 = J1;
if(i1==0) /* sta 0 作为参考测站 */
{ for(k2=k1+1;k2<odat_n[k];k2++)
{ id_i12j12(&id_odat[k2+I10*Ista*k],&I1,&I2,&J1,&J2);
i2 = I1-1; j2 = J1;
if(i2!=0 && j1==j2) /* 非参考测站和同样卫星 */
{ id_sdat[i+k*I10*Ista]=(i1+1)*10000+(i2+1)*100+j2;
for(j=0;j<ob_num;j++)
{ if(iodat1[j+k1*Isato+Isato*I10*Ista*k]==1 &&
iodat1[j+k2*Isato+Isato*I10*Ista*k]==1)
{ sdat[j+i*Isato+Isato*I10*Ista*k]=
odat1[j+k2*Isato+Isato*I10*Ista*k]-
odat1[j+k1*Isato+Isato*I10*Ista*k];
isdat[j+i*Isato+Isato*I10*Ista*k]=1;}
else
{ sdat[j+i*Isato+Isato*I10*Ista*k]=0.;
isdat[j+i*Isato+Isato*I10*Ista*k]=0;}}
i += 1;}}}}
sdat_n[k]=i;
}
/* 为检查输出单差数据 */
if(O_sdtdat==0)
{ for(k=0;k<Iepoches;k++)
{ 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");}}}
/* 组成双差数据 ( 径向基线情况) */
for(k=0;k<Iepoches;k++) /* 每一历元 */
{ for(i=0,k1=0;k1<sdat_n[k];k1++) /* 每一单差 sdat */
{ id_i12j12(&id_sdat[k1+I10*Ista*k],&I1,&I2,&J1,&J2);
/*printf("%8d %2d %2d %2d %2d\n",
id_sdat[k1+I10*Ista*k],I1,I2,J1,J2);*/
i1 = I1*100+I2; j1 = J1;
if(j1==ref_sv) /* 如果相同卫星是参考卫星 */
{ for(k2=0;k2<sdat_n[k];k2++)
{ id_i12j12(&id_sdat[k2+I10*Ista*k],&I1,&I2,&J1,&J2);
i2 = I1*100+I2; j2 = J1;
if(i1==i2 && j2!=j1) /* 非参考卫星和相同基线 */
{ id_ddat[i+k*I10*Ista]=i1*10000+j1*100+j2;
for(j=0;j<ob_num;j++)
{ if(isdat[j+k1*Isato+Isato*I10*Ista*k]==1 &&
isdat[j+k2*Isato+Isato*I10*Ista*k]==1)
{ ddat[j+i*Isato+Isato*I10*Ista*k]=
sdat[j+k2*Isato+Isato*I10*Ista*k]-
sdat[j+k1*Isato+Isato*I10*Ista*k];
iddat[j+i*Isato+Isato*I10*Ista*k]=1;}
else
{ ddat[j+i*Isato+Isato*I10*Ista*k]=0.;
iddat[j+i*Isato+Isato*I10*Ista*k]=0;}}
i += 1;}}}}
ddat_n[k]=i;
}
/* 为检查输出双差数据 */
if(O_sdtdat==0)
{ for(k=0;k<Iepoches;k++)
{ 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");}}}
/* 组成三差数据 ( 每两个连续的双差数据) */
tdat_n[0] = 0;
for(k=0;k<Iepoches-1;k++) /* 每一历元 */
{ i=0; tdat_n[k+1] = 0; kp1 = k+1;
if(ddat_n[k]!=0 & ddat_n[k+1]!=0)
{ for(k1=0;k1<ddat_n[k];k1++) /* k 历元双差 ddat */
{ i1 = id_ddat[k1+I10*Ista*k];
for(k2=0;k2<ddat_n[k+1];k2++) /* k+1 历元双差 ddat */
{ i2 = id_ddat[k2+I10*Ista*kp1];
if(i1==i2) /* 如果两个 id_ddat 相同 */
{ id_tdat[i+kp1*I10*Ista]=id_ddat[k1+I10*Ista*k];
for(j=0;j<ob_num;j++)
{ if(iddat[j+k1*Isato+Isato*I10*Ista*k]==1 &&
iddat[j+k2*Isato+Isato*I10*Ista*kp1]==1)
{ tdat[j+i*Isato+Isato*I10*Ista*kp1]=
ddat[j+k2*Isato+Isato*I10*Ista*kp1]-
ddat[j+k1*Isato+Isato*I10*Ista*k];
itdat[j+i*Isato+Isato*I10*Ista*kp1]=1;}
else
{ tdat[j+i*Isato+Isato*I10*Ista*kp1]=0.;
itdat[j+i*Isato+Isato*I10*Ista*kp1]=0;}}
i += 1;}}}}
tdat_n[kp1]=i;
}
/* 为检查输出三差数据 */
if(O_sdtdat==0)
{ for(k=0;k<Iepoches;k++)
{ 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"); }}}}
void odat_tc(ista,xhidat,odat_t,odat1,odat_dt,odat_n,odat,id_odat,iodat,iodat1
,singlet)
/* 此函数内插 GPS RINEX 数据，将所有的观测值归算为同一时间 odat_t，删去零向量 odat_dt1 */
long ista;
long xhidat[Ista*Isat*2];
long odat_n[Iepoches];
double odat[Isato*I10*Ista*Iepoches],odat1[Isato*I10*Ista*Iepoches];
double odat_dt[Ista*Iepoches],odat_t[Iepoches];
long id_odat[I10*Ista*Iepoches];
unsigned int iodat[Isato*I10*Ista*Iepoches],iodat1[Isato*I10*Ista*Iepoches];
double singlet[Ista*Iepoches];
{
long i,i1,i2,i3,i4,j,j1,j2,j3,k,km1,kp1,k1,k2,ob_num,id;
int I1,I2,J1,J2;
double rit,sec,t0,t1,t2,y0,y1,y2,var;
double dct,dct1,dct2,dct3;
int hour,minute;
ob_num = (int)xhidat[0];
/* 将原始数据 odat 和 iodat 放入 odat1 和 iodat1 */
for(k=0;k<Iepoches;k++)
{ if(odat_n[k]!=0)
{ for(j1=0;j1<odat_n[k];j1++)
{ for(j=0;j<ob_num;j++)
{ odat1[j+j1*Isato+Isato*I10*Ista*k] =
odat[j+j1*Isato+Isato*I10*Ista*k];
iodat1[j+j1*Isato+Isato*I10*Ista*k]=
iodat[j+j1*Isato+Isato*I10*Ista*k];}}}}
if(ReceivC==0) goto function_end;
/* 内插数据，使新数据无时间改正 (odat_dt==0) */
for(k=0;k<Iepoches-1;k++) /* 每一历元 */
{ if(odat_n[k]==0) goto function_end;
for(j1=0;j1<odat_n[k];j1++) /* 每一观测量 */
{ id = id_odat[j1+I10*Ista*k];
id_i12j12(&id,&I1,&I2,&J1,&J2);
/*printf("%8d %2d %2d %2d %2d %lf\n",
id_odat[j1+I10*Ista*k],I1,I2,J1,J2,
odat_dt[I1-1+k*Ista]);*/
dct=odat_dt[I1-1+k*Ista];
if(ReceivC==2) dct=singlet[I1-1+k*Ista]*1.e-6;/* 这里只用原始数据 I1,J1 */
if(fabs(dct)>=0.000001) /* dt != 0 以相同id求最后历元数据 */
{ i1 = -1; i2 = -1;
if(k!=0)
{ km1 = k-1;
for(j2=0;j2<odat_n[km1];j2++)
if(id==id_odat[j2+I10*Ista*km1]) i1=j2;
t_hms(odat_t[km1],&hour,&minute,&sec);
t0 = hour*3600.+minute*60.+sec;
dct1 = odat_dt[i1+km1*Ista];
if(ReceivC==2) dct1=singlet[i1+km1*Ista]*1.e-6;
if(i1!=-1) t0 += dct1; }
if(odat_n[k+1]!=0) /* 以相同id求下一个历元数据 */
{ kp1 = k+1;
for(j2=0;j2<odat_n[kp1];j2++)
if(id==id_odat[j2+I10*Ista*kp1]) i2=j2;
t_hms(odat_t[kp1],&hour,&minute,&sec);
t2 = hour*3600.+minute*60.+sec;
dct2 = odat_dt[i2+kp1*Ista];
if(ReceivC==2) dct2=singlet[i2+kp1*Ista]*1.e-6;
if(i2!=-1) t2 += dct2; }
if(i1==-1 && i2==-1) /* 数据是单个，不能内插，设置相关状态 iodat=0 */
{ for(j=0;j<ob_num;j++) iodat1[j+j1*Isato+Isato*I10*Ista*k]=0; }
t_hms(odat_t[k],&hour,&minute,&sec); /* 内插时间 rit */
rit = hour*3600.+minute*60.+sec;
t1 = rit + dct;
for(j=0;j<ob_num;j++) /* 每一观测值 */
{ if(iodat[j+j1*Isato+Isato*I10*Ista*k]==1)
{ y1 = odat[j+j1*Isato+Isato*I10*Ista*k];
if(i1!=-1)
{ if(odat[j+i1*Isato+Isato*I10*Ista*km1]==0) i1=-1;
else y0 = odat[j+i1*Isato+Isato*I10*Ista*km1]; }
if(i2!=-1)
{ if(odat[j+i2*Isato+Isato*I10*Ista*kp1]==0) i2=-1;
else y2 = odat[j+i2*Isato+Isato*I10*Ista*kp1]; }
if(i1!=-1 && i2!=-1)
var=line2_interpo(t0,t1,t2,y0,y1,y2,rit);
if(i1!=-1 && i2==-1) var = line_interpo(t0,t1,y0,y1,rit);
if(i1==-1 && i2!=-1) var = line_interpo(t1,t2,y1,y2,rit);
if(i1==-1 && i2==-1)
{ var=0.; iodat1[j+j1*Isato+Isato*I10*Ista*k]=0; }
odat1[j+j1*Isato+Isato*I10*Ista*k] = var;}}}
if(DopplerI==1) /* 应用多普勒数据内插，内插的时间 rit */
{ t_hms(odat_t[k],&hour,&minute,&sec);
rit = hour*3600.+minute*60.+sec;
t1 = rit + dct;
for(j=0;j<ob_num;j++) /* 伪距和相位 */
{ if(iodat[j+j1*Isato+Isato*I10*Ista*k]==1)
{ y1 = odat[j+j1*Isato+Isato*I10*Ista*k];
y2 = odat1[PD1+j1*Isato+Isato*I10*Ista*k];
if(j==PC1||j==PC2)
{ var=y1-y2*Lambda1*dct*Doppler_sign;
odat1[j+j1*Isato+Isato*I10*Ista*k] = var;}
if(j==L1)
{ var=y1-y2*dct*Doppler_sign;
odat1[j+j1*Isato+Isato*I10*Ista*k] = var;}
if(j==L2) { rit=1./Lambda2; t1=Lambda1;
var=y1-y2*dct*Doppler_sign*(60./77.);
odat1[j+j1*Isato+Isato*I10*Ista*k] = var;}}}}/*若(DopplerI==1)结束*/
}}
function_end:; /* 为检查输出数据 */
}
double line2_interpo(t0,t1,t2,y0,y1,y2,t)
/* f(t) = a+b(t-t1)+c*(t-t1)*(t-t1) 对时间t 内插 */
double t0,t1,t2,y0,y1,y2,t;
{ double b,c,var;
c = ((y2-y1)/(t2-t1) - (y1-y0)/(t1-t0))/(t2-t0);
b = (y2-y1)/(t2-t1) - c*(t2-t1);
var = y1 + (t-t1)*b + c*(t-t1)*(t-t1);
return(var);}
double line_interpo(t1,t2,y1,y2,t)
/* 线性内插 f(t) = a + b(t-t1) */
double t1,t2,y1,y2,t;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -