基于Matlab的TEQC绘图程序代码

%在matlab下新建一个m文件，将以下代码干脆拷贝进去，即可执行。%须要一个TEQC生成的plot文件作为参数functionout=teqcplot3(files);%读取TEQC生成的Plot文件，绘制数据图表，支持Copmact、Compact2、Compact3格式%选取一个TEQC的Plot文件%格式说明%*.sn1载波L1的信噪比Signaltonoiseratio(S/N)%*.sn2载波L2的信噪比Signaltonoiseratio(S/N)CarrierL2%*.iod*.d12*.d21电离层延迟观测值改变率(米/秒)Derivativeofionosphericdelayobservable(m/s)%*.ion*.i12*.i21电离层延迟观测值(米)Ionosphericdelayobservable(m)%*.mp1*.m12载波L1的多路径误差MultipathCarrierL1%*.mp2*.m21载波L2的多路径误差MultipathCarrierL2%*.azi卫星方位角（°）Satelliteazimuthaldata(degrees)%*.ele卫星高度角（°）Satelliteelevationdata(degrees)ifnargin==0[]=uigetfile('*.sn1;*.sn2;*.iod;*.ion;*.mp1;*.mp2;*.azi;*.ele;*.i12;*i21;*.m12;m21;*.d12;*d21',...'请选择TEQC报告文件：');else[path,]=(files);path=[path'\'];filen={[filenext]};end%读取这个文件(filen);%按行读取文件至数组A[A]=importdata([pathfile],'\t');%定义SAT，存放卫星数据%GPS有32颗卫星，存放序号1-32，%GLONASS有32颗卫星，存放序号33-64，%BEIDOU有35颗卫星，存放序号65-99SAT(1:length(A),1:99)=NaN;%sats(1:length(A),1:99)=NaN;%存放采样时间，单位秒tsec(1:length(A))=NaN;%读取文件的第一行{1};%推断是哪种格式switchfilelxcase'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));%转成时间序列数字，dateserialnumber，从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));%转成时间序列数字，dateserialnumber，从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,55t_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;otherwisedisp('数据格式存在问题');returnend%n=1;i=3;%sats=str2num(A{i});snygg(filen,'_','-');%生成进度条h=waitbar(0,['正在读取数据，请稍等……'char(snyggfilen)]);whilei<length(A);waitbar(i/length(A),h);drawnow%读取卫星编号数据行采样时间卫星数量卫星编号COMAPCT3satbhstr=char(A(i));s=splitstr(satbhstr,'**',32);%假如卫星编号数据为0，说明数据缺失，退出本函数ifstrtrim(satbhstr)=='0'disp('数据为空！');close(h);returnendswitchfilelxcase'COMPACT3'%记录采样时间，单位秒tsec(n)=str2num(s{1});ifs{2}=='-1'%假如为-1，运用上次的卫星编号行satbhstr=oldsatbh;s=splitstr(satbhstr,'**',32);elseoldsatbh=satbhstr;end%获得当前卫星数量satcount=str2num(s{2});case{'COMPACT','COMPACT2'}%记录采样时间，单位秒tsec(n)=(n-1)*T_SAMP;ifs{1}=='-1'%假如为-1，运用上次的卫星编号行satbhstr=oldsatbh;s=splitstr(satbhstr,'**',32);elseoldsatbh=satbhstr;end%获得当前卫星数量satcount=str2num(s{1});end%读取卫星数据行satdatastr=char(A(i+1));sdata=splitstr(satdatastr,'**',32);fork=3:satcount+2switchfilelxcase'COMPACT'%假如是COMPACT格式，卫星编号前面没有字母，默认是GPS卫星，在前面增加字符Gsatbh=['G'char(s{k-1})];case'COMPACT2'%假如是COMPACT2格式，卫星编号行第2个数据起先是卫星编号satbh=char(s{k-1});case'COMPACT3'%假如是COMPACT3格式，卫星编号行第3个数据起先是卫星编号satbh=char(s{k});endswitchsatbh(1);case'G'%假如是GPS卫星，获得与编号对应的存储序号，范围1-32index=str2num(satbh(2:3));case'R'%假如是GLONASS卫星，获得与编号对应的存储序号，范围33-64index=32+str2num(satbh(2:3));case'C'%假如是BEIDOU卫星，获得与编号对应的存储序号，范围65-99index=64+str2num(satbh(2:3));otherwise%其他状况，获得与编号对应的存储序号，范围1-32index=str2num(satbh(2:3));end%switch%获得对应的卫星数据sdatastr=char(sdata{k-2});%假如卫星数据不是数值型，用0替代ifisempty(str2num(sdatastr))==1sdatastr='0.000';end%将卫星数据存入SAT数组的对应位置SAT(n,index)=str2num(sdatastr);end%forn=n+1;%依据数据文件结构，一行是卫星编号，下一行就是对应的卫星数据%程序一次读取2行数据i=i+2;end%while%数据读取完毕，关闭进度条close(h);%获得已运用的卫星编号，返回给sat_bhsat_bh=getsatbh(SAT);%增加最终一个结束字符sat_bh_temp=[sat_bh'{'End'}];sat_bh_end=sat_bh_temp';%删去全部都为NaN值的列，返回给sat_datasat_data=delnancol(SAT);%删去全部NaN值的采样时间列t_seconds=delnancol(tsec);%删去全部都为NaN值的行，返回sat_datassat_datas=delnanrow(sat_data);%增加一个nan列[row,col]=size(sat_datas);bb(1:row,1)=NaN;sat_datas=[sat_datasbb];%将含有NaN值的数据替换为0，返回给s_data%s_data=repnan20(sat_data);%将采样时间转成依次日期serialdatenumberswitchfilelxcase{'COMPACT','COMPACT2'}t_s_jd_time=MJD_START;case'COMPACT3't_s_jd_time=datenum(t_s_time);endfori=1:length(t_seconds)%将每次采样时间转换成相应的依次日期t_jd_time(i)=t_s_jd_time+t_seconds(i)/60/60/24;end%绘制图表++++++++++++++++++++++++++%依据不同的文件类型，设置坐标轴的范围[type,maxy,miny]=get_(file);figure;boxon;holdon[row,col]=size(sat_datas);%绘制渐变彩色图pcolor(t_jd_time',[1:col],sat_datas');set(gcf,'renderer','zbuffer');shadingflatset(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',[minymaxy],type);set(gca,'ylim',[1col+1])set(gca,'ytick',[1.5:1:col+0.5])set(gca,'yticklabel',sat_bh_end);set(gca,'fontsize',7);colormap(flipud(jet));caxis([minymaxy]);%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.())=sat;%out.T_samp=T_SAMP;out.Start=datestr(t_jd_time(1));out.Stop=datestr(t_jd_time(end));%+++++++++++++++++++++++++++++++++++++++%以下是本函数中所运用的一些相关函数%++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++++++++++++++++++++++++functions=splitstr(str,split,num);%用指定分隔符分隔字符串函数%str='name02010-05-0633QtBKB';%'**'表示特别分隔符回车、空格、换行、制表符等%num表示分隔字符串的个数string=cell(num,1);index=1;temp=[];ifsplit=='**'%特别分隔符fori=1:length(str)ifisspace(str(i))ifisempty(temp)==falsestring{index}=temp;index=index+1;temp=[];endelsetemp=[temp,str(i)];endendelse%常规分隔符fori=1:length(str)ifstr(i)==splitifisempty(temp)==falsestring{index}=temp;index=index+1;temp=[];endelsetemp=[temp,str(i)];endendendifisempty(temp)==falsestring{index}=temp;ends=string;%+++++++++++++++++++++++++++++++++++functionB=delnancol(A);%删除矩阵中全部值都是nan的列[xy]=size(A);fori=y:-1:1ifall(isnan(A(:,i)))==1A(:,i)=[];endendB=A;%+++++++++++++++++++++++++++++++++++functionB=delnanrow(A);%删除矩阵中全部值都是nan的行[xy]=size(A);fori=x:-1:1ifall(isnan(A(i,:)))==1A(i,:)=[];endendB=A;%+++++++++++++++++++++++++++++++++++functionB=repnan20(data);%将矩阵中的NAN值用0替代[datas,features]=size(data);fork=1:featuresfori=1:datasifisnan(data(i,k))==1data(i,k)=0;endendendB=data;%++++++++++++++++++++++++++++++++functionsatbh=getsatbh(satdata);%获得卫星数据中的卫星编号，返回一个字符串数组s=cell(32,1);temp=[];index=1;fori=1:99ifall(isnan(satdata(:,i)))~=1switchi>32case0temp=['G',num2str(i,'%02d')];case1ifi>64temp=['C',num2str(i-64,'%02d')];elsetemp=['R',num2str(i-32,'%02d')];endends{index}=temp;index=index+1;endendifindex>1satbh=s(1:index-1);elsesatbh=[];end%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++%FUNCTION:GETINFO+++++++++++++++++++++++++++++++++++++function[out,maxy,miny]=get_(teqfile);%依据文件类型，设置坐标轴的范围[path,name,ext,ver]=(teqfile);switchextcase'.sn1'%out='SignaltonoiseratioS/NL1';out='L1载波上的信噪比(dBHz)';maxy=60;miny=15;case'.sn2'%out='SignaltonoiseratioS/NL2';out='L2载波上的信噪比(dBHz)';maxy=60;miny=15;case{'.mp1','.m12'}%out='MultipathL1';out='L1多路径观测值';maxy=1;miny=-1;case{'.mp2','.m21'}%out='MultipathL2';out='L2多路径观测值';maxy=1;miny=-1;case{'.iod','.d12','d21'}%out='Derivativeofionosphericdelayobservable(m/s)';out='电离层延迟观测值改变率（米/秒）';maxy=1;miny=-1;case{'.ion','.i12','i21'}maxy=2;miny=-2;%out='Ionosphericdelayobservable(m)';out='电离层延迟观测值(米)';case'.ele'maxy=90;miny=0;%out='Satelliteelevationdata';out='卫星高度角(°)';case'.azi'maxy=180;miny=-180;%out='Satelliteazimuthaldata';out='卫星方位角(°)';otherwisedisp('Somethingswrong..!')end%FUNCTION:PLACEAMODIFIEDCOLORBAR++++++++++++++++++++++++++++++functionCB=cbar(loc,range,label);%.............................................................%CB=cbar(loc,range,label)%placesacolorbarat:%loc='v'inverticalor'h'inhorizontal%positionincurrentfigurescaledbetween:%range=[minmax]witha:%label='string'.%%fontsizeisreducedto10andwidthofbarishalfdefault.%%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)]);switchloccase'v'CB=colorbar('vertical');set(CB,'ylim',[range(1)range(2)]);POS=get(CB,'position');set(CB,'position',[POS(1)POS(2)0.03POS(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%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++functionjd=cal2jd(cal)%cal2jd1 将公历年月日时分秒转换到儒略日。%jd=cal2jd(cal)返回儒略日%cal：1x6矩阵，6列分别为年月日时分秒。构造cal时可以省略末尾的0%%公元1582年10月4日24:00点之前运用儒略历，公元1582年10月15日00:00点之后运用公历iflength(cal)<6 cal(6)=0;endyear=cal(1);month=cal(2);day=cal(3)+(cal(4)*3600+cal(5)*60+cal(6))/86400;y=year+4800;%4801B.C.isacenturyyearandalsoaleapyear.if(year<0) y=y+1; %Pleasenotethatthereisnoyear0A.D.endm=month;if(m<=2) %JanuaryandFebruarycomeafterDecember. m=m+12; y=y-1;ende=floor(30.6*(m+1));a=floor(y/100); %numberofcenturies%教皇格雷戈里十三世于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);%numberofcenturyyearsthatarenotleapyearsendc=floor(365.25*y);%Juliancalendaryearsandleapyearsjd=b+c+e+day-32167.5;%+++++++++++++++++++++++++++++++++++++++++++++++++functioncal=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%startfromJulianMarch1,4801B.C.if(J<2299160.5) %before1582.10.4.24:00isJuliancalender j0=floor(J+0.5); dd=J+0.5-j0;else %after1582.10.15.00:00isGregoriancalender %numberofcerturyyearsthatarenotleapyear j0=n1+n2+n3+J+10; dd=j0+0.5-floor(j0+0.5); j0=floor(j0+0.5);endj0=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);ifmonth>12 month=month-12; year=year+1;endyear=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];%+++++++++++++++++++++++++++++++++++++++++++++++++functionmjd=cal2mjd(cal)%cal2mjd 将公历年月日时分秒转换到简化儒略日。%mjd=cal2mjd(cal)返回简化儒略日%cal：1x6矩阵，6列分别为年月日时分秒。构造cal时可以省略末尾的0jd=cal2jd(cal);mjd=jd-2400000.5;%FUNCTION:SECO