Discretising a size distribution function and area under the curve

I have some parameters for some data and I am plotting a size distribution function. The probability function I am using is given below:

Currently I have the following code which plots the function:

clear
clc
close all

mu = 0.015;  % geometric mean radius
sigma = 1.6; % geometric standard deviation
Ntot = 850;  % total number concentration
nbins = 200; % number of bins

rMin = (mu/(10*sigma)); % determine the minimum radius
rMax = (mu*10*sigma);   % determine the maximum radius
rs = logspace(log10(rMin),log10(rMax),nbins+1); % vector containing rs

for i = 1:length(rs)
    sizedist = Ntot/((sqrt(2*pi))*log(sigma))*exp(-(log(rs./mu).^2)./(2*log(sigma^2)));
end

% plot the size distribution function 
semilogx(rs,sizedist,'Linewidth',2);
xlim([10e-4 10e1]);
xlabel('Particle Radius, \mum');
ylabel('Number Concentration, cm^3');

The first problem I have occurs here. Shouldn't the area under the curve be equal to the total number concentration specified? As this doesn't seem to be the case.

For the model I am coding it is necessary for the aerosol population to be discretised into bins.

Now suppose I wish to discretise this function into bins, currently the vector rs is the bin edges. I would like to know how to find the number concentration in each bin, but as stated above the area under the curve doesn't seem to equal Ntot. Can anyone see whether i'm doing anything wrong?

What I would like to do is to determine the Number Concentration in each bin, which when added together should total the value of Ntot.

Hopefully someone can help!!!

Answer 1

[Alan Stevens]

1. You should have log(sigma)^2, not log(sigma^2).
2. Don't forget the "dx" part when integrating the curve.

mu = 0.015;  % geometric mean radius
sigma = 1.6; % geometric standard deviation
Ntot = 850;  % total number concentration
nbins = 200; % number of bins

rMin = (mu/(10*sigma)); % determine the minimum radius
rMax = (mu*10*sigma);   % determine the maximum radius
rs = logspace(log10(rMin),log10(rMax),nbins+1); % vector containing rs

dx = (log(rMax)-log(rMin))/nbins;
sizedist = Ntot/(sqrt(2*pi)*log(sigma))*exp(-log(rs./mu).^2./(2*log(sigma)^2));
    
disp(sum(sizedist)*dx)   % This should display the area under the curve

% plot the size distribution function 
semilogx(rs,sizedist,'Linewidth',2);
xlim([10e-4 10e1]);
xlabel('Particle Radius, \mum');
ylabel('Number Concentration, cm^3');

[Liam Holbeche-Smith]

Thankyou Alan I hadn't realised the typo.

Suppose I wished to split the area under the curve into the number of bins given by nbins. Is there away to find the area in each bin to determine the number concentration within it? Then assign this value to the midpoint in each bin?

I normally assume the bin centres would be given by:

a = rs(1:end-1);
b = rs(2:end);
rmid = (b+a)/2; % value in the middle of each bin

But i'm not sure if this is still valid for the lognormal spacing?

[Alan Stevens]

Your "bin" values are exactly at the values of rs, because they are where you calculate the values of sizedist. It would be reasoable to assume that they occur in the centre of bins of width dx (see my listing for dx).

[Liam Holbeche-Smith]

Thanks Alan,

Do you know if there a way to visualise this as a histogram, or would it not work due to the lognormal axis?

[Alan Stevens]

You could do the following

mu = 0.015;  % geometric mean radius
sigma = 1.6; % geometric standard deviation
Ntot = 850;  % total number concentration
nbins = 200; % number of bins
rMin = (mu/(10*sigma)); % determine the minimum radius
rMax = (mu*10*sigma);   % determine the maximum radius
rs = logspace(log10(rMin),log10(rMax),nbins+1); % vector containing rs
dx = (log(rMax)-log(rMin))/nbins;
sizedist = Ntot/(sqrt(2*pi)*log(sigma))*exp(-log(rs./mu).^2./(2*log(sigma)^2));
    
disp(sum(sizedist)*dx)   % This should display the area under the curve

nbins2 = 21;
rs2 = logspace(log10(rMin),log10(rMax),nbins2+1); % vector containing rs2
sizedist2 = Ntot/(sqrt(2*pi)*log(sigma))*exp(-log(rs2./mu).^2./(2*log(sigma)^2));
x = []; y = []; bw = rs2(1)/2;

for i = 1:nbins2
   x = [x rs2(i) rs2(i) rs2(i+1) rs2(i+1)];
   yav = (sizedist2(i+1)+sizedist2(i))/2;
   y = [y 0 yav yav 0];
end

% plot the size distribution function 
semilogx(rs,sizedist,'Linewidth',2);
hold on
semilogx(x,y)
xlabel('Particle Radius, \mum');
ylabel('Number Concentration, cm^3');

which produces

Answer 2

[Liam Holbeche-Smith]

As an update after some trial and error and some reading up on the subject I have now come up with some code that creates the correct plot and returns the correct vectors:

mu = 0.015;  % this is the geometric mean radius
sigma = 1.6; % this is the geometric standard deviation
Ntot = 850;  % this is the total number concentration of the aerosol 
nr = 200;    % this is the number of radii being tracked

rMin = (mu/(10*sigma)); % minimum dry radius
rMax = (mu*10*sigma);   % maximum dry radius
rvec = logspace(log10(rMin),log10(rMax),nr+1); % these are the bin edges

a = rvec(1:end-1); % lower bound in each bin
b = rvec(2:end);   % upper bound in each bin

pdfa = (Ntot/(sqrt(2.0*pi)*log(sigma))./a).*exp(-log(a./mu).^2./(2*log(sigma)^2));
pdfb = (Ntot/(sqrt(2.0*pi)*log(sigma))./b).*exp(-log(b./mu).^2./(2*log(sigma)^2));

Nconcs = 0.5*(b-a).*(pdfa+pdfb); % number concentration in each bin
rds = sqrt(a.*b);                % radius value in centre of each bin

semilogx(rds,Nconcs,'k'); % plot the aerosol population
xlim([10e-4 10e1]);
xlabel('Aerosol Dry Radius, \mum');
ylabel('Aerosol Number Concentration, cm^3');
hold on;
area(rds,Nconcs,'Facecolor','#CC0066');
legend('Aerosol1','Fontsize',12);

Suppose I had another aerosol population that I wish to plot on the same graph, with different characteristics. If I made mu, sigma, Ntot and nr into arrays is there a way to achieve this? As an example lets say my starting arrays are:

mus = [0.015 0.85];
sigmas = [1.6 1.2];
Ntots = [850 10];
nr = [200 40];

I imagine using a for loop is best, i'm quite new to MATLAB so i'm unsure if the complexity of the calculations will make this too difficult to implement.

Once again, thankyou for your help

[Alan Stevens]

The following loop structure does it:

mus = [0.015 0.85];
sigmas = [1.6 1.2];
Ntots = [850 10];
nrs = [200 40];

for i = 1:numel(mus)

mu = mus(i);  % this is the geometric mean radius
sigma = sigmas(i); % this is the geometric standard deviation
Ntot = Ntots(i);  % this is the total number concentration of the aerosol 
nr = nrs(i);    % this is the number of radii being tracked

rMin = (mu/(10*sigma)); % minimum dry radius
rMax = (mu*10*sigma);   % maximum dry radius
rvec = logspace(log10(rMin),log10(rMax),nr+1); % these are the bin edges

a = rvec(1:end-1); % lower bound in each bin
b = rvec(2:end);   % upper bound in each bin

pdfa = (Ntot/(sqrt(2.0*pi)*log(sigma))./a).*exp(-log(a./mu).^2./(2*log(sigma)^2));
pdfb = (Ntot/(sqrt(2.0*pi)*log(sigma))./b).*exp(-log(b./mu).^2./(2*log(sigma)^2));

Nconcs = 0.5*(b-a).*(pdfa+pdfb); % number concentration in each bin
rds = sqrt(a.*b);                % radius value in centre of each bin

figure
semilogx(rds,Nconcs,'k'); % plot the aerosol population
xlim([10e-4 10e1]);
xlabel('Aerosol Dry Radius, \mum');
ylabel('Aerosol Number Concentration, cm^3');
hold on;
area(rds,Nconcs,'Facecolor','#CC0066');
legend('Aerosol1','Fontsize',12);
hold off

end

[Liam Holbeche-Smith]

This seems to produce two seperate figures, I was hoping to be able to have both plots on the same graph if at all possible?

[Alan Stevens]

In that case remove the words "figure" and "hold off".