% input information
nav_file = 'epgga2.010';
sp3_file = 'igs11484.sp3';
dow = 4;
sv = 31;
doy = 010;
year = 2002;

% read sp3 file
[Xs,Ys,Zs,dT,Ts] = read_sp3(sp3_file,sv);

% convert to meters
SP3 = [Xs'; Ys'; Zs'].*1e3;

% read rinex navigation file
eph = read_rinexn(nav_file);

% start at first epoch found in sp3 file (in sec of GPS week)
t = Ts(1) + dow*24*3600;
% every fifteen minutes
dt = 15*60;
EPH = [];
T = [];
% compute position of satellite sv
for i=t:dt:t+(24*3600-15*60)
  satpos = get_satpos(i,sv,eph,0); % use closest toe
  %satpos = get_satpos(i,sv,eph,1); % use last toe (as GPS receiver would do)
  EPH = [EPH  satpos];
end

% length of residual vector
R = sqrt((EPH(1,:)-SP3(1,:)).^2 + (EPH(2,:)-SP3(2,:)).^2 + (EPH(3,:)-SP3(3,:)).^2);

% range of SP3 and EPH
Rsp3 = sqrt(SP3(1,:).^2+SP3(2,:).^2+SP3(3,:).^2);
Reph = sqrt(EPH(1,:).^2+EPH(2,:).^2+EPH(3,:).^2);

% plot
figure(); hold on;
plot((Ts)./3600,Reph-Rsp3,'-ob');
xlabel(['time (hours)']);
ylabel(['range difference (m)']);
title(['range difference, PRN' num2str(sv) ', day ' num2str(doy) '/' num2str(year)]);

