function [zwd_o] = gzwd_bp(doy,dlat,dlon,h_ell,mod_bp,i_scale,o_scale,mod_type,...
    surf_e,surf_t)
%
% This script is used to determine a priori values for tropospheric zenith
% wet delays at the vacinity of the Earth's surface based on the GZWD_BP model. 
% 
% input parameters:
% doy           I   day of year, scalar
% dlat          I   latitude in [degree], vector
% dlon          I   longitude in [degree], vector
% h_ell         I   elliposidal height in [meter], vector
% mod_bp        I   coefficients of the BPNN-based ZWD model, structure
% i_scale       I   scale for the normalization of the input parameters, structure
% o_scale       I   scale for the normalization of the output parameters, structure
% mod_type      I   model type (1: gztd_bp_a; 2: gztd_bp_b; 3: gztd_bp_c; 4: gztd_rf_d)
% surf_e        I   surface water vapor pressure in [hPa], vector
% surf_t        I   surface temperature in [Kelvin], vector
% zwd_init      I   a priori values for zenith wet delays in [m], vector
%
% output parameters:
% zwd_o         O   zenith wet delays in [meter], vector
%
% (c) Faculty of Geosciences and Environmental Engineering, Southwest
% Jiaotong University, 611756, P.R.China
% (c) Department of Higher Geodesy, Vienna University of Techonology (TU Wien),
% Wien 1040, Austria
%
% Qinzheng Li
% version 1, February 1, 2022
% version 2, July 28, 2022, (Global ZWD modeling based on the combination of 
% radiosonde measurements and GPS radio occultations)

if mod_type==1
    narginchk(8,10);
elseif mod_type==2
    narginchk(9,10);
elseif mod_type==3
    narginchk(10,10);
elseif mod_type==4
    narginchk(9,10);
end

doy_i = rem(doy,365.25);
if length(doy_i)==1
    if doy_i<0 || 366<doy_i % DOY setting
        error('Day of year should be at the range of [0 366]');
    end
elseif 1<length(doy_i)
    if doy_i<0 | 366<doy_i % DOY setting
        error('Day of year should be at the range of [0 366]');
    end
end

sn = length(dlat);
if length(dlon)~=sn || length(h_ell)~=sn 
    error('The length of dlat, dlon, h_ell should be consistent!');
end

if mod_type==2 && isempty(surf_e)
    error('Site-wise surface water vapor pressure data are required:');
elseif mod_type==3 && (isempty(surf_e) || isempty(surf_t))
    error('Site-wise surface water vapor pressure and surface temperature data are required:');
end

if mod_type==2 && length(surf_e)~=sn
    error('The length of dlat and surf_e should be consistent!');
elseif mod_type==3 && (length(surf_e)~=sn || length(surf_t)~=sn)
    error('The length of dlat, surf_e, and surf_t should be consistent!');
end

% Initialization
zwd_o = NaN(sn,1);
if size(dlat,1)==1
    dlat = dlat';
end
if size(dlon,1)==1
    dlon = dlon';
end
if size(h_ell,1)==1
    h_ell = h_ell';
end

% Constrcut input parameters matrix
ind_lon = find(-180<dlon & dlon<0);
if ~isempty(ind_lon)
    dlon(ind_dlon) = dlon(ind_dlon) + 360;
end

switch (mod_type)
    case {1}
        if length(doy_i)==1
            inp_data = [ones(sn,1)*doy_i,dlat,dlon,h_ell];
        elseif 1<length(doy_i)
            inp_data = [doy_i,dlat,dlon,h_ell];
        end
    case {2,4}
        if size(surf_e,1)==1
            surf_e = surf_e';
        end
        if length(doy_i)==1
            inp_data = [ones(sn,1)*doy_i,dlat,dlon,h_ell,surf_e];
        elseif 1<length(doy_i)
            inp_data = [doy_i,dlat,dlon,h_ell,surf_e];
        end
    case {3}
        if size(surf_e,1)==1
            surf_e = surf_e';
        end
        if size(surf_t,1)==1
            surf_t = surf_t';
        end
        
        if length(doy_i)==1
            inp_data = [ones(sn,1)*doy_i,dlat,dlon,h_ell,surf_e,surf_t];
        elseif 1<length(doy_i)
            inp_data = [doy_i,dlat,dlon,h_ell,surf_e,surf_t];
        end
end

% Normalization
[paras,~] = mapstd('apply',inp_data',i_scale);

% Predict the results using trained network
zwd_o(1:sn,1) = sim(mod_bp,paras);

% Normalization reverse
zwd_o = mapstd('reverse',zwd_o,o_scale);

clear inp_data;
end



        
    



