%% 获取大气质量相关信息
function wind = GetWindData_ERA5(file_wind, Facility)
wind.u = ncread(file_wind,'u');  % in m s-1  东风
wind.v = ncread(file_wind,'v');  % in m s-1  北风

wind.lon = ncread(file_wind,'longitude');
wind.lat = ncread(file_wind,'latitude');
wind.time = ncread(file_wind,'time');
wind.time = datetime('1900-01-01 00:00:00') + hours(wind.time);    % 将小时数转换为日期时间
% 寻找最近邻元
% lon_ind = find_nearest_index(wind.lon, Facility.lon);
% lat_ind = find_nearest_index(wind.lat, Facility.lat);
% time_ind = find_nearest_index(wind.time, Facility.satellite_time);
% 插值
X = double(wind.lon);
Y = double(wind.lat);
Z = datenum(wind.time);
U = wind.u;
V = wind.v;
Xq = Facility.lon;
Yq = Facility.lat;
Zq = datenum(Facility.time); % satellite_time是卫星时间（部分为UTC+8）
if(length(Z)>1)
    wind.Uq = interp3(X,Y,Z,U,Xq,Yq,Zq,'nearest');
    wind.Vq = interp3(X,Y,Z,V,Xq,Yq,Zq,'nearest');
    wind.Uq2 = interp3(X,Y,Z,U,Xq,Yq,Zq,'linear');
    wind.Vq2 = interp3(X,Y,Z,V,Xq,Yq,Zq,'linear');
    wind.Uq3 = interp3(X,Y,Z,U,Xq,Yq,Zq,'spline');
    wind.Vq3 = interp3(X,Y,Z,V,Xq,Yq,Zq,'spline');
    wind.Uq4 = interp3(X,Y,Z,U,Xq,Yq,Zq,'makima');
    wind.Vq4 = interp3(X,Y,Z,V,Xq,Yq,Zq,'makima');
    wind.Uq5 = interp3(X,Y,Z,U,Xq,Yq,Zq,'cubic');
    wind.Vq5 = interp3(X,Y,Z,V,Xq,Yq,Zq,'cubic');
else
    wind.Uq = NaN;
    wind.Vq = NaN;
end
% figure
% [X,Y,Z] = meshgrid(X,Y,Z)
% slice(X,Y,Z,V,Xq,Yq,Zq);
% shading flat

if(isnan(wind.Vq) || isnan(wind.Uq)) %
    disp('请检查VMA数据一致性')
    wind.interpdata = wind.data;
else
    % 计算实际风速和有效风速
    wind.speed = sqrt(wind.Uq^2 + wind.Vq^2); % 风速
    wind.speed2 = sqrt(wind.Uq2^2 + wind.Vq2^2); % 风速
    wind.speed3 = sqrt(wind.Uq3^2 + wind.Vq3^2); % 风速
    wind.speed4 = sqrt(wind.Uq4^2 + wind.Vq4^2); % 风速
    wind.direction = atan2d(wind.Vq, wind.Uq); % 风向
    wind.direction = mod((wind.direction + 360), 360);
    % Cusworth,2023
    q1 = 1.1;
    q2 = 0.6;
    wind.Ueff1 = q1 * log(wind.speed) + q2; % 单位：m/s
    wind.Ueff2 = q1 * log(wind.speed2) + q2; % 单位：m/s
    wind.Ueff3 = q1 * log(wind.speed3) + q2; % 单位：m/s
    wind.Ueff4 = q1 * log(wind.speed4) + q2; % 单位：m/s
    Ueff_list = [wind.Ueff1, wind.Ueff2, wind.Ueff3, wind.Ueff4];
    wind.Ueff = mean(Ueff_list);
    wind.Ueff_uncertainty = std(Ueff_list); % 有效风速不确定度
end
end