eqdd_s.c

GPS应用程序设计.程序源代码.GPS应用程序设计.程序源代码

eqdd_s.c

void eq_dd_s(icon,x0,Height,xhidat,In0,odat_t,odat_n,id_odat,iodat1,odat1,
             sdat,id_sdat,isdat,sdat_n,ddat,id_ddat,iddat,ddat_n,
             tdat,id_tdat,itdat,tdat_n,
             slips_t,id_slips,
borb,itborb,Tr1,
Fi0,N0,idN0,N0t,N0te,am,lo,idlo) 

int icon[I10*Isat];
double x0[3*Ista],Height[Ista];
long xhidat[2*Isat*Ista],*In0;
double odat1[Isato*10*Ista*Iepoches],odat_t[Iepoches];
long odat_n[Iepoches],id_odat[I10*Ista*Iepoches];
unsigned int iodat1[Isato*10*Ista*Iepoches];
double sdat[Isato*10*Ista*Iepoches],ddat[Isato*10*Ista*Iepoches],
       tdat[Isato*10*Ista*Iepoches]; /* sdtdat */
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*10*Ista*Iepoches],iddat[Isato*10*Ista*Iepoches],
             itdat[Isato*10*Ista*Iepoches];
double slips_t[Islips]; /* cycle_slips_test */
long id_slips[Islips];
long itborb[8]; /* orbit */
double borb[Isat*I96*8];
double Tr1[I3*Ista*I96]; /* tide */

long Fi0[Iunknown],N0[Iunknown],idN0[Iunknown],N0t[Iunknown],idlo[Iunknown];
long N0te[Iunknown];
double am[Isato*Ista*Iunknown],lo[Isato*Ista];
/* Function eq_dd_s using double difference data dd at epoch In0 
   for all stations and all satellites to form the observations
   equation system in static case
*/
{
     int i,i1,i2,j,j1,j2,k,k1,k2,k3,k4,k5,I1,I2,J1,J2,K=0,PC,hour,minute;
     double T=20.,P=1013.,rh=0.5,fc,Lambdafc,sum,sumc,sec,time;
     double zdis,z1,z2,z3,z4,z11,z12,z13,z14,f1,f2,f3,f4,height;
     double xj1i1,yj1i1,zj1i1,xj1i2,yj1i2,zj1i2,xj2i1,yj2i1,zj2i1,
            xj2i2,yj2i2,zj2i2;
     double p0j1i2,p0j2i2,p0j1i1,p0j2i1,dtj1i1,dtj1i2,dtj2i1,dtj2i2;
     double xi1,yi1,zi1,xi2,yi2,zi2,xj1,yj1,zj1,xj2,yj2,zj2;
     double dti1,dti2,dti3,dti4;
     double xj1dot,yj1dot,zj1dot,xj2dot,yj2dot,zj2dot,dcj1,dcj2;
     double orbj1[8],orbj2[8],pj[6],pi[3],Tr[3*Ista],dd_rela,dd_trop,dd_cloc;
     double txi1,tyi1,tzi1,txi2,tyi2,tzi2,dt,ff;
double xj1doti1,yj1doti1,zj1doti1,xj1doti2,yj1doti2,zj1doti2;
double xj2doti1,yj2doti1,zj2doti1,xj2doti2,yj2doti2,zj2doti2;
/* Set frequency for combined phase obs. */ 
     if (icon[20]==2) fc = FL1 - FL2*(FL2/FL1);
     if (icon[20]==3) fc = 0.77*FL1 - 0.6*FL2;
     if (icon[20]==4) fc = FL1 - FL2;
     if (icon[20]==5) fc = icon[21]*FL1/icon[22]+icon[23]*FL2/icon[24];
     if (icon[20]==0) fc = FL1;
     if (icon[20]==1) fc = FL2;
     if (icon[20]>=6) fc = FL1;
     Lambdafc = C / fc;
/* For epoch In0, every observation ddat */
     k=*In0;
     t_hms(odat_t[k],&hour,&minute,&sec); dt=odat_t[k+1]-odat_t[k];
     if(ddat_n[k]==0) icon[200]=0;
     for (j1=0;j1<ddat_n[k];j1++) {
      /* set am = 0. */
        for (k1=0;k1<Isato*Ista;k1++) am[k1+j1*Iunknown] = 0.;
      /* find out the id_ddat[] idex I1,I2,J1,J2, set idlo[] */
        id_i12j12(&id_ddat[j1+I10*Ista*k],&I1,&I2,&J1,&J2);
        idlo[K]=id_ddat[j1+I10*Ista*k];
      /* decide if to cancel this obs. */
        if(Run_n>=2) {
           for(k5=0;k5<Ncancel_n;k5++) {
              if(idlo[K]==N00id[k5]) {
              if(N00tb[k5]/10.-0.2<=odat_t[k] && odat_t[k]<=N00te[k5]/10.-0.2)
              goto nextj1;
        }}}
      /* find out coord. for sta. I1 and I2, sat. J1 and J2,
         get clock corrections of sat. J1 and J2 (unit s/1000000) */
        xi1=x0[0+(I1-1)*3]; yi1=x0[1+(I1-1)*3]; zi1=x0[2+(I1-1)*3];
        xi2=x0[0+(I2-1)*3]; yi2=x0[1+(I2-1)*3]; zi2=x0[2+(I2-1)*3];
        k5=getorb(itborb,borb,J1,odat_t[*In0],icon[25],orbj1);/*sat. J1*/
        if(k5==0) {printf("no sat_%d data available\n",J1); exit(0);}
        xj1=orbj1[0]; yj1=orbj1[1]; zj1=orbj1[2]; 
        xj1dot=orbj1[3]; yj1dot=orbj1[4]; zj1dot=orbj1[5]; dcj1=orbj1[6]; 
        k5=getorb(itborb,borb,J2,odat_t[*In0],icon[25],orbj2);/*sat. J2*/
        if(k5==0) {printf("no sat_%d data available\n",J2); exit(0);}
        xj2=orbj2[0]; yj2=orbj2[1]; zj2=orbj2[2]; 
        xj2dot=orbj2[3]; yj2dot=orbj2[4]; zj2dot=orbj2[5]; dcj2=orbj2[6]; 
      /* tropospheric effects f1,f2,f3,f4, dd_trop,
        z1,z2,z3,z4 are zenith distances between j2i2, j2i1, j1i2, j1i1.*/
        /* canceled here */
        f1=0.;f2=0.;f3=0.;f4=0.;
        /* because : odat[], i1_ref.sta. j1_ref.sat. 
        sdat(i1,i2,j) = odat(i2,j) - odat(i1,j)
        ddat(i1,i2,j1,j2) = sdat(i1,i2,j2) - sdat(i1,i2,j1)
        therefore : dd = od(i2,j2) - od(i1,j2) - od(i2,j1) + od(i1,j1) */ 
        dd_trop = (f1 - f2 - f3 + f4); /* in meter */
      /* relative correction : z1,z2,z3,z4 (canceled here)*/
        z1=0.; z2=0.; z3=0.; z4=0.;
        dd_rela = z1-z2-z3+z4;
      /* sat. clock correction dcj1,dcj2 : dd_cloc here == 0 */
        dd_cloc = 0.;
      /* tide correction (canceled): one epoch, every stations */
      /* time transmitting correction xyzj2i2,xyzj2i1,xyzj1i2,xyzj1i1 
        (rotation correction already considered in orbit data) */
        xj1i1=0.;yj1i1=0.;zj1i1=0.;xj1i2=0.;yj1i2=0.;zj1i2=0.;
        xj2i1=0.;yj2i1=0.;zj2i1=0.;xj2i2=0.;yj2i2=0.;zj2i2=0.;
        z11=0.; z12=0.; z13=0.; z14=0.;
        k5 = 0;
        transtime:;
        /* Distance between sat. and sta. p0j2i2,p0j2i1,p0j1i2,p0j1i1 */
if(Ixyz>=2) { /* transmitting and rotation */
        xj1i1=xj1;yj1i1=yj1;zj1i1=zj1;xj1i2=xj1;yj1i2=yj1;zj1i2=zj1;
        xj2i1=xj2;yj2i1=yj2;zj2i1=zj2;xj2i2=xj2;yj2i2=yj2;zj2i2=zj2;
        xj1doti1=xj1dot;yj1doti1=yj1dot;zj1doti1=zj1dot;
        xj1doti2=xj1dot;yj1doti2=yj1dot;zj1doti2=zj1dot;
        xj2doti1=xj2dot;yj2doti1=yj2dot;zj2doti1=zj2dot;
        xj2doti2=xj2dot;yj2doti2=yj2dot;zj2doti2=zj2dot;

        transrot1(&xj1i1,&yj1i1,&zj1i1,&xj1doti1,&yj1doti1,&zj1doti1,
                 xi1,yi1,zi1,&p0j1i1,&dtj1i1);
        xj1i1-=xj1;yj1i1-=yj1;zj1i1-=zj1;
        transrot1(&xj1i2,&yj1i2,&zj1i2,&xj1doti2,&yj1doti2,&zj1doti2,
                 xi2,yi2,zi2,&p0j1i2,&dtj1i2);
        xj1i2-=xj1;yj1i2-=yj1;zj1i2-=zj1;

        transrot1(&xj2i1,&yj2i1,&zj2i1,&xj2doti1,&yj2doti1,&zj2doti1,
                 xi1,yi1,zi1,&p0j2i1,&dtj2i1);
        xj2i1-=xj2;yj2i1-=yj2;zj2i1-=zj2;
        transrot1(&xj2i2,&yj2i2,&zj2i2,&xj2doti2,&yj2doti2,&zj2doti2,
                 xi2,yi2,zi2,&p0j2i2,&dtj2i2);
        xj2i2-=xj2;yj2i2-=yj2;zj2i2-=zj2;
}
        /* check the transmitting time */
        /* correction of receiver clock error related to satellite movement */
        /* calculate receiver clock error tci1,tci2 */
        /* ref. sat. clock error correction is dcj1 */
        /* find out code measurements related to ref. sat. of both sta. */
        if(Code==0) {if(P1<100) PC=P1; else {
                        if(P2<100) PC=P2; else {
                           if(C1<100) PC=C1; else PC=C2;
        }}} else {
           if(Code==1) PC=P1; if(Code==2) PC=P2;
           if(Code==3) PC=C1; if(Code==4) PC=C2;
        }
        ff = (FL1*FL1-FL2*FL2);
        for(k5=0;k5<odat_n[k];k5++) {
           if(id_odat[k5+I10*Ista*k]==I1*100+J1) {
              z11=odat1[PC+k5*Isato+Isato*I10*Ista*k];
              z11=z11*(FL1*FL1)/ff;
              z11=z11-odat1[PC2+k5*Isato+Isato*I10*Ista*k]*(FL2*FL2)/ff;
           }
           if(id_odat[k5+I10*Ista*k]==I2*100+J1) {           
              z12=odat1[PC+k5*Isato+Isato*I10*Ista*k];
              z12=z12*(FL1*FL1)/ff;
              z12=z12-odat1[PC2+k5*Isato+Isato*I10*Ista*k]*(FL2*FL2)/ff;
           }
           if(id_odat[k5+I10*Ista*k]==I1*100+J2) {
              z13=odat1[PC+k5*Isato+Isato*I10*Ista*k];
              z13=z13*(FL1*FL1)/ff;
              z13=z13-odat1[PC2+k5*Isato+Isato*I10*Ista*k]*(FL2*FL2)/ff;
           }
           if(id_odat[k5+I10*Ista*k]==I2*100+J2) {           
              z14=odat1[PC+k5*Isato+Isato*I10*Ista*k];
              z14=z14*(FL1*FL1)/ff;
              z14=z14-odat1[PC2+k5*Isato+Isato*I10*Ista*k]*(FL2*FL2)/ff;
           }
        }
        dti1 = (z11-p0j1i1)/C ; 
        dti2 = (z12-p0j1i2)/C ; 
        dti3 = (z13-p0j2i1)/C ; 
        dti4 = (z14-p0j2i2)/C ; 
        dti1=-dti1; dti2=-dti2; 
        dti3=-dti3; dti4=-dti4;
xj2doti2=xj2doti2*dti4;yj2doti2=yj2doti2*dti4;zj2doti2=zj2doti2*dti4;
if(Ixyz>=3) earthrot(dtj2i2,&xj2doti2,&yj2doti2,&zj2doti2);
xj2i2 += xj2doti2; yj2i2 += yj2doti2; zj2i2 += zj2doti2;
        txi1 = xj2 + xj2i2 - xi2 ;
        tyi1 = yj2 + yj2i2 - yi2 ;
        tzi1 = zj2 + zj2i2 - zi2 ;
        p0j2i2 = sqrt(txi1*txi1 + tyi1*tyi1 + tzi1*tzi1);
xj2doti1=xj2doti1*dti3;yj2doti1=yj2doti1*dti3;zj2doti1=zj2doti1*dti3;
if(Ixyz>=3) earthrot(dtj2i1,&xj2doti1,&yj2doti1,&zj2doti1);
xj2i1 += xj2doti1; yj2i1 += yj2doti1; zj2i1 += zj2doti1;
        txi1 = xj2 + xj2i1 - xi1 ;
        tyi1 = yj2 + yj2i1 - yi1 ;
        tzi1 = zj2 + zj2i1 - zi1 ;
        p0j2i1 = sqrt(txi1*txi1 + tyi1*tyi1 + tzi1*tzi1);
xj1doti2=xj1doti2*dti2;yj1doti2=yj1doti2*dti2;zj1doti2=zj1doti2*dti2;
if(Ixyz>=3) earthrot(dtj1i2,&xj1doti2,&yj1doti2,&zj1doti2);
xj1i2 += xj1doti2; yj1i2 += yj1doti2; zj1i2 += zj1doti2;
        txi1 = xj1 + xj1i2 - xi2 ;
        tyi1 = yj1 + yj1i2 - yi2 ;
        tzi1 = zj1 + zj1i2 - zi2 ;
        p0j1i2 = sqrt(txi1*txi1 + tyi1*tyi1 + tzi1*tzi1);
xj1doti1=xj1doti1*dti1;yj1doti1=yj1doti1*dti1;zj1doti1=zj1doti1*dti1;
if(Ixyz>=3) earthrot(dtj1i1,&xj1doti1,&yj1doti1,&zj1doti1);
xj1i1 += xj1doti1; yj1i1 += yj1doti1; zj1i1 += zj1doti1;