function testinterpeph(td,sats,nmax,opt)
if nargin<1, td=caltomjd([2004,10,3]); end
if nargin<2, sats={}; end
if nargin<3, nmax=10; end
if nargin<4, opt=0; end

sp3dir='k:\gps\sp3\%Y%m';
%ephdir='k:\gps\eph\%Y%m';
ephdir='k:\gps\eph\%Y%m\cod';
t=-43200:900:86400+43200;
time=0:30:86400;

if isempty(sats)
    sats={}; for n=1:31, sats={sats{:},sprintf('GPS%02d',n)}; end
end

%ephs=readeph(td,t,sats,sp3dir,'igs');     % original ephemeris
ephs=readeph(td,t,sats,sp3dir,'cod');
ephe=readeph(td,time,sats,ephdir,'ephs'); % interpolated by orbit-model

disp(['            3D(m)     X(m)     Y(m)     Z(m)     VX(m/s)  VY(m/s)  VZ(m/s) n=',...
     num2str(nmax)]);
for n=1:length(sats)
    [pos,vel]=interplag(t,ephs(:,1:3,n),time,nmax); % lagrange interpolation
    err=[pos,vel]-ephe(:,:,n);
    rms(n,:)=sqrt(mean(err.^2));
    disp(sprintf('%-7s : %8.4f %8.4f %8.4f %8.4f  %8.4f %8.4f %8.4f',sats{n},...
         sqrt(sum(rms(n,1:3).^2)),rms(n,:)))
    if opt, ploterr(td,time,err,sats{n},nmax); end
end
rms=sqrt(mean(rms(~isnan(rms(:,1)),:).^2,1));
disp(sprintf('TOTAL   : %8.4f %8.4f %8.4f %8.4f  %8.4f %8.4f %8.4f',sqrt(sum(rms.^2)),rms))

%-------------------------------------------------------------------------------
function ploterr(td,time,err,sat,nmax)
margin=[0.08,0.03,0.04,0.012];
gut('newfigg','','',[600,400]);
ggt('subplotv','',1,2,'taxis',td,'topts','nolabel','xlim',[time(1),time(end)]/3600,...
    'ylim',[-0.02,0.02],'margin',margin);
plot(time/3600,err(:,1:3));
rms=sqrt(mean(err(:,1:3).^2));
ggt('mtext','_1',sprintf('RMS 3D:%.4fm X:%.4fm Y:%.4fm Z:%.4fm',sqrt(sum(rms.^2)),rms));
ylabel('Satellite Position Error (m)');

title(['Interpolated Ephemeris Error : ',sat,' ',tstr(td,time(1)),'-',tstr(td,time(end)),...
       ' n=',num2str(nmax)])

ggt('subplotv','',2,2,'taxis',td,'xlim',[time(1),time(end)]/3600,...
    'ylim',[-0.02,0.02],'margin',margin);
plot(time/3600,err(:,4:6));
rms=sqrt(mean(err(:,4:6).^2));
ggt('mtext','_2',sprintf('RMS 3D:%.4fm/s X:%.4fm/s Y:%.4fm/s Z:%.4fm/s',sqrt(sum(rms.^2)),rms));
ylabel('Satellite Velocity Error (m/s)');

function s=tstr(td,time), s=sprintf('%04d/%02d/%02d %02d:%02d:%02d',mjdtocal(td,time));
