This repository has been archived by the owner on Nov 22, 2021. It is now read-only.
-
-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathb_point.m
59 lines (55 loc) · 1.91 KB
/
b_point.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
function pos = b_point(Eph,pr,sv,time) %%%%%
%B_POINT Prepares input to the Bancroft algorithm for finding
% a preliminary position of a receiver. The input is
% name of ephemeris file, four or more pseudoranges
% and the coordinates of the satellites.
%Kai Borre and C.C. Goad 11-24-96
%Copyright (c) by Kai Borre
%$Revision: 1.0 $ $Date: 1997/09/26 $
v_light = 299792458;
dtr = pi/180;
m = size(pr,1);
%%%%Eph = get_eph(efile);
for t = 1:m
col_Eph(t) = find_eph(Eph,sv(t),time);
end
pos = zeros(4,1); % First guess Earth center
for l = 1:2
B = [];
for jsat = 1:m
k = col_Eph(jsat);
tx_RAW = time - pr(jsat)/v_light;
toc = Eph(21,k);
dt = check_t(tx_RAW-toc);
tcorr = (Eph(2,k)*dt + Eph(20,k))*dt + Eph(19,k);
tx_GPS = tx_RAW-tcorr;
dt = check_t(tx_GPS-toc);
tcorr = (Eph(2,k)*dt + Eph(20,k))*dt + Eph(19,k);
tx_GPS = tx_RAW-tcorr;
X = satpos(tx_GPS, Eph(:,k));
if l ~= 1
rhosq = (X(1)-pos(1))^2+(X(2)-pos(2))^2+(X(3)-pos(3))^2;
traveltime = sqrt(rhosq)/v_light;
Rot_X = e_r_corr(traveltime,X);
rhosq = (Rot_X(1)-pos(1))^2 + ...
(Rot_X(2)-pos(2))^2 + ...
(Rot_X(3)-pos(3))^2;
traveltime = sqrt(rhosq)/v_light;
[az,el,dist] = topocent(Rot_X, Rot_X-pos(1:3,:));
trop = tropo(sin(el*dtr),h/1000,1013.0,293.0,50.0,...
0.0,0.0,0.0);
else trop = 0;
end
corrected_pseudorange = pr(jsat)+v_light*tcorr-trop;
B(jsat,1) = X(1);
B(jsat,2) = X(2);
B(jsat,3) = X(3);
B(jsat,4) = corrected_pseudorange;
end; % jsat-loop
pos = bancroft(B);
[phi,lambda,h] = togeod(6378137, 298.257223563, ...
pos(1), pos(2), pos(3));
end; % l-loop
clock_offset = pos(4)/v_light;
pos(4) = clock_offset*1.e+9; % ns
%%%%%%%%%%% b_point.m %%%%%%%%%%%%%%%%%%%%%%%%%