📄 rinex.c
字号:
{ double var,dydt;
dydt = (y2-y1)/(t2-t1);
var = y1 + (t-t1)*dydt;
return(var);}
void id_i12j12(id,i1,i2,j1,j2)
/* 改变观测值的 id 放入 sta_i1, sta_i2, sat_j1, sat_j2 */
long *id;
int *i1,*i2,*j1,*j2;
{ if(*id>=1000000) /* 双差 id */
{ *i1 = (int)(*id/1000000);
*i2 = (int)((*id-*i1*1000000)/10000);
*j1 = (int)((*id-*i1*1000000-*i2*10000)/100);
*j2 = (int)(*id-*i1*1000000-*i2*10000-*j1*100); }
if(10000<*id && *id<1000000) /* 单差 id */
{ *i1 = (int)(*id/10000);
*i2 = (int)((*id-*i1*10000)/100);
*j1 = (int)(*id-*i1*10000-*i2*100);
*j2 = 0; }
if(*id<10000) /* 原始差 id */
{ *i1 = (int)(*id/100);
*j1 = (int)(*id-*i1*100);
*i2 = 0;
*j2 = 0; }}
void getdat_o(fp,ista,xhidat,odat_t,odat_dt,odat_n,In0,rt,dt,odat,id_odat,iodat)
/* 此函数得到所有测站在Iepoches的 GPS RINEX 数据 。若 In0==0，由 epoch rt 开始读数据，直到下一个 Iepoches。否则将数据从( In0--Iepoches)移到( 0--Iepoches-In0)，继续读数据 */
long ista;
FILE *fp[Ista];
long xhidat[Ista*Isat*2];
long odat_n[Iepoches],*In0;
double odat[Isato*I10*Ista*Iepoches],odat_t[Iepoches],rt,dt;
double odat_dt[Ista*Iepoches];
long id_odat[I10*Ista*Iepoches];
unsigned int iodat[Isato*I10*Ista*Iepoches];
{ long i,i1,i2,j,j1,j2,j3,k,k1,k2,ob_num;
long xdidat1[Isat*Ista],io[Ista];
double xddat1[I10*Isato*Ista],rit,sec;
int hour,minute;
unsigned int ixddat1[I10*Isato*Ista];
ob_num = (int)xhidat[0];
/* 移动数据 */
i2=Iepoches;
if(*In0==0) {i1=0; rit = rt-dt;}
else
{ for(j2=0,k=*In0-1;k<i2;k++)
{ k1 = k - (*In0-1);
for(i=0;i<ista;i++)
{ odat_dt[i+Ista*k1] = odat_dt[i+Ista*k];}
odat_n[k1] = odat_n[k];
odat_t[k1] = odat_t[k];
j2+=odat_n[k];
for (j1=0;j1<odat_n[k];j1++)
{ id_odat[j1+I10*Ista*k1] = id_odat[j1+I10*Ista*k];
for(j=0;j<ob_num;j++)
odat[j+j1*Isato+Isato*I10*Ista*k1]= odat[j+j1*Isato+Isato*I10*Ista*k];
for(j=0;j<ob_num;j++)
iodat[j+j1*Isato+Isato*I10*Ista*k1]=iodat[j+j1*Isato+Isato*I10*Ista*k];}}
rit = odat_t[Iepoches-1];
i1 = Iepoches - *In0 + 1; *In0 = 1; }
/* 读取新值 */
if(dt<=0.49) dt = xhidat[14]+xhidat[15]/10.;
for(k=i1;k<i2;k++)
{ rit += dt;
getdat(fp,ista,&rit,xhidat,xdidat1,io,xddat1,ixddat1);
odat_t[k] = rit; odat_n[k] = 0;
for(j2=0,i=0;i<ista;i++)
{ odat_dt[i+k*Ista] = -rit + xdidat1[3+Isat*i]*3600. + xdidat1[4+Isat*i]*60.
+ xdidat1[5+Isat*i] + xdidat1[6+Isat*i]/10000000.;
if(io[i]==2) xdidat1[9+Isat*i]=0;
odat_n[k] += xdidat1[9+Isat*i];
for (j1=0;j1<xdidat1[9+Isat*i];j1++)
{ id_odat[j2+I10*Ista*k] = (i+1)*100+xdidat1[9+Isat*i+1+j1];
for(j=0;j<ob_num;j++)
{odat[j2*Isato+j+k*Ista*Isato*I10] = xddat1[j+j1*ob_num+Isato*I10*i];
iodat[j2*Isato+j+k*Ista*Isato*I10] = ixddat1[j+j1*ob_num+Isato*I10*i]; }
j2 += 1;}}} /* 为检查输出原始数据 */
}
void t_hms(t,hour,minute,sec)
/* 换算时间 t 为时、分、秒 */
double t,*sec;
int *hour,*minute;
{ *hour = (int)(t/3600);
*minute = (int)((t-*hour*3600)/60);
*sec = t - *hour*3600 - *minute*60;}
void getdat(fp,ista,rt,xhidat,xdidat1,io,xddat1,ixddat1)
long ista;
FILE *fp[Ista];
long xhidat[Ista*Isat*2];
long xdidat1[Isat*Ista],io[Ista];
double xddat1[I10*Isato*Ista];
unsigned int ixddat1[I10*Isato*Ista];
double *rt;
/* 函数 getdat 应用函数 r_rinexd 得到所有测站在时间 rt的gps数据，若某测站无数据，则取下一个历元的数据，并返回到时间rt */
{ int i,i1,i2,j,j1,j2,k;
long xdidat[Isat],i3[Ista],ob_num;
double xddat[I10*Isato],f3[Ista],f1,f2,f,rit;
unsigned int ixddat[I10*Isato];
/* 设置初始值 */
rit = *rt; for(i=0;i<ista;i++) i3[i] = ista;
repeat:;
/* 读每一测站数据并放在一块 */
for (i=0;i<ista;i++)
{ if(i!=i3[i])
{ r_rinexd(fp[i],xhidat,rit,xdidat,xddat,ixddat,&i2);
io[i] = i2;
if(xdidat[7]==0) {io[i]=2; }
ob_num = (int)xhidat[0]; k=ob_num;
for(j=0;j<10+xdidat[9];j++) xdidat1[j+Isat*i]=xdidat[j];
for (j1=0;j1<xdidat1[9+Isat*i];j1++)
{ for(j=0;j<ob_num;j++) xddat1[j+j1*k+Isato*I10*i]=xddat[j+j1*k];
for(j=0;j<ob_num;j++)ixddat1[j+j1*k+Isato*I10*i]=ixddat[j+j1*k]; }}
else io[i]=1; }
/* printf("xx %f %d %d\n",rit,io[0],io[1]);*/
/* 如果某数据有缺，读下一个历元数据 */
for(i1=0,i=0;i<ista;i++)
{ if(io[i]==0)
{ i1+=1;
f = xdidat1[3+Isat*i]*3600.+xdidat1[4+Isat*i]*60+xdidat1[5+Isat*i];
f1 = xdidat1[6+Isat*i]/10000000.;
if(0.4<f1 && f1<0.6) f += 0.5;
if(0.9<f1 && f1<1.1) f += 1.;
f3[i] = f;}
else f3[i] = rit; }
if(i1>=1)
{ f = f3[0];
for(i=1;i<ista;i++) if(f<f3[i]) f=f3[i];
for(i=0;i<ista;i++) if(fabs(f-f3[i])<0.1) i3[i] = i;
/* printf("%f %d \n",f,i1); */
*rt = f; rit = f; goto repeat; }
/* 检查双频相位观测值，如果丢失，则删除该观测值 */
if(Frequency_n==2)
{ once1:;
for (i=0;i<ista;i++)
{for (j1=0;j1<xdidat1[9+Isat*i];j1++)
{if(ixddat1[L1+j1*k+Isato*I10*i]==0||ixddat1[L2+j1*k+Isato*I10*i]==0)
{for (j2=j1;j2<xdidat1[9+Isat*i]-1;j2++)
{for(j=0;j<ob_num;j++)
{ xddat1[j+j2*k+Isato*I10*i]=xddat1[j+(j2+1)*k+Isato*I10*i];
ixddat1[j+j2*k+Isato*I10*i]=ixddat1[j+(j2+1)*k+Isato*I10*i];}
xdidat1[j2+10+Isat*i]=xdidat1[j2+10+1+Isat*i];}
xdidat1[9+Isat*i]-=1;
goto once1; }}}}
/* 为检查输出数据 */
}
void r_rinexd(fp,xhidat,rit,xdidat,xddat,ixddat,Id)
/* 函数 r_rinexd 读取rinex 文件的数据部分 */
FILE *fp;
long *Id;
long xdidat[Isat],xhidat[Ista*Isat*2];
double xddat[I10*Isato];
unsigned int ixddat[I10*Isato];
double rit;
{ int i,j,j1,j2,k,k1=0,c,c1,c2,c3,c4,ob_num;
double f,f1;
char str[Iline],str1[17],str2[4];
*Id = 0; ob_num = (int)xhidat[0];
for(i=0;i<I10*Isato;i++) ixddat[i]=0; /* 读数据部分 */
loop:; k1+=1; /* 第一行 */
xdidat[7]=1;
fgets(str,16,fp);
if(str[2]!=' ' && str[3]==' ')
{ sscanf(str,"%d %d %d %d %d", &c,&c1,&c2,&c3,&c4);
xdidat[0]=(long)c; xdidat[1]=(long)c1;
xdidat[2]=(long)c2; xdidat[3]=(long)c3;
xdidat[4]=(long)c4;}
else
{ xdidat[7] = 0; goto finished; }
fgets(str1,12,fp); f = atof(str1);
xdidat[5] = (long)f;
xdidat[6] = (long)((f-xdidat[5])*10000000);
if (xdidat[6]>=9000000) /* 使 xdidat[6] 变小 */
{ xdidat[6] -= 10000000;
xdidat[5] += 1; }
if (xdidat[5]>=60) /* 避免 xdidat[5] >= 60 */
{ xdidat[5] -= 60;
xdidat[4] += 1; }
if (xdidat[4]>=60) /* 避免 xdidat[4] >= 60 */
{ xdidat[4] -= 60;
xdidat[3] += 1; }
for (i=8;i<10;i++)
{ fgets(str2,4,fp);xdidat[i]=atoi(str2);}
for (i=10;i<10+xdidat[9];i++)
{ fgets(str2,4,fp);xdidat[i]=atoi(str2); }
c = getline(fp,str);
/* 时间检查 */
f = xdidat[3]*3600.+xdidat[4]*60.+xdidat[5]+xdidat[6]/10000000.;
f1 = rit;
if (f-0.02 < f1 && f1 < f+0.02)
{ *Id = 1; }
else
{ if (f1 >= f+0.2)
{ ob_num = (int)xhidat[0];
for (i=0;i<xdidat[9];i++)
c = getline(fp,str); /* 观测值第一行 */
if(ob_num>=6)
for (i=0;i<xdidat[9];i++)
c = getline(fp,str); /* 观测值第二行 */
goto loop; }}
ob_num = (int)xhidat[0]; /* 每一卫星 */
k = ob_num;
for (i=0;i<xdidat[9];i++)
{ c = getline(fp,str); /* 观测值第一行 */
j2=ob_num; if(j2>=6) j2=5;
for (j=0;j<j2;j++)
{ if(str[10+j*16]==' ') xddat[j+i*k] = 0.;
else
{ for (j1=0;j1<14;j1++) str1[j1]=str[j1+j*16];
xddat[j+i*k] = atof(str1);
ixddat[j+i*k]=1;} }
if(ob_num>=6)
{ c = getline(fp,str); /* 观测值第二行 */
for (j=5;j<ob_num;j++)
{ if(str[10+(j-5)*16]==' ') xddat[j+i*k] = 0.;
else
{ for (j1=0;j1<14;j1++) str1[j1]=str[j1+(j-5)*16];
xddat[j+i*k] = atof(str1);
ixddat[j+i*k]=1;}}}}
finished:;
}
void r_rinexh(fp,ifp,xhidat)
/* 函数r_rinexh 读取rinex 文件的字头部分 */
FILE *fp[Ista];
long ifp,xhidat[2*Isat*Ista];
{
int i,i1,j,k,k1,c,c1,c2,c3,c4,c5;
double f,f1,f2,d,dtstep;
char str[Iline];
for(i1=0;i1<ifp;i1++)
{ for (i=0;i<2*Isat;i++) xhidat[i+i1*2*Isat] = 0;
nextline:; /* 读字头部分 */
k1 = getline(fp[i1],str);
if(O_rinexh==0)
fprintf(stdout,"%s",str);
if(strstr(str,"VERSION"))
{ sscanf(str,"%d",&c); xhidat[11+i1*2*Isat]=(long)c;
goto nextline; }
if(strstr(str,"APPROX"))
{ sscanf(str,"%lf %lf %lf",&f,&f1,&f2);
xhidat[23+i1*2*Isat]=(long)f;
xhidat[24+i1*2*Isat]=(long)((f-xhidat[23+i1*2*Isat])*1000+0.5);
xhidat[25+i1*2*Isat]=(long)f1;
xhidat[26+i1*2*Isat]=(long)((f1-xhidat[25+i1*2*Isat])*1000+0.5);
xhidat[27+i1*2*Isat]=(long)f2;
xhidat[28+i1*2*Isat]=(long)((f2-xhidat[27+i1*2*Isat])*1000+0.5);
goto nextline;}
if(strstr(str,"ANTENNA"))
{ sscanf(str,"%lf %lf %lf",&f,&f1,&f2);
xhidat[29+i1*2*Isat]=(long)f;
xhidat[30+i1*2*Isat]=(long)((f-xhidat[29+i1*2*Isat])*1000)+0.5;
xhidat[31+i1*2*Isat]=(long)f1;
xhidat[32+i1*2*Isat]=(long)((f1-xhidat[31+i1*2*Isat])*1000+0.5);
xhidat[33+i1*2*Isat]=(long)f2;
xhidat[34+i1*2*Isat]=(long)((f2-xhidat[33+i1*2*Isat])*1000+0.5);
goto nextline;}
if(strstr(str,"WAVELENGTH"))
{ sscanf(str,"%d %d",&c,&c1);
xhidat[12+i1*2*Isat]=(long)c; xhidat[13+i1*2*Isat]=(long)c1;
goto nextline; }
if(strstr(str,"TYPES"))
{ sscanf(str,"%d",&c);xhidat[0+i1*2*Isat]=(long)c;
for (i=1;i<=xhidat[0+i1*2*Isat];i++)
{ c = str[4+i*6]; k = str[5+i*6];
switch(c)
{ case 'L': {xhidat[i+i1*2*Isat] = 1;break;}
case 'C': {xhidat[i+i1*2*Isat] = 3;break;}
case 'P': {xhidat[i+i1*2*Isat] = 5;break;}
case 'D': {xhidat[i+i1*2*Isat] = 7;break;}
case 'T': {xhidat[i+i1*2*Isat] = 9;break;} }
if (k=='2') xhidat[i+i1*2*Isat] += 1; }
goto nextline; }
if(strstr(str,"INTERVAL"))
{ for(k=0,i=0;i<16;i++) if(str[i]=='.') k=1;
if(k==0) {sscanf(str,"%d",&c); xhidat[14+i1*2*Isat]=(long)c;}
if(k==1) {sscanf(str,"%lf",&d); xhidat[14+i1*2*Isat]=(long)d;
xhidat[15+i1*2*Isat]=(long)((d-xhidat[14+i1*2*Isat])*10.);}
goto nextline; }
if(strstr(str,"FIRST"))
{ sscanf(str,"%d %d %d %d %d %lf",&c,&c1,&c2,&c3,&c4,&f);
xhidat[16+i1*2*Isat]=(long)c; xhidat[17+i1*2*Isat]=(long)c1;
xhidat[18+i1*2*Isat]=(long)c2; xhidat[19+i1*2*Isat]=(long)c3;
xhidat[20+i1*2*Isat]=(long)c4; xhidat[21+i1*2*Isat]=(long)f;
xhidat[22+i1*2*Isat]=(long)((f-xhidat[21+i1*2*Isat])*1000000);
goto nextline; }
if(strstr(str,"HEADER")) goto dataend;
if(k1>10) goto nextline;;
dataend:; /* 输出 xhidat[] */
}
/* 检查是否数据有相同的结构 */
for (i=0,i1=1;i1<ifp;i1++)
if(xhidat[0]!=xhidat[0+i1*2*Isat]) i+=1 ;
if(i==0)
{ for(j=0,k=1;k<=xhidat[0];k++)
{ for(i1=1;i1<ifp;i1++) if(xhidat[k]!=xhidat[k+i1*2*Isat]) j+=1 ;}
if(j==0) printf("All rinex data files have the same structure\n");
else printf("Rinex data files have not the same structure\n"); }
else printf("Rinex data files have different obs_number\n");}
int getline(fp,s)
/* 从文件fp得到行写入s，返回s的长度。s最大长度为 Mline，现定义为100*/
FILE *fp;
char s[Iline];
{ int c,i=0;
while ((c=fgetc(fp)) != EOF)
{ s[i] = c; i += 1;
if (c=='\n') { s[i]='\0'; return(i); } }
return(c);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -