這是我在國外論壇找到的一個半導體模擬， 但一開始該程式作者提到aj和eta定義，到底function要怎打? 以下是程式: function xx=fermi(aj,eta) % Matlab m file to evaluate the Half-order Fermi-Dirac integral % For semiconductor/solid state applications % Half-order implies that j (= aj)= -1/2 or 1/2 or 3/2 or 5/2 % Created on 5-Jan-2007 by Natarajan and Mohankumar %======================================================== % Based on the Published work below % Title: The accurate numerical evaluation of half-order % Fermi-Dirac integrals % Authors: N.Mohankumar & A.Natarajan % Journal: Physica Status Solidi(b) vol.188, 1995, pp. 635-644 %============================================================ % xx is a dummy output % eta and aj are input values to be given by the user. % Printed output is the integral value. % Input variable aj is same as j in the definition of F(j,eta) given below % Accuracy: For solid state physics applications, eta values % typically lie in the interval [-5,25]; For eta in [-5,25], % the code below is guaranteed to yield 14 digit accuracy. %============================================================ % Integral definition is given below % % Integral is defined by F(j,eta) % Integrand is defined by (x^j)/[Gamma(1+j)* (exp(x-eta)+1)] % Lower limit=0; Upper limit=infinity % Integrated with respect to x % Please check whether your definition of the integral % needs the Gamma(1+j) factor or not. %=================================================================== % Trapezoidal Integration in y after the transformation % from the original integration variable x to y % where x= y^2 % Residue correction for the poles of the transformed integrand is % added to the trapezoidal integration sum to expedite convergence %========================================================= %Program begins format long e; %============================================================== % Evaluation of Trapezoidal sum begins range=8.; if eta > 0. range=sqrt(eta+64.);end; h=0.5; nmax=range/h; sum=0.; if aj== (-0.5) sum=1./(1.+exp(-eta));end; for i=1:nmax u=i*h; ff=2.*(u^(2.*aj+1))/(1.+exp(u*u-eta)); sum=sum+ff;end; %Trapezoidal Summation ends %============================================================== % Pole correction for trapezoidal sum begins pol=0.; npole=0; % Fix the starting value of BK1 to start while loop bk1=0; while bk1 <= 14.*pi npole=npole+1; bk=(2*npole -1)*pi; rho=sqrt(eta*eta+bk*bk); t1=1; t2=0; if eta < 0; tk=- aj*(atan(-bk/eta)+pi); elseif eta ==0; tk=0.5*pi*aj; else; eta > 0; tk=aj*atan(bk/eta); end; rk=- (rho^aj); tk=tk+0.5*atan(t2/t1); if eta < 0 rk= -rk; end; ak=(2.*pi/h)*sqrt(0.5*(rho+eta)); bk1=(2.*pi/h)*sqrt(0.5*(rho-eta)); if bk1 <= (14.*pi) gama=exp(bk1); t1=gama*sin(ak+tk)-sin(tk); t2=1.-2.*gama*cos(ak)+gama*gama; pol=pol+4.*pi*rk*t1/t2; end; %ends if loop above end; % Top while loop ends npole=npole-1; fdp=sum*h+pol; % Program ends with the following output %=============================================================================== disp('Fermi-Dirac Integral Value'); disp(fdp/gamma(1+aj)); disp('Number of trapezoidal points & number of poles'); disp([round(nmax),npole]); --

※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 104.38.74.64 ※ 文章網址: http://webptt.com/m.aspx?n=bbs/MATLAB/M.1410069913.A.26D.html