igs.c

用C语言写的用来读IGS 精密星历子程序！

igs.c

void igs(fp,orb,igs_h,I_la) 
FILE *fp;
long igs_h[8],I_la;
double orb[8*Isat*I96];
{ igsdat(fp,stdout,igs_h,orb,I_la);}
void igsdat(fp,fp1,igs_h,orb,I_la) 
FILE *fp,*fp1;
long igs_h[8],I_la;
double orb[8*Isat*I96];
/* 函数 igsdat 读取 igs gps 精密星历 */
{       long i,i1,i2,j,k,l,l1=5,deltat=900,igs_sat_num,itb;
        long idat[8],iidat[Isat];
    	double ddat[8*Isat],dat[I_LA],rit,rr,rit1,rit2,f1,f2,f=0.;;
/* 读取数据 */
        r_igs_h(fp,fp1,igs_h); 
        igs_sat_num = igs_h[0]; igs_h[7] = I96;
        for(j=0;j<I96;j++) 
         { r_igs_d(fp,fp1,igs_sat_num,idat,ddat,&l);
           if(l>70) l1=7;
           if(l1==5) 
            { for(k=0;k<igs_sat_num;k++) 
               { orb[0+k*8+j*8*Isat]=ddat[0+k*8];
                 for(i=1;i<4;i++) orb[i+k*8+j*8*Isat]=ddat[i+k*8]*1000; 
                 if(Igs_v<2)
                 for(i=4;i<7;i++) orb[i+k*8+j*8*Isat]=0.;
                 if(Igs_v==2)
                 for(i=4;i<7;i++) orb[i+k*8+j*8*Isat]=ddat[i+1+k*8]*1000;
                 orb[7+k*8+j*8*Isat]=ddat[4+k*8]; }} 
            else 
            { for(k=0;k<igs_sat_num;k++) 
               { orb[0+k*8+j*8*Isat]=ddat[0+k*8];
                 for(i=1;i<7;i++) orb[i+k*8+j*8*Isat]=ddat[i+k*8]*1000; 
                 orb[7+k*8+j*8*Isat]=0.;}}}
/* 如果需要内查速度 */
     if(l1==5 && Igs_v!=2) 
     { for(j=0;j<I96;j++) 
        {  rit = igs_h[4]*3600.+j*deltat+1;
           rit1=rit-1-0.5; rit2=rit-0.5;
           i1 = I_la/2; if(j<=i1) i2=0; else i2+=1;
           if(j>=I96-i1-1) i2=I96-I_la-1;
           itb = igs_h[4]*3600+i2*deltat; 
           for (k=0;k<igs_sat_num;k++) 
            { for(i=1;i<=3;i++) 
              {  for(i1=i2;i1<=i2+I_la;i1++) dat[i1-i2]=orb[i+k*8+i1*8*Isat];
                 f1=lagrange(dat,I_la,deltat,rit1,itb);
                 f2=lagrange(dat,I_la,deltat,rit2,itb);
                 orb[i+3+k*8+j*8*Isat]=f2-f1; } } } }
/* 在 l1=5 的情况,存在未定义的钟差，此数据不能采用朗格日法内查，采用简单平均 */
     if(l1==5) 
     { for (k=0;k<igs_sat_num;k++) {ddat[k] = 0.;iidat[k]=0;}
     for(j=0;j<I96;j++) 
     { for (k=0;k<igs_sat_num;k++) 
       { if(orb[7+k*8+j*8*Isat]<999.) 
         { iidat[k]+=1; ddat[k]+=orb[7+k*8+j*8*Isat]; }}}
     for (k=0;k<igs_sat_num;k++) 
     { if(iidat[k]!=0) ddat[k]=ddat[k]/iidat[k];}
     for(j=0;j<I96;j++) 
     { for (k=0;k<igs_sat_num;k++) 
       { if(orb[7+k*8+j*8*Isat]>999.) 
         { orb[7+k*8+j*8*Isat]= ddat[k]; }}}}
/* 为控制输出数据 */
        for (i=1;i<7;i++) idat[i]=igs_h[i];
}
double lagrange(dat,I_la,deltat,rit,itb) 
long deltat,itb,I_la;
double dat[I_LA],rit;
/* 函数 lagrange 应用拉格朗日函数 I_la 按时间次序内查 */
{   long i,j,it0;
    double sum,f1,f2,rk;
    it0 = itb;
    if (it0 < 0) it0 = 0;
    for (i=0;i<=I_la;i++) if(rit==it0+i*deltat) return(dat[i]);
    sum = 0.0;
    rk = rit-it0;
for (i=0;i<=I_la;i++) 
{ f1 = 1.0;
      for (j=0;j<=I_la;j++)
      if(i!=j) { f2 = (rk-j*deltat)/((double)(deltat*(i-j)));
      f1 *= f2;  }
    sum += f1*dat[i]; }
    return(sum); }
 void r_igs_h(fp,fp1,igs_h) 
/*函数 r_igs_h 读取标准 igs 文件字头数据 */  
FILE *fp,*fp1;
long igs_h[8];
{ long i,j,k;
  char c,c1,c2,s[Iline];
  i = getline(fp,s); /* 第一行 */   
  sscanf(s,"%1c%1c%1c%d%d%d%d%d",&c,&c1,&c2,
  &igs_h[1],&igs_h[2],&igs_h[3],&igs_h[4],&igs_h[5]);
  igs_h[6] = 0;
  if(Igs_v==2) 
  { igs_h[3]+=1; igs_h[4]=0; igs_h[5]=0; }
  i = getline(fp,s); /* 第2--3 行 */
  i = getline(fp,s); 
  sscanf(s,"%1c%1c%d",&c,&c1,&igs_h[0]);
  for (k=0;k<19;k++) /* 下面的行数 */
  { i = getline(fp,s); }
/*  对新版本文件重读数据的处理 */
        if(Igs_v==2)
        {  j = (igs_h[0]*2+1)*6;
           for (k=0;k<j;k++) 
            { i = getline(fp,s); }}}
void r_igs_d(fp,fp1,igs_sat_num,igs_idat,igs_ddat,id) 
/* 函数 r_igs_d 读取标准 igs 数据文件 */  
FILE *fp,*fp1;
long igs_sat_num,igs_idat[8];
double igs_ddat[Isat*8];
int *id;
{ long i,j,k,l;
  char c,c1,s[Iline]; double f;
/* 读取 igs 文件的数据部分 */
  repeat:; /* 第一行 */
  i=getline(fp,s);  
  if(s[0]!='*') 
  { if(s[0]==' '&&s[1]=='*') {;}
    else{ goto repeat;} }
   sscanf(s,"%1c%1c%d %d %d %d %d",&c,&c1,&igs_idat[0],&igs_idat[1],
            &igs_idat[2],&igs_idat[3],&igs_idat[4]); 
   igs_idat[5] = 0;
/* 对每个卫星的一行 */
     if(Igs_v<2) /* 对老版本 igs 数据格式 */
     { for (k=0;k<igs_sat_num;k++) 
         { l=getline(fp,s); if(k==0) *id = l;
           if(l<70) 
           {  sscanf(s,"%1c%d%lf%lf%lf%lf",&c,&j,&igs_ddat[1+k*8],
                  &igs_ddat[2+k*8],&igs_ddat[3+k*8],&igs_ddat[4+k*8]);
              igs_ddat[0+k*8]=(double)j;}
            else 
           {  sscanf(s,"%1c%d%lf%lf%lf%lf%lf%lf", &c,&j,&igs_ddat[1+k*8],
                      &igs_ddat[2+k*8],&igs_ddat[3+k*8],
                      &igs_ddat[4+k*8],&igs_ddat[5+k*8],&igs_ddat[6+k*8]);
              igs_ddat[0+k*8]=(double)j;}}}
             else			/* 对新版本 igs 数据格式 */
            { for (k=0;k<igs_sat_num;k++) 
               { l=getline(fp,s); if(k==0) *id = l;
                 sscanf(s,"%1c%d%lf%lf%lf%lf",&c,&j,&igs_ddat[1+k*8],
                        &igs_ddat[2+k*8],&igs_ddat[3+k*8],&igs_ddat[4+k*8]);
                 igs_ddat[0+k*8]=(double)j;
                 l=getline(fp,s); if(k==0) *id = l;
                 sscanf(s,"%1c%d%lf%lf%lf%lf",&c,&j,&igs_ddat[5+k*8],
                        &igs_ddat[6+k*8],&igs_ddat[7+k*8],&f);}}
}