orbit.c

int refsat(itborb,borb,ista,x0,rtb,rte)
/* 函數refsat 應用軌道數據 itborb, borb 和測站坐標為rtb-rte觀測時間段確定參考衛星*/
long itborb[8],ista;
double borb[Isat*I96*8],rtb,rte,x0[3*Ista];
{    int i,i1,i2,j,j1,j2,k,k1,k2;
     long deltat=900,ite,itb,sat_num,it,It0,idz[I96*Isat];
     double xj,yj,zj,xi,yi,zi,z,sum,zmark,zij[Isat*Ista*I96],zmean[I96*Isat];
/* 輸出軌道數據 */
     sat_num = itborb[0];
     It0 = itborb[4]*3600+itborb[5]*60+itborb[6];
     itb = (long)((rtb-It0)/deltat);
     if(rtb-(itb*deltat+It0)>=600.) itb += 1;
     ite = (long)((rte-It0)/deltat);
     if(rte-(ite*deltat+It0)>=300.) ite += 1;
     if(ite>=I96) ite=I96;
/* 輸出每個測站的坐標，計算每個歷元衛星到測站的天頂距，存儲并輸出 */
     for (i=itb;i<=ite;i++) 
      { it = It0 + i*deltat;
        for (i1=0;i1<sat_num;i1++) 
         { xj = borb[1+i1*8+i*8*Isat]; 
           yj = borb[2+i1*8+i*8*Isat]; 
           zj = borb[3+i1*8+i*8*Isat];
           for (i2=0;i2<ista;i2++) 
            { xi = x0[0+i2*3]; yi = x0[1+i2*3]; zi = x0[2+i2*3]; 
              z = zenith1(xj,yj,zj,xi,yi,zi);
              zij[i2+i1*Ista+i*Ista*Isat] = z*90/asin(1.0);}}}
/* 計算全部衛星的天頂距，設置 向量jn[]-1為有效，天頂距設為90-15=75度 */
     zmark = 75.0; 
     for (i=0;i<Isat;i++) {zmean[i] = 0.0; idz[i]=0;}
 for (i=itb;i<=ite;i++) 
 { for (i1=0;i1<sat_num;i1++) 
   { for (i2=0;i2<ista;i2++) 
      {   sum = zij[i2+i1*Ista+i*Isat*Ista];
              zmean[i1] += sum/(ista*(ite-itb+1)); 
              if(sum<=zmark) idz[i1] += 1;}}}
/* 輸出平均天頂距和有效向量，決定參考衛星 */
     for(i=0,i2=0;i2<sat_num;i2++) if(i<idz[i2]) i=idz[i2];
     sum=90.;
     for(i1=0;i1<sat_num;i1++) 
     { if(i==idz[i1])
       {if(zmean[i1]<sum) {sum=zmean[i1]; j=i1;}}}
     i = (long)borb[j*8];
     return(i); }
int getorb(itborb,borb,isat_id,rt,I_la,orb)
/* 函數getorb 應用輸入 isat_id 和時間 rt計算有關的衛星軌道 orb[] */
long itborb[8],isat_id,I_la;
double borb[Isat*I96*8],rt,orb[8];
{    int k,j,i,i1,i2,i3=1;
     long deltat=900,itb;
     double dat[I_LA];
     i=isat_id;				 /* 由 sat_id 求k */
     for(k=-1,j=0;j<itborb[0];j++) 
     {  if(i==(int)borb[0+j*8+0*8*Isat]) 
       { k=j; goto finditb;}}
     finditb:; 				/* 由rt 求 itb  */
     if(k!=-1) 
     {  j = (itborb[4]*3600 + itborb[5]*60)/900;
        i2 = ((int)rt)/900 - I_la/2 -j;
        if(i2>itborb[7]-1-I_la) i2=itborb[7]-1-I_la;
        if(i2<0) i2 = 0;
        itb = itborb[4]*3600 + i2*deltat;
        for(j=1;j<=7;j++) 
         { for(i1=i2;i1<=i2+I_la;i1++) dat[i1-i2]=borb[j+k*8+i1*8*Isat];
           orb[j-1]=lagrange(dat,I_la,deltat,rt,itb);}}
      else 
      { i3 = 0;}  
     return(i3);
}