function [ldistvec, lEnergymeds,  matched_Bmeds] = quantiles_B_d_lE(ldist,matched_B,lEnergy,index,a,b,c)
% Function to identify the magnetic field data points within each distance
% bin associated with four different lightning energies, lE1 < quantile
% 0.25, quantile 0.25 < lE2 < median, median < lE3 < quantile 0.75, and lE4
% > quantile 0.75
%lEnergy{index}(:,2)=nan(length(lEnergy{index}),1);
%figure;
hold on;
ldistb=discretize(ldist{index},a:b:c);
for i = 1:length([a:b:c])-1
%ax(i)=subplot(3,4,i);
ldistb_indices = find(ldistb==i);
if isempty(ldistb_indices)
    lEnergymeds(1:4,i) = nan(4,1);
    matched_Bmeds(1:4,i) = nan(4,1);
    continue
end
lEnergyd = lEnergy{index}(ldistb_indices);
lEnergyd_quantindex = discretize(lEnergyd,[0 quantile(lEnergyd,0.25) quantile(lEnergyd,0.5) quantile(lEnergyd,0.75) max(lEnergyd)]);
scatcols = ['b','r','y','m'];
for j = 1:length([0 quantile(lEnergyd,0.25) quantile(lEnergyd,0.5) quantile(lEnergyd,0.75) max(lEnergyd)])-1
% find the median lEnergy value and associated B and ldist point for each
% of the four lEnergy bins in each ldist bin
lEnergymeds(j,i) = median(lEnergy{index}(ldistb_indices(lEnergyd_quantindex == j)));
matched_Bmeds(j,i) = median(abs(matched_B{index}(ldistb_indices(lEnergyd_quantindex == j))));
hs=scatter(ldist{index}(ldistb_indices(lEnergyd_quantindex == j)),abs(matched_B{index}(ldistb_indices(lEnergyd_quantindex == j))),30,scatcols(j),'.');
end

end
ldistvec = a+b/2:b:c-b/2;
end