function Fig9(wp,radar,disdro,startDate,endDate)
% wp: Windprofiler.mat
% radar: Remote_sensing.nc (output from readRemoteSensing.m)
% disdro: Disdrometer.mat
% startDate/endDate: datenum(yyyy,mm,dd,HH,MM,SS)

% time index
idx_radar = radar.time >= startDate & radar.time <= endDate;
idx_disdro = disdro.time>= startDate & disdro.time<= endDate;
idx_wp = wp.t >= startDate & wp.t <= endDate;
time_radar = radar.time(idx_radar);
time_disdro = disdro.time(idx_disdro);
time_wp = reshape(wp.t(idx_wp),[length(unique(wp.z)) length(unique(wp.t(idx_wp)))]);
speed_wp = transpose(reshape(wp.speed(idx_wp),[length(unique(wp.z)) length(unique(wp.t(idx_wp)))]));
z_wp = unique(wp.z);

% shear
for i = 2:8
    shear_wp(:,i-1) = (speed_wp(:,i) - speed_wp(:,i-1))./((z_wp(i) - z_wp(i-1)).*1e-3);
    h_wp(i-1) = (z_wp(i) + z_wp(i-1))./2;
end

% radar
Z = radar.data.Z(:,idx_radar);
Z(Z<-100) = NaN;

% disdro
precip = disdro.data.rain_intensity(idx_disdro);
n = 3;
disdro_mean = arrayfun(@(i) nanmean(precip(i:i+n-1)),1:n:length(precip)-n+1)'; % the averaged vector
time_disdro_mean = arrayfun(@(i) nanmean(time_disdro(i:i+n-1)),1:n:length(time_disdro)-n+1)'; % the averaged vector

% figure
fig = figure;
fig.Units = 'centimeters';
fig.Position = [0 0 12.5 7.2];
fig.PaperPositionMode = 'auto';
pos = fig.Position;

set(fig,'PaperPositionMode','Auto','PaperUnits','Centimeters','PaperSize',[pos(3), pos(4)])
set(gca,'color','none');
ax1 = gca; % current axes
ax1.Position = [0.1 0.13 0.68 0.85];
FS = 8;

yyaxis left;
Z(Z<-60) = NaN;
p = pcolor(time_radar,radar.height,Z);
p.EdgeColor = 'None';
set(p,'facealpha',0.7)
cmap = viridis(20);
colormap(cmap);
hold on;
[X1,Y1] = meshgrid(time_wp(1,:),h_wp);
Z1 = transpose(shear_wp);
[M,p1] = contour(X1,Y1,Z1,'ShowText','on','Color','k','LineWidth',0.8);
clabel(M,p1,'FontSize',FS);
xlim([startDate endDate]);
set(gca,'Xtick',startDate:datenum(0,0,0,1,0,0):endDate);
datetick('x','HH:MM','keeplimits','keepticks');
hold on;
cbar = colorbar;
x1 = get(gca,'Position');
cbar.Position(3) = 0.5.*cbar.Position(3);
cbar.Position(1) = x1(3) + 0.21;
set(gca,'Position',x1);
caxis([-40 10]);
ylabel(cbar','Reflectivity (dBZ)','Fontsize',FS);
ylim([1600 3600]);
set(gca,'YTick',1600:400:3600);
ax1 = gca;
ax1.YAxis(1).Color = 'k';
ax1.XAxis(1).Color = 'k';
xlabel('Time (UTC)');
ylabel('Altitude (m)');
hold on;
yyaxis right;
c3 = [0.05 0.17 0.55];
p4 = bar(time_disdro_mean,disdro_mean,'FaceColor',c3);
ax2 = gca;
ax2.YAxis(2).Color = c3;
ylim([0 1]);
hold on;
ylabel('Precipitation rate (mm h^{-1})');
box on;
set(gca,'Fontsize',8,'YTick',0:0.1:0.2);
end