%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pypdf.m %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main program for the calculation of Percus Yevick pair distribution
% function and structure factor for a medium with spherical particles
% -- Part of the Electromagnetic Wave MATLAB Library (EWML)--
% <http://www.emwave.com/>
% Original: by K.H. Ding, November 1998
% input parameters
fv=input('Enter the volume fraction of particles : ');
dia=input('Enter the diameter of particles : ');
tau= 999999999.
tauc=(2-sqrt(2))/6;
xic=(3*sqrt(2)-4)/2;
rm=5;
nk=512;
dk=pi/(rm*dia);
dr=rm*dia/(nk+1);
hk=zeros(2*nk+2,1);
vol=pi*dia^3/6;
rho=fv/vol;
xi=fv;
xi1=1-xi;
nu=tau+xi/xi1;
gama=xi*(1+xi/2)/(3*xi1^2);
lamd=6*(nu-sqrt(nu^2-gama))/xi;
mu=lamd*xi*xi1;
if mu > 1+2*xi
break;
end
cst1=xi/xi1;
cst2=1-lamd*xi+3*cst1;
cst3=3-lamd*xi1;
fpk=fopen('pysf.dat','w+');
pk=0;
pp1=cst1*(4-lamd+3*cst1)+1;
pp2=0;
pyhk=(1/(pp1^2+pp2^2)-1)/rho;
pysf=1+pyhk*rho;
fprintf(fpk,'%6u %14.9f \n',pk,pysf);
for ik=1:nk
pk=ik*dk;
x=pk*dia/2;
snx=sin(x);
csx=cos(x);
psix=snx/x;
phix=3*(snx-x*csx)/(x^3);
pp1=cst1*(cst2*phix+cst3*psix)+csx;
pp2=cst1*x*phix+snx;
pyhk=(1/(pp1^2+pp2^2)-1)/rho;
pysf=1+pyhk*rho;
fprintf(fpk,'%6u %14.9f \n',pk,pysf);
hk(ik+1)=pk*pyhk;
end
fclose(fpk);
hw=-2*fft(hk);
hr=imag(hw(2:nk+1));
fpr=fopen('pypdf.dat','w+');
for ir=1:nk
r=ir*dr/dia;
if r >= 1
g=1+hr(ir)/(rm*4*pi*r*dia^2);
fprintf(fpr,'%14.9f %14.9f \n',r,g);
end
end
fclose(fpr);
% plot pair distribution function
clear;
load pdfpy.dat -ascii;
load pysf.dat -ascii;
figure;
plot(pypdf(:,1),pypdf(:,2),'r:');
axis([0.0,5.0,0.0,3.0]);
title('PY pair distribution function');
xlabel('r/(2a) ');
ylabel('g(r) ');
legend('fv=0.3');
figure;
plot(pysf(:,1),pysf(:,2),'r:');
axis([0.0,150.0,0.0,2.0]);
title('PY structure factor');
xlabel('k ');
ylabel('S(k) ');
legend('fv=0.3');
clear;
pypdf.rar_theory
版权申诉
94 浏览量
2022-09-24
02:11:47
上传
评论
收藏 1KB RAR 举报
朱moyimi
- 粉丝: 61
- 资源: 1万+