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

37 views (last 30 days)

Show older comments

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!!!

Alan Stevens
on 26 Feb 2021

- You should have log(sigma)^2, not log(sigma^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');

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

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

Start Hunting!