Calculate structure factor using fft function

22 次查看（过去 30 天）
Hao Hoang
Hao Hoang 2018-11-29
回答： Sayyed Ahmad Khadem 2019-5-16
Hi guys, I'm coding a program on matlab to find and graph the structure factor of materials by doing a fourier transform of the radial distribution functions.
My code for the radial distribution is as follow:
function res = radialDistribution2D(coords,Lx,Ly,dr)
%calculate rad Distribution of 1 frame with coordinates of particles in set
%COORDS and frame boundaries Lx Ly, R STEP dr
rangeC = min(Lx,Ly)/4;
centerParts = cFilter(coords,Lx/2,Ly/2,rangeC);
nCenter = length(centerParts);
raw_pair = zeros(nCenter,ceil(rangeC/dr)+1);
for i = 1:nCenter
temp_ref = cFilter(coords,centerParts(i,1),centerParts(i,2),rangeC);
vec_dist = temp_ref - centerParts(i,:);
dist = sqrt(sum((vec_dist).^2,2));
temp = dist((dist>dr) & (dist<rangeC));
temp = round(temp/dr);
for j=1:length(temp)
val = temp(j)+1;
raw_pair(i,val) = raw_pair(i,val)+1;
end
normCoeff = rangeC^2/2/length(temp_ref)/nCenter/dr^2;
for k = 1:size(raw_pair,2)
raw_pair(i,k) = raw_pair(i,k)*normCoeff/k;
end
end
res = sum(raw_pair);
end
where the cFilter function is:
function centerPart = cFilter(coords,cLx,cLy,range)
%get circle of particles radius RANGE with CENTER(cLx,cLy)in picture frame Lx Ly
nPart = length(coords);
c_vect = [cLx cLy];
temp = zeros(nPart,2);
j = 1;
for i = 1:nPart
rad_vect = coords(i,:) - c_vect;
if (sum(rad_vect.*rad_vect)<range^2)
temp(j,:) = coords(i,:);
j = j+1;
end
end
sizeCenter = j - 1;
centerPart = zeros(sizeCenter,2);
for i = 1:sizeCenter
centerPart (i,:) = temp(i,:);
end
end
The COORDINATEs of particles is given by a data set that I acquires from the images of my experiments (with colloids).
After this I get the matrix avr_res containing the radial distribution values, then, according to this book Introduction to Liquids State Physics, I should apply the FFT on the function (avr_res -1) then normalize the received matrix as the code below:
L = length(avr_res); %length of radial distribution signal
Fs = 1/dr; %sampling frequency
nfft= 2^nextpow2(L);
temp = fft((avr_res-1),nfft);
temp = - imag(temp);
test = temp;
for i = 1:length(temp)
temp(i) = 1 + 2*3.14*rho*temp(i)/(i*Fs);
end
temp = temp(1:nfft/2+1);
f = Fs*(0:(nfft/2))/nfft;
plot(f,temp);
This should works fine, theoretically, but still, I get very strange results compared to the known structure factor graphs of typical liquids, here is my result for a set of liquid state particles, but normally, the very first values of the structure factors should be around 0.
Could anyone points out where was I wrong in my code or method ?
P/s: As I have read from certain sources then the built-in matlab fft give poor results at low spatial frequency? Is this correct?
Thanks for any help.
S graph.png

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

回答（1 个）

Sayyed Ahmad Khadem
Hi Hao,
I have similar problem, but I am seeing that no has answered yet. Did you realize where the problem is? I was wondering if you could share your experince with us.
Thanks,

类别

在 Help Center 和 File Exchange 中查找有关 Statics and Dynamics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by