|
|
function [dtsvS] = GpsEph2Dtsv(gpsEph,tS) |
|
|
%[dtsvS] = GpsEph2Dtsv(gpsEph,tS) |
|
|
% |
|
|
% Calculate satellite clock bias |
|
|
% Inputs: |
|
|
% gpsEph = vector of GPS ephemeris structures defined in ReadRinexNav |
|
|
% tS GPS time of week (secs) at time of transmission, calculate this as: |
|
|
% tS = trx - PR/c, |
|
|
% where trx = time of reception (receiver clock) |
|
|
% PR = pseudorange |
|
|
% |
|
|
% tS may be a vector, gpsEph is a structure or vector of structures |
|
|
% If gpsEph is a vector then t must have the same number of rows as gpsEph |
|
|
% tS(i) is interpreted as being associated with gpsEph(i,:) |
|
|
% If gpsEph is a single structure then tS may be any length |
|
|
% |
|
|
% Outputs: |
|
|
% dtsvS = sat clock bias (seconds). GPS time = satellite time - dtsvS |
|
|
% length(dtsvS) = length(tS) |
|
|
% |
|
|
% See IS-GPS-200 and RTCM Paper 99-92 SC104-88 for details of data |
|
|
|
|
|
%Author: Frank van Diggelen |
|
|
%Open Source code for processing Android GNSS Measurements |
|
|
|
|
|
%% Check size of inputs |
|
|
if min(size(tS))>1 |
|
|
error('tS must be a vector or a scalar, not a matrix') |
|
|
end |
|
|
tS=tS(:)';%make t a row vector, to match [gpsEph.*] vectors below |
|
|
pt=length(tS); |
|
|
[p]=length(gpsEph); |
|
|
if (p>1 && pt~=p), |
|
|
error('If gpsEph is a vector tS must be a vector with #rows = length(gpsEph),\n') |
|
|
end |
|
|
%% |
|
|
%%Extract the necessary variables from gpsEph |
|
|
TGD = [gpsEph.TGD]; |
|
|
Toc = [gpsEph.Toc]; |
|
|
af2 = [gpsEph.af2]; |
|
|
af1 = [gpsEph.af1]; |
|
|
af0 = [gpsEph.af0]; |
|
|
Delta_n = [gpsEph.Delta_n]; |
|
|
M0 = [gpsEph.M0]; |
|
|
e = [gpsEph.e]; |
|
|
Asqrt = [gpsEph.Asqrt]; |
|
|
Toe = [gpsEph.Toe]; |
|
|
%% |
|
|
%%Calculate dependent variables ------------------------------------------------ |
|
|
|
|
|
tk = tS - Toe; %time since time of applicability |
|
|
I = find(tk > 302400.0); |
|
|
if any(I), tk(I) = tk(I)-GpsConstants.WEEKSEC; end, |
|
|
I = find(tk < -302400.0); |
|
|
if (I), tk(I) = tk(I)+GpsConstants.WEEKSEC; end, |
|
|
|
|
|
A = Asqrt.^2; %semi-major axis of orbit |
|
|
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec) |
|
|
n=n0+Delta_n; %Corrected Mean Motion |
|
|
Mk=M0+n.*tk; %Mean Anomaly |
|
|
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly |
|
|
%% |
|
|
%%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1) ---------------- |
|
|
dt = tS - Toc; |
|
|
I = find(dt > 302400.0); |
|
|
if any(I), dt(I) = dt(I)-GpsConstants.WEEKSEC; end, |
|
|
I = find(dt < -302400.0); |
|
|
if (I), dt(I) = dt(I)+GpsConstants.WEEKSEC; end, |
|
|
|
|
|
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ... |
|
|
GpsConstants.FREL.*e.*Asqrt.*sin(Ek) -TGD; |
|
|
|
|
|
end % end of function GpsEph2Dtsv |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|
|
|
|
|
% Copyright 2016 Google Inc. |
|
|
% |
|
|
% Licensed under the Apache License, Version 2.0 (the "License"); |
|
|
% you may not use this file except in compliance with the License. |
|
|
% You may obtain a copy of the License at |
|
|
% |
|
|
% http://www.apache.org/licenses/LICENSE-2.0 |
|
|
% |
|
|
% Unless required by applicable law or agreed to in writing, software |
|
|
% distributed under the License is distributed on an "AS IS" BASIS, |
|
|
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
|
|
% See the License for the specific language governing permissions and |
|
|
% limitations under the License. |
|
|
|
|
|
|
