function zd = SHAtrop1(fn, blh, mjd)
%  zd = SHAtrop1(fn, blh, mjd)  % subroutine for SHAtrop1 ZTD calculation
% 
% (c) GNSS Analysis Center at Shanghai Astronomical Observatory, 2017
%
% The copyright in this document is vested in theGNSS Analysis Center at 
%  Shanghai Astronomical Observatory (SHA). 
% 
% Contacts:
%   Professor. Dr.  Junping Chen 
%   http://202.127.29.4/shao_gnss_ac/
%
%
%   Input/Output
%   ----------------
%
%   Input:
%           fn       string
%                    file name for saved SHAtrop1 grid data
%           blh      double
%                    input site's latitude-longitude and ellipsoidal 
%                    height, in degree, degree and meter
%           mjd      Modified Julian Date
%                    could be a 1-D matrix, for several days
%
%   Output
%           zd        zenith total delay, in meter (-999 for error)
%
%   Created/Modified
%   ----------------
%
%      DATE          WHO                 WHAT
%      ----          ---                 ----	
%      Nov. 2017     Jungang Wang        Subroutine  created
% 
%
%   General Comments:
%   -----------------
%    1. For convenience, if the input coordinate is exactly at one grid
%       point, we make tiny modification and ensure the siplicity of the
%       program
%    2. If any error happens the return ZTD value will be -999 as default
%    3. The grid data file should be give, and usually should be put in the
%       same directory with this program
%    4. The SHAtrop grid data is declared as persistent parameter
%
%   References: 
%   ---------- 
%
% Test cases
% Example 1
% input:
% fn = 'SHAtrop1.mat'
% blh = [ 37.4720  101.4009    2.4506 ]
% mjd = 57754
%
% output:
%     2.3789
%  


zd = -999; 

persistent SHAtrop

if isempty( SHAtrop )  
    try
        load( fn )
    catch
        ['Input data file do not exist, ERROR!'] 
        return
    end
end
 
% check if the input coordinate is legal
if size(blh,1)*size(blh,2) ~= 3
    ['Input BLH is not for one sites, ERROR!'] 
    return
end

blhi = blh;


% Height correction coefficient, depending on the latitude
i1 = max( find( SHAtrop.hgt_interval <= blhi(1) ) );
coefi =  exp(  SHAtrop.hgt_coef(i1) * blhi(3) ) ;

% to avoid the case where the input site is exactly on of the grid point,
% and therefore make the following calculation complicated
if mod(blhi(1), SHAtrop.dlat) == 0; blhi(1) = blhi(1) + 1d-9; end
if mod(blhi(2), SHAtrop.dlon) == 0; blhi(2) = blhi(2) + 1d-9; end

% find the grid where the input site locating in
i1 = find( abs( SHAtrop.pts(:,1) - blhi(:,1) ) < SHAtrop.dlat ) ;
i2 = find( abs( SHAtrop.pts(:,2) - blhi(:,2) ) < SHAtrop.dlon ) ;
i0 = intersect( i1, i2 ) ;

% if the input site is not in one of the grid of this model, we do not
% ouput any ZTD value
if length(i0) ~= 4
    ['Site not in mainland China!']
    blh
    return
end
 
ndim = size( SHAtrop.grd,2 ) ;

% search for the right for the input site
dbtmp = SHAtrop.pts(i0,1) - blhi(:,1) ;
dltmp = SHAtrop.pts(i0,2) - blhi(:,2) ;
dx = abs( dbtmp .* dltmp ) * ones(1, ndim)  ;
dy = SHAtrop.grd( i0, : ) ;
ipt = [4 3 2 1];
dy = sum(dx .* dy(ipt,:)) / ( SHAtrop.dlat * SHAtrop.dlon ) ;

w1 = 2*pi/365.25;   w2 = 4*pi/365.25;

% parameters 
a = dy(1:5) ;
zd = a(1) + a(2).* cos( w1.*(mjd-a(3)) ) + a(4).* cos( w2.*(mjd-a(5)) ) ;


% height correction 
zd = zd * coefi ;




return 





end