{  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);
}