project2.dpr

Mie散射数值计算的Delphi代码,用来计算粒子的光散射

program Project1;

{$APPTYPE CONSOLE}

///uses
///SysUtils;

Const
nfimax = 2001;
nHfimax = (nfimax+1) shr 1;

Type
AngleType = array[0..nfimax-1] of double;
SAngleType = array[0..nHfimax-1] of double;
PTType = array[0..(nHfimax shl 1) - 1] of double;
FourType = array[1..4] of double;
EOType = array[0..nHfimax-1] of FourType;

VAR
s1r,s1i,
s2r,s2i : ^AngleType;
nfi,j: integer;
//n : longint;
s11,s12,s33,s34,
qg,qb,t,
qs,qe,
x,m1,m2 : double;
ft : text;

{--------------------------------------------------------------}
Procedure NEXTPT8(nhfi:integer; n :longint;
var mu,p0 : SAngleType; Var pt:PTType);

Var
n2,nm1,np1 : longint;
i,j : integer;
y : double;
Begin
j:= 0;
n2:= (n shl 1) - 1;
nm1:= n - 1;
np1:= n + 1;
for i:= 0 to nhfi-1 do
begin
y:= (n2 * mu[i] * PT[j] - n * p0[i]) /nm1;
p0[i]:= PT[j];
PT[j]:= y;
inc(j);
PT[j]:= n * mu[i] * y - np1 * p0[i];
inc(j)
end
End;


{------------------------------------------------}
Procedure NextBR(twonm1:longint; xinv:extended;
var brf1,brf2:extended);
Var
y : extended;
Begin
y:= brf2;
brf2:= twonm1 * xinv * y - brf1;
brf1:= y
End;


{-----------------------------------}
Procedure NextCDn(n:longint; rxr,rxi:extended;
var dnr,dni : extended);
Var
y1,y2,z : extended;
Begin
y1:= n * rxr;
y2:= n * rxi;
dnr:= y1 - dnr;
dni:= y2 - dni;
z:= 1/( sqr(dni) + sqr(dnr));
dnr:= dnr * z - y1;
dni:= - dni * z - y2
End;


{--------------------------------------}
Procedure Svert4( nh:integer; ar,ai,br,bi : extended;
var b:integer; var Pt:PTType; var Ev,Od : EOType);
Var
i,j : integer;
Begin
j:= 0;
if b > 0
then for i:= 0 to nh-1 do
begin
Ev[i][1]:= Ev[i][1] + ar * PT[j];
Ev[i][2]:= Ev[i][2] + ai * PT[j];
Ev[i][3]:= Ev[i][3] + br * PT[j];
Ev[i][4]:= Ev[i][4] + bi * PT[j];
inc(j);
Od[i][1]:= Od[i][1] + br * PT[j];
Od[i][2]:= Od[i][2] + bi * PT[j];
Od[i][3]:= Od[i][3] + ar * PT[j];
Od[i][4]:= Od[i][4] + ai * PT[j];
inc(j)
end
else for i:= 0 to nh-1 do
begin
Od[i][1]:= Od[i][1] + ar * PT[j];
Od[i][2]:= Od[i][2] + ai * PT[j];
Od[i][3]:= Od[i][3] + br * PT[j];
Od[i][4]:= Od[i][4] + bi * PT[j];
inc(j);
Ev[i][1]:= Ev[i][1] + br * PT[j];
Ev[i][2]:= Ev[i][2] + bi * PT[j];
Ev[i][3]:= Ev[i][3] + ar * PT[j];
Ev[i][4]:= Ev[i][4] + ai * PT[j];
inc(j)
end;
b:= -b;
End;


{--------------------------------}
Procedure GetMod(br,bi,ps0,ps1,x0,x1:extended; var cr,ci:extended);
Var
rch,ich,
rzn,izn,z : extended;
Begin
rch:= br * ps1 - ps0;
ich:= bi * ps1;
rzn:= rch + bi * x1;
izn:= ich - br * x1 + x0;
z:= sqr(rzn) + sqr(izn);
cr:= ( rch * rzn + ich * izn) /z;
ci:= ( rzn * ich - rch * izn) /z
End;


{-------------------------------}
//m1相对折射率实部; m2相对折射率虚部(虚部应该取正号，我想才正常。);x无因次粒径；nfi,0~pi节点数;
//qe,qs,qg,qb分别为Qe,Qs,g,Qb
Procedure BigP2(m1,m2,x:extended;nfi:integer; var qe,qs,qg,qb:double;
var b1,b2,b3,b4 :AngleType);
Const
Eps = 1e-10;
Var
Pe,Po : ^EOType;
PiTau : PTType;
mua,PiPrev : SAngleType;
hnr,hni,
nix,cf,
ani,anr,
bni,bnr,
ix,y,mi,mr,
imi,imr,
imxi,imxr,
xi1,xi0,
psi1,psi0,
d1i,d1r : extended;
i,j,k,nh : integer;
nbrmax,
tnp1,np1,L,nm1 : longint;
sbi,sbr,sgqs,
anim1,anrm1,
bnim1,bnrm1,
SQE,SQS,
maxerror,dop : double;
b : integer;
Begin

sbi:= 0;
sbr:= 0;
sgqs:= 0;
anim1:= 0;
anrm1:= 0;
bnim1:= 0;
bnrm1:= 0;
SQE:= 0;
SQS:= 0;
if nfi = 1
then begin
nfi:= 2;
nh:= 0
end
else nh:= (nfi+1) shr 1;
new(pe);
new(po);

y:= Pi/(nfi-1);
for i:= 0 to nh-1 do
mua[i]:= cos(i*y);
j:= 0;
for i:= 0 to nh-1 do
begin
pitau[j]:= 1;
inc(j);
pitau[j]:= mua[i];
inc(j);
for k:= 1 to 4 do
begin
Pe^[i][k]:= 0;
Po^[i][k]:= 0
end;
PiPrev[i]:= 0
end;

nbrmax:= 2+round(x + 4*exp(ln(x)/3));

maxerror := EPS/nbrmax;
ix:= 1/x;
mr:= m1;
mi:= m2;

y:= 1/( sqr(mr) + sqr(mi) );
imr:= mr * y;
imi:= - mi * y;
imxr:= imr*ix; { Re(1/(m*x)) }
imxi:= imi*ix; { Im }
y:= x * mr;
y:= y + y;
xi0:= sin(y);
psi0:= cos(y);
y:= x * mi;

if y < 5677 then
begin
y:= exp(y + y);
xi1:= 0.5*(1/y - y);
psi0:= 0.5*(1/y + y) - psi0;
psi1:= 0;
if psi0 > 1e-10
then begin
d1r:= xi0/Psi0;
d1i:= xi1/psi0
end
else begin
y:= 2*x*mi;
if y < 1e-2400 then
y:= 1e-2400;
psi0:= 2/(sqr(xi0)+sqr(y));
d1r:= xi0*psi0;
d1i:= -y*psi0
end
end
else begin
d1r:= 0;
d1i:= -1
end;

psi0:= cos(x);
psi1:= sin(x);
xi0:= -psi1;
xi1:= psi0;

tnp1:= 1; { 2*n+1}
np1:= 1; { n+1 }
L:= 0;
b:= 1; { b = true }
repeat
nm1:= L;
L:= np1;
inc(np1);
NextBR(tnp1,ix,psi0,psi1);
NextBR(tnp1,ix,xi0,xi1);
NextCDN(l,imxr,imxi,d1r,d1i);
nix:= L * ix;
hnr:= nix + d1r*imr - d1i*imi;
hni:= d1r*imi + d1i*imr;
GetMod(hnr,hni,psi0,psi1,xi0,xi1,anr,ani);
hnr:= nix + d1r*mr - d1i*mi;
hni:= d1r*mi + d1i*mr;
GetMod(hnr,hni,psi0,psi1,xi0,xi1,bnr,bni);
if nm1<>0 then
NEXTPT8(nh,L,mua,piprev,PiTau);
tnp1:= tnp1 + 2;
cf:= tnp1/L/np1;

dop:= (anr+bnr)*tnp1;
sqe:= sqe + dop;
sqs:= sqs + (sqr(anr)+sqr(ani)+sqr(bnr)+sqr(bni) )*tnp1;
sbr:= (anr - bnr)*tnp1 - sbr;
sbi:= (ani - bni)*tnp1 - sbi;
sgqs:= sgqs +
(L - 1/L) * ( anr * anrm1 + ani * anim1 +
bnr * bnrm1 + bni * bnim1 ) +
cf* ( anr * bnr + ani * bni );
dop:= dop/sqs;
anim1:= ani;
anrm1:= anr;
bnim1:= bni;
bnrm1:= bnr;
anr:= anr * cf;
bnr:= bnr * cf;
ani:= ani * cf;
bni:= bni * cf;
Svert4(nh,anr,ani,bnr,bni,b,PiTau,Pe^,Po^);
until (L=nbrmax) or (dop<maxerror);

ix:= sqr(ix);
QE:= SQE*2*ix;
QS:= SQS*2*ix;
QG:= SGQS*4*ix/QS;
QB:= ix*(sqr(sbr)+sqr(sbi));

j:= nfi;
for i:= 0 to nh-1 do
begin
dec(j);
b1[i]:= Pe^[i][1] + Po^[i][1];
b1[j]:= Pe^[i][1] - Po^[i][1];
b2[i]:= Pe^[i][2] + Po^[i][2];
b2[j]:= Pe^[i][2] - Po^[i][2];
b3[i]:= Pe^[i][3] + Po^[i][3];
b3[j]:= Pe^[i][3] - Po^[i][3];
b4[i]:= Pe^[i][4] + Po^[i][4];
b4[j]:= Pe^[i][4] - Po^[i][4]
end;
dispose(po);
dispose(pe);

writeln(' N_max=',nbrmax:8,' N= ',L:8);

End;

{---------main program --------}

Begin
write(' m1 =');
readln(m1);

write(' m2 =');
readln(m2);

repeat
write(' nfi[2..',nfimax,'] ='); //取值2到nfimax=20001
readln(nfi);
until (nfi>1) and (nfi<=nfimax);
write(' x = ');
readln(x);

new(s1r);
new(s1i);
new(s2r);
new(s2i);
writeln('x=',x:14:3);
BigP2(m1,m2,x,nfi,qe,qs,qg,qb,s1r^,s1i^,s2r^,s2i^);

assign(ft,'out2.txt');
{$I-}
//append(ft);
{$I+}
if IOResult<>0 then
rewrite(ft);
writeln(ft);

writeln(ft,' x=',x:17,' Re(m)=',m1:17,' Im(m)=',m2:17);
writeln(ft,' qext=',qe:17,' qsca=',qs:17);
writeln(ft,' g=',qg:17);
writeln(ft,' Angles   s11                s12              s33              s34');

if nfi > 1 then
for j:= 0 to nfi-1 do
begin
t:= j*180.0/(nfi-1);

s11:=(s2r^[j]*s2r^[j]+s2i^[j]*s2i^[j]+s1r^[j]*s1r^[j]+s1i^[j]*s1i^[j])/2;
s12:=(s2r^[j]*s2r^[j]+s2i^[j]*s2i^[j]-s1r^[j]*s1r^[j]-s1i^[j]*s1i^[j])/2;
s33:=(s1r^[j]*s2r^[j]+s1i^[j]*s2i^[j]);
s34:=(s1r^[j]*s2i^[j]-s2r^[j]*s1i^[j]);

writeln(ft,t:7:2,' ',
' ',s11:16,' ',s12:16,
' ',s33:16,' ',s34:16)
end;
close(ft);
End.