Dear Colleagues,
I am trying to repeat a series termination effect calculation displayed as
figure 2 in a publihsed paper
(http://www.ncbi.nlm.nih.gov/pubmed/12215645). Formula
(1) was used to implement this calculation. Since f(s) is not defined in
detail in this paper, I used formula and parameters listed in another
paper (http://scripts.iucr.org/cgi-bin/paper?a05896) to calculate it.
However, the result I got is not consistent with figure 2 of the first
paper. I am not sure if the formulas I used are right or not. Or if there
is any problem in the MatLab code, which I list below:
###########
clear all;clc;format compact;format long;
% matrix of a, b, c coefficients:
% rows: Fe, S, Fe1, Mo
% columns: A1; B1; A2; B2; A3; B3; A4; B4; C
fM = ...
[11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339 1.40818;...
7.18742 1.43280 5.88671 0.02865 5.15858 22.1101 1.64403 55.4651
-3.87732;...
11.9185 4.87394 7.04848 0.34023 3.34326 15.9330 2.27228 79.0339 1.40818;...
19.3885 0.97877 11.8308 10.0885 3.75919 31.9738 1.46772 117.932 5.55047];
%%% store radius data:
% distance from: origin
% columns: Fe, S, Fe, Mo
R_el = [2.0 3.3 3.5 3.5];
RHO_t = zeros(4,400);
for numel = 1:4
EL = numel;
RHO = zeros(1,400);
dmax = zeros(1,400);
for iter = 1:400
dmax(iter) = iter/100; % in angstroms
% numerical integration
int_fun = @(s) 4*pi*(s.^2).* ...
(fM(EL,1).*exp(-fM(EL,2).*(s.^2)*0.25) + ...
fM(EL,3).*exp(-fM(EL,4).*(s.^2)*0.25) + ...
fM(EL,5).*exp(-fM(EL,6).*(s.^2)*0.25) + ...
fM(EL,7).*exp(-fM(EL,8).*(s.^2)*0.25) + fM(EL,9)).* ...
sin(2*pi*s*R_el(EL))./(2*pi*s*R_el(EL));
RHO(iter) = quad(int_fun,0,1/dmax(iter));
clc;display(iter);display(numel);
end
RHO_t(numel,:) = RHO;
end
RHO_t(1,:)= 6*RHO_t(1,:);
RHO_t(2,:)= 9*RHO_t(2,:);
figure;
axis([0.5 3.5 -10 10]); hold on;
plot(dmax,RHO_t(1,:),...
dmax,RHO_t(2,:),...
dmax,RHO_t(3,:),...
dmax,RHO_t(4,:),...
dmax,sum(RHO_t,1));
title('Electron Density Profile');
legend('Fe','S','Fe1','Mo','Sum');
xlabel('d_m_a_x'); ylabel('Rho(r)');
set(gca,'XDir','reverse');
##############
Any suggestions will be appreciated. Thanks!
Niu
