PPP
1.读取文件
1.1RINEX观测文件(.**o)
观测文件这部分在SPP里面已经进行详细的介绍
1.2精密星历文件(.sp3)
调用readsp3()函数
function nav=readsp3(nav,fname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% only surpport to read ".sp3" or ".eph" file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
%首先是查找文件名,是否是空的,如果没有发现文件名的格式为.sp3、eph等,提示错误,文件格式不支持
if ~isempty(strfind(fname,'.'))
if isempty(strfind(fname,'sp3'))&&isempty(strfind(fname,'SP3'))&&...
isempty(strfind(fname,'eph'))&&isempty(strfind(fname,'EPH'))
error('The file %s is not surpported!!!',fname);
end
end
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);%这段代码的作用是从字符串 fname 中截取最后一个分隔符 glc.sep 后面的部分,并将结果保存到变量 fname0 中。这可能用于文件名的处理,提取文件名中的文件扩展名或最后一个目录的名称。
fprintf('Info:reading sp3 file %s...',fname0);%读取后缀为.sp3文件的名字,最终在命令行窗口输出Info:reading obs file wum20464.sp3...over
fid=fopen(fname);%打开sp3文件
% read sp3 header
[headinfo,fid]=decode_sp3h(fid);%读取精密星历的头文件
% read sp3 body
nav=decode_sp3b(nav,headinfo,fid);%读取精密星历的主体文件
fprintf('over\n');
return
在读取精密星历文件的时候分为两部分,先通过调用decode_sp3h()函数来读取精密星历的头文件,之后通过调用decode_sp3b()函数来读取精密星历的主体文件。
调用decode_sp3h()函数
function [headinfo,fid]=decode_sp3h(fid)%一个函数 decode_sp3h(fid),用于解析 SP3 文件的头部信息。
%函数的输入参数是文件句柄 fid,表示要解析的 SP3 文件。函数通过逐行读取文件,提取文件头部信息,并将其保存在一个 headinfo 结构体中。
%函数首先定义了一些变量 type、time、ns、sats、tsys 和 bfact,用来存储头部信息的各个字段。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read file header information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc gls
type=' '; time=gls.gtime; ns=0;
sats=zeros(1,glc.MAXSAT); tsys=''; bfact=zeros(1,2);
n=0; idxns=0; idxc=0; idxf=0; idx_break=0;
while ~feof(fid)%然后,函数通过循环读取文件的每一行,并根据行的不同内容,执行不同的操作来解析头部信息。
line=fgetl(fid);n=n+1;
if n==1%当 n 等于 1 时,解析第一行,获取文件类型 type 和时间信息 time。
type=line(3);
time=str2time(line(4:31));
elseif strcmp(line(1),'+')%当行以 + 开头时,表示卫星编号信息,解析卫星编号,并存储在 sats 数组中。
if n==3,ns=str2double(line(3:6));end
if idxns<ns
j=0;
while j<17
sys=code2sys(line(10+3*j));%在每次迭代中,首先使用函数 code2sys(line(10+3*j)) 将字符串中的字符转换为卫星系统代码(如 GPS、GLO、GAL 等),存储在变量 sys 中。
prn=str2double(line(11+3*j:12+3*j));%使用 str2double(line(11+3*j:12+3*j)) 解析出卫星的 PRN 号码,存储在变量 prn 中。
if sys==glc.SYS_QZS,prn=prn+192;end%如果卫星系统为 QZSS (glc.SYS_QZS),则将 PRN 号码加上 192,因为 QZSS 的 PRN 号码是 GPS PRN 号码加上 192。
if idxns<=glc.MAXSAT
sats(idxns+1)=satno(sys,prn);
idxns=idxns+1;
end%最后,如果索引 idxns 小于等于 glc.MAXSAT,则将卫星的编号存储在数组 sats 的对应位置上,并将索引 idxns 加 1。
j=j+1;
end
end
elseif strcmp(line(1:2),'%c')&&idxc==0%当行以 %c 开头时,表示时钟系统信息,解析并存储在 tsys 变量中。
tsys=line(10:12);
idxc=idxc+1;
elseif strcmp(line(1:2),'%f')&&idxf==0 %当行以 %f 开头时,表示钟差比例因子信息,解析并存储在 bfact 数组中。
bfact(1)=str2double(line(4:13));
bfact(2)=str2double(line(15:26));
idxf=idxf+1;
elseif strcmp(line(1:2),'/*')%当行以 /* 开头时,表示头部信息解析结束,跳出循环。
idx_break=idx_break+1;
if idx_break==4
break;
end
end
end
%最后,将解析得到的头部信息存储在 headinfo 结构体中,并返回该结构体作为函数的输出。
headinfo.type=type;
headinfo.time=time;
headinfo.ns=ns;
headinfo.sats=sats;
headinfo.tsys=tsys;
headinfo.bfact=bfact;
return
调用decode_sp3b()函数
function nav=decode_sp3b(nav,headinfo,fid)
global glc gls
peph=gls.peph;%精密星历
NMAX=10000; nav.peph=repmat(gls.peph,NMAX,1);
if headinfo.type=='P'
ns=headinfo.ns;
else
ns=2*headinfo.ns;
end
%这段代码的作用是初始化导航结构体 nav 中的精密星历部分,并根据头部信息结构体 headinfo 中的数据确定导航结构体中精密星历的数量 ns。
%按照时间步进遍历)是一种处理数据或执行操作的方法。它通常用于按照时间顺序逐个处理一系列数据点或执行一系列操作。在时间序列数据处理中,"traversal by epoch"指的是按照时间的顺序,从最早的数据点开始,逐个处理每个数据点,并按照一定的步进(如固定时间间隔)依次处理下一个数据点,直到遍历完所有的数据点。
while ~feof(fid)
line=fgetl(fid);
if strcmp(line,'EOF'),continue;end
if line(1)~='*'
fprintf('sp3 invalid epoch\n');
continue;
end
time=str2time(line(4:31));
if strcmp(headinfo.tsys,'UTC'),time=utc2gpst(time);end
peph.time=time;
%这段代码的作用是解析 SP3 文件中的每个时刻的数据,并将时间信息存储在精密星历结构体 peph 中的 time 字段中。该代码还进行了一些错误处理,例如检查是否到达文件末尾、验证时间信息的有效性等。
%初始化精密星历的结构体
for i=1:glc.MAXSAT
for j=1:4
peph.pos(i,j)=0;%位置
peph.std(i,j)=0;%标准差
peph.vel(i,j)=0;%速度
peph.vst(i,j)=0;%速度标准差
end
for j=1:3
peph.cov(i,j)=0;%协方差
peph.vco(i,j)=0;%速度协方差
end
end
%这段代码的作用是将精密星历结构体 peph 中每个卫星的位置、标准差、速度、速度标准差、协方差矩阵和速度协方差矩阵等字段初始化为 0,为后续的数据填充和存储做好准备。
for i=1:ns
line=fgetl(fid);
if size(line,2)<4||(line(1)~='P'&&line(1)~='V'),continue;end
if line(2)==' '
sys=glc.SYS_GPS;
else
sys=code2sys(line(2));
end
prn=str2double(line(3:4));
if sys==glc.SYS_QZS,prn=prn+192;end
sat=satno(sys,prn);
if ~sat,continue;end
%此代码段的作用是解析 SP3 文件中每个卫星的信息,并通过解析卫星系统类型和 PRN 号将其转换为卫星号。同时,此代码段还进行了一些错误处理,例如检查行的长度是否符合要求、验证卫星系统类型和 PRN 号的有效性等。
for j=1:4%对于前三个 j 值,执行以下赋值操作:m=2:将变量 m 的值设置为 2。k1=1000:将变量 k1 的值设置为 1000。%base=headinfo.bfact(1):将变量 base 的值设置为 headinfo.bfact(1) 的值,headinfo.bfact 是一个结构体或数组,表示其中的第一个元素。%base=headinfo.bfact(1):将变量 base 的值设置为 headinfo.bfact(1) 的值,headinfo.bfact 是一个结构体或数组,表示其中的第一个元素。%k3=0.1:将变量 k3 的值设置为 0.1。k4=10^-7:将变量 k4 的值设置为 0.0000001。
if j<=3
m=2;k1=1000;base=headinfo.bfact(1);k2=10^-3;k3=0.1;k4=10^-7;
else
m=3;k1=10^-6;base=headinfo.bfact(2);k2=10^-12;k3=10^-10;k4=10^-16;
end%m=3:将变量 m 的值设置为 3。k1=10^-6:将变量 k1 的值设置为 0.000001。base=headinfo.bfact(2):将变量 base 的值设置为 headinfo.bfact(2) 的值,headinfo.bfact 是一个结构体或数组,表示其中的第二个元素。k2=10^-12:将变量 k2 的值设置为 0.000000000001。k3=10^-10:将变量 k3 的值设置为 0.0000000001。
val=str2double(line(4+14*(j-1)+1:4+14*(j-1)+14));
if size(line,2)>61
std=str2double(line(61+3*(j-1)+1:61+3*(j-1)+m));
else
std=0;
end%段代码的作用是解析文本行 line 中特定位置的字符串,并将其转换为数字,分别存储在变量 val 和 std 中。具体的位置计算是根据变量 j 和之前的条件判断得出的。
if line(1)=='P'
if val~=0&&abs(val-999999.999999)>=10^-6
peph.pos(sat,j)=val*k1; v=1;%这两行代码的作用是对变量 peph.pos 和 v 进行赋值操作。根据给定的索引 sat 和 j,将 val 与 k1 相乘的结果赋值给 peph.pos。同时,将变量 v 的值设置为 1。
end
if base>0&&std>0
peph.std(sat,j)=base^std*k2;
end%这段代码的作用是根据不同的条件,对 peph 数据结构进行赋值操作。如果 line 的第一个字符为 'P',并且满足一定条件,将 val 乘以特定的系数并赋值给 peph.pos,将 base 的 std 次方乘以另一个系数并赋值给 peph.std。同时,在满足特定条件时,设置变量 v 的值为 1。
elseif v
if val~=0&&abs(val-999999.999999)>=10^-6
peph.vel(sat,j)=val*k3;
end
if base>0&&std>0
peph.vst(sat,j)=base^std*k4;
end%这段代码是在之前的代码的if line(1)=='P'判断条件不成立的情况下执行的,它执行了以下操作:首先,通过判断变量 v 是否为真(非零)来执行接下来的代码块。
%判断 val 不等于 0,并且 val 的绝对值与 999999.999999 的差值大于等于 10 的负 6 次方(即绝对值大于等于 0.000001)。如果条件满足,执行以下操作:
%将变量 val 乘以 k3,并将结果赋值给 peph.vel(sat, j)。peph.vel 是一个数组或矩阵,并且 sat 和 j 用作索引。
%如果 base 大于 0 且 std 大于 0,则执行以下操作:
%将 base 的 std 次方乘以 k4,并将结果赋值给 peph.vst(sat, j)。peph.vst 是一个数组或矩阵,并且 sat 和 j 用作索引。
%简而言之,这段代码的作用是在 v 为真的情况下,根据一些条件对 peph 数据结构的 vel 和 vst 进行赋值操作。如果条件满足,将 val 乘以特定的系数,并赋值给 peph.vel 或 peph.vst。
end
end
end
if v
if nav.np>size(nav.peph,1)
nav.peph(nav.np+1:nav.np+NMAX,1)=repmat(gls.peph,NMAX,1);
end%nav.np:这是一个变量,表示当前位置;nav.peph:这是一个数组或矩阵,存储了位置信息。size(nav.peph,1):这是一个函数,返回 nav.peph 的行数。nav.np < size(nav.peph,1):这是一个条件判断,检查 nav.np 是否小于 nav.peph 的行数。
nav.peph(nav.np+1)=peph;
nav.np=nav.np+1;
end%这段代码的作用是在 v 为真的情况下,根据一些条件将变量 peph 添加到 nav.peph 中的特定位置。如果条件满足,则根据需要扩展 nav.peph 的大小,并将 peph 添加到 nav.peph 中的下一个空闲位置,同时更新 nav.np 的值。
end
fclose(fid);
if nav.np<size(nav.peph,1)
nav.peph(nav.np+1:end,:)=[];
end%这段代码的目的是将 nav.peph 中超出当前索引范围的部分数据移除。如果 nav.np 较小,导致 nav.peph 中存在多余的行数,则通过删除超出范围的行来保持数据的正确性。
return
1.3精密钟差文件(.clk)
调用readclk()函数
function nav=readclk(nav,opt,fname)
global glc gls
NMAX=10000; nav.pclk=repmat(gls.pclk,NMAX,1);
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);%这段代码主要是对变量进行赋值操作。它将 gls.pclk 在行方向上复制了 NMAX 次,并将结果赋值给 nav.pclk。然后,通过 strfind 函数找到 fname 中最后一个斜杠的位置,并将斜杠后面的部分赋值给 fname0。
fprintf('Info:reading clk file %s...',fname0);%读取后缀为.clk文件的名字,最终在命令行窗口输出Info:reading obs file wum20464.clk...over
fid=fopen(fname);%代开精密钟差文件
% 读取精密钟差的头文件
while ~feof(fid)
line=fgetl(fid);
label=line(61:end);
if size(line,2)<=60
continue;
elseif ~isempty(strfind(label,'RINEX VERSION / TYPE'))
ver=str2double(line(1:9)); %#示例数据版本为3.0,已经经过对比源文件进行查看
type=line(21);%根据示例数据类型为'C'
elseif ~isempty(strfind(label,'END OF HEADER'))
break;
end
end
if type~='C'%type为C的话,那么表示下载的是正确的clk文件的数据,如果type不为C的话,提示错误“该文件不是clock file”
error('The file %s is not clock file!!!',fname0);
end
mask=set_sysmask(opt.navsys);%这个函数可以用于根据输入参数的不同设置系统屏蔽状态,进而控制系统的行为。
% 读取精密钟差的主体文件
while ~feof(fid)
buff(1:glc.MAXRNXLEN)=' ';
line=fgets(fid);
buff(1:size(line,2))=line;
if size(buff,2)<8,continue; end
if ~strcmp(buff(1:2),'AS'),continue;end
%这段代码的作用是按行读取文件内容,并对读取到的每行进行一些条件判断。它首先将 buff 数组初始化为空格字符,然后使用 fgets 函数逐行读取文件内容,并将其复制到 buff 数组中。接着,它进行一些条件判断,例如判断行的长度是否足够、是否以特定的字符开头等,从而决定是否继续处理该行。
satid=buff(4:8);
sat=satid2no(satid);
[sys,~]=satsys(sat);
if sys==0,continue;end
if mask(sys)==0,continue;end
%这段代码提取了行内容中的卫星标识符,并根据卫星标识符获取卫星编号和导航系统类型。然后,根据导航系统类型进行一些条件判断,例如判断导航系统是否被禁用,从而决定是否继续处理该行。
time=str2time(buff(8:34));%time=str2time(buff(8:34));:将 buff 数组中第 8 到第 34 个元素(即行内容的第 8 到第 34 个字符)作为字符串参数,调用函数 str2time 进行解析和转换,并将转换后的时间结果赋值给变量 time。该行代码的目的是将字符子串解析为时间值。
clk=str2double(buff(41:59));%clk=str2double(buff(41:59));:将 buff 数组中第 41 到第 59 个元素(即行内容的第 41 到第 59 个字符)作为字符串参数,调用函数 str2double 进行解析和转换,并将转换后的数值结果赋值给变量 clk。该行代码的目的是将字符子串解析为双精度浮点数。
std=str2double(buff(61:79));%std=str2double(buff(61:79));:将 buff 数组中第 61 到第 79 个元素(即行内容的第 61 到第 79 个字符)作为字符串参数,调用函数 str2double 进行解析和转换,并将转换后的数值结果赋值给变量 std。该行代码的目的是将字符子串解析为双精度浮点数。
if isnan(clk),clk=0;end
if isnan(std),std=0;end
%这段代码的作用是在变量 clk 和 std 的值为 NaN 时,将它们置为 0,以避免计算或使用这些值时出现错误。NaN 是一种特殊的浮点数值,表示不是一个数字。
if nav.nc<=0%这段代码是用于更新卫星钟差信息的逻辑判断和操作。下面是每行代码的含义:
nav.nc=nav.nc+1;%if nav.nc<=0:如果变量 nav.nc 小于或等于 0,即当前没有卫星钟差信息,执行下面的操作。
%nav.nc=nav.nc+1;:将变量 nav.nc 增加 1,表示有一个卫星钟差信息。
nav.pclk(nav.nc).time=time;%nav.pclk(nav.nc).time=time;:将当前时间 time 赋值给 nav.pclk 数组中的 nav.nc 索引的 time 字段。
elseif abs(timediff(time,nav.pclk(nav.nc).time))>1e-9
if nav.nc+1>size(nav.pclk,1)
nav.pclk(nav.nc+1:nav.nc+NMAX)=repmat(gls.pclk,NMAX,1);
end
%%elseif abs(timediff(time,nav.pclk(nav.nc).time))>1e-9:如果当前时间 time 与 nav.pclk 数组中的 nav.nc 索引的 time 字段之间的时间差值大于 1e-9 秒(即差值超过 1 纳秒),执行下面的操作。
%if nav.nc+1>size(nav.pclk,1):如果 nav.nc+1 大于 nav.pclk 数组的行数(即超过数组的大小),执行下面的操作。
%nav.pclk(nav.nc+1:nav.nc+NMAX)=repmat(gls.pclk,NMAX,1);:将 gls.pclk 数组重复 NMAX 次,并赋值给 nav.pclk 数组中从 nav.nc+1 到 nav.nc+NMAX 索引的位置,即扩展 nav.pclk 数组的大小以容纳更多的卫星钟差信息。
%nav.nc=nav.nc+1;:将 nav.nc 增加 1,表示有一个新的卫星钟差信息。
%nav.pclk(nav.nc).time=time;:将当前时间 time 赋值给 nav.pclk 数组中的 nav.nc 索引的 time 字段。
nav.nc=nav.nc+1;
nav.pclk(nav.nc).time=time;%该代码段的作用是更新卫星钟差信息,在 nav.pclk 数组中加入一个新的时间值,并将 nav.nc 的值增加 1。这样就可以在 nav.pclk 数组中持续追踪和存储卫星钟差信息。
end
nav.pclk(nav.nc).clk(sat,1)=clk;
nav.pclk(nav.nc).std(sat,1)=std;%这段代码的作用是将输入的时钟值和标准差赋值给卫星钟差信息数组 nav.pclk 中对应位置的元素。通过 nav.nc 索引获取当前的卫星钟差信息,然后根据 sat 索引指定具体的卫星,将对应的时钟值 clk 和标准差 std 赋值给相应的字段。这样就更新了卫星钟差信息的时钟值和标准差。
end
if nav.nc<size(nav.pclk,1)
nav.pclk(nav.nc+1:end,:)=[];
end
%该代码段的作用是检查 nav.pclk 数组的大小是否超过 nav.nc 的值(即卫星钟差信息的数量),如果超出,则删除多余的行,以确保数组的大小与 nav.nc 的值相匹配。通过这样的操作,只保留最新的卫星钟差信息,并删除旧的无效信息。
fclose(fid);
fprintf('over\n');
return
在读取精密星历文件的时候主要处理两部分,一部分是精密星历头文件,另外一部分是主体文件。
1.4天线校正文件(.atx)
调用readatx()函数
function nav=readatx(nav,fname)
global glc gls
freq=0; pcv0=gls.pcv; pcvs=gls.pcvs;
NMAX=10000; pcvs.pcv=repmat(gls.pcv,NMAX,1);
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading atx file %s...',fname0);
fid=fopen(fname);%打开天线文件
%这段代码通过设置变量、重复数组、提取字符串以及打印消息等操作,为读取 ATX 文件做好了准备。最后,使用 fopen 函数打开了文件,并返回文件句柄 fid,可以用于后续的文件读取操作。
while ~feof(fid)
line=fgets(fid);
if size(line,2)<=60 || ~isempty(strfind(line(61:end),'COMMENT')),continue;end%找到相关的标签是“COMMENT”。
if ~isempty(strfind(line(61:end),'START OF ANTENNA'))
pcv=pcv0;%如果第61-最后一列中发现有标签“START OF ANTENNA”,将pcv0的值赋给pcv。
end
if ~isempty(strfind(line(61:end),'END OF ANTENNA'))
if pcvs.n+1>size(pcvs.pcv,1)
pcvs.pcv(pcvs.n+1:pcvs.n+NMAX)=repmat(gls.pcv,NMAX,1);
end
pcvs.pcv(pcvs.n+1)=pcv;
pcvs.n=pcvs.n+1;
continue;
end
%代码段会检查读取到的 ATX 文件中的行,如果某行包含字符串 “END OF ANTENNA”,则在特定条件下进行数据处理和更新。这些操作包括扩展数组大小、赋值操作和变量的更新。最后,使用 continue 语句跳过当前循环,继续执行下一轮循环。
%这段代码用于检测并处理读取到的 ATX 文件中的行,如果某行包含字符串 “END OF ANTENNA”,则执行一系列的操作。下面是代码的含义:
if ~isempty(strfind(line(61:end),'TYPE / SERIAL NO'))%查看标签“TYPE / SERIAL NO”
pcv.type=line(1:20);%pcv.type 是一个变量,表示天线相位中心变量的类型
pcv.code=line(21:40);%pcv.code 是一个变量,表示天线的编码(code)
prn=fix(str2num(pcv.code(2:3)));
%这行代码对字符串 pcv.code 进行处理,提取其中的第2个和第3个字符,并将其转换为数值类型,然后使用 fix 函数对结果进行向下取整。最后将得到的结果赋值给变量 prn。
%fix() 是一个MATLAB(或类似语言)中的函数,用于将输入的数值向下取整到最接近的较小整数。
if ~prn,continue;end
if strcmp(pcv.code(4:11),' ')%pcv.code(4:11):这个索引操作提取了 pcv.code 字符串的第4个到第11个字符,即一个8个字符长的子字符串。' ':这个字符串是由8个空格组成的,用于与 pcv.code(4:11) 的子字符串进行比较。
pcv.sat=satid2no(pcv.code);%如果 pcv.code(4:11) 的子字符串与8个空格相等,则调用 satid2no() 函数,并将 pcv.code 作为参数传入。satid2no() 函数可能用于根据 pcv.code 获取卫星的编号或其他相关信息,并将其赋值给 pcv.sat 变量
end
elseif ~isempty(strfind(line(61:end),'VALID FROM'))%有效起始日期(值得是)
time=str2time(line(1:43));%将含有“有效起始日期”标签的第一列到43列存储到time中。
pcv.ts=time; %将time值赋给pcv.ts
continue;
elseif ~isempty(strfind(line(61:end),'VALID UNTIL'))%有效截止日期
time=str2time(line(1:43));%将含有“有效截止日期”标签的第一列到43列存储到time中。
pcv.te=time; %将time值赋给pcv.te
continue;
elseif ~isempty(strfind(line(61:end),'DAZI'))%“DAZI” 代表着差分方位角(Differential Azimuth)。方位角是指某一点相对于参考点的方向角度。在差分GPS中,使用DAZI来表示接收机与参考站之间的方位角差异。
pcv.dazi=str2num(line(3:8)); %#ok将第3-8列的字符以str2num的形式赋值给pcv.dazi
continue;
elseif ~isempty(strfind(line(61:end),'ZEN1 / ZEN2 / DZEN'))%“ZEN1” 和 “ZEN2” 代表测站天顶角(Zenith Angle),即接收机天线指向接收到信号的角度。这些参数用于估计接收机的卫星可见性和观测品质,对于实时运动定位和静态测量都很重要。
%“DZEN” 指的是测站天顶角的变化。在RTKLIB的解算过程中,由于卫星位置的变化或者天线运动引起的变化,接收机的测站天顶角也会发生变化。“DZEN” 用于表示测站天顶角的变化量。
pcv.zen1=str2num(line(3:8)); %#ok第3-8列的值以str2num的形式赋值给pcv.zen1
pcv.zen2=str2num(line(9:14)); %#ok第9-14列的值以str2num的形式赋值给pcv.zen2
pcv.dzen=str2num(line(15:20)); %#ok第15-20列的值以str2num的形式赋值给pcv.zen3
continue;
elseif ~isempty(strfind(line(61:end),'START OF FREQUENCY'))%表示的是开始的频率
%它检查文本行中是否包含"START OF FREQUENCY"的标识,如果是的话,就表示接下来要处理的是频率信息。
[f,count]=sscanf(line(5:6),'%d');
if count<1,continue;end
[csys,count]=sscanf(line(4),'%c');%“csys” 是一个表示坐标系统的变量
if count<1,continue;end
if csys=='G'
%如果坐标系统是 ‘G’,则将频率值直接赋给 freq。
freq=f;
elseif csys=='R'
%如果坐标系统是 ‘R’,则将频率值加上glc.NFREQ(频率数量)后赋给 freq。
freq=f+glc.NFREQ;
elseif csys=='E'%主要是根据频率的特性进行计算
%如果坐标系统是 ‘E’,根据频率值的不同情况进行赋值:
if f==1,freq=f+2*glc.NFREQ;
%如果频率值为 1,将其加上 2 乘以 glc.NFREQ 后赋给 freq。
elseif f==5,freq=2+2*glc.NFREQ;
%如果频率值为 5,将其加上 2 乘以 glc.NFREQ 后再加 2 后赋给 freq。
elseif f==6,freq=3+2*glc.NFREQ;
%如果频率值为 6,将其加上 2 乘以 glc.NFREQ 后再加 3 后赋给 freq。
else ,freq=0;
%其他情况下,将 freq 赋为 0。
end
elseif csys=='C'
%如果坐标系统是 ‘C’,则将频率值加上 3 乘以 glc.NFREQ 后赋给 freq。
freq=f+3*glc.NFREQ;
elseif csys=='J'
%如果坐标系统是 ‘J’,根据频率值的不同情况进行赋值:
if f<5,freq=f+4*glc.NFREQ;
%如果频率值小于 5,将其加上 4 乘以 glc.NFREQ 后赋给 freq。
elseif f==5,freq=3+4*glc.NFREQ;
%如果频率值为 5,将其加上 4 乘以 glc.NFREQ 后再加 3 后赋给 freq。
else ,freq=0;
%其他情况下,将 freq 赋为 0。
end
else
freq=0;
%如果坐标系统不属于上述情况,则将 freq 赋为 0。
end%最终,根据不同的坐标系统和频率值,通过计算得到对应的频率索引值 freq。
elseif ~isempty(strfind(line(61:end),'END OF FREQUENCY'))%表示的是结束的频率
freq=0; %此时将频率为0
elseif ~isempty(strfind(line(61:end),'NORTH / EAST / UP'))%北、东与垂直方向
%if freq<1||freq>glc.NFREQ,continue;end
[neu,count]=decodef(line,3);%调用decodef()函数
if count<3,continue;end
if freq<1,continue;end
if pcv.sat~=0%如果卫星不为0的时候
pcv.off(freq,:)=neu;%将neu值赋给pcv.off()
else
pcv.off(freq,:)=[neu(2),neu(1),neu(3)];%其他的话将neu(2)、neu(1)与neu(3)的值赋给pcv.off()
end
elseif ~isempty(strfind(line,'NOAZI'))%在句子里面查找NOAZI,找到的话执行下面的操作
if freq<1,continue;end
dd=(pcv.zen2-pcv.zen1)/pcv.dzen+1;%其中,pcv.zen2 和 pcv.zen1 是探测器的 Zenith 角度,pcv.dzen 表示 Zenith 角度的间隔。
if dd~=myRound(dd)||dd<=1
%fprintf('error');
continue;
end
if pcv.dazi==0
[val,count]=decodef(line(9:end),fix(dd));%调用decodef()函数
%decodef() 是一个用于解码字符串并返回浮点数数组的函数。
if count<=0
%fprintf('error');
continue
elseif count~=fix(dd)
%fprintf('error');
continue
end
pcv.var(freq,1:count)=val;
else
id=fix(360/pcv.dazi)+1;
%计算 id 的值,该值为 360/pcv.dazi 的取整值加 1。
for i=1:id
line=fgets(fid);
[val,j]=decodef(line(9:end),fix(dd));
if j<=0,continue
elseif j~=fix(dd),continue%fix()函数用于取整数部分
end
pcv.var(freq,(i-1)*fix(dd)+1:(i-1)*fix(dd)+j)=val;
end%复制数组:可以使用这个语句将一个数组的值复制到另一个数组的特定位置。
end%pcv.var(1, 6:9) = val; % 将数组 val 的值复制到 pcv.var 的第一个行中的索引 6 到 9 的元素
end%更新部分数据:可以使用这个语句将特定位置的元素更新为新的值。
%pcv.var(2, 3:6) = val; % 将数组 val 的值更新到 pcv.var 的第二行中的索引 3 到 6 的元素
end
fclose(fid);
fprintf('over\n');
if pcvs.n<size(pcvs.pcv,1)%如果 pcvs.n 的值小于 pcvs.pcv 数组的行数(即数组的大小),则执行条件判断内的语句。
pcvs.pcv(pcvs.n+1:end,:)=[];%在条件判断内部,pcvs.pcv(pcvs.n+1:end,:)=[]; 用于删除 pcvs.pcv 数组中从 pcvs.n+1 行到末尾的所有行。
end
nav.ant_para=pcvs;%pcvs的结果赋值给nav.ant_para里面
return
1.5差分码偏差文件(.DCB)
调用readdcb()函数
function nav=readdcb(nav,obs,fname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%only read dcb for GPS or GLONASS(CODE procuct)
%convert CODE product to consistent with CAS product
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
type=0; DCBPAIR=glc.DCBPAIR;%导入dcbpair的参数
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading DCB file %s...',fname0);%读到该文件的时候,在命令行输出“Info:reading DCB file %s...”
fid=fopen(fname);
%nav.cbias
%双差码偏差参数通常用于补偿不同频率的伪距观测之间的差异。将这个差异转化为以米为单位的数值,可以更准确地对卫星导航系统进行校正。
while ~feof(fid)
line=fgets(fid);
if ~isempty(strfind(line,'DIFFERENTIAL (P1-P2) CODE BIASES')),type=1;
elseif ~isempty(strfind(line,'DIFFERENTIAL (P1-C1) CODE BIASES')),type=2;
elseif ~isempty(strfind(line,'DIFFERENTIAL (P2-C2) CODE BIASES')),type=3;
end%将三个文件p1p2、p1c1与p2c2文件,将三种文件对应不同的类型,第一种类型是type=1;第二种类型是type=2,第三种类型是type=3.
if type==0,continue;end
if isempty(line),continue;end
if line(1)~='G'&&line(1)~='R',continue;end
buff=strsplit(line);%strsplit() 是一个字符串操作函数,通常用于将一个字符串分割成多个子字符串,并将这些子字符串存储在一个数组中。
%buff 一般指的是用于存储或处理数据的缓冲区或缓冲器。
cbias=str2double(line(27:35));%将第27-35列的值以str2double的形式赋值给cbias(abias表示的是载波偏差)。
%对于接收机来说
if strcmp(buff(1),'G')||strcmp(buff(1),'R')%如果存储的系统的类型为G或者R的时候。
if ~isfield(obs,'sta'),continue;end
if strcmp(buff(1),'G'),j=1;else,j=2;end
if strcmpi(obs.sta.name,char(buff(2)))
nav.rbias(j,type)=cbias*1E-9*glc.CLIGHT;%nav.rbias(j,type) = cbias * 1E-9 * glc.CLIGHT:如果上述条件满足,将计算得到的值 cbias 乘以 1E-9(表示1倍 1−910 −9,即纳秒转为秒的缩放因子),再乘以 glc.CLIGHT(光速)并赋值给 nav.rbias(j,type)。
end%rbias 是一个术语,表示接收机的接收机钟差(Receiver Clock Bias)。
%这段代码片段的功能是,如果 obs 中有 sta 字段,然后根据 buff 字符串中的第一个字符确定变量 j 的值,最后通过比较字符串 obs.sta.name 与 buff(2) 来决定是否给 nav.rbias(j,type) 赋值。
else
%对于卫星来说
sat=satid2no(char(buff(1)));%调用的是satid2no()函数,该函数用于将卫星的标识符(satid)转换为卫星编号(sats)。
if sat~=0
[sys,~]=satsys(sat);%调用的是satsys()函数,输入参数主要包括sat。
if type==1%如果是第一类型,p1p2的文件。
dcb_pair='C1W-C2W'; %P1-P2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1W-C2W’。
%在这个字符串中,‘C1W’ 表示的是 P1 码,‘C2W’ 表示的是 P2 码。P1 码和 P2 码都是 GPS 卫星导航的伪距测量中的两个频率。
if sys==glc.SYS_GPS%如果系统是GPS系统的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统为GPS的卫星的双差码偏差。
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;%对于给定的卫星编号sat和双差码偏差参数索引i,将卫星卫星导航系统的dcbs数组中的元素设置为cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT。
elseif sys==glc.SYS_GLO%如果系统是GLO系统话
dcb_pair='C1P-C2P'; %P1-P2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1P-C2P’。
for i=1:glc.MAXDCBPAIR%%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统GLO的卫星的双差码偏差
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;%;%对于给定的卫星编号sat和双差码偏差参数索引i,将卫星卫星导航系统的dcbs数组中的元素设置为cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT。
end
elseif type==2%如果是p1c1类型的文件的话
dcb_pair='C1C-C1W'; %C1-P1
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1C-C1W’。
if sys==glc.SYS_GPS%如果系统是GPS系统的话
for i=1:glc.MAXDCBPAIR%%%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统为GPS的卫星的双差码偏差。
nav.cbias(sat,i)=-cbias*1E-9*glc.CLIGHT;
%nav.cbias(sat, i) = -cbias * 1E-9 * glc.CLIGHT; 这一行代码的作用是将导航系统为GPS的卫星的双差码偏差的负值赋值给nav.cbias(sat, i)。
elseif sys==glc.SYS_GLO%如果系统是GLO系统的话
dcb_pair='C1C-C1P'; %C1-P1
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1C-C1P’。
for i=1:glc.MAXDCBPAIR%%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统GLO的卫星的双差码偏差
nav.cbias(sat,i)=-cbias*1E-9*glc.CLIGHT; %note the minus sign !
%nav.cbias(sat, i) = -cbias * 1E-9 * glc.CLIGHT; 这一行代码的作用是将导航系统为GLO的卫星的双差码偏差的负值赋值给nav.cbias(sat, i)。
end
elseif type==3%如果是p2c2类型的文件的话
dcb_pair='C2W-C2L'; %P2-C2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C2W-C2L’。
if sys==glc.SYS_GPS%如果卫星系统为GPS的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算卫星系统GPS的双差码偏差
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;
%在这个代码中,将导航系统为GPS的卫星的双差码偏差赋值给nav.cbias(sat, i)。具体而言,将cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,然后将结果赋给nav.cbias(sat, i)。
elseif sys==glc.SYS_GLO%如果卫星系统为GLO的话
dcb_pair='C2C-C2P'; %C2-P2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C2C-C2P’。
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统GLO的卫星的双差码偏差
nav.cbias(sat,i)=-cbias*1E-9*glc.CLIGHT; %note the minus sign !
%将cbias(双差码偏差值)乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,并将结果赋值为nav.cbias(sat, i)的相反数。这样做的目的可能是根据特定的需求或公式要求,将双差码偏差的值取负。
end
end
end
end
end
fclose(fid);
fprintf('over\n');
return