基于Matlab的TEQC绘图程序代码.doc

%在 matlab下新建一个m文件,将以下代码直接拷贝进去,即可执行。
％需要一个TEQC生成的plot文件作为参数
function out=teqcplot3(files）;
％读取TEQC生成的Plot文件，绘制数据图表，支持Copmact、Compact2、Compact3格式
%选取一个TEQC的Plot文件
％ 格式说明
％ ＊。sn1 载波 L1的信噪比 Signal to noise ratio （S/N）
％ *。sn2 载波 L2的信噪比 Signal to noise ratio （S/N） Carrier L2
％ ＊。iod ＊。d12 ＊。d21 电离层延迟观测值变化率(米/秒） Derivative of ionospheric delay observable （m/s）
％ ＊。ion ＊。i12 *.i21 电离层延迟观测值(米）Ionospheric delay observable (m）
％ *。mp1 *.m12 载波L1的多路径误差 Multipath Carrier L1
% ＊。mp2 ＊。m21 载波L2的多路径误差 Multipath Carrier L2
％ ＊.azi 卫星方位角(°） Satellite azimuthal data (degrees）
% ＊.ele 卫星高度角（°) Satellite elevation data （degrees）
if nargin==0
[filen,path］=uigetfile（’＊。sn1；＊.sn2；＊.iod；＊。ion；*.mp1;＊。mp2；＊.azi;*.ele;*。i12；＊i21;＊.m12；m21；＊。d12；＊d21’,.。.
’请选择TEQC报告文件：')；
else
［path，filen，ext]=fileparts（files)；
path=[path ’\’］;
filen=｛［filen ext］};
end
%读取这个文件
file=char（filen）；
%按行读取文件至数组A
［A]=importdata([path file］,’\t’)；
%定义SAT，存放卫星数据
%GPS有32颗卫星，存放序号1-32，
%GLONASS有32颗卫星，存放序号33—64,
％BEIDOU有35颗卫星，存放序号65—99
SAT（1:length（A），1:99)=NaN;
％sats（1:length（A）,1：99)=NaN；
％存放采样时间，单位秒
tsec（1：length(A))=NaN；
％读取文件的第一行
filelx=A{1}；
%判断是哪种格式
switch filelx
case ’COMPACT'
%读取数据采样间隔
t_samp=char(A（3)）;
％读取开始时间
mjl=char（A（4)）；
%读取数据采样间隔
T_SAMP=str2num（t_samp（max（find(t_samp==’ '））：end）)；
％读取数据采样开始时间
MJL_START=str2num（mjl(max（find(mjl==’ ’)）：end）);
％转成时间序列数字,date serial number，从0000年1月1日0时0分0秒开始计算的十进制天数
MJD_START=MJL_START+678941。999999741；
%i为行号
n=1;i=5;
case ’COMPACT2’
%读取数据采用间隔
t_samp=char(A（2))；
%读取开始时间
mjl=char（A（3））；
％读取数据采样间隔
T_SAMP=str2num（t_samp（max（find（t_samp==' '））:end）);
％读取数据采样开始时间
MJL_START=str2num(mjl（max（find（mjl==’ '））:end））；
％转成时间序列数字，date serial number，从0000年1月1日0时0分0秒开始计算的十进制天数
MJD_START=MJL_START+678941.999999741；
n=1；i=4;
case 'COMPACT3’
％读取开始时间
t_start=char(A（2））；
t_start = deblank（t_start）；
s = splitstr（t_start,'*＊’，6）；
%t_start_time=［char（s｛2｝) ’年' char(s｛3｝) ’月’ char（s{4}) ’日’ char(s{5｝） '时' char（s｛6}） '分' num2str（str2num（char（s{7｝）），'%02d') ’秒'] ;
%获取采样的开始时间,2013，12，7，03，05,55
t_s_time=［str2num（char（s{2}））， str2num(char（s{3｝）)，str2num(char（s｛4｝）)，str2num(char（s｛5｝)),str2num（char(s{6｝)）,str2num（char（s｛7｝))］;
n=1；i=3;
otherwise
disp（’数据格式存在问题')；
return
end
％n=1;i=3;
％sats=str2num（A{i}）；
snyggfilen=strrep（filen，’_’，’-')；
％生成进度条
h=waitbar（0，［'正在读取数据,请稍等…… ’ char（snyggfilen）］）；
while i<length（A)；
waitbar(i/length(A）,h）;drawnow
%读取卫星编号数据行 采样时间 卫星数量 卫星编号 COMAPCT3
satbhstr=char(A（i)）；
s = splitstr(satbhstr，’＊＊’,32）；
％如果卫星编号数据为0,说明数据缺失，退出本函数
if strtrim（satbhstr）==’0’
disp (’数据为空!’）;
close（h）；
return
end
switch filelx
case ’COMPACT3’
%记录采样时间，单位秒
tsec（n)=str2num（s{1｝)；
if s｛2｝==’—1’
％如果为—1,使用上次的卫星编号行
satbhstr=oldsatbh;
s = splitstr（satbhstr，'*＊’，32)；
else
oldsatbh=satbhstr；
end
％获取当前卫星数量
satcount=str2num（s｛2｝）;
case ｛'COMPACT',’COMPACT2’｝
%记录采样时间，单位秒
tsec(n)=(n—1)*T_SAMP；
if s{1}=='-1’
％如果为-1,使用上次的卫星编号行
satbhstr=oldsatbh；
s = splitstr（satbhstr，'**’,32）；
else
oldsatbh=satbhstr;
end
％获取当前卫星数量
satcount=str2num（s｛1｝）;
end
％读取卫星数据行
satdatastr=char（A(i+1)）；
sdata=splitstr(satdatastr,'＊＊',32)；
for k=3：satcount+2
switch filelx
case 'COMPACT'
%如果是COMPACT格式，卫星编号前面没有字母，默认是GPS卫星,在前面增加字符G
satbh=[’G’ char(s{k—1｝)]；
case ’COMPACT2'
%如果是COMPACT2格式，卫星编号行第2个数据开始是卫星编号
satbh=char(s｛k-1}）；
case 'COMPACT3'
％如果是COMPACT3格式，卫星编号行第3个数据开始是卫星编号
satbh=char（s{k}）;
end
switch satbh(1）；
case 'G'
％如果是GPS卫星,获得与编号对应的存储序号,范围1—32
index=str2num（satbh（2:3）);
case ’R'
％如果是GLONASS卫星,获得与编号对应的存储序号,范围33—64
index=32+str2num（satbh（2:3));
case 'C'
％如果是BEIDOU卫星，获得与编号对应的存储序号,范围65—99
index=64+str2num（satbh(2:3)）;
otherwise
%其他情况，获得与编号对应的存储序号，范围1—32
index=str2num（satbh（2：3）);
end ％switch
%获得对应的卫星数据
sdatastr=char（sdata｛k—2}）;
％如果卫星数据不是数值型,用0替代
if isempty（str2num(sdatastr)）==1
sdatastr=’0。000’；
end
％将卫星数据存入SAT数组的对应位置
SAT(n,index）=str2num(sdatastr）;
end ％for
n=n+1;
%根据数据文件结构，一行是卫星编号，下一行就是对应的卫星数据
%程序一次读取2行数据
i=i+2;
end %while
%数据读取完毕,关闭进度条
close（h）;
％获取已使用的卫星编号，返回给sat_bh
sat_bh=getsatbh（SAT);
％增加最后一个结束字符
sat_bh_temp=［sat_bh’ {’End’｝］；
sat_bh_end=sat_bh_temp’；
％删去所有都为NaN值的列,返回给sat_data
sat_data=delnancol（SAT）；
％删去所有NaN值的采样时间列
t_seconds=delnancol（tsec);
%删去所有都为NaN值的行，返回sat_datas
sat_datas=delnanrow（sat_data);
％增加一个nan列
［row，col]=size(sat_datas）;
bb（1：row，1)=NaN；
sat_datas=［sat_datas bb］;
%将含有NaN值的数据替换为0，返回给s_data
%s_data=repnan20（sat_data)；
％将采样时间转成顺序日期serial date number
switch filelx
case ｛’COMPACT’,’COMPACT2’}
t_s_jd_time=MJD_START;
case ’COMPACT3'
t_s_jd_time=datenum（t_s_time)；
end
for i=1：length（t_seconds）
％将每次采样时间转换成相应的顺序日期
t_jd_time（i)=t_s_jd_time+t_seconds(i）/60/60/24；
end
％绘制图表++++++++++++++++++++++++++
%根据不同的文件类型，设置坐标轴的范围
[type,maxy,miny］=get_filetype(file）;
figure;box on；hold on
[row,col］=size(sat_datas）；
％绘制渐变彩色图
pcolor（t_jd_time’,［1：col］，sat_datas’)；
set(gcf，'renderer'，'zbuffer')；
shading flat
set（gca，’xticklabel’，[t_jd_time］）；
set（gca，’xlim’，[t_jd_time(1） t_jd_time（end)］)
%t_start_times=［t_start_h ’：’ t_start_m ’：’ t_start_s］；
dateaxis（'X’,13）
cbar（'v'，［miny maxy］，type）;
set（gca,’ylim'，［1 col+1]）
set（gca,’ytick'，［1。5：1：col+0。5］)
set(gca,'yticklabel’,sat_bh_end）；
set（gca,’fontsize’，7);
colormap(flipud（jet))；
caxis（［miny maxy］)；
％t_end_time=cal2et（t_s_time，t_seconds(end)）；
xlabel（［datestr（t_jd_time(1）) ’ ｜——-——-—— 采样间隔： ’ num2str（t_seconds（2）—t_seconds（1)) ' 秒 —-—-———-| ' datestr（t_jd_time（end）)])
ylabel（'卫星编号’）
％timestr=secs2hms（length（sat）)；
T=title（［’TEQC 报告文件： ' strrep（file，’_’，’—’）］）；
set（T,’fontsize'，8）
％out。(file（end-2:end））=sat;
％out.T_samp=T_SAMP；
out。Start=datestr（t_jd_time(1));
out。Stop=datestr(t_jd_time（end））；
％+++++++++++++++++++++++++++++++++++++++
%以下是本函数中所使用的一些相关函数
％++++++++++++++++++++++++++++++++++++++++++++++++++
％+++++++++++++++++++++++++++++++++++++++++++++++++
function s=splitstr（str，split,num);
%用指定分隔符分隔字符串函数
%str='name 0 2010-05—06 33 Qt BKB';
%'＊＊'表示特殊分隔符回车、空格、换行、制表符等
%num表示分隔字符串的个数
string=cell（num,1)；
index=1;
temp=［］;
if split==’＊＊’
％特殊分隔符
for i=1:length(str）
if isspace（str（i)）
if isempty（temp）==false
string{index}=temp;
index=index+1；
temp=［];
end
else
temp=［temp，str（i)］；
end
end
else
％常规分隔符
for i=1:length（str)
if str（i）==split
if isempty（temp）==false
string｛index}=temp；
index=index+1；
temp=［]；
end
else
temp=［temp,str(i）］；
end
end
end
if isempty(temp）==false
string{index}=temp;
end
s=string；
%+++++++++++++++++++++++++++++++++++
function B=delnancol（A);
％删除矩阵中所有值都是nan的列
［x y]=size（A)；
for i=y:—1：1
if all（isnan(A（：，i)))==1
A（：，i）=［］；
end
end
B=A；
%+++++++++++++++++++++++++++++++++++
function B=delnanrow（A);
%删除矩阵中所有值都是nan的行
［x y]=size（A)；
for i=x:-1：1
if all（isnan(A（i，：)）)==1
A（i，:）=［］;
end
end
B=A;
%+++++++++++++++++++++++++++++++++++
function B=repnan20（data）;
%将矩阵中的NAN值用0替代
［datas,features]=size(data);
for k=1：features
for i=1：datas
if isnan（data（i，k))==1
data（i，k)=0;
end
end
end
B=data；
%++++++++++++++++++++++++++++++++
function satbh=getsatbh(satdata);
％获取卫星数据中的卫星编号，返回一个字符串数组
s=cell(32，1）;
temp=[］；
index=1；
for i=1:99
if all(isnan(satdata（：，i)）)~=1
switch i〉32
case 0
temp=[’G'，num2str(i，’％02d'）］;
case 1
if i〉64
temp=['C’，num2str（i—64，'％02d’)］；
else
temp=［’R’，num2str（i—32，’％02d’)];
end
end
s｛index}=temp；
index=index+1；
end
end
if index〉1
satbh=s（1：index—1）；
else
satbh=[］；
end
％+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
％ FUNCTION: GET FILETYPE INFO +++++++++++++++++++++++++++++++++++++
function ［out，maxy,miny］=get_filetype（teqfile）；
%根据文件类型，设置坐标轴的范围
［path，name，ext,ver］=fileparts(teqfile）；
switch ext
case ’.sn1'
％out='Signal to noise ratio S/N L1’；
out=’L1载波上的信噪比(dBHz）’；
maxy=60;
miny=15；
case ’。sn2’
%out=’Signal to noise ratio S/N L2’；
out='L2载波上的信噪比（dBHz）’；
maxy=60；
miny=15；
case ｛’。mp1'，’。m12'｝
％out=’Multipath L1’;
out=’L1多路径观测值'；
maxy=1；
miny=—1;
case ｛'。mp2'，’。m21'｝
％out='Multipath L2’；
out=’L2多路径观测值’;
maxy=1;
miny=-1；
case ｛'。iod’,’。d12’,'d21’}
%out='Derivative of ionospheric delay observable (m/s)’；
out=’电离层延迟观测值变化率（米/秒)’；
maxy=1；
miny=—1;
case {’.ion’，’。i12'，’i21’｝
maxy=2；
miny=—2；
%out=’Ionospheric delay observable (m）’；
out=’电离层延迟观测值(米）'；
case ’.ele'
maxy=90；
miny=0；
%out=’Satellite elevation data’；
out=’卫星高度角(°）'；
case ’。azi’
maxy=180；
miny=—180;
％out=’Satellite azimuthal data’；
out=’卫星方位角（°）’；
otherwise
disp（'Somethings wrong。。!’）
end
％ FUNCTION: PLACE A MODIFIED COLORBAR ++++++++++++++++++++++++++++++
function CB=cbar(loc，range,label）;
％ 。.。。。..。。。..。.。。.。。。..。..。..。.。。。。.。。。。。。..。.。。。.。.。。.。。.。。..
% CB = cbar（loc,range,label）
% places a colorbar at:
％ loc = 'v' in vertical or ’h’ in horizontal
% position in current figure scaled between：
% range = ［min max］ with a:
% label = ’string’。
％
％ fontsize is reduced to 10 and width of bar is half default.
％
% Example: ［X，Y，Z］=peaks(25）；
％ range=[min（min（Z）） max（max（Z））］;
％ pcolor(X,Y,Z）;
% cbar(’v’,range,’Elevation （m)'）
% 。。.。。。。.。..。.....。。。。。..。。。。。。.。...。。。。。.。。。..。。。.。。..。。。.。。。
caxis(［range（1) range(2）]）;
switch loc
case ’v’
CB=colorbar（’vertical’）；
set（CB，'ylim’,［range（1） range(2)]）；
POS=get(CB，’position’）;
set(CB，’position'，[POS（1) POS（2） 0.03 POS（4）]）；
set(CB,’fontsize',8）;
set（get（CB，’ylabel'），'string'，label）;
set（CB，’box’,’on’）
case ’h’
CB=colorbar（'horizontal’）；
set(CB，’xlim',［range（1） range(2)］）；
POS=get（CB,’position’）；
set（CB，’position'，[POS(1） POS（2） POS（3) 0。03])；
set（CB,'fontsize’，8）；
set（get（CB，’xlabel’）,’string',label）
end
％++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function jd=cal2jd（cal）
％ cal2jd1 将公历年月日时分秒转换到儒略日。
% jd=cal2jd(cal） 返回儒略日
％ cal：1x6矩阵，6列分别为年月日时分秒。构造cal时可以省略末尾的0
％
％ 公元1582年10月4日24：00点之前使用儒略历,公元1582年10月15日00：00点之后使用公历
if length（cal） 〈 6
cal（6)=0;
end
year=cal(1);
month=cal（2）;
day=cal（3）+（cal(4)＊3600+cal（5)*60+cal(6）)/86400；
y = year + 4800; ％4801 B.C。 is a century year and also a leap year。
if（ year 〈 0 )
y =y+ 1； ％ Please note that there is no year 0 A.D。
end
m=month;
if( m 〈= 2 ) ％ January and February come after December。
m = m+12;
y = y — 1；
end
e=floor（30.6 * （m+1））；
a=floor（y/100)； % number of centuries
% 教皇格雷戈里十三世于1582年2月24日以教皇训令颁布，将1582年10月5日至14抹掉。1582年10月4日过完后第二天是10月15日
if( year <1582 ）|(year==1582&month〈10）｜（year==1582& month==10 ＆day〈15)
b = —38；
else
b = floor((a/4） — a); ％ number of century years that are not leap years
end
c=floor(365.25＊ y）； ％ Julian calendar years and leap years
jd= b + c + e + day — 32167。5；
％+++++++++++++++++++++++++++++++++++++++++++++++++
function cal=jd2cal（J）
% 从儒略日计算公历年月日时分秒
% cal=jd2cal(J）
％ 返回的cal是1x6矩阵，6列分别为年月日时分秒
%
％ 公元1582年10月4日24：00点之前使用儒略历，公元1582年10月15日00:00点之后使用公历
if （J 〈 1721423。5） % 公元1月1日0时
BC = 1;
else
BC = 0；
end
％start from Julian March 1, 4801 B。C。
if（ J 〈 2299160。5 ) ％ before 1582。10。4。 24：00 is Julian calender
j0=floor(J+0。5）；
dd=J+0。5-j0；
else % after 1582.10。15. 00:00 is Gregorian calender
％number of certury years that are not leap year
n1=floor((J—2342031.5）/36524。25/4）+1; %1700。3。1。0
n2=floor(（J—2378555。5)/36524。25/4）+1; ％1800.3。1。0
n3=floor（（J-2415079。5)/36524.25/4)+1； %1900.3.1.0
j0=n1+n2+n3+J+10;
dd=j0+0.5—floor（j0+0。5）；
j0=floor（j0+0.5）；
end
j0=j0+32083;
year0=ceil(j0/365。25）—1;
year=year0—4800；
day=j0—floor（year0＊365。25）；
month=floor（（day—0。6)/30。6)+3;
day=day—round（（month—3)*30.6）；
if month>12
month=month—12；
year=year+1；
end
year=year—BC;
sec=round（dd＊86400)；
hour=floor(sec/3600）;
sec=sec-hour＊3600；
min=floor（sec/60）;
sec=sec-min＊60；
cal=[year, month, day， hour, min， sec］;
％+++++++++++++++++++++++++++++++++++++++++++++++++
function mjd=cal2mjd(cal）
% cal2mjd 将公历年月日时分秒转换到简化儒略日。
% mjd=cal2mjd（cal） 返回简化儒略日
% cal：1x6矩阵,6列分别为年月日时分秒.构造cal时可以省略末尾的0
jd = cal2jd（cal）；
mjd = jd — 2400000。5；
% FUNCTION： SECONDS TO HOURS， MINUTES and SECONDS ++++++++++++++++
function timestr=secs2hms（SECS）
HOURS=SECS/60/60；
hours=floor(HOURS)；
MINUTES=(HOURS—hours）＊60；
minutes=floor(MINUTES）；
seconds=（MINUTES—minutes）*60;
HH=num2str（hours）；
MM=num2str(minutes）；
SS=num2str（seconds）；
if seconds<10；
SS=['0' num2str(SS）］；
else
SS=num2str（SS)；
end
if minutes〈10;
MM=［'0’ num2str(MM）］；
else
MM=num2str(MM）；
end
timestr=[HH ’：’ MM ’:' SS];
％+++++++++++++++++++++++++++++++++++++
function timestr=cal2et(t_start，secs)
%根据开始时间,时长（秒），计算结束时间
％t_start 1×6矩阵，年，月，日,时,分,秒
％secs 时长,单位秒
t_s_days=t_start（4)/24+t_start(5）/60/24+t_start（6）/60/60/24；
t_e_days=t_s_days+secs/60/60/24；
HOURS=t_e_days＊24；
hours=floor（HOURS）;
MINUTES=(HOURS-hours)*60;
minutes=floor（MINUTES);
seconds=round（(MINUTES—minutes）*60）;
HH=num2str（hours,’%02d’）;
MM