-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Dear Niu,
as far as I can tell, all your parameters are correct and the scattering term for f(s) you use is also correct. f(s) furthermore matches very closely those tabulated in the Intl. Tables C Tab. 6.1.1.1. My reproduction of the mentioned formula with r=2.0 using MAXIMA also shows quite a different graph. The graph does not make much sense: as d_max -> 0 1/d_max -> infinity and the integrand goes to infinity, because f(s) contains a constant positive term. Hence the integral should oscilatingly approach infinity and not stabilise as the upper integral limit approaches infinity. Best wishes, Tim On 09/13/2012 11:09 PM, Niu Tou wrote: > 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 > - -- - -- Dr Tim Gruene Institut fuer anorganische Chemie Tammannstr. 4 D-37077 Goettingen GPG Key ID = A46BEE1A -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.12 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/ iD8DBQFQUx2lUxlJ7aRr7hoRAi6IAJ9cU5etz4nb8Y7ti9zSjjL5P3ptJgCgtsl4 RM8Zn+l70MLPXoO1MpGKYUI= =Lq9Z -----END PGP SIGNATURE-----
