examplecr_gps.m

此功能包用于各种GPS坐标和时间的转换

% ex_gps.m
%
% Example script for generating GPS measurements. 
% This example is for a static receiver fixed to the surface of the Earth.

% Written by: Jimmy LaMance 
% Copyright (c) 1998 by Constell, Inc.
% Mods. C.Rizos, UNSW 2002

% functions called: UTC2GPS, GPST2SEC, FIND_ALM, READYUMA, ALM2GEPH,
%                   PROPGEPH, ECEF2NED, PSEUDO_R, ECEF2LLA,
%                   NED2AZEL, VIS_DATA, SEC2GPST, LLA2ECEF,
%                   LOS, PLOTPASS, VIS_DATA
% ------------------------------------------------------------------------
% Variables to be used in this example:
% A BETTER MEANS OF INPUTTING THIS INFORMATION?
start_time = [2003 03 20 9 0 0];        % UTC start time
duration = 1*3600;                     % 1 hour duration
base_station = [-33*pi/180 151*pi/180 80];       % ground station coordinates

sample_rate = 10;                      % data sampling every 20 seconds

% Conversion factor from degree to radians3
d2r = pi / 180;

% Compute the stop time based on the start time and duration:
start_gps = UTC2GPS(start_time);
start_gps_sec = GPST2SEC(start_gps);
stop_gps_sec = start_gps_sec + duration;
stop_gps = SEC2GPST(stop_gps_sec);

% Find the almanac that is most recent to the start time:
% A BETTER MEANS OF INPUTTING THIS INFORMATION? 
% WHERE DOES IT GET THE ALMANAC FILE?
alm_file = FIND_ALM(start_gps(1));

% Read in the almanac:
alm = READYUMA(alm_file);

% Sort out the unhealthy satellites:
I_healthy = find(alm(:,2) == 0);
alm = alm(I_healthy,:);

% Convert from almanac to GPS/GLONASS ephemeris format:
gps_ephem = ALM2GEPH(alm);

% Compute satellite positions in ECEF system:
[t_gps,prn_gps,x_gps,v_gps] = PROPGEPH(gps_ephem,start_gps,stop_gps,sample_rate);

t_base = t_gps(1,:);

% Convert the station coordinates to ECEF:
x_base = ones(size(t_base,1),1) * LLA2ECEF(base_station);  
v_base = zeros(size(x_base));

% MASK ANGLE, HERE 10 DEGREES
% A BETTER MEANS OF INPUTTING THIS INFORMATION?
base_mask = 10*d2r;

% Model tropo, iono, and receiver clock.  Do not model satellite motion, SA,
% Earth rotation, satellite clocks, line baises, or relativity:
% A BETTER MEANS OF INPUTTING THIS INFORMATION?
base_model = [0 0 1 1 1 1 0 0 1 0 0];          % turn on error models
base_code_noise = .2;
base_carrier_noise = .01;
base_seed = 0; 

[t_pr_base,prn_base,pr_base,pr_errors_base,base_ndex] = ...
               PSEUDO_R(t_base,x_base,v_base,t_gps,[prn_gps x_gps],...
                        v_gps,base_model,base_seed, ...
                        base_code_noise, base_carrier_noise);

% Compute masking for the base station.  Start by compute the LOS from
% the base to the satellites in ECEF.  This is easily done using the indices
% returned from PSEUDO_R where the LOS has been vectorized.  Masking is
% done external to PSEUDO_R so that it can be performed in any coordinate
% system.  This example is for a station on the surface of the Earth with
% the antenna boresight pointed up.  Therefore, the NED system is used.
los_base_ecef = x_gps(base_ndex(:,2),:) - x_base(base_ndex(:,1),:);

% Rotate the LOS to NED, using the station as the reference for the NED
% coordinate system and rotation:
ref_lla = ECEF2LLA(x_base(base_ndex(:,1),:));
los_base_ned = ECEF2NED(los_base_ecef,ref_lla);

% Compute the azimuth/elevation of the NED vectors:
[az, el] = NED2AZEL(los_base_ned);

% Apply the masking in the NED coordinate system:
[visible_data, I_pass] = VIS_DATA(base_mask,[az el]);

% Remove all of the station data that did not pass the elevation mask test:
t_pr_base = t_pr_base(I_pass,:);
prn_base = prn_base(I_pass,:);
pr_base = pr_base(I_pass,:);   
pr_errors_base = pr_errors_base(I_pass,:);
base_ndex = base_ndex(I_pass,:);
x_ned = los_base_ned(I_pass,:);

% Rename the pr_errors for easier use later:
pr_sa_errb = pr_errors_base(:,1) + pr_errors_base(:,2);
trop_base = pr_errors_base(:,3) + pr_errors_base(:,4);
iono_base = pr_errors_base(:,5);
clk_biasb = pr_errors_base(:,6);
clk_driftb = pr_errors_base(:,7);

% plot troposphere contribution to the PR error
fh(1) = PLOTPASS(t_pr_base,trop_base,prn_base,...
         'Example of Troposphere Effects on Pseudorange Measurements',...
         'Tropo (m)');

% plot ionosphere contribution to the PR error
fh(2) = PLOTPASS(t_pr_base,iono_base,prn_base,...
         'Example of Ionosphere Effects on Pseudorange Measurements',...
         'Iono (m)');

% plot the receiver clock bias for a single satellite 
% the receiver clock bias is common to all satellites
I = find(prn_base == prn_base(2));
fh(3) = PLOTPASS(t_pr_base(I,:),clk_biasb(I),prn_base(I),...
         'Example of Receiver Clock Bias on Pseudorange Measurements',...
         'Clock Bias (m)');

fh(4) = PLOTPASS(t_pr_base(I,:),clk_driftb(I),prn_base(I),...
         'Example of Receiver Clock Drift on Pseudorange Measurements',...
         'Clock Drift (m/2)');