Discretising a size distribution function and area under the curve

8 ビュー (過去 30 日間)
Liam Holbeche-Smith
Liam Holbeche-Smith 2021 年 2 月 26 日
コメント済み: 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 件のコメント
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
Liam Holbeche-Smith 2021 年 2 月 28 日
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 件のコメント
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?
Alan Stevens
Alan Stevens 2021 年 2 月 28 日
In that case remove the words "figure" and "hold off".

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeMATLAB についてさらに検索

製品


リリース

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by