%                              uclock.m
%  Scope:   This MATLAB macro computes user clock bias and drift by using a 
%           simplified second order model [1]. WGS-84 constants are used.
%  Usage:   [cbias,cdrift] = uclock(nsteps,sb,sf,dt,iseed)
%           [cbias,cdrift] = uclock(nsteps,sb,sf,dt) when seed is not reset 
%           [cbias,cdrift] = uclock(nsteps,sb,sf)    when seed is not reset and 
%                                                    dt = 1 as default values 
%  Description of parameters:
%           nsteps  - input, number of steps to be executed
%           sb      - input, white noise spectral density amplitude for the first 
%                     component corresponding to clock bias, in seconds
%           sf      - input, white noise spectral density amplitude for the second 
%                     component corresponding to clock drift, in second/second
%           dt      - input (optional), time interval between two consecutive steps,
%                     in seconds; the default value is 1 second
%           iseed   - input (optional), seed value
%           cbias   - output, clock bias vector for nsteps, in meters
%           cdrift  - output, clock drift vector for nsteps, in meters/sec
%  External Matlab macros used:  wgs84con
%  Remark:  The initial clock bias and drift values can be changed; the default 
%           values are 0.
%  References:  [1] Brown, Grover R., Hwang, P.Y.C., Introduction to random signals 
%                   and applied Kalman filtering. Third Edition. John Wiley $ Sons,
%                   New York, 1997, pp. 428-431
%  Last update:  10/08/00
%  Copyright (C) 1999-00 by LL Consulting. All Rights Reserved.

function    [cbias,cdrift] = uclock(nsteps,sb,sf,dt,iseed)

wgs84con;
% global constants used: c_speed

if nargin == 5
   randn('seed',iseed);
end
if nargin < 4
   dt = 1.;
end
cbias = [];
cdrift = [];

sbias = sqrt(sb*dt);
sdrift = sqrt(sf*dt);

wnbias = sbias * randn(nsteps-1,1);
wndrift = sdrift * randn(nsteps-1,1);

dt2 = 0.5*dt;
cbias(1) = 0.;                      % initial clock bias  
cdrift(1) = 0.;                     % initial clock drift 
for k = 1:nsteps-1
   cbias(k+1,1) = cbias(k,1) + dt*cdrift(k,1) + wnbias(k);
   cdrift(k+1,1) = cdrift(k,1) + dt2*wnbias(k) + wndrift(k);
end   

cbias = c_speed * cbias;            % clock bias in meters
cdrift = c_speed * cdrift;          % clock drift in meters/second