tropunb1.m

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

function [trop_dry, trop_wet] = tropunb1(lla,t_gps,elev);

% [trop_dry, trop_wet] = tropunb1(lla,t_gps,elev);
%
% Computes the wet and dry troposphere delays using
% the University of New Brunswick (UNB1), altitude dependent, 
% troposphere model.
%
% The UNB1 model is a composite model that uses the explicit 
% forms of Saastamoinen's delay algorithms combined with 
% Niell's mapping functions. The surface and lapse rate parameters 
% are from the U.S. 1976 Standard Atmosphere. There will be a
% several centimetre bias at the poles and the equator unless 
% surface met data is used.  This model was developed at the
% University of New Brunswick under contract to Nav Canada.
%
% Input:
%   lla        - matrix of geodetic position (1x3 or nx3) [lat, lon, height]
%                 lat and lon are in rad, height is in m
%                 valid latitudes are -pi/2 -> pi/2
%                 valid longitudes are -pi -> 2*pi
%   t_gps      - GPS time vector for az/el pairs [GPS_week, GPS_sec] (nx2)
%                 valid GPS_week values are 1-3640 (years 1980-2050)
%                 valid GPS_sec values are 0-604799
%   elev       - elevation angle to GPS satellites (rad) (nx1)
%                 Valid elevations are between -pi/2 and pi/2.
%						Elevations below .1 degree will have the a 
%        			delay mapped to .1 degree.  No mapping below
% 						zero degree elevation is modeled.
% Output:
%   trop_dry   - Dry troposphere component (m) (nx1)
%   trop_wet   - Wet troposphere component (m) (nx1)
%
% See also PSEUDO_R 

% Written in MATLAB by: Jimmy LaMance 4/12/99
% Copyright (c) 1999 by Constell, Inc.
%
% Original FORTRAN by: Paul Collins, University of New Brunswick
% UNB intellectual property rights reserved.
% For more information on UNB visit:
% http://www.unb.ca/GGE/Personnel/Langley/Langley.html

% Reference for UNB1 Troposphere Model: 
% Collins, J.P. and R.B. Langley. "A Tropospheric Delay Model for the User of
% the Wide Area Augmentation System." Final contract report prepared for Nav
% Canada, Department of Geodesy and Geomatics Engineering Technical Report
% No. 187, University of New Brunswick, Fredericton, N.B., Canada.
%
% Available in PDF format from
% http://gauss.gge.unb.ca/papers.pdf/waas.tropo.oct96.pdf

%%%%% BEGIN VARIABLE CHECKING CODE %%%%%
% declare the global debug mode
global DEBUG_MODE

% Initialize the output variables
trop_dry=[]; 
trop_wet=[]; 

% Check the number of input arguments and issues a message if invalid
msg = nargchk(3,3,nargin);
if ~isempty(msg)
  fprintf('\n');
  dbstack
  fprintf('\n');
  fprintf('%s  See help on TROPUNB1 for details.\n',msg);
  fprintf('Returning with empty outputs.\n\n');
  return
end

estruct.func_name = 'TROPUNB1';

% Develop the error checking structure with required dimension, matching
% dimension flags, and input dimensions.
estruct.variable(1).name = 't_gps';
estruct.variable(1).req_dim = [901 2];
estruct.variable(1).var = t_gps;
estruct.variable(1).type = 'GPS_TIME';
  
estruct.variable(2).name = 'lla';
estruct.variable(2).req_dim = [901 3; 901 2; 1 3; 1 2];
estruct.variable(2).var = lla;
estruct.variable(2).type = 'LLA_RAD';
  
estruct.variable(3).name = 'elev';
estruct.variable(3).req_dim = [901 1];
estruct.variable(3).var = elev;
estruct.variable(3).type = 'ELEVATION_RAD';
  
% Call the error checking function
stop_flag = err_chk(estruct);
  
if stop_flag == 1           
  fprintf('Invalid inputs to %s.  Returning with empty outputs.\n\n', ...
           estruct.func_name);
  return
end % if stop_flag == 1

%%%%% END VARIABLE CHECKING CODE %%%%%

%%%%% BEGIN ALGORITHM CODE %%%%%

% Set constants for T, P, and e. 
tempkelv = 288.15 - 0.00650 * lla(:,3);
I_zero = find(tempkelv < .1);
if ~isempty(I_zero)
   tempkelv(I_zero) = ones(size(tempkelv(I_zero)))*.1;
end

presmbar = 1013.25 *(tempkelv ./ 288.15).^5.256;
wvpmbar  = 11.691 *(tempkelv ./ 288.15).^21.024;
      
% Convert the elevation angle to degrees for future use
I_zero = find(elev < .0017);
if ~isempty(I_zero)
   elev(I_zero) = ones(size(elev(I_zero)))*.0017;
end
elev = elev * 180/pi;

% Compute the vertical delays for the wet and dry tropo
%   P,T,e equivalent to 324.8 N units at surface/msl
%   Uses Essen & Froome constants for consistency with Saastamoinen
[hzd,wad] = tzdds87(presmbar,tempkelv,wvpmbar,0.0065,3.0,lla);
     
% Compute the day of year from the current GPS time 
[ut, ls, day_of_year] = gps2utc(t_gps);   						% integer day
day_of_year = day_of_year + rem(t_gps(:,2),86400) ./ 86400;	% int + fractional day

% Compute the mapping function for the dry tropo
latdeg = lla(:,1)*180/pi;
hmf = nhmf2(day_of_year,latdeg,lla(:,3),elev);

% Compute the mapping function for the wet tropo
wmf = nwmf2(latdeg,elev);
     
% Compute the total troposphere delay
trop_dry = hzd.*hmf(:,1);
trop_wet = wad.*wmf(:,1);

% end of tropunb1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Built-in subroutines for UNB tropo model.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tzd1, tzd2] = tzdds87(P, T, E, t_lapse, lambda, lla);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Comments from the original FORTRAN code.
%
%     FUNCTION: Computes the hydrostatic zenith delay using the Davis,
%               et.al. formulation of Saastamoinen's model and the wet
%               delay based on an updated 'Smith & Weintraub' expression
%               and incorporating the two-parameter formula first derived
%               by Saastamoinen and then Askne & Nordius.
%
%     REFERENCES: (1) Davis,J.L., T.A. Herring, I.I. Shapiro, A.E.E. Rogers
%                 and G. Elgered. (1985). "Geodesy by radio interferometry:
%                 Effects of atmospheric modelling errors on estimates of
%                 baseline length." Radio Science, Vol. 20, No. 6, pp. 1593
%                  - 1607.
%
%                 (2) Askne,J. and H. Nordius (1987). "Estimation of
%                 tropospheric delay for microwaves from surface weather
%                 data." Radio Science, Vol. 22, No. 3, pp. 379 - 386.
%
%                 (3) Saastamoinen, J. (1973). "Contributions to the theory
%                 of atmospheric refraction." In three parts. Bulletin
%                 Geodesique, No. 105, pp.279-298; No. 106, pp.383-397;
%                 No. 107, pp.13-34.
%
%     HISTORY: 
%               ORIGINAL CODE:    Paul Collins        19.05.96
%   
%     PARAMETERS AND VARIABLES:
%
%         INPUT:
%                P        - Total Atmospheric Pressure          [mbar]
%                T        - surface temperature                 [K]
%                E        - surface water vapor pressure        [mbar]
%                t_lapse   - temperature lapse rate              [K/m]
%                lambda   - the "lambda" parameter         
%                height       - height of the site                  [m]
%                latitude      - geodetic latitude                   [rad]
%
%         LOCAL:
%                rd       - molar gas constant of dry air       [J/(K.kg)]
%                k1,K3'   - Refractivity constants:
%                  k1     - constant for dry/moist gas          [K/mbar]
%                  K3'    - derived constant for water vapour   [K^2/mbar]
%
%                Values are from Thayer (1974) quoted in Davis.
%                k1                = 77.604 +/- 0.014
%                k3' = k3 + k2'.Tm = 382000 +/- 4000
%
%                c1       - constant for Rd, g0, k1
%                c2       - constant for Rd, k3'
%
%        OUTPUT:
%                tzd1     - hydrostatic zenith delay             [m]
%                tzd2     - wet zenith delay                     [m]
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% Get height and latitude from lla
height = lla(:,3);
latitude = lla(:,1);

% Set up constants
k1 = 77.604;
k3prime = 3.82d5;
g0 = 9.784;
rd = 287.054;
c1 = 1e-6 * rd * k1 / g0;
c2 = 1e-6 * rd * k3prime;

%-----------------------------------------------------------------------
%        compute the acceleration at the mass center
%        of a vertical column of the atmosphere
%-----------------------------------------------------------------------
dgref = 1.0 - 2.66d-03*cos(2.0 *latitude) - 2.8e-07 * height;

%-----------------------------------------------------------------------
%               compute zenith hydrostatic delay
%-----------------------------------------------------------------------
tzd1 = c1 ./ dgref .* P;

%-----------------------------------------------------------------------
%               compute zenith wet delay
%-----------------------------------------------------------------------
tzd2 = c2 ./ (g0 .* dgref *(lambda + 1.0) - rd * t_lapse) .* E ./ T;

% end of tzdds87

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin nwmf2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function wmf = nwmf2(latitude, elev)
 
% new aen 930517 Routine to compute the new wmf2.0 mapping function which