tropo_mh.m

这是国外关于卫星导航方面一书的源代码

function [dr]=tropo_mh(kosinus_z,t,p,h)
% -----------------------------------------------------------------------------
% GPSLab:   TROPO_MH.M
%
% Correction of the tropospheric refraktion following the model of
% Modified Hopfield (Goad & Goodman)
%
% How: [dr]=tropo_mh(kosinus_z,t,p,h);
% Input:   kosinus_z Cosine (zenithal angle) from the direction cosines Geocenter to Topocenter
%                    and Topocenter to Satellite, corrected by aberration (earth rotation)
%          t         Temperature in 癈 following the Standard Atmosphere
%          p         Absolute Pressure in mbar following the Standard Atmosphere
%          h         Relative Humidity in Percent following the Standard Atmosphere
%
% Output:  dr [m]    Tropospheric Correction of the measurement
%                    between receiver and satellite
%
% Standard atmosphere:   like given in the calling program
%
% History:
% -----------
% Since 1969 Helen S. Hopfield worked on the basis of the deterministic
% approach of Smith & Weintraub to derive the formulae, which describe the
% tropospheric distance error (zenithal part) due to refraktion of radio waves.
% Succeeding models are almost just modifications or simlifications of 
% her approach, or refined it - as done by Goad & Goodman in 1974 resulting in 
% the set of formulae implemented here.
%
% Besides the standard atmosphere (T,p,H) one has also to set the altitudes, where the
% impact of each the dry and the wet component (hd, hw) gets zero.
% H. Hopfield only gave course values as recommendations.
% Improved values presented Kaniuth (DGFI Munich) in the 80ies. 
% As a standard one should use the values given in the common publications.
% A selection can be found in the source code of this M-file.
%
% Literature: Like mentioned in the source code below.
%
% -----------------------------------------------------------------------------
% (c) iapg, tu m黱chen, zebhauser, 1996-07-01, 1996-07-23, 1999-10-11

%
% not yet (!) vectorized regarding the parameter kosinus_z !

% kosinus_z: Scalar or Vector, do not input a Matrix

if min(size(kosinus_z))>1
error('ERROR T1: Input parameter of function TROPO_MH is no scalar or vector');
end

% kosinus_z has to be a lying Vector
kosinus_z=kosinus_z(:)';

% Earth radius following WGS84
RE=6378137;

% Partial pressure of the steam
e=h./100.*exp(-37.2465+0.213166.*(t+273.15)-0.000256908.*(t+273.15).^2);

% Dry part refractivity index
Nd=77.64.*p./(t+273.15);

% Wet part refractivity index
Nw=-12.96.*e./(t+273.15)+3.718e5.*e./(t+273.15)./(t+273.15);

% Height of null impact of the dry component [m]
% hd=40152+147.95.*t;           % following Helen Hopfield
% hd=41640+148.00.*t;           % following Kaniuth
% hd=40082+148.98.*t;           % following Hopfield in Akhundov & Stotskii (1992)
% hd=40126+148.72.*t;           % following Hopfield in Klein (1990)
% hd=5e6*0.002277.*p./Nd;       % following Hopfield in Leick (1990)
%                                 Bei 20/1013/50: 40000+148*t m
hd=40136+148.72.*t;		% following Hopfield in Hofmann-Wellenhof et al. (1992)

% Height of null impact of the wet component [m]
% hw=11462+62.*t;			% following Kaniuth
% B=48/180*pi;                  % .........for the model:
% hw=10226+78.*t+94.*cos(B);    % following Kaniuth & Tremel (1992), only if latitude B will also be provided
% hw=5e6*0.002277.*e./Nw.*(1255./(t+273.15)+0.5);       % following Hopfield in Leick (1990)
%                                 Bei 20/1013/50: 11464+62*t m
hw=11000;			% following Hopfield in Hofmann-Wellenhof et al. (1992)

% Calculation of cosine(E) and sine(E)
coseq=1-kosinus_z.^2;
cose =sqrt(coseq);
sine =kosinus_z;

% Model coefficients
ad=-sine./hd;
bd=-coseq./(2.*RE.*hd);
aw=-sine./hw;
bw=-coseq./(2.*RE.*hw);

% ############################# from here on: to be vectorized (matrices!) #######################
aad=ones(size(kosinus_z));
aaw=aad;

aad=[aad;4.*ad;6.*ad.*ad+4.*bd;4.*ad.*(ad.*ad+3.*bd);ad.^4+12.*ad.^2.*bd+6.*bd.^2];
aad=[aad;4.*ad.*bd.*(ad.^2+3.*bd);bd.^2.*(6.*ad.^2+4.*bd);4.*ad.*bd.^3;bd.^4];

aaw=[aaw;4.*aw;6.*aw.*aw+4.*bw;4.*aw.*(aw.*aw+3.*bw);aw.^4+12.*aw.^2.*bw+6.*bw.^2];
aaw=[aaw;4.*aw.*bw.*(aw.^2+3.*bw);bw.^2.*(6.*aw.^2+4.*bw);4.*aw.*bw.^3;bw.^4];

% Distance to the locations, where Nd or Nw becomes 0 
rd=sqrt((RE+hd).^2-(RE^2).*coseq)-RE.*sine;
rw=sqrt((RE+hw).^2-(RE^2).*coseq)-RE.*sine;

k0=(1:9)';
k = k0*ones(size(kosinus_z));

rdk = ones(9,1)*rd;
rwk = ones(9,1)*rw;

% Dry part
delta_d=1e-6.*Nd.*sum(aad./k.*rdk.^k);

% Wet part
delta_w=1e-6.*Nw.*sum(aaw./k.*rwk.^k);

% The total tropospheric correktion at the actual zenithal distance
dr=delta_d+delta_w;

% delta_d'
% delta_w'

% -----------------------------------------------------------------------------
% (c) iapg, zebhauser, 1996-07-01, 1996-07-23, 1999-10-11
% -----------------------------------------------------------------------------