Calculate structure factor using fft function

22 Ansichten (letzte 30 Tage)
Hao Hoang
Hao Hoang am 29 Nov. 2018
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

Antworten (1)

Sayyed Ahmad Khadem
Sayyed Ahmad Khadem am 16 Mai 2019
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,

Kategorien

Mehr zu Statics and Dynamics finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by