Discretising a size distribution function and area under the curve

4 次查看（过去 30 天）
Liam Holbeche-Smith
评论： Alan Stevens 2021-2-28
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!!!

请先登录，再回答此问题。

回答（2 个）

Alan Stevens
Alan Stevens 2021-2-26
  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 个评论
显示 2更早的评论
Liam Holbeche-Smith
Liam Holbeche-Smith 2021-2-27
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
Alan Stevens 2021-2-27
编辑：Alan Stevens 2021-2-27
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

请先登录，再进行评论。

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
  3 个评论
显示 1更早的评论
Liam Holbeche-Smith
Liam Holbeche-Smith 2021-2-28
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by