function testtec
% test driver for TECDEMO
%
%          Copyright (C) 2006 by T.TAKASU, All rights reserved.
%
% version: $Revision: 2 $ $Date: 06/02/24 18:36 $
% history: 2006/02/24 1.1 new

t0=datenum(2004,10,1);     % date number by datenum.m (days)
%obsp='l:\gps\obs_gsi\20041001\';
obsp='k:\gps\obs_igs\200410\';
navp='k:\gps\nav\200410\';

% stations list
rcvs={...
'93017', '93022', '93044', '93076', '93101', '940001','940002','940003',...
'940004','940006','940007','940010','940011','940013','940018','940019',...
'940020','940024','940026','940027','940029','940030','940032','940038',...
'940040','940042','940044','940047','940051','940053','940055','940057',...
'940062','940065','940069','940070','940073','940075','940078','940079',...
'940080','940082','940085','940088','940093','940096','940099','940100',...
'95113'};
%rcvs={...
%'ALGO','ALIC','AMC2','AREQ','AUCK','BAHR','BRUS','CAS1','CEDU','CHAT',...
%'CHUR','COCO','DARW','DAV1','DRAO','FAIR','FLIN','GODE','GOLD','GOUG',...
%'GRAZ','GUAM','HOB2','IISC','IRKT','KARR','KELY','KERG','KIRU','KOUR',...
%'LHAS','LPGS','MAC1','MAG0','MAS1','MATE','MCM4','MKEA','NKLG','NLIB',...
%'NRC1','NRIL','OHI2','ONSA','OUS2','PERT','PIMO','POL2','POTS','REUN',...
%'RIOG','STJO','SYOG','THTI','TIDB','TROM','TSKB','USNO','VILL','YELL'};

% read rinex observation data/navigation messages
if exist('obs1.mat')~=2
    fs={}; dv=datevec(t0);
    s=sprintf('%03d0.%02d',t0-datenum(dv(1),1,1)+1,mod(dv(1),100));
    for r=rcvs, fs={fs{:},[obsp,r{1}(end-3:end),s,'o']}; end
    [obs,iobs,nav,inav]=readrnx(t0,0,86400,{fs{:},[navp,'brdc',s,'n']});
    save('obs1.mat','t0','obs','iobs','nav','inav');
else
    load('obs1.mat');
end

time=(1:2:23)*3600;        % estimation time (sec)
ts=-3600:300:3600;         % data sampling time relative to time (sec)
model=0;                   % tec model (0:local,1:global)
nmax=8;                   % max degree of tec model
pos0=[38,138]*pi/180;      % origin latitude/longitude of local tec model (rad)

% read station positions
posr=readrpos(rcvs);

for n=1:length(time)
    % tec demo
    [enm,dcb(n,:)]=tecdemo(t0,time(n)+ts,obs,iobs,nav,inav,posr,model,nmax,pos0);
    
    % plot tec map
    plottec(t0,time(n),enm,model,nmax,pos0,posr,[0,60]);
end
% plot dcbs
plotdcb(t0,time,dcb(:,1:31),'Satellite');
plotdcb(t0,time,dcb(:,32:end),'Receiver');

% plot tec map -----------------------------------------------------------------
function plottec(t0,time,enm,model,nmax,pos0,posr,range)
if ~model, lats=pos0(1)*180/pi+(-10:10); lons=pos0(2)*180/pi+(-10:10);
else lats=-90:5:90; lons=-180:5:180; end
for n=1:length(lats)
    for m=1:length(lons)
        tec(n,m)=tecmodel(time,[lats(n),lons(m)]*pi/180,enm,model,pos0);
    end
end
tec(tec<0)=0;
figure, hold on, box on, axis equal, caxis(range);
[c,h]=contourf(lons,lats,tec,range(1):range(2)); set(h,'edgecolor','none');
xlim(lons([1,end])), ylim(lats([1,end]))
if model, colorbar('horizontal'); else colorbar; end
plot(posr(:,2)*180/pi,posr(:,1)*180/pi,'r.'); plotmap('w');
dv=datevec(t0+time/86400); 
title(sprintf('TEC MAP (TECU) : %04d/%02d/%02d %02d:%02d',dv(1:5)))

% plot dcbs --------------------------------------------------------------------
function plotdcb(t0,time,dcb,ti)
figure, plot(time/3600,dcb,'.-'), xlabel('time (H)'), ylabel('DCB (m)')
dv=datevec(t0); title(sprintf('%s DCB : %04d/%02d/%02d',ti,dv(1:3)));

% plot map ---------------------------------------------------------------------
function plotmap(s)
persistent p, if isempty(p), load('gshhs_l'); end
for q=p, lon=q.lon; lon(abs(diff(lon))>180)=nan; plot(lon,q.lat,s); end
