rel_pos.c

国内外搜集的GPS程序代码

ssgsoft.c

/* Main program of kinematic/static GPS data processing software */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "xu.h"

main(argc,argv)
int argc;
char *argv[];
{
     long ifp,In0,I,Itend,iter=0,Iequa_n,Iequa_n1,Iequa_n2,Iequa_n3;
     char conf_name[I50],str[Iline],docuname[2][I50];
     char docux[I50];
     FILE *fp_rinex[Ista],*fp_orbit[Ista],*fp_docu[6],*fp,*fpd,*fpd1;
     long xhidat[2*Isat*Ista]; /* r_rinexh */
     double rt,dt; /* getdat, odat_tc */
     double odat[Isato*10*Ista*Iepoches],odat1[Isato*10*Ista*Iepoches];
     double odat_t[Iepoches],odat_dt[Ista*Iepoches];
     long odat_n[Iepoches],id_odat[I10*Ista*Iepoches];
     unsigned int iodat[Isato*10*Ista*Iepoches],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],ion_r[Isato*10*Ista*Iepoches];/*cycle slips_test*/
     long id_slips[Islips];
     unsigned int iion_r[Isato*10*Ista*Iepoches];
     int hour,minute; /* t_hms */
     double sec;
     int icon[I10*Isat],ici[Iunknown],icon20,icon201,icon231; /* r_con */
     double x0[3*Ista];  
     long itborb[8],isat_id,I_la,deltat=900,itb,ista; /* orbit */
     double borb[Isat*I96*8],rtb,rte,dt_interval;
     double Tr1[I3*Ista*I96]; /* tide */
     long Fi0[Iunknown],N0[Iunknown],idN0[Iunknown],N0t[Iunknown],
          idlo[Iunknown],idN00[Iunknown],N0tb[Iunknown];/*eq*/
     long N0te1[Iunknown],N0te2[Iunknown];/*add to eq*/
     double am[Isato*Ista*Iunknown],lo[Isato*Ista];
     double an_s[Nlow],b1_s[N],an[Nlow],b1[N],c_dd[Nlow],ltpl_s;/*norm.c*/
     double N00[Iunknown],m00,ltpl,accur[N]
            ,accur_[N],N00ac[Iunknown],dx01ac,b1_[N];
     double Nfixed1[N],Nfixed2[N],Nfixed3[N];

     double an_k[Nlow],b1_k[N],an_s1[Nlow],b1_s1[N],ltpl_k,ltpl_s1;
     double ltpl1,lctplc1,lptplp1,lctplc2,lptplp2;
     long i,i1,i2,i3,i4,i5,j,j1,j2,j3,j4,j5,k,k1,k2,k3,k4,k5;/* local */
     double phi,lambda,height,xi,yi,zi,xj,yj,zj;
     double pi[3],pj[6],m0,m0_,t00,t00_,t00_1;
     double Fi[Ista],La[Ista],Height[Ista],x01[3*Ista];
     double T=18.,P=1013.,rh=0.5,zdis,trop,Hour,phi0,lambda0;
     double z01[3*Ista];/*phase*/
     double dx01[3*Ista];/*phase*/
     long year,month,day; 
     double min,max,ave,var,min1,max1,ave1,var1;
     double mi1[Isat],ma1[Isat],av1[Isat],va1[Isat];
     double Dx,Dy,Dz,Cfi,Cla,Sfi,Sla,Dfi,Dla,A1,A2,A3,B1,B2,B3; char c2;
     double disx,wn,dttt;
     double singlex[3*Ista*Iepoches],singlet[Ista*Iepoches];
     double singlx[3*Ista],singlt[Ista];

/* 输入文件名 */
     printf("Command Line should be: command <conf_name> \n");
     printf("used command components : %d\n",argc);
     Run_n4:;
     for(i=0;i<argc;i++) { printf("%s\n",argv[i]);
     if(i==1) strcpy(conf_name,argv[i]); 
     }
/* 读控制文件和打开文件 */
     for(i=0;i<Ista;i++) TropCoef[i]=1.0;

     r_con(fp_rinex,fp_orbit,fp_docu,conf_name,icon,argc,x0,
           idN00,N00,N0tb,&m00,ici,docux);

     ifp=icon[0]; 
/* 检查定义 */
/* 运行run_n>=3， 第二次运行， 读模糊度向量 */
     if(icon[12]!=1) { /* 动态情况 */
       ;
     } else { /* 精态情况 */
        /* 如果 Run_n==4，读数据 */
        if(Run_n==4){
           if((fp_docu[3] = fopen("nfiles","r")) == NULL)
           { printf("Cannot open the file : nfiles\n");
           exit(0); }
           i=0;
           onelinemores:;   
           k1 = getline(fp_docu[3],str); if(k1<20) goto nomores;
           sscanf(str,"%ld %ld %ld %lf %ld %ld %ld",
           &N0xid[i],&N0xtb[i],&N0xte[i],&N00f[i],&N00n[i],&N0x[i],&Fi0x[i]); 
           printf("N0xidtbtevnN0Fi0= %8d %8d %8d %10.4f %6d %8d %9d\n",
           N0xid[i],N0xtb[i],N0xte[i],N00f[i],N00n[i],N0x[i],Fi0x[i]);
           i+=1; if(k1>=20) goto onelinemores;
           nomores:; Nfix=i; printf("Nfix=%d\n",Nfix);
           fclose(fp_docu[3]);
        }
     }
/* 读GPS数据字头 */
     r_rinexh(fp_rinex,ifp,xhidat); Ob_num = (int)xhidat[0];
     for(i1=0;i1<ifp;i1++) {
        for (i=0;i<=34;i++) fprintf(stdout,"%ld ",xhidat[i+i1*2*Isat]);
        fprintf(stdout,"\n"); 
     }
     if(xhidat[14]==0 && xhidat[15]==0) {
        printf("input data interval\n");
        scanf("%lf",&dt); printf("input: %lf\n",dt);
        xhidat[14]=(int)dt; xhidat[15]=(int)(10*(dt-xhidat[14]));
     }
     id_obs(xhidat,&C1,&C2,&L1,&L2,&P1,&P2,&D1,&D2);
     /* 确定 PC1 和 PC2 */
     PC12(&PC1,&PC2,&PD1,&PD2);
/* 确定所采用的坐标 */
     if(icon[3]==2) {
        printf("coordinates in RINEX files used:\n");
        fprintf(fp_docu[0],"coordinates in RINEX files used:\n");
        for(i=0;i<icon[0];i++) {
           x0[0+i*3]=xhidat[23+i*2*Isat]+xhidat[24+i*2*Isat]/1000.;
           x0[1+i*3]=xhidat[25+i*2*Isat]+xhidat[26+i*2*Isat]/1000.;
           x0[2+i*3]=xhidat[27+i*2*Isat]+xhidat[28+i*2*Isat]/1000.;
           printf("%14.3f %14.3f %14.3f\n",x0[0+i*3],x0[1+i*3],x0[2+i*3]);
           fprintf(fp_docu[0],"%14.3f %14.3f %14.3f\n",
                   x0[0+i*3],x0[1+i*3],x0[2+i*3]);
        }
     }
/* 迭代计算初始化 */
     for(i=0;i<icon[0];i++) { /* 为了迭代 */
        x01[0+i*3]=x0[0+i*3];x01[1+i*3]=x0[1+i*3]; x01[2+i*3]=x0[2+i*3];
     }
/* 计算每个测站的相对高程 */
     fprintf(fp_docu[0],"coordinates in WGS84 geodetic ellipsoid:\n");
     for(i=0;i<icon[0];i++) {
        xi=x0[0+i*3]; yi=x0[1+i*3]; zi=x0[2+i*3];
        if(O_main==0) printf("xyz=%lf %lf %lf\n",xi,yi,zi);
        TRANF(2,&phi,&lambda,&height,&xi,&yi,&zi,0.,0.);
        Fi[i]=phi; La[i]=lambda; Height[i]=height;
        fprintf(fp_docu[0],"%14.8f %14.8f %14.3f\n",
                Fi[i]*180/PI,La[i]*180/PI,Height[i]);
     }
     for(i=0;i<Kstation;i++) Height[i]=0.;/*set kinematic Height as zero*/
     if(icon[12]==2) lambda0=La[0];
/* 时间 */
     rtb = icon[4]*3600 + icon[5]*60. + icon[6]; 
     rte = icon[7]*3600 + icon[8]*60. + icon[9]; 
      dt = xhidat[14]+xhidat[15]/10.; /*  RINEX 数据中的dt  */
      dt = icon[10]*1.0/icon[11]; /* 想要的 dt */
     rt=rtb; dt_interval=dt; Itend = (long)((rte-rtb)/dt)+1; 
     printf("rtb,rte,dt,Itend= %7.1f %7.1f %4.1f %7d\n",rtb,rte,dt,Itend);
/* 得到轨道数据并决定参考卫星 Ref_sv */
     I_la = icon[25]; 
     if(icon[1]==1) igs(fp_orbit[0],borb,itborb,I_la);
     else broadcast(fp_orbit,icon[2],borb,itborb);
     Ref_sv=icon[14];
     if(icon[14]<=0 || icon[14]>31) 
     Ref_sv=refsat(itborb,borb,icon[0],x0,rtb,rte);
     printf("refsat= %d\n",Ref_sv);
     for(i=0;i<8;i++) printf("%d ",itborb[i]); printf("\n");
     i1=icon[0];
     for(i=0;i<I96;i++) for(j=0;j<3*Kstation;j++) Tr1[j+i*3*i1] = 0.;
/* 为得到GPS数据的初始化 */
     In0=-1;
/* 为计算矩阵 an_s[]，向量 b1_s[]、 ltpl_s 的初始化 */
     for (i=0;i<Iunknown;i++) { 
        for(j=0;j<=i;j++) {an_s[j+i*(i+1)/2] = 0.; }
        b1_s[i] = 0.; accur[i]=0.; N00_n[i]=0; Nk[i]=0.;
     } ltpl_s = 0.; icon[202]=0; icon[207]=0; icon20=icon[20];
/* 每个历元数据进行循环 */ 
     for(I=1;I<Itend;I++) { /* 历元间数据循环，求得GPS数据，数据标号为In0 */
        In0+=1; 
        icon[206]=0; if(In0==Iswitch) icon[206]=1;
        gps_dat(fp_rinex,ifp,xhidat,odat_t,odat_dt,odat_n,&In0,rt,dt,odat,
                id_odat,iodat,odat1,iodat1,
                sdat,id_sdat,isdat,sdat_n,ddat,id_ddat,iddat,ddat_n,
                tdat,id_tdat,itdat,tdat_n,Ref_sv,
                slips_t,id_slips,ion_r,iion_r
                ,icon,x0,Height,borb,itborb,Tr1,singlex,singlet,singlx,singlt
        );
        t_hms(odat_t[In0],&hour,&minute,&sec);
        if(I==1) t00_1=odat_t[In0-1]; else t00_1=t00_;