基于Matlab的TEQC繪圖程序代碼

%在matlab下新建一個(gè)m文件，將以下代碼干脆拷貝進(jìn)去，即可執(zhí)行。%須要一個(gè)TEQC生成的plot文件作為參數(shù)functionout=teqcplot3(files);%讀取TEQC生成的Plot文件，繪制數(shù)據(jù)圖表，支持Copmact、Compact2、Compact3格式%選取一個(gè)TEQC的Plot文件%格式說(shuō)明%*.sn1載波L1的信噪比Signaltonoiseratio(S/N)%*.sn2載波L2的信噪比Signaltonoiseratio(S/N)CarrierL2%*.iod*.d12*.d21電離層延遲觀測(cè)值改變率(米/秒)Derivativeofionosphericdelayobservable(m/s)%*.ion*.i12*.i21電離層延遲觀測(cè)值(米)Ionosphericdelayobservable(m)%*.mp1*.m12載波L1的多路徑誤差MultipathCarrierL1%*.mp2*.m21載波L2的多路徑誤差MultipathCarrierL2%*.azi衛(wèi)星方位角（°）Satelliteazimuthaldata(degrees)%*.ele衛(wèi)星高度角（°）Satelliteelevationdata(degrees)ifnargin==0[]=uigetfile('*.sn1;*.sn2;*.iod;*.ion;*.mp1;*.mp2;*.azi;*.ele;*.i12;*i21;*.m12;m21;*.d12;*d21',...'請(qǐng)選擇TEQC報(bào)告文件：');else[path,]=(files);path=[path'\'];filen={[filenext]};end%讀取這個(gè)文件(filen);%按行讀取文件至數(shù)組A[A]=importdata([pathfile],'\t');%定義SAT，存放衛(wèi)星數(shù)據(jù)%GPS有32顆衛(wèi)星，存放序號(hào)1-32，%GLONASS有32顆衛(wèi)星，存放序號(hào)33-64，%BEIDOU有35顆衛(wèi)星，存放序號(hào)65-99SAT(1:length(A),1:99)=NaN;%sats(1:length(A),1:99)=NaN;%存放采樣時(shí)間，單位秒tsec(1:length(A))=NaN;%讀取文件的第一行{1};%推斷是哪種格式switchfilelxcase'COMPACT'%讀取數(shù)據(jù)采樣間隔t_samp=char(A(3));%讀取起先時(shí)間mjl=char(A(4));%讀取數(shù)據(jù)采樣間隔T_SAMP=str2num(t_samp(max(find(t_samp=='')):end));%讀取數(shù)據(jù)采樣起先時(shí)間MJL_START=str2num(mjl(max(find(mjl=='')):end));%轉(zhuǎn)成時(shí)間序列數(shù)字，dateserialnumber，從0000年1月1日0時(shí)0分0秒起先計(jì)算的十進(jìn)制天數(shù)MJD_START=MJL_START+678941.999999741;%i為行號(hào)n=1;i=5;case'COMPACT2'%讀取數(shù)據(jù)采納間隔t_samp=char(A(2));%讀取起先時(shí)間mjl=char(A(3));%讀取數(shù)據(jù)采樣間隔T_SAMP=str2num(t_samp(max(find(t_samp=='')):end));%讀取數(shù)據(jù)采樣起先時(shí)間MJL_START=str2num(mjl(max(find(mjl=='')):end));%轉(zhuǎn)成時(shí)間序列數(shù)字，dateserialnumber，從0000年1月1日0時(shí)0分0秒起先計(jì)算的十進(jìn)制天數(shù)MJD_START=MJL_START+678941.999999741;n=1;i=4;case'COMPACT3'%讀取起先時(shí)間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})'時(shí)'char(s{6})'分'num2str(str2num(char(s{7})),'%02d')'秒'];%獲得采樣的起先時(shí)間，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('數(shù)據(jù)格式存在問(wèn)題');returnend%n=1;i=3;%sats=str2num(A{i});snygg(filen,'_','-');%生成進(jìn)度條h=waitbar(0,['正在讀取數(shù)據(jù)，請(qǐng)稍等……'char(snyggfilen)]);whilei<length(A);waitbar(i/length(A),h);drawnow%讀取衛(wèi)星編號(hào)數(shù)據(jù)行采樣時(shí)間衛(wèi)星數(shù)量衛(wèi)星編號(hào)COMAPCT3satbhstr=char(A(i));s=splitstr(satbhstr,'**',32);%假如衛(wèi)星編號(hào)數(shù)據(jù)為0，說(shuō)明數(shù)據(jù)缺失，退出本函數(shù)ifstrtrim(satbhstr)=='0'disp('數(shù)據(jù)為空！');close(h);returnendswitchfilelxcase'COMPACT3'%記錄采樣時(shí)間，單位秒tsec(n)=str2num(s{1});ifs{2}=='-1'%假如為-1，運(yùn)用上次的衛(wèi)星編號(hào)行satbhstr=oldsatbh;s=splitstr(satbhstr,'**',32);elseoldsatbh=satbhstr;end%獲得當(dāng)前衛(wèi)星數(shù)量satcount=str2num(s{2});case{'COMPACT','COMPACT2'}%記錄采樣時(shí)間，單位秒tsec(n)=(n-1)*T_SAMP;ifs{1}=='-1'%假如為-1，運(yùn)用上次的衛(wèi)星編號(hào)行satbhstr=oldsatbh;s=splitstr(satbhstr,'**',32);elseoldsatbh=satbhstr;end%獲得當(dāng)前衛(wèi)星數(shù)量satcount=str2num(s{1});end%讀取衛(wèi)星數(shù)據(jù)行satdatastr=char(A(i+1));sdata=splitstr(satdatastr,'**',32);fork=3:satcount+2switchfilelxcase'COMPACT'%假如是COMPACT格式，衛(wèi)星編號(hào)前面沒(méi)有字母，默認(rèn)是GPS衛(wèi)星，在前面增加字符Gsatbh=['G'char(s{k-1})];case'COMPACT2'%假如是COMPACT2格式，衛(wèi)星編號(hào)行第2個(gè)數(shù)據(jù)起先是衛(wèi)星編號(hào)satbh=char(s{k-1});case'COMPACT3'%假如是COMPACT3格式，衛(wèi)星編號(hào)行第3個(gè)數(shù)據(jù)起先是衛(wèi)星編號(hào)satbh=char(s{k});endswitchsatbh(1);case'G'%假如是GPS衛(wèi)星，獲得與編號(hào)對(duì)應(yīng)的存儲(chǔ)序號(hào)，范圍1-32index=str2num(satbh(2:3));case'R'%假如是GLONASS衛(wèi)星，獲得與編號(hào)對(duì)應(yīng)的存儲(chǔ)序號(hào)，范圍33-64index=32+str2num(satbh(2:3));case'C'%假如是BEIDOU衛(wèi)星，獲得與編號(hào)對(duì)應(yīng)的存儲(chǔ)序號(hào)，范圍65-99index=64+str2num(satbh(2:3));otherwise%其他狀況，獲得與編號(hào)對(duì)應(yīng)的存儲(chǔ)序號(hào)，范圍1-32index=str2num(satbh(2:3));end%switch%獲得對(duì)應(yīng)的衛(wèi)星數(shù)據(jù)sdatastr=char(sdata{k-2});%假如衛(wèi)星數(shù)據(jù)不是數(shù)值型，用0替代ifisempty(str2num(sdatastr))==1sdatastr='0.000';end%將衛(wèi)星數(shù)據(jù)存入SAT數(shù)組的對(duì)應(yīng)位置SAT(n,index)=str2num(sdatastr);end%forn=n+1;%依據(jù)數(shù)據(jù)文件結(jié)構(gòu)，一行是衛(wèi)星編號(hào)，下一行就是對(duì)應(yīng)的衛(wèi)星數(shù)據(jù)%程序一次讀取2行數(shù)據(jù)i=i+2;end%while%數(shù)據(jù)讀取完畢，關(guān)閉進(jìn)度條close(h);%獲得已運(yùn)用的衛(wèi)星編號(hào)，返回給sat_bhsat_bh=getsatbh(SAT);%增加最終一個(gè)結(jié)束字符sat_bh_temp=[sat_bh'{'End'}];sat_bh_end=sat_bh_temp';%刪去全部都為NaN值的列，返回給sat_datasat_data=delnancol(SAT);%刪去全部NaN值的采樣時(shí)間列t_seconds=delnancol(tsec);%刪去全部都為NaN值的行，返回sat_datassat_datas=delnanrow(sat_data);%增加一個(gè)nan列[row,col]=size(sat_datas);bb(1:row,1)=NaN;sat_datas=[sat_datasbb];%將含有NaN值的數(shù)據(jù)替換為0，返回給s_data%s_data=repnan20(sat_data);%將采樣時(shí)間轉(zhuǎn)成依次日期serialdatenumberswitchfilelxcase{'COMPACT','COMPACT2'}t_s_jd_time=MJD_START;case'COMPACT3't_s_jd_time=datenum(t_s_time);endfori=1:length(t_seconds)%將每次采樣時(shí)間轉(zhuǎn)換成相應(yīng)的依次日期t_jd_time(i)=t_s_jd_time+t_seconds(i)/60/60/24;end%繪制圖表++++++++++++++++++++++++++%依據(jù)不同的文件類型，設(shè)置坐標(biāo)軸的范圍[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('衛(wèi)星編號(hào)')%timestr=secs2hms(length(sat));T=title(['TEQC報(bào)告文件:'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));%+++++++++++++++++++++++++++++++++++++++%以下是本函數(shù)中所運(yùn)用的一些相關(guān)函數(shù)%++++++++++++++++++++++++++++++++++++++++++++++++++%+++++++++++++++++++++++++++++++++++++++++++++++++functions=splitstr(str,split,num);%用指定分隔符分隔字符串函數(shù)%str='name02010-05-0633QtBKB';%'**'表示特別分隔符回車、空格、換行、制表符等%num表示分隔字符串的個(gè)數(shù)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%常規(guī)分隔符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);%獲得衛(wèi)星數(shù)據(jù)中的衛(wèi)星編號(hào)，返回一個(gè)字符串?dāng)?shù)組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);%依據(jù)文件類型，設(shè)置坐標(biāo)軸的范圍[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多路徑觀測(cè)值';maxy=1;miny=-1;case{'.mp2','.m21'}%out='MultipathL2';out='L2多路徑觀測(cè)值';maxy=1;miny=-1;case{'.iod','.d12','d21'}%out='Derivativeofionosphericdelayobservable(m/s)';out='電離層延遲觀測(cè)值改變率（米/秒）';maxy=1;miny=-1;case{'.ion','.i12','i21'}maxy=2;miny=-2;%out='Ionosphericdelayobservable(m)';out='電離層延遲觀測(cè)值(米)';case'.ele'maxy=90;miny=0;%out='Satelliteelevationdata';out='衛(wèi)星高度角(°)';case'.azi'maxy=180;miny=-180;%out='Satelliteazimuthaldata';out='衛(wèi)星方位角(°)';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 將公歷年月日時(shí)分秒轉(zhuǎn)換到儒略日。%jd=cal2jd(cal)返回儒略日%cal：1x6矩陣，6列分別為年月日時(shí)分秒。構(gòu)造cal時(shí)可以省略末尾的0%%公元1582年10月4日24:00點(diǎn)之前運(yùn)用儒略歷，公元1582年10月15日00:00點(diǎn)之后運(yùn)用公歷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日以教皇訓(xùn)令頒布，將1582年10月5日至14抹掉。1582年10月4日過(guò)完后其次天是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)%從儒略日計(jì)算公歷年月日時(shí)分秒%cal=jd2cal(J)%返回的cal是1x6矩陣，6列分別為年月日時(shí)分秒%%公元1582年10月4日24:00點(diǎn)之前運(yùn)用儒略歷，公元1582年10月15日00:00點(diǎn)之后運(yùn)用公歷if(J<1721423.5) %公元1月1日0時(shí) 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 將公歷年月日時(shí)分秒轉(zhuǎn)換到簡(jiǎn)化儒略日。%mjd=cal2mjd(cal)返回簡(jiǎn)化儒略日%cal：1x6矩陣，6列分別為年月日時(shí)分秒。構(gòu)造cal時(shí)可以省略末尾的0jd=cal2jd(cal);mjd=jd-2400000.5;%FUNCTION:SECO