unb3mr.m

it is a source code from mathworks

%
function [DRATE HZD DHMFDEL WZD DWMFDEL]=UNB3MR(LATRAD,HEIGHTM,DAYOYEAR,ELEVRAD)
%
%---|---1|0---|---2|0---|---3|0---|---4|0---|---5|0---|---6|0---|---7|0-
%=======================================================================
%
%  Subroutine UNB3MR
%  =================
%
%     This subroutine uses UNB3m model to calculate the slant
%     neutral atmospheric delay rate for a given latitude, height, 
%     day of year (DOY) and elevation angle. 
%
%   How this function works:
%   ------------------------
%
%     1. A look-up table is used to calculate mean sea level (MSL)
%        values for pressure, temperature, relative humidity (RH) 
%        and lapse rate parameters (for given latitude and DOY).
%
%     2. MSL water vapor pressure (WVP) is computed according to
%        IERS Conventions 2003.
%
%     3. Pressure, temperature and WVP values are computed for 
%        station height.
%
%     4. Zenith hydrostatic and non-hydrostatic delays are 
%        computed using Saastamoinen formulas as modified by
%        Davis et. al (1985).
%
%     5. Computation of Hydrostatic and Non-hydrostatic Niell 
%        mapping function rates. 
%
%     6. Total slant delay rate is determined using (4) and (5).
%
%   General Comments:
%   -----------------
%
%     1. This subroutine was adapted from subroutine UNB3M.
%
%     2. The internal variables for RH are called E0, EAVG and EAMP.
%
%   References: 
%   ----------
%
%     Leandro R.F., Santos, M.C. and Langley R.B., (2006). 
%     UNB Neutral Atmosphere Models: Development and Performance. 
%     Proceedings of ION NTM 2006, pp. 564-573, Monterey, California, 
%     January, 2006.
%
%     Davis, J.L., T.A. Herring, I.I. Shapiro, A.E.E. Rogers, 
%     and G. Elgered (1985). Geodesy by Radio Interfereometry: 
%     Effects of Atmospheric Modeling Errors on Estimates of 
%     Baseline Length. Radio Science, Vol. 20, No. 6, pp. 1593-1607.
%
%     IERS (2004). IERS Conventions (2003), eds. D.D. McCarthy and 
%     G. Petit, IERS Technical Note No. 32, International Earth Rotation 
%     and Reference Systems Service, Verlag des Bundesmates f黵 Kartographie 
%     und Geodasie, Frankfurt am Main.
%
%   Created/Modified
%   ----------------
%
%   DATE           WHO			   WHAT
%   ----           ---			   ----	
%   2006/01        Rodrigo Leandro          SUBROUTINE UNB3M.f created
%   2006/06/22     Rodrigo Leandro          SUBROUTINE UNB3MR.f created 
%                                           (addapted from UNB3M.f)
%   2006/11/23     Rodrigo Leandro          Code re-formatted. Capital case
%                                           should be always used in code 
%                                           from now on.
%
%
%   Input/Output
%   ------------
%
%     PARAMETER   I/O   TYPE       DESCRIPTION
%     ---------   ---   ----       -----------
%     LATRAD      IN    DOUBLE     Station geodetic latitude (radians)
%     HEIGHTM     IN    DOUBLE     Station orthometric height (m)
%     DAYOYEAR    IN    DOUBLE     Day of year
%     ELEVRAD     IN    DOUBLE     Elevation angle (radians)
%     HZD         OUT   DOUBLE     Hydrostatic zenith delay (m)
%     DHMFDEL     OUT   DOUBLE     Hydrostatic niell mapping function rate (rad^-1)
%     WZD         OUT   DOUBLE     Non-hyd. zenith delay (m)
%     DWMFDEL     OUT   DOUBLE     Non-hyd. Niell mapping function rate (rad^-1)
%     DRATE       OUT   DOUBLE     Total slant delay rate (m/rad)
%
%=======================================================================
%
%
%
%=======================================================================
% Initialize UNB3m look-up table
%-----------------------------------------------------------------------
AVG=[  15.0  1013.25  299.65  75.00  6.30  2.77
       30.0  1017.25  294.15  80.00  6.05  3.15
       45.0  1015.75  283.15  76.00  5.58  2.57
       60.0  1011.75  272.15  77.50  5.39  1.81
       75.0  1013.00  263.65  82.50  4.53  1.55];
%      lat      P       T      RH    beta lambda
%
AMP=[  15.0   0.00   0.00   0.00  0.00  0.00
       30.0  -3.75   7.00   0.00  0.25  0.33
       45.0  -2.25  11.00  -1.00  0.32  0.46
       60.0  -1.75  15.00  -2.50  0.81  0.74
       75.0  -0.50  14.50   2.50  0.62  0.30];
%      lat     P      T      RH    beta lambda
%=======================================================================
%
EXCEN2 = 6.6943799901413e-03;
MD     = 28.9644;
MW     = 18.0152;
K1     = 77.604;
K2     = 64.79;
K3     = 3.776e5;
R      = 8314.34;
C1     = 2.2768e-03;
K2PRIM = K2 - K1*(MW/MD);
RD     = R / MD;
DTR = 1.745329251994329e-02;
DOY2RAD=(0.31415926535897935601e01)*2/365.25;

%---|---1|0---|---2|0---|---3|0---|---4|0---|---5|0---|---6|0---|---7|0-

%=======================================================================
% Initialize NMF tables
%-----------------------------------------------------------------------
ABC_AVG=[ 15.0 1.2769934e-3 2.9153695e-3 62.610505e-3
          30.0 1.2683230e-3 2.9152299e-3 62.837393e-3
          45.0 1.2465397e-3 2.9288445e-3 63.721774e-3
          60.0 1.2196049e-3 2.9022565e-3 63.824265e-3
          75.0 1.2045996e-3 2.9024912e-3 64.258455e-3];
     
ABC_AMP=[ 15.0 0.0          0.0          0.0   
          30.0 1.2709626e-5 2.1414979e-5 9.0128400e-5 
          45.0 2.6523662e-5 3.0160779e-5 4.3497037e-5 
          60.0 3.4000452e-5 7.2562722e-5 84.795348e-5 
          75.0 4.1202191e-5 11.723375e-5 170.37206e-5];

A_HT = 2.53e-5;
B_HT= 5.49e-3;
C_HT = 1.14e-3;

HT_TOPCON = 1 + A_HT/(1 + B_HT/(1 + C_HT));

ABC_W2P0=[ 15.0 5.8021897e-4 1.4275268e-3 4.3472961e-2
           30.0 5.6794847e-4 1.5138625e-3 4.6729510e-2
           45.0 5.8118019e-4 1.4572752e-3 4.3908931e-2
           60.0 5.9727542e-4 1.5007428e-3 4.4626982e-2
           75.0 6.1641693e-4 1.7599082e-3 5.4736038e-2];
%=======================================================================
%
%
%=======================================================================
% Transform latitude from radians to decimal degrees
%-----------------------------------------------------------------------
LATDEG = LATRAD * 180.0 / pi;
%=======================================================================
%
%
%=======================================================================
% Deal with southern hemisphere and yearly variation
%-----------------------------------------------------------------------
TD_O_Y = DAYOYEAR;
if LATDEG<0
    TD_O_Y = TD_O_Y + 182.625;
end
COSPHS = cos((TD_O_Y - 28) * DOY2RAD );
%=======================================================================
%
%
%=======================================================================
% Initialize pointers to lookup table
%-----------------------------------------------------------------------
LAT = abs( LATDEG );
if LAT>=75
    P1 = 5;
    P2 = 5;
    M = 0;
elseif LAT<=15
    P1 = 1;
    P2 = 1;
    M = 0;
else
    P1 = fix((LAT - 15)/15.D0) + 1;
    P2 = P1 + 1;
    M = (LAT - AVG(P1,1) ) / ( AVG(P2,1) - AVG(P1,1) );