MATLAB Answers

Discretising a size distribution function and area under the curve

6 views (last 30 days)
Liam Holbeche-Smith
Liam Holbeche-Smith on 26 Feb 2021
Commented: Alan Stevens on 28 Feb 2021
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!!!

Answers (2)

Alan Stevens
Alan Stevens on 26 Feb 2021
  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');
  4 Comments
Alan Stevens
Alan Stevens on 27 Feb 2021
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

Sign in to comment.


Liam Holbeche-Smith
Liam Holbeche-Smith on 28 Feb 2021
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
  3 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!