在MATLAB中复制一篇文章的图

``````% implementing equation 9 and figure 4
step   = 0.01;    t = 1:step:3600;
d      = 3;     % dimension
N      = 8000;  % number of molecules
H      = 0.01;  % H = [0.01,0.1,1] is in mol/micrometer^3
H      = H*6.02214078^5; % hence I scaled the Avogadro's number (right or wrong?)
D      = 10;    % diffusion coefficient in micrometer^2/sec

u(1)   =  1./(1.^(d/2)); % inner function in equation 9; first pulse

for i = 2:numel(t)/1000
u(i)     =  u(i-1)+(1./(i.^(d/2))); %  u-> the pulse number
lmda(i)  = (1/(4*pi*D))*((N/(H)).*sum(u)).^(2/d);
end

figure;plot(lmda)
``````

ps：这不是作业。

``````H = H*6.02214078^5;
``````

``````H = H * 6.02214078e23 % = 6.02214078 * 10^23 : the Avogadro number
``````

``````for i = 2:numel(t)
% ...
end
% Then plot
plot(t, lmda)
``````

``````% implementing equation 9 and figure 4
d = 3;     % dimension
N = 8000;  % number of molecules
D = 100;   % diffusion coefficient in micrometer^2/sec

dtPulses = 10; % Seconds between pulses
tPulses = 1:dtPulses:3600; % Time array to plot against
nt = numel(tPulses);
i = 1:nt;  % pulse numbers
u = 1 ./ (i.^(d/2)); % inner function in equation 9: individual pulse
for k = 2:nt % Running sum
u(k) = u(k-1)+u(k);
end
% Now plot for different H (mol/micrometer^3)
H = [0.01, 0.1, 1];
figure; hold on; linestyles = {':k', '--k', '-k'};
for nH = 1:3
lmda = ((1/(4*pi*D))*(N/H(nH)).*u).^(2/d);
plot(tPulses, lmda, linestyles{nH}, 'linewidth', 2)
end
``````

0 条评论