RTKLIB源码学习之后处理定位postpos.c

本博客为转载学习之用，感谢：http://t.csdnimg.cn/JyzsD

postpos函数为后处理定位的主入口函数，根据tu分计算时间段，调用execses_b()进行下一步解算，输入文件包括观测文件、导航文件、精密星历文件等，postpos在处理输入文件时有两种方法：一种是输入文件可以只包含关键词，然后通过函数处理，将关键词用时间、基准站编号、流动站编号等代替；另一种是直接调用输入文件的文件名，postpos主要是用来判断是哪一种输入方式，然后调用相应函数

postpos流程图

1. 主程序

/* post-processing positioning -------------------------------------------------
* post-processing positioning
*args    ：gtime_t ts       I   处理的起始时间，写0表示不限制
*        ：gtime_t te       I   处理的起始时间，写0表示不限制
*          double ti        I   处理的间隔时间 (s)，写0表示不限制，全处理
*          double tu        I   处理的单元时间（s)，写0表示全部做一个单元处理
*          prcopt_t *popt   I   处理选项结构体
*          solopt_t *sopt   I   结果选项结构体
*          filopt_t *fopt   I   文件选项结构体
*          char   **infile  I   传入文件路径数组首地址
*          int    n         I   传入文件数量
*          char   *outfile  I   输出文件的路径，写0表示stdout终端
*          char   *rov      I   流动站ID列表，空格隔开
*          char   *base     I   基准站ID列表，空格隔开
* return : status (0:ok,0>:error,1:aborted)
* notes  : input files should contain observation data, navigation data, precise 
*          ephemeris/clock (optional), sbas log file (optional), ssr message
*          log file (optional) and tec grid file (optional). only the first 
*          observation data file in the input files is recognized as the rover
*          data.
*
*          the type of an input file is recognized by the file extention as ]
*          follows:
*              .sp3,.SP3,.eph*,.EPH*: precise ephemeris (sp3c)
*              .sbs,.SBS,.ems,.EMS  : sbas message log files (rtklib or ems)
*              .lex,.LEX            : qzss lex message log files
*              .rtcm3,.RTCM3        : ssr message log files (rtcm3)
*              .*i,.*I              : tec grid files (ionex)
*              others               : rinex obs, nav, gnav, hnav, qnav or clock
*
*          inputs files can include wild-cards (*). if an file includes
*          wild-cards, the wild-card expanded multiple files are used.
*
*          inputs files can include keywords. if an file includes keywords,
*          the keywords are replaced by date, time, rover id and base station
*          id and multiple session analyses run. refer reppath() for the
*          keywords.
*
*          the output file can also include keywords. if the output file does
*          not include keywords. the results of all multiple session analyses
*          are output to a single output file.
*
*          ssr corrections are valid only for forward estimation.
*-----------------------------------------------------------------------------*/
extern int postpos(gtime_t ts, gtime_t te, double ti, double tu,
                   const prcopt_t *popt, const solopt_t *sopt,
                   const filopt_t *fopt, char **infile, int n, char *outfile,
                   const char *rov, const char *base)

 1.2 处理流程

 变量定义，stat默认为0，flag默认为1

    gtime_t tts,    //解算单元的开始时间
            tte,    //解算单元的结束时间
            ttte;   //读取星历文件的结束时间
    double tunit,   //
           tss;     //
    int i,j,k,      //循环和数组下标控制
            nf,     //文件路径数组下标控制
        stat=0,     //接收返回状态值，为1
        week,       //用于存GPST的周
        flag=1,
        index[MAXINFILE]={0};
    char *ifile[MAXINFILE],
          ofile[1024],
          *ext;

 调用openses函数读取天线，大地水准面文件

 /* open processing session */   //开始处理,文件读取，赋值navs、pcvs、pcvsr
    if (!openses(popt,sopt,fopt,&navs,&pcvss,&pcvsr)) return -1;

判断起始解算时间ts、结束结算时间te、结算时间单元tu的关系是否正确，一共有三种情况：

    if (ts.time!=0&&te.time!=0&&tu>=0.0) {  //判断起始时间ts、te、处理单位时间是否大于0
        if (timediff(te,ts)<0.0) {  //结束时间早于开始时间
            showmsg("error : no period");
            closeses(&navs,&pcvss,&pcvsr);  //不合理则关闭处理，释放navs、pcvs、pcvsr
            return 0;
        }
         for (i=0;i<MAXINFILE;i++) {
            if (!(ifile[i]=(char *)malloc(1024))) { //为infile数组malloc开辟空间
                for (;i>=0;i--) free(ifile[i]);     //开辟失败则释放已开辟的空间，关闭处理释放navs、pcvs、pcvsr
                closeses(&navs,&pcvss,&pcvsr);
                return -1;
            }
        }
        if (tu==0.0||tu>86400.0*MAXPRCDAYS) tu=86400.0*MAXPRCDAYS;  //解算处理时间单元处理，0或者时间大于100天，设为100天
        settspan(ts,te);    //设置时间跨度，好像是空函数，需要自己实现
        
        tunit=tu<86400.0?tu:86400.0;    //tunit：如果tu小于一天就为tu；否则为一天
        tss=tunit*(int)floor(time2gpst(ts,&week)/tunit);   //
        
        //根据解算时间单元，分时间段循环处理，算出来tts>te或过程有错误，结束循环
        //很多时候解算单元时间直接设0.0，只循环一次，tts=ts，tte=te
        for (i=0;;i++) { /* for each periods */
            tts=gpst2time(week,tss+i*tu);       //解算单元开始时间，每次循环加上一个i个tu？
            tte=timeadd(tts,tu-DTTOL);          //解算结束时间tte=tu-DTTOL
            if (timediff(tts,te)>0.0) break;   //算出来tts>te结束循环
            if (timediff(tts,ts)<0.0) tts=ts;   //分时间段后tts若早于ts，设为ts
            if (timediff(tte,te)>0.0) tte=te;   //分时间段后tte若早于te，设为te
            
            strcpy(proc_rov ,"");   //流动站、基准站值赋空
            strcpy(proc_base,"");   
            if (checkbrk("reading    : %s",time_str(tts,0))) {
                stat=1;
                break;
            }
            for (j=k=nf=0;j<n;j++) {    //遍历infile[]，根据后缀名
                
                ext=strrchr(infile[j],'.'); //ext：文件路径中.后缀开始的地址
                
                if (ext&&(!strcmp(ext,".rtcm3")||!strcmp(ext,".RTCM3"))) {  //rtcm3文件
                    strcpy(ifile[nf++],infile[j]);
                }
                else {      //星历文件，包括精密星历和广播星历
                    /* include next day precise ephemeris or rinex brdc nav */
                    ttte=tte;
                    if (ext&&(!strcmp(ext,".sp3")||!strcmp(ext,".SP3")||
                              !strcmp(ext,".eph")||!strcmp(ext,".EPH"))) {
                        ttte=timeadd(ttte,3600.0);  //精密星历加一小时
                    }
                    else if (strstr(infile[j],"brdc")) {
                        ttte=timeadd(ttte,7200.0);  //广播星历加两小时
                    }
                    nf+=reppaths(infile[j],ifile+nf,MAXINFILE-nf,tts,ttte,"","");
                }
                while (k<nf) index[k++]=j;
                
                if (nf>=MAXINFILE) {
                    trace(2,"too many input files. trancated\n");
                    break;
                }
            }
            if (!reppath(outfile,ofile,tts,"","")&&i>0) flag=0;
            
            /* execute processing session */
            stat=execses_b(tts,tte,ti,popt,sopt,fopt,flag,ifile,index,nf,ofile,
                           rov,base);
            
            if (stat==1) break;
        }

        for (i=0;i<MAXINFILE;i++) free(ifile[i]);
    }

①：若ts、te不为0，且tu大于等于0

·若te的时间早于ts，return0

·反之，为ifile[]数组开辟空间，若开辟失败则释放已开辟的空间，释放navs、pcvs、pcvsr，return0

·开辟成功后，解算处理时间单元tu，等于0或大于100天，设为100天

·给tunit赋值：如果tu小于一天就为tu；否则为一天

·循环处理每个时间单元tts到tte：
·计算解算时间单元的开始tts、结束tte，判断若tts<ts则tts设为ts，若tte>te则tte设为te

·给流动站、基准站名赋空值
·遍历infile[]，strrchr函数找文件后缀名所在地址，strcmp函数判断后缀名，strstr函数在infile[]中查找是否含有广播星历标识 ：

        rtcm3：直接把infile[j]中路径赋值到ifile[]中
        星历文件：精密星历ttte=tte+一小时、广播星历ttte=tte+两小时，根据tts、ttte调用reppaths函数将infile[j]中路径展开到ifile[nf]中
·把infile[]的下标j存到index[]中

·调用reppath函数将outfile中的通配符进行替换，存到ofile中

·调用execses_b函数对基站进行处理

else if (ts.time!=0) {  //如果起始时间不为0，结束时间为0或处理单元时间小于0
        for (i=0;i<n&&i<MAXINFILE;i++) {
            if (!(ifile[i]=(char *)malloc(1024))) {
                for (;i>=0;i--) free(ifile[i]);
                return -1;
            }
            reppath(infile[i],ifile[i],ts,"","");
            index[i]=i;
        }
        reppath(outfile,ofile,ts,"","");
        
        /* execute processing session */
        stat=execses_b(ts,te,ti,popt,sopt,fopt,1,ifile,index,n,ofile,rov,
                       base);
        
        for (i=0;i<n&&i<MAXINFILE;i++) free(ifile[i]);
    }

②：若①不成立且ts不为0（不考虑te和tu）

·为ifile[]开辟空间，循环替换infile[i]的路径展开到ifile[i]中

·调用reppath函数将outfile中的通配符进行替换，存到ofile中

·调用execses_b函数对基站进行处理

    else {  //如果起始时间为0
        for (i=0;i<n;i++) index[i]=i;
        
        /* execute processing session */
        stat=execses_b(ts,te,ti,popt,sopt,fopt,1,infile,index,n,outfile,rov,
                       base);
    }

③：ts为0

直接把infile[]的下标j存到index[]中，调用execses_b函数对基站进行处理

 思考

         为什么要判断时间关系？ ·拆分时间段解算需要tu值有效、调用reppath需要ts有效，调用reppaths需要ts和te有效

         ifile[]、ofile[]作用：·infile[]、ofile[]里的路径替换处理后存到ifile[]、ofile[]，传入execses_b()进行之后的解算。

        index[]的作用： ·会传给execses_b()，再传给execses_r()，再传给execses()，再传给readobsnav() 。如果不需要根据tu分时间段解算，index存的是0~n，如果需要分时间段解算，index存的是对应时间段内文件的下标

调用closeses()，释放openses()开辟的内存，返回stat

 /* close processing session */
    closeses(&navs,&pcvss,&pcvsr);
    
    return stat;

 1.3 返回值

处理一切正常会接收execxes_b()的返回值，失败返回0，内存失败返回-1

execses_b()正常会接收execses_b()的返回值，失败返回0

execses_r()正常会接收execses()的返回值，失败返回0，aborts返回1

2. 调用函数

2.1 openses函数：读取导航信息

openses函数实现读取filopt_t fopt中的参数

filopt_t结构体

typedef struct {        /* file options type */
    char satantp[MAXSTRPATH]; /* satellite antenna parameters file */
    char rcvantp[MAXSTRPATH]; /* receiver antenna parameters file */
    char stapos [MAXSTRPATH]; /* station positions file */
    char geoid  [MAXSTRPATH]; /* external geoid data file */
    char iono   [MAXSTRPATH]; /* ionosphere data file */
    char dcb    [MAXSTRPATH]; /* dcb data file */
    char eop    [MAXSTRPATH]; /* eop data file */
    char blq    [MAXSTRPATH]; /* ocean tide loading blq file */
    char tempdir[MAXSTRPATH]; /* ftp/http temporaly directory */
    char geexe  [MAXSTRPATH]; /* google earth exec file */
    char solstat[MAXSTRPATH]; /* solution statistics file */
    char trace  [MAXSTRPATH]; /* debug trace file */
} filopt_t;

/* open procssing session ----------------------------------------------------*/
static int openses(const prcopt_t *popt, const solopt_t *sopt,
                   const filopt_t *fopt, nav_t *nav, pcvs_t *pcvs, pcvs_t *pcvr)

 *popt为处理方案选项类型

*sopt为解决方案选项类型

*fopt为文件处理选项类型，可以为{0}

*nav为导航数据类型

*pcvs和*pcer分别为卫星和接收机的天线参数类型

处理过程

    /* read satellite antenna parameters */
    if (*fopt->satantp&&!(readpcv(fopt->satantp,pcvs))) {
        showmsg("error : no sat ant pcv in %s",fopt->satantp);
        trace(1,"sat antenna pcv read error: %s\n",fopt->satantp);
        return 0;
    }
    /* read receiver antenna parameters */
    if (*fopt->rcvantp&&!(readpcv(fopt->rcvantp,pcvr))) {
        showmsg("error : no rec ant pcv in %s",fopt->rcvantp);
        trace(1,"rec antenna pcv read error: %s\n",fopt->rcvantp);
        return 0;
    }
    /* open geoid data */
    if (sopt->geoid>0&&*fopt->geoid) {
        if (!opengeoid(sopt->geoid,fopt->geoid)) {
            showmsg("error : no geoid data %s",fopt->geoid);
            trace(2,"no geoid data %s\n",fopt->geoid);
        }
    }
    return 1;

 openses函数读取*fopt中的卫星天线参数、dcb参数、电离层文件、接收机天线参数、大地水准面数据和地球自转参数

以读取卫星天线参数流程为例：

·readpcv函数用来读取天线参数并存储到pcvs中

·若没有卫星天线参数，showmsg函数显示报错信息，并返回0

·完成所有文件读取后返回1

2.2  execses_b函数：处理基站信息

 execses_b函数处理基站信息

/* execute processing session for each base station -------------------------
gtime_t ts              I   处理的起始时间，写0表示不限制
gtime_t te              I   处理的起始时间，写0表示不限制
double ti               I   处理的间隔时间 (s)，写0表示不限制，全处理
const prcopt_t *popt    I   处理选项结构体
const solopt_t *sopt    I   结果选项结构体
const filopt_t *fopt    I   文件选项结构体
int flag                I   用于控制输出
char **infile           I   传入文件路径数组首地址
const int *index        I   传入文件路径数组首地址
int n                   I   传入文件数量
char *outfile           I   输出文件的路径，写0表示stdout终端
const char *rov         I   流动站ID列表，空格隔开
const char *base        I   基准站ID列表，空格隔开
---------------------------------------------------------------------------*/
static int execses_b(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt,
                     const solopt_t *sopt, const filopt_t *fopt, int flag,
                     char **infile, const int *index, int n, char *outfile,
                     const char *rov, const char *base)

 调用readpreceph函数读取精密星历和SBAS数据

/* read prec ephemeris and sbas data */
    readpreceph(infile,n,popt,&navs,&sbss); //读取精密星历和SBAS数据
    

 遍历infile[]，寻找基准站替换符%b：

    //%b：基准站ID的替换符
    for (i=0;i<n;i++) if (strstr(infile[i],"%b")) break;

       · 找不到基准站ID的替换符，直接调用execses_r()进行下一步解算 

else {  //infile[i]都没有有基准站ID的替换符，直接调用execses_r()进行下一步解算
        stat=execses_r(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile,rov);
    }

        ·找到了infile[i]含有基准站ID的替换符，遍历基准站：
                1、将基准站ID赋值给proc_base

    //如果某个infile[i]含有基准站ID的替换符
    if (i<n) { /* include base station keywords */
        //为base_开辟空间，将base赋值给base_
        if (!(base_=(char *)malloc(strlen(base)+1))) {  
            freepreceph(&navs,&sbss);
            return 0;
        }
        strcpy(base_,base); 
        
        for (i=0;i<n;i++) {     //为ifile[]开辟空间
            if (!(ifile[i]=(char *)malloc(1024))) { 
                free(base_); for (;i>=0;i--) free(ifile[i]);
                freepreceph(&navs,&sbss);
                return 0;
            }
        }
        //遍历base_基准站字符串
        for (p=base_;;p=q+1) { /* for each base station */
            if ((q=strchr(p,' '))) *q='\0'; //拆出一个基准站
            
            if (*p) {   
                strcpy(proc_base,p);    //把基准站名赋值给proc_base
                if (ts.time) time2str(ts,s,0); else *s='\0';
                if (checkbrk("reading    : %s",s)) {
                    stat=1;
                    break;

                2、循环替换infile[i]里的基准站ID的替换符到ifile[i] 

                3、替换outfile里的基准站ID替换符到ofile

                4、调用execses_r函数进行下一步解算 

                //循环替换infile[i]里的基准站ID的替换符到ifile[i]
                for (i=0;i<n;i++) reppath(infile[i],ifile[i],t0,"",p);
                //替换outfile里的基准站ID替换符到ofile
                reppath(outfile,ofile,t0,"",p); 
                //调用execses_r()进行下一步解算
                stat=execses_r(ts,te,ti,popt,sopt,fopt,flag,ifile,index,n,ofile,rov);
            }
            if (stat==1||!q) break;
        }
        free(base_); for (i=0;i<n;i++) free(ifile[i]);
    }

 调用freepreceph函数，释放readpreceph函数开辟的空间，返回stat

    /* free prec ephemeris and sbas data */
    freepreceph(&navs,&sbss);
    
    return stat;

2.2.1 execses_r函数：替换流动站替换符

 execses_r函数功能实现与execses_b函数相似，故不再赘述，仅展示代码

static int execses_r(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt,
                     const solopt_t *sopt, const filopt_t *fopt, int flag,
                     char **infile, const int *index, int n, char *outfile,
                     const char *rov)
{
    gtime_t t0={0};
    int i,stat=0;
    char *ifile[MAXINFILE],ofile[1024],*rov_,*p,*q,s[64]="";
    
    trace(3,"execses_r: n=%d outfile=%s\n",n,outfile);

    //%r：流动站ID的替换符
    for (i=0;i<n;i++) if (strstr(infile[i],"%r")) break;
    
    //如果某个infile[i]含有流动站ID的替换符
    if (i<n) { /* include rover keywords */
        if (!(rov_=(char *)malloc(strlen(rov)+1))) return 0;
        strcpy(rov_,rov);
        
        for (i=0;i<n;i++) {
            if (!(ifile[i]=(char *)malloc(1024))) {
                free(rov_); for (;i>=0;i--) free(ifile[i]);
                return 0;
            }
        }
        for (p=rov_;;p=q+1) { /* for each rover */
            if ((q=strchr(p,' '))) *q='\0';
            
            if (*p) {
                strcpy(proc_rov,p);
                if (ts.time) time2str(ts,s,0); else *s='\0';
                if (checkbrk("reading    : %s",s)) {
                    stat=1;
                    break;
                }
                for (i=0;i<n;i++) reppath(infile[i],ifile[i],t0,p,"");
                reppath(outfile,ofile,t0,p,"");
                
                /* execute processing session */
                stat=execses(ts,te,ti,popt,sopt,fopt,flag,ifile,index,n,ofile);
            }
            if (stat==1||!q) break;
        }
        free(rov_); for (i=0;i<n;i++) free(ifile[i]);
    }
    else {
        /* execute processing session */
        stat=execses(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile);
    }
    return stat;
}

 2.2.1.1 execses函数：读取各种文件并获取基准站的位置

execses函数读取各种文件，并将文件中的内容赋值到程序的结构体内，获取基准站的位置，根据滤波方向调用procpos函数进行下一步解算

.trace文件的生成、文件读取相关trace文件内容的生成，均在execses中

/* execute processing session --------------------------------------------
gtime_t ts              I   处理的起始时间，写0表示不限制
gtime_t te              I   处理的起始时间，写0表示不限制
double ti               I   处理的间隔时间 (s)，写0表示不限制，全处理
const prcopt_t *popt    I   处理选项结构体
const solopt_t *sopt    I   结果选项结构体
const filopt_t *fopt    I   文件选项结构体
int flag                I   用于控制输出
char **infile           I   传入文件路径数组首地址
const int *index        I   传入文件路径数组首地址
int n                   I   传入文件数量
char *outfile           I   输出文件的路径，写0表示stdout终端
const char *rov         I   流动站ID列表，空格隔开
const char *base        I   基准站ID列表，空格隔开
------------------------------------------------------------------------*/
static int execses(gtime_t ts, gtime_t te, double ti, const prcopt_t *popt,
                   const solopt_t *sopt, const filopt_t *fopt, int flag,
                   char **infile, const int *index, int n, char *outfile)

 处理过程

 初始化定义

{
    FILE *fp;
    prcopt_t popt_=*popt;
    char tracefile[1024],statfile[1024],path[1024],*ext;
    
    trace(3,"execses : n=%d outfile=%s\n",n,outfile);
    

 调用traceclose函数 、traceopen函数 、tracelevel函数，先关闭原有trace，打开trace文件，并设置trace等级

    /* open debug trace */  //打开trace文件，并设置trace等级
    if (flag&&sopt->trace>0) {
        if (*outfile) {
            strcpy(tracefile,outfile);
            strcat(tracefile,".trace");
        }
        else {
            strcpy(tracefile,fopt->trace);
        }
        traceclose();
        traceopen(tracefile);
        tracelevel(sopt->trace);
    }

调用readtec函数 ，读取电离层TEC文件，TEC:Total electronic content 总电子含量 

 /* read ionosphere data file */ //读取电离层TEC文件
    if (*fopt->iono&&(ext=strrchr(fopt->iono,'.'))) {
        if (strlen(ext)==4&&(ext[3]=='i'||ext[3]=='I')) {
            reppath(fopt->iono,path,ts,"","");
            readtec(path,&navs,1);  //TEC:Total electronic content 总电子含量
        }
    }

调用readerp函数，读取地球自转参数ERP文件

    /* read erp data */ //读取地球自转参数ERP文件
    if (*fopt->eop) {
        free(navs.erp.data); navs.erp.data=NULL; navs.erp.n=navs.erp.nmax=0;
        reppath(fopt->eop,path,ts,"","");
        if (!readerp(path,&navs.erp)) {
            showmsg("error : no erp data %s",path);
            trace(2,"no erp data %s\n",path);
        }
    }

调用readobsnav函数，读取OBS和NAV文件 

    /* read obs and nav data */ //读取OBS和NAV文件
    if (!readobsnav(ts,te,ti,infile,index,n,&popt_,&obss,&navs,stas)) return 0;

调用readdcb函数，读取差分码偏差DCB参数，一种硬件误差 

    /* read dcb parameters */   //读取差分码偏差DCB参数，一种硬件误差
    if (*fopt->dcb) {
        reppath(fopt->dcb,path,ts,"","");
        readdcb(path,&navs,stas);
    }

调用setpcv函数，读取天线参数，PCV：天线相位中心变化 

 /* set antenna paramters */ //读取天线参数，PCV：天线相位中心变化
    if (popt_.mode!=PMODE_SINGLE) {
        setpcv(obss.n>0?obss.data[0].time:timeget(),&popt_,&navs,&pcvss,&pcvsr,
               stas);
    }

调用readotl函数，读取潮汐参数 

   /* read ocean tide loading parameters */    //读取潮汐参数
    if (popt_.mode>PMODE_SINGLE&&*fopt->blq) {
        readotl(&popt_,fopt->blq,stas);
    }
 

①FIXED模式，调用antpos函数得到流动站坐标 

②DGPS、KINEMA、STATIC模式，调用antpos函数得到基准站坐标 

   /* rover/reference fixed position */    //FIXED模式，调用antpos()得到流动站坐标
    if (popt_.mode==PMODE_FIXED) {
        if (!antpos(&popt_,1,&obss,&navs,stas,fopt->stapos)) {
            freeobsnav(&obss,&navs);
            return 0;
        }
    }
    else if (PMODE_DGPS<=popt_.mode&&popt_.mode<=PMODE_STATIC) {    //DGPS、KINEMA、STATIC模式，调用antpos()得到基准站坐标
        if (!antpos(&popt_,2,&obss,&navs,stas,fopt->stapos)) {
            freeobsnav(&obss,&navs);
            return 0;
        }
    }
   

调用rtkclosestat函数 、rtkopenstat函数，打开结果统计文件 

    /* open solution statistics */  //打开结果统计文件
    if (flag&&sopt->sstat>0) {
        strcpy(statfile,outfile);
        strcat(statfile,".stat");
        rtkclosestat();
        rtkopenstat(statfile,sopt->sstat);
    }

调用outhead函数，写输出结果文件的文件头 。结果文件的文件尾在procpos函数内调用outsol函数输出

    /* write header to output file */   //写输出结果文件的文件头
    if (flag&&!outhead(outfile,infile,n,&popt_,sopt)) {
        freeobsnav(&obss,&navs);
        return 0;
    }
    iobsu=iobsr=isbs=revs=aborts=0;

判断滤波类型，用不同的方式调用procpos函数进行下一步解算：

①forward 前向滤波：iobsu=iobsr=isbs=revs=0，直接调用procpos函数
②backward 后向滤波：revs=1，iobsu=iobsr=obss.n-1 ，isbs=sbss.n-1 ，再调用procpos函数
③combined ：先算前向滤波的结果，设置revs、iobsu、iobsr、isbs值之后再算后向滤波的结果，最后调用combress函数结合

注：

        前向滤波和后向滤波调用procpos函数传参相同，两者区别在于procpos函数内会调用inputobs函数。针对不同的滤波解算类型，inputobs函数内读取文件数据的顺序不同

        revs：0：forward；1：backward

        iobsu：当前流动站观测数据下标

        iobsr：当前参考站观测数据下标        

        isbs：当前sbas数据下标

    if (popt_.mode==PMODE_SINGLE||popt_.soltype==0) {
        if ((fp=openfile(outfile))) {
            procpos(fp,&popt_,sopt,0); /* forward */    //前向滤波
            fclose(fp);
        }
    }
    else if (popt_.soltype==1) {
        if ((fp=openfile(outfile))) {
            revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;
            procpos(fp,&popt_,sopt,0); /* backward */   //后向滤波
            fclose(fp);
        }
    }
    else { /* combined */
        //开辟内存空间
        solf=(sol_t *)malloc(sizeof(sol_t)*nepoch);     //前向结果
        solb=(sol_t *)malloc(sizeof(sol_t)*nepoch);     //后向结果
        rbf=(double *)malloc(sizeof(double)*nepoch*3);  //前向基准站坐标
        rbb=(double *)malloc(sizeof(double)*nepoch*3);  //后向基准站坐标
        
        if (solf&&solb) {   //判断内存开辟成功
            isolf=isolb=0;
            procpos(NULL,&popt_,sopt,1); /* forward */      //前向滤波
            revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;
            procpos(NULL,&popt_,sopt,1); /* backward */     //后向滤波

            //虽然前向滤波和后向滤波调用procpos函数的源代码相同（如下所示），
            //但是两者最主要的一个区别就是由于procpos函数内会调用inputobs函数，
            //然而针对不同的滤波解算类型，inputobs函数内读取文件数据的顺序不同
            /* combine forward/backward solutions */
            if (!aborts&&(fp=openfile(outfile))) {
                combres(fp,&popt_,sopt);
                fclose(fp);
            }
        }
        else showmsg("error : memory allocation");
        free(solf);
        free(solb);
        free(rbf);
        free(rbb);
    }

调用freeobsnav函数释放obs->data 、nav->eph 、nav->geph 、nav->seph

    /* free obs and nav data */
    freeobsnav(&obss,&navs);
    
    return aborts?1:0;
}

 2.2.1.1.1 antpos函数：获取流动站或基站坐标

得到坐标，参数rcvno传1得到流动站坐标，传0得到基准站坐标

static int antpos(prcopt_t *opt, int rcvno, const obs_t *obs, const nav_t *nav,
                  const sta_t *sta, const char *posfile)

 初始化定义

{
    double *rr=rcvno==1?opt->ru:opt->rb,
            del[3],pos[3],dr[3]={0};
    int i,
    postype=rcvno==1?opt->rovpos:opt->refpos;
    char *name;
    
    trace(3,"antpos  : rcvno=%d\n",rcvno);

①postype=POSOPT_SINGLE ：调用avepos函数利用基准站的观测文件计算其SPP定位结果作为基准站的坐标

avepos函数：通过nav和多个obs单点定位计算位置，存到rr中 

    if (postype==POSOPT_SINGLE) { /* average of single position */  //利用基准站的观测文件计算其SPP定位结果作为基准站的坐标
        if (!avepos(rr,rcvno,obs,nav,opt)) {
            showmsg("error : station pos computation");
            return 0;
        }
    }

②postype=POSOPT_FILE ：调用getstapos函数从pos文件读取基准站坐标 

getstapos函数：从posfile中读取基准站坐标

    else if (postype==POSOPT_FILE) { /* read from position file */  //从pos文件读取基准站坐标
        name=stas[rcvno==1?0:1].name;
        if (!getstapos(posfile,name,rr)) {
            showmsg("error : no position of %s in %s",name,posfile);
            return 0;
        }
    }

③postype=POSOPT_RINEX ：从rinex头文件中获取测站经过相位中心改正的位置数据。头文件中的测站数据经过读取后已存到stas中

    else if (postype==POSOPT_RINEX) { /* get from rinex header */   //从基准站的OBS观测文件的文件头部分读取基准站坐标
        if (norm(stas[rcvno==1?0:1].pos,3)<=0.0) {      //如果没有坐标数据，报错
            showmsg("error : no position in rinex header");
            trace(1,"no position position in rinex header\n");
            return 0;
        }
        //天线相位中心偏差改正
        /* antenna delta */
        if (stas[rcvno==1?0:1].deltype==0) { /* enu */
            for (i=0;i<3;i++) del[i]=stas[rcvno==1?0:1].del[i];
            del[2]+=stas[rcvno==1?0:1].hgt;
            ecef2pos(stas[rcvno==1?0:1].pos,pos);
            enu2ecef(pos,del,dr);
        }
        else { /* xyz */
            for (i=0;i<3;i++) dr[i]=stas[rcvno==1?0:1].del[i];
        }
        for (i=0;i<3;i++) rr[i]=stas[rcvno==1?0:1].pos[i]+dr[i];
    }
    return 1;
}

sta_t结构体

static sta_t stas[MAXRCV];      /* station infomation */

typedef struct {        /* station parameter type */
    char name   [MAXANT]; /* marker name */
    char marker [MAXANT]; /* marker number */
    char antdes [MAXANT]; /* antenna descriptor */
    char antsno [MAXANT]; /* antenna serial number */
    char rectype[MAXANT]; /* receiver type descriptor */
    char recver [MAXANT]; /* receiver firmware version */
    char recsno [MAXANT]; /* receiver serial number */
    int antsetup;       /* antenna setup id */
    int itrf;           /* ITRF realization year */
    int deltype;        /* antenna delta type (0:enu,1:xyz) */
    double pos[3];      /* station position (ecef) (m) */
    double del[3];      /* antenna position delta (e/n/u or x/y/z) (m) */
    double hgt;         /* antenna height (m) */
    int glo_cp_align;   /* GLONASS code-phase alignment (0:no,1:yes) */
    double glo_cp_bias[4]; /* GLONASS code-phase biases {1C,1P,2C,2P} (m) */
} sta_t;

  2.2.1.1.2 outhead函数：创建输出结果文件，写入文件头

static int outhead(const char *outfile, char **infile, int n,
                   const prcopt_t *popt, const solopt_t *sopt)
{
    FILE *fp=stdout;    //fp默认初始为stdout
    
    trace(3,"outhead: outfile=%s n=%d\n",outfile,n);
    
    if (*outfile) {
        createdir(outfile); //递归的创建文件夹
        
        if (!(fp=fopen(outfile,"wb"))) {    //wb：以写的方式打开二进制文件
            showmsg("error : open output file %s",outfile);
            return 0;
        }
    }
    /* output header */
    outheader(fp,infile,n,popt,sopt);
    
    if (*outfile) fclose(fp);
    
    return 1;
}

  2.2.1.1.3 openfile函数：以追加的方式打开结果文件

static FILE *openfile(const char *outfile)
{
    trace(3,"openfile: outfile=%s\n",outfile);
    
    return !*outfile?stdout:fopen(outfile,"ab");    //ab：以追加的方式打开二进制文件
}

  2.2.1.1.4 procpos函数：逐历元解算流动站和基准站数据

procpos函数开始正式整个流动站和基准站逐历元处理。每次循环都通过inputobs函数读取一个历元的数据，并调用rtkpos函数对该历元的数据进行解算

/* process positioning -------------------------------------------
FILE *fp	   		   I/O 输出结果文件指针  
const prcopt_t *popt    I   处理选项结构体
const solopt_t *sopt    I   结果选项结构体
const filopt_t *fopt    I   文件选项结构体
int mode			   I   0：forward/backward、1：combined
-----------------------------------------------------------------*/
static void procpos(FILE *fp, const prcopt_t *popt, const solopt_t *sopt,
                    int mode)

 处理过程

初始化定义

{
    gtime_t time={0};
    sol_t sol={{0}};
    rtk_t rtk;
    obsd_t obs[MAXOBS*2]; /* for rover and base */
    double rb[3]={0};
    int i,
    nobs,
    n,
    solstatic,
    pri[]={6,1,2,3,4,5,1,6};

判断结果是否为静态,处理选项和结果选项都为静态才算静态

    solstatic=sopt->solstatic&&     //先判断结果是否为静态,处理选项和结果选项都为静态才算静态
              (popt->mode==PMODE_STATIC||popt->mode==PMODE_PPP_STATIC); 

调用rtkinit函数初始化rtk_t ，将popt结构体赋值给rtk的部分成员

    rtkinit(&rtk,popt);    //初始化rtk_t，主要将popt结构体赋值给rtk的部分成员
    rtcm_path[0]='\0';

while大循环，调用inputobs函数，每次取一个历元的观测数据obs[]

·排除禁用卫星的观测值
·PPP模式设置需要，调用corr_phase_bias_ssr函数进行相位的小数轴偏差改正
·调用rtkpos函数对当前历元进行解算

    //对每一个历元进行遍历求解和输出
    //获取当前历元观测值数nobs以及当前历元各观测记录obs[MAXOBS*2]
    while ((nobs=inputobs(obs,rtk.sol.stat,popt))>=0) { 
        /* exclude satellites */
        for (i=n=0;i<nobs;i++) {
            //satsys:传入satellite number，返回卫星系统(SYS_GPS,SYS_GLO,...) ，通过传入的指针prn传出PRN码。
            if ((satsys(obs[i].sat,NULL)&popt->navsys)&&
                popt->exsats[obs[i].sat-1]!=1) obs[n++]=obs[i]; //排除禁用卫星的观测值
        }
        if (n<=0) continue;
        //如果ppp模式设置了fractional cycle bias相位的小数轴偏差
        /* carrier-phase bias correction */
        if (!strstr(popt->pppopt,"-ENA_FCB")) {     
            corr_phase_bias_ssr(obs,n,&navs);   
        }
        //调用rtkpos()进行解算
        if (!rtkpos(&rtk,obs,n,&navs)) continue;
        
        //单forward/backward模式
        if (mode==0) { /* forward/backward */
            if (!solstatic) {   //不是静态模式就直接输出结果
                outsol(fp,&rtk.sol,rtk.rb,sopt);
            }
            else if (time.time==0||pri[rtk.sol.stat]<=pri[sol.stat]) {
                sol=rtk.sol;    
                for (i=0;i<3;i++) rb[i]=rtk.rb[i];
                if (time.time==0||timediff(rtk.sol.time,time)<0.0) {
                    time=rtk.sol.time;      //记录上一历元的时间
                }
            }
        }
        else if (!revs) { /* combined-forward */
            if (isolf>=nepoch) return;
            solf[isolf]=rtk.sol;
            for (i=0;i<3;i++) rbf[i+isolf*3]=rtk.rb[i]; 
            isolf++;
        }
        else { /* combined-backward */
            if (isolb>=nepoch) return;
            solb[isolb]=rtk.sol;
            for (i=0;i<3;i++) rbb[i+isolb*3]=rtk.rb[i];
            isolb++;
        }
    }

根据模式，输出结果，记录当前历元时间

    if (mode==0&&solstatic&&time.time!=0.0) {
        sol.time=time;
        outsol(fp,&sol,rb,sopt);
    }
    rtkfree(&rtk);
}

  2.2.1.1.5 combres函数：结合前后向滤波结果

调用smoother函数结合前后向滤波的结果

static void combres(FILE *fp, const prcopt_t *popt, const solopt_t *sopt)

处理过程

初始化定义

{
    gtime_t time={0};
    sol_t sols={{0}},sol={{0}};
    double tt,Qf[9],Qb[9],Qs[9],rbs[3]={0},rb[3]={0},rr_f[3],rr_b[3],rr_s[3];
    int i,j,k,solstatic,pri[]={0,1,2,3,4,5,1,6};
    
    trace(3,"combres : isolf=%d isolb=%d\n",isolf,isolb);

判断静态模式，处理选项和结果选项都得为静态

    //判断静态模式，处理选项和结果选项都得为静态
    solstatic=sopt->solstatic&&
              (popt->mode==PMODE_STATIC||popt->mode==PMODE_PPP_STATIC);

开始大循环，i:从前到后，取前向滤波的结果 ，j:从后到前，取后向滤波的结果 ，判断前后向滤波结果的时间差 tt

    //i:从前到后，取前向滤波的结果
    //j:从后到前，取后向滤波的结果
    for (i=0,j=isolb-1;i<isolf&&j>=0;i++,j--) {

①时间差大于DTTOL ，sols、rbs取时间早的结果，另一个结果的下标不变，进行下一次循环的判断

        //判断前后向滤波结果的时间差，时间差大于DTTOL，
        //sols、rbs取时间早的结果，另一个结果的下标不变，进行下一次循环的判断
        if ((tt=timediff(solf[i].time,solb[j].time))<-DTTOL) {  //如果前向时间迟于后向时间
            sols=solf[i];                                       
            for (k=0;k<3;k++) rbs[k]=rbf[k+i*3];                //把前向基站坐标赋值给rbs[]
            j++;            //j不变
        }
        else if (tt>DTTOL) {                        //如果前向时间早于后向时间
            sols=solb[j];                           
            for (k=0;k<3;k++) rbs[k]=rbb[k+j*3];    //把后向基站坐标赋值给rbs[]
            i--;            //i不变
        }

②时间差很小，solution status不同，sols、rbs取solution status小的结果

 //时间差很小，solution status不同，sols、rbs取solution status小的结果
        else if (solf[i].stat<solb[j].stat) {
            sols=solf[i];
            for (k=0;k<3;k++) rbs[k]=rbf[k+i*3];
        }
        else if (solf[i].stat>solb[j].stat) {
            sols=solb[j];
            for (k=0;k<3;k++) rbs[k]=rbb[k+j*3];
        }

③时间差很小，solution status相同，进行结合

        sols取前向滤波结果 ，时间取前后向时间的平均
        相对定位模式，若结果为固定解，调用valcomb()检验，如果失败将fix降级为float
        赋值前后向协方差给Qf、Qb ，调用smoother()进行前后向滤波结果结合，位置存在sols.rr[]，方差存在sols.qr[]
        同样的方式，对速度进行结合，位置存在sols.rr[]，方差存在sols.qv[]

        //时间差很小，solution status相同
        else {
            sols=solf[i];   //sols取前向滤波结果
            sols.time=timeadd(sols.time,-tt/2.0);   //时间取前后向时间的平均
            //相对定位模式，若结果为固定解，调用valcomb()检验，如果失败将fix降级为float
            if ((popt->mode==PMODE_KINEMA||popt->mode==PMODE_MOVEB)&&
                sols.stat==SOLQ_FIX) {
                /* degrade fix to float if validation failed */
                if (!valcomb(solf+i,solb+j)) sols.stat=SOLQ_FLOAT;
            }
//赋值前后向协方差给Qf、Qb，
            for (k=0;k<3;k++) {     //k+k*3是取对角线元素
                Qf[k+k*3]=solf[i].qr[k];
                Qb[k+k*3]=solb[j].qr[k];
            }
            Qf[1]=Qf[3]=solf[i].qr[3];  //赋值非对角线元素
            Qf[5]=Qf[7]=solf[i].qr[4];
            Qf[2]=Qf[6]=solf[i].qr[5];
            Qb[1]=Qb[3]=solb[j].qr[3];
            Qb[5]=Qb[7]=solb[j].qr[4];
            Qb[2]=Qb[6]=solb[j].qr[5];
            
            //调用smoother()进行前后向滤波结果结合，位置存在sols.rr[]，方差存在sols.qr[]
            if (popt->mode==PMODE_MOVEB) {  //如果是移动基线模式
                for (k=0;k<3;k++) rr_f[k]=solf[i].rr[k]-rbf[k+i*3]; //流动站坐标-基准站坐标得到基线
                for (k=0;k<3;k++) rr_b[k]=solb[j].rr[k]-rbb[k+j*3];
                if (smoother(rr_f,Qf,rr_b,Qb,3,rr_s,Qs)) continue;
                for (k=0;k<3;k++) sols.rr[k]=rbs[k]+rr_s[k];
            }
            else {
                if (smoother(solf[i].rr,Qf,solb[j].rr,Qb,3,sols.rr,Qs)) continue;
            }
            sols.qr[0]=(float)Qs[0];
            sols.qr[1]=(float)Qs[4];
            sols.qr[2]=(float)Qs[8];
            sols.qr[3]=(float)Qs[1];
            sols.qr[4]=(float)Qs[5];
            sols.qr[5]=(float)Qs[2];
            
            /* smoother for velocity solution */
            if (popt->dynamics) {
                for (k=0;k<3;k++) {
                    Qf[k+k*3]=solf[i].qv[k];
                    Qb[k+k*3]=solb[j].qv[k];
                }
                Qf[1]=Qf[3]=solf[i].qv[3];
                Qf[5]=Qf[7]=solf[i].qv[4];
                Qf[2]=Qf[6]=solf[i].qv[5];
                Qb[1]=Qb[3]=solb[j].qv[3];
                Qb[5]=Qb[7]=solb[j].qv[4];
                Qb[2]=Qb[6]=solb[j].qv[5];
                if (smoother(solf[i].rr+3,Qf,solb[j].rr+3,Qb,3,sols.rr+3,Qs)) continue;
                sols.qv[0]=(float)Qs[0];
                sols.qv[1]=(float)Qs[4];
                sols.qv[2]=(float)Qs[8];
                sols.qv[3]=(float)Qs[1];
                sols.qv[4]=(float)Qs[5];
                sols.qv[5]=(float)Qs[2];
            }
        }
        if (!solstatic) {
            outsol(fp,&sols,rbs,sopt);
        }
        else if (time.time==0||pri[sols.stat]<=pri[sol.stat]) {
            sol=sols;
            for (k=0;k<3;k++) rb[k]=rbs[k];
            if (time.time==0||timediff(sols.time,time)<0.0) {
                time=sols.time;
            }
        }
    }
        

 循环处理完之后，如果是静态模式且时间存在，调用outsol函数输出结果

    //循环处理完之后，如果是静态模式且时间存在，调用outsol()输出结果
    if (solstatic&&time.time!=0.0) {    
        sol.time=time;
        outsol(fp,&sol,rb,sopt);
    }
}

 输出结果

结果状态的宏定义：

#define SOLQ_NONE   0                   /* solution status: no solution */
#define SOLQ_FIX    1                   /* solution status: fix */
#define SOLQ_FLOAT  2                   /* solution status: float */
#define SOLQ_SBAS   3                   /* solution status: SBAS */
#define SOLQ_DGPS   4                   /* solution status: DGPS/DGNSS */
#define SOLQ_SINGLE 5                   /* solution status: single */
#define SOLQ_PPP    6                   /* solution status: PPP */
#define SOLQ_DR     7                   /* solution status: dead reconing */
#define MAXSOLQ     7                   /* max number of solution status */

• sol_t结构体：

因为协方差矩阵是对称的，qr、qv都只用6个元素就可存协方差矩阵，但计算的时候得转成3*3矩阵才行 

typedef struct {        /* solution type */
    gtime_t time;       /* time (GPST) */
    double rr[6];       /* position/velocity (m|m/s) */
                        /* {x,y,z,vx,vy,vz} or {e,n,u,ve,vn,vu} */
    float  qr[6];       /* position variance/covariance (m^2) */
                        /* {c_xx,c_yy,c_zz,c_xy,c_yz,c_zx} or */
                        /* {c_ee,c_nn,c_uu,c_en,c_nu,c_ue} */
    float  qv[6];       /* velocity variance/covariance (m^2/s^2) */
    double dtr[6];      /* receiver clock bias to time systems (s) */
    uint8_t type;       /* type (0:xyz-ecef,1:enu-baseline) */
    uint8_t stat;       /* solution status (SOLQ_???) */
    uint8_t ns;         /* number of valid satellites */
    float age;          /* age of differential (s) */
    float ratio;        /* AR ratio factor for valiation */
    float thres;        /* AR ratio threshold for valiation */
} sol_t;

2.2.2 readpreceph函数：

readpreceph函数遍历infile[]，判断，调用readsp3函数读取精密星历、调用readrnxc函数读取精密钟差，调用sbsreadmsg函数读取sbas文件，将RCTM的路径赋值给rtcm_file，调用init_rtcm函数初始化rtcm控制结构体

static void readpreceph(char **infile, int n, const prcopt_t *prcopt,
                        nav_t *nav, sbs_t *sbs)
{
    seph_t seph0={0};
    int i;
    char *ext;
    
    trace(2,"readpreceph: n=%d\n",n);
    
    nav->ne=nav->nemax=0;
    nav->nc=nav->ncmax=0;
    sbs->n =sbs->nmax =0;
    
    /* read precise ephemeris files */  //读精密星历sp3
    for (i=0;i<n;i++) {
        if (strstr(infile[i],"%r")||strstr(infile[i],"%b")) continue;
        readsp3(infile[i],nav,0);
    }
    /* read precise clock files */  //读精密钟差
    for (i=0;i<n;i++) {
        if (strstr(infile[i],"%r")||strstr(infile[i],"%b")) continue;
        readrnxc(infile[i],nav);
    }
    /* read sbas message files */   //读sbas文件
    for (i=0;i<n;i++) {
        if (strstr(infile[i],"%r")||strstr(infile[i],"%b")) continue;
        sbsreadmsg(infile[i],prcopt->sbassatsel,sbs);
    }
    /* allocate sbas ephemeris */   //为nav->seph开辟空间
    nav->ns=nav->nsmax=NSATSBS*2;   
    if (!(nav->seph=(seph_t *)malloc(sizeof(seph_t)*nav->ns))) {
         showmsg("error : sbas ephem memory allocation");
         trace(1,"error : sbas ephem memory allocation");
         return;
    }
    for (i=0;i<nav->ns;i++) nav->seph[i]=seph0; 
    
    /* set rtcm file and initialize rtcm struct */
    rtcm_file[0]=rtcm_path[0]='\0'; fp_rtcm=NULL;
    
    //遍历ifile，将后缀为RTCM3的路径赋值到rtcm_file,初始化rtcm控制结构体
    for (i=0;i<n;i++) {
        if ((ext=strrchr(infile[i],'.'))&&      
            (!strcmp(ext,".rtcm3")||!strcmp(ext,".RTCM3"))) {
            strcpy(rtcm_file,infile[i]);
            init_rtcm(&rtcm);
            break;
        }
    }
}

2.3  closeses函数：结束解算程序

结束解算程序，释放天线、geoid、erp、trace、fp_stat 。会调用closegeoid函数、rtkclosestat函数 、traceclose函数 

static void closeses(nav_t *nav, pcvs_t *pcvs, pcvs_t *pcvr)
{
    trace(3,"closeses:\n");
    
    /* free antenna parameters */
    free(pcvs->pcv); pcvs->pcv=NULL; pcvs->n=pcvs->nmax=0;
    free(pcvr->pcv); pcvr->pcv=NULL; pcvr->n=pcvr->nmax=0;
    
    /* close geoid data */
    closegeoid();
    
    /* free erp data */
    free(nav->erp.data); nav->erp.data=NULL; nav->erp.n=nav->erp.nmax=0;
    
    /* close solution statistics and debug trace */
    rtkclosestat();
    traceclose();
}

2.4 reppaths函数：替换文件路径的替换符并生成多个路径

根据ts、te分时间段，循环调用reppath函数，替换path[]中的替换符，存到repath[]中，返回文件数量

/* replace keywords in file path and generate multiple paths -------------------
* replace keywords in file path with date, time, rover and base station id
* generate multiple keywords-replaced paths
* args   : char   *path     I   file path (see below)
*          char   *rpath[]  O   file paths in which keywords replaced
*          int    nmax      I   max number of output file paths
*          gtime_t ts       I   time start (gpst)
*          gtime_t te       I   time end   (gpst)
*          char   *rov      I   rover id string        ("": not replaced)
*          char   *base     I   base station id string ("": not replaced)
* return : number of replaced file paths
* notes  : see reppath() for replacements of keywords.
*          minimum interval of time replaced is 900s.
*-----------------------------------------------------------------------------*/
extern int reppaths(const char *path, char *rpath[], int nmax, gtime_t ts,
                    gtime_t te, const char *rov, const char *base)
{
    gtime_t time;
    double tow,tint=86400.0;
    int i,n=0,week;
    
    trace(3,"reppaths: path =%s nmax=%d rov=%s base=%s\n",path,nmax,rov,base);
    
    if (ts.time==0||te.time==0||timediff(ts,te)>0.0) return 0; //如果起止时间为0，或ts>te，直接return
    
    if (strstr(path,"%S")||strstr(path,"%M")||strstr(path,"%t")) tint=900.0;    //15分钟
    else if (strstr(path,"%h")||strstr(path,"%H")) tint=3600.0;     //一小时
    
    tow=time2gpst(ts,&week);
    time=gpst2time(week,floor(tow/tint)*tint);  
    
    while (timediff(time,te)<=0.0&&n<nmax) {
        reppath(path,rpath[n],time,rov,base);
        if (n==0||strcmp(rpath[n],rpath[n-1])) n++;
        time=timeadd(time,tint);
    }
    for (i=0;i<n;i++) trace(3,"reppaths: rpath=%s\n",rpath[i]);
    return n;
}

2.4.1 reppath函数：替换文件路径的替换符

如果输入文件（file）中，含有替换符，则 reppath函数将文件名中的替换符调用repstr函数进行替换，保存到rpath中。替换符如下：

/* replace keywords in file path -----------------------------------------------
* replace keywords in file path with date, time, rover and base station id
* args   : char   *path     I   file path (see below)
*          char   *rpath    O   file path in which keywords replaced (see below)
*          gtime_t time     I   time (gpst)  (time.time==0: not replaced)
*          char   *rov      I   rover id string        ("": not replaced)
*          char   *base     I   base station id string ("": not replaced)
* return : status (1:keywords replaced, 0:no valid keyword in the path,
*                  -1:no valid time)
* notes  : the following keywords in path are replaced by date, time and name
*              %Y -> yyyy : year (4 digits) (1900-2099)
*              %y -> yy   : year (2 digits) (00-99)
*              %m -> mm   : month           (01-12)
*              %d -> dd   : day of month    (01-31)
*              %h -> hh   : hours           (00-23)
*              %M -> mm   : minutes         (00-59)
*              %S -> ss   : seconds         (00-59)
*              %n -> ddd  : day of year     (001-366)
*              %W -> wwww : gps week        (0001-9999)
*              %D -> d    : day of gps week (0-6)
*              %H -> h    : hour code       (a=0,b=1,c=2,...,x=23)
*              %ha-> hh   : 3 hours         (00,03,06,...,21)
*              %hb-> hh   : 6 hours         (00,06,12,18)
*              %hc-> hh   : 12 hours        (00,12)
*              %t -> mm   : 15 minutes      (00,15,30,45)
*              %r -> rrrr : rover id
*              %b -> bbbb : base station id
*-----------------------------------------------------------------------------*/
extern int reppath(const char *path, char *rpath, gtime_t time, const char *rov,
                   const char *base)

reppaths()需要ts和te、而reppath只用ts 

extern int reppath(const char *path, char *rpath, gtime_t time, const char *rov,
                   const char *base)
{
    double ep[6],ep0[6]={2000,1,1,0,0,0};
    int week,dow,doy,stat=0;
    char rep[64];
    
    strcpy(rpath,path);
    
    if (!strstr(rpath,"%")) return 0;           //找不到%号直接结束
    if (*rov ) stat|=repstr(rpath,"%r",rov );   //如果有，替换基准站、流动站名
    if (*base) stat|=repstr(rpath,"%b",base);
    if (time.time!=0) {
        //把时间从gtime_t转为ep数组、DOW、DOY
        time2epoch(time,ep);    
        ep0[0]=ep[0];
        dow=(int)floor(time2gpst(time,&week)/86400.0);
        doy=(int)floor(timediff(time,epoch2time(ep0))/86400.0)+1;
        //把要替换的内容存到rep中，再用rep替换
        sprintf(rep,"%02d",  ((int)ep[3]/3)*3);   stat|=repstr(rpath,"%ha",rep);
        sprintf(rep,"%02d",  ((int)ep[3]/6)*6);   stat|=repstr(rpath,"%hb",rep);
        sprintf(rep,"%02d",  ((int)ep[3]/12)*12); stat|=repstr(rpath,"%hc",rep);
        sprintf(rep,"%04.0f",ep[0]);              stat|=repstr(rpath,"%Y",rep);
        sprintf(rep,"%02.0f",fmod(ep[0],100.0));  stat|=repstr(rpath,"%y",rep);
        sprintf(rep,"%02.0f",ep[1]);              stat|=repstr(rpath,"%m",rep);
        sprintf(rep,"%02.0f",ep[2]);              stat|=repstr(rpath,"%d",rep);
        sprintf(rep,"%02.0f",ep[3]);              stat|=repstr(rpath,"%h",rep);
        sprintf(rep,"%02.0f",ep[4]);              stat|=repstr(rpath,"%M",rep);
        sprintf(rep,"%02.0f",floor(ep[5]));       stat|=repstr(rpath,"%S",rep);
        sprintf(rep,"%03d",  doy);                stat|=repstr(rpath,"%n",rep);
        sprintf(rep,"%04d",  week);               stat|=repstr(rpath,"%W",rep);
        sprintf(rep,"%d",    dow);                stat|=repstr(rpath,"%D",rep);
        sprintf(rep,"%c",    'a'+(int)ep[3]);     stat|=repstr(rpath,"%H",rep);
        sprintf(rep,"%02d",  ((int)ep[4]/15)*15); stat|=repstr(rpath,"%t",rep);
    }
    else if (strstr(rpath,"%ha")||strstr(rpath,"%hb")||strstr(rpath,"%hc")||
             strstr(rpath,"%Y" )||strstr(rpath,"%y" )||strstr(rpath,"%m" )||
             strstr(rpath,"%d" )||strstr(rpath,"%h" )||strstr(rpath,"%M" )||
             strstr(rpath,"%S" )||strstr(rpath,"%n" )||strstr(rpath,"%W" )||
             strstr(rpath,"%D" )||strstr(rpath,"%H" )||strstr(rpath,"%t" )) {
        return -1; /* no valid time */
    }
    return stat;
}

3、postpos的main函数实现

 在rnx2rtkp函数中修改main函数为如下即可（在运行之前需要在vs中配置环境，不赘述）

rnx2rtkp函数为rtklib自带的函数

观测值文件和导航文件均为rtklib自带的测试数据

/* rnx2rtkp main -------------------------------------------------------------*/
int main(int argc, char** argv)
{
    
        gtime_t ts = { 0 }, te = { 0 };
        prcopt_t opt = prcopt_default;
        solopt_t sopt = solopt_default;
        filopt_t filopt = { 0 };
        char* infile[3];
        infile[0] = "‘输入rtklib的存储路径’  RTKLIB-rtklib_2.4.3\\test\\data\\rinex\\07590920.05o";
        infile[1] = "‘输入rtklib的存储路径’  RTKLIB-rtklib_2.4.3\\test\\data\\rinex\\30400920.05o";
        infile[2] = "‘输入rtklib的存储路径’  RTKLIB-rtklib_2.4.3\\test\\data\\rinex\\07590920.05n";
        char* outfile = "C:\\Users\\Administrator\\Desktop\\mypos.pos";
        opt.navsys = SYS_GPS;
        opt.mode = PMODE_MOVEB;
        sopt.posf = SOLF_ENU;//输出的坐标为enu形式
        postpos(ts, te, 0.0, 0.0, &opt, &sopt, &filopt, infile, 3, outfile, "", "");

}

4、rtkpos函数：计算接收机的位置、速度和钟差

rtkpos函数实际在procpos函数下运行，因其比较重要，故后续单独编写rtkpos.c的内容