unb3mr.m

it is a source code from mathworks

end
%=======================================================================
%
%
%=======================================================================
% Compute average surface tropo values by interpolation
%-----------------------------------------------------------------------
PAVG = M * ( AVG(P2,2) - AVG(P1,2) ) + AVG(P1,2);
TAVG = M * ( AVG(P2,3) - AVG(P1,3) ) + AVG(P1,3);
EAVG = M * ( AVG(P2,4) - AVG(P1,4) ) + AVG(P1,4);
BETAAVG   = M * ( AVG(P2,5) - AVG(P1,5) ) + AVG(P1,5);
LAMBDAAVG = M * ( AVG(P2,6) - AVG(P1,6) ) + AVG(P1,6);
%=======================================================================
%
%
%=======================================================================
% Compute variation of average surface tropo values
%-----------------------------------------------------------------------
PAMP = M * ( AMP(P2,2) - AMP(P1,2) ) + AMP(P1,2);
TAMP = M * ( AMP(P2,3) - AMP(P1,3) ) + AMP(P1,3);
EAMP = M * ( AMP(P2,4) - AMP(P1,4) ) + AMP(P1,4);
BETAAMP   = M * ( AMP(P2,5) - AMP(P1,5) ) + AMP(P1,5);
LAMBDAAMP = M * ( AMP(P2,6) - AMP(P1,6) ) + AMP(P1,6);
%=======================================================================
%
%
%=======================================================================
% Compute surface tropo values
%-----------------------------------------------------------------------
P0 = PAVG - PAMP * COSPHS;
T0 = TAVG - TAMP * COSPHS;
E0 = EAVG - EAMP * COSPHS;
BETA = BETAAVG - BETAAMP * COSPHS;
BETA   = BETA / 1000.D0;
LAMBDA = LAMBDAAVG - LAMBDAAMP * COSPHS;
%=======================================================================
%
%
%=======================================================================
% Transform from relative humidity to WVP (IERS Conventions 2003)
%-----------------------------------------------------------------------
ES = 0.01 * exp(1.2378847e-5 * (T0 ^ 2) - 1.9121316e-2 * ...
    T0 + 3.393711047e1 - 6.3431645e3 * (T0 ^ -1));
FW = 1.00062 + 3.14e-6 * P0 + 5.6e-7 * ((T0 - 273.15) ^ 2);
E0 = (E0 / 1.00e2) * ES * FW;
%=======================================================================
%
%
%=======================================================================
% Compute power value for pressure & water vapour
%-----------------------------------------------------------------------
EP = 9.80665 / 287.054 / BETA;
%=======================================================================
%
%
%=======================================================================
% Scale surface values to required height
%-----------------------------------------------------------------------
T = T0 - BETA * HEIGHTM;
P = P0 * ( T / T0 ) ^ EP;
E = E0 * ( T / T0 ) ^ ( EP * (LAMBDA+1) );
%=======================================================================
%
%
%=======================================================================
% Compute the acceleration at the mass center
%  of a vertical column of the atmosphere
%-----------------------------------------------------------------------
GEOLAT = atan((1.0-EXCEN2)*tan(LATRAD));
DGREF = 1.0 - 2.66e-03*cos(2.0*GEOLAT) - 2.8e-07*HEIGHTM;
GM    = 9.784 * DGREF;
DEN   = ( LAMBDA + 1.0 ) * GM;
%=======================================================================
%
%
%=======================================================================
% Compute mean temperature of the water vapor
%-----------------------------------------------------------------------
TM  = T * (1 - BETA * RD / DEN);
%=======================================================================
%
%
%=======================================================================
% Compute zenith hydrostatic delay
%-----------------------------------------------------------------------
HZD = C1 / DGREF * P;
%=======================================================================
%
%
%=======================================================================
% Compute zenith wet delay
%-----------------------------------------------------------------------
WZD = 1.0e-6 * ( K2PRIM + K3/TM) * RD * E/DEN;
%=======================================================================
%
%
%=======================================================================
% Compute average NMF(H) coefficient values by interpolation
%-----------------------------------------------------------------------
A_AVG = M * ( ABC_AVG(P2,2) - ABC_AVG(P1,2) ) + ABC_AVG(P1,2);
B_AVG = M * ( ABC_AVG(P2,3) - ABC_AVG(P1,3) ) + ABC_AVG(P1,3);
C_AVG = M * ( ABC_AVG(P2,4) - ABC_AVG(P1,4) ) + ABC_AVG(P1,4);
%=======================================================================
%
%
%=======================================================================
% Compute variation of average NMF(H) coefficient values
%-----------------------------------------------------------------------
A_AMP = M * ( ABC_AMP(P2,2) - ABC_AMP(P1,2) ) + ABC_AMP(P1,2);
B_AMP = M * ( ABC_AMP(P2,3) - ABC_AMP(P1,3) ) + ABC_AMP(P1,3);
C_AMP = M * ( ABC_AMP(P2,4) - ABC_AMP(P1,4) ) + ABC_AMP(P1,4);
%=======================================================================
%
%
%=======================================================================
% Compute NMF(H) coefficient values
%-----------------------------------------------------------------------
A = A_AVG - A_AMP * COSPHS;
B = B_AVG - B_AMP * COSPHS;
C = C_AVG - C_AMP * COSPHS;
%=======================================================================
%
%
%=======================================================================
% Compute sine of elevation angle
%-----------------------------------------------------------------------
SINE = sin(ELEVRAD);
DSINEDEL = cos(ELEVRAD);
%=======================================================================
%
%
%=======================================================================
% Compute NMF(H) value
%-----------------------------------------------------------------------
ALPHA  = B/(SINE + C );
DALPHADEL = (-B/((SINE+C)^2))*DSINEDEL;
GAMMA  = A/(SINE + ALPHA);
DGAMMADEL = (-A/((SINE+ALPHA)^2))*(DSINEDEL+DALPHADEL);
TOPCON = (1 + A/(1 + B/(1 + C)));
DHMFDEL = (-TOPCON/((SINE+GAMMA)^2))*(DSINEDEL+DGAMMADEL);
%=======================================================================
%
%
%=======================================================================
% Compute and apply height correction
%-----------------------------------------------------------------------
ALPHA  = B_HT/( SINE + C_HT );
DALPHADEL = (-B_HT/((SINE+C_HT)^2))*DSINEDEL;
GAMMA  = A_HT/( SINE + ALPHA);
DGAMMADEL = (-A_HT/((SINE+ALPHA)^2))*(DSINEDEL+DALPHADEL);
DHCCDEL = ((HT_TOPCON/((SINE+GAMMA)^2))-1/(SINE^2))*DSINEDEL+ ...
    (HT_TOPCON/((SINE+GAMMA)^2))*DGAMMADEL;
DHT_CORRDEL = (HEIGHTM / 1000)*DHCCDEL;
DHMFDEL = DHMFDEL+DHT_CORRDEL;
%=======================================================================
%
%
%=======================================================================
% Comput average NMF(W) coefficient values by interpolation
%-----------------------------------------------------------------------
A = M * ( ABC_W2P0(P2,2) - ABC_W2P0(P1,2) ) + ABC_W2P0(P1,2);
B = M * ( ABC_W2P0(P2,3) - ABC_W2P0(P1,3) ) + ABC_W2P0(P1,3);
C = M * ( ABC_W2P0(P2,4) - ABC_W2P0(P1,4) ) + ABC_W2P0(P1,4);
%=======================================================================
%
%
%=======================================================================
% Compute NMF(W) value
%-----------------------------------------------------------------------
ALPHA = B/( SINE + C );
DALPHADEL=(-B/((SINE+C)^2))*DSINEDEL;
GAMMA = A/( SINE + ALPHA);
DGAMMADEL=(-A/((SINE+ALPHA)^2))*(DSINEDEL+DALPHADEL);
TOPCON = (1 + A/(1 + B/(1 + C)));
DWMFDEL=(-TOPCON/((SINE+GAMMA)^2))*(DSINEDEL+DGAMMADEL);
%=======================================================================
%
%
%=======================================================================
% Compute total slant delay rate
%-----------------------------------------------------------------------
DRATE=HZD*DHMFDEL+WZD*DWMFDEL;
%=======================================================================
%
%
%=======================================================================
% End of function
%=======================================================================