Program Shmie;

{ Written by

Eugeny B. Shybanov,
Department of Optics,
Marine Hydrophysical Institute,
Kapitanskaya 2, Sevastopol,
Crimea, Ukraine, 99000
e-mail: <lee@alpha.mhi.iuf.net>

This code is based on:

Comments:
This program computes complex amplitudes and effectiveness factors for Mie scattering
of light by spherical particles with size parameter 'x' up to 1 000 000 000,
or, more precisely, a number of series members should be less than 2**30 - 1.
Dlya chastits s pogloz'eniem tochnost'
rezul'tatov ne garantiruetsya, poskol'ku algoritm osnovan na pryamoy
rekursii dlya funktsiy Besselya-Rikatti.
( Drugaya moya programma dlya proizvol'nyh m1,m2 schitaet do 'x' =16 000 000,
ne ispol'zuya pri etom extended memory )
V to vremya kak, razrabotchikami sovremennyh kompilyatorov prilagayutsya
bol'shie usiliya k povysheniyu bystrodeystviya na desyatki protsentov,
nami uzhe razrabotany algoritmy raschetov napisannye na Assemblere,
kotorye dayut uvelichenie skorosti v tri raza po sraneniyu so
standartnymi kompilyatorami, a dlya komp'yuterov tipa Pentium-2
(naprimer INTEL(R) CELERON(TM)-MMX ) v 5 - 6 raz. Aktual'nost'
uvelicheniya bystrodeystviya raschetov dlya krupnyh chastits sostoit v
sleduyuz'em:
1) chislo ekstremumov indikatrisy rasseyaniya na monodispersnyh
chastitsah poryadka parametra difraktsii. Poetomu dlya polidispersnoy
vzvesi nuzhno integrirovat' po chislu tochek poryadka 'x' _srednee.
( Vozmozhno tak zhe predvaritel'no sglazhivat' indikatrisu, a dlya
etogo zhelatel'no raschityvat' uglovye funktsii dlya bol'shego chisla
tochek po uglu. Nizhe budet priveden timing, dlya vybora optimal'nogo
rezhima raschetov ).
2) Suz'estvuet metod ucheta nesferichnosti, kogda nesfericheskie chastitsy
odnogo razmera zamenyayutsya nekotorym raspredeleniem po razmeram
sfericheskih chastits. I,takim obrazom, model'noe raspredelenie stanovitsya
shire.

Neobhodimo tak zhe otmetit', chto vvidu nalichiya rezonansnyh chastot
elektricheskih i magnitnyh mod ( a_n, b_n kak funktsii 'x' ), luchshim
metodom integrirovaniya po raspredeleniyu budet metod pryamougol'nikov !!
Dlya shirokih raspredeleniy sleduet zadavat' neskol'ko intervalov s
razlichnym shagom po 'x', shag v kazhdom intervale i ih polozhenie
opredelyaetsya funktsiey raspredeleniya, a tak zhe tem, chto bolee vazhno
dlya zadachi: polnoe rasseyanie, obratnoe rasseyanie ili rasseyanie v
opredelennom intervale uglov.
Dannaya programma pervonachal'no byla napisana na vstroennom
Assemlere Turbo Pascal-ya 6.0, i poetomu chisto paskalevskaya versiya
mozhet byt' dlya cheloveka ne ochen' ponyatna. Tak naprimer, v privedennom
kode otsutstvuet opisanie tipa komleksnyh chisel.
Hotya dlya bol'shih pogloz'ayuz'ih chastits programma nepravil'no rabotaet,
v programme realizovan algoritm dlya pogloz'ayuz'ih chastits. Esli kogo-
nibud' interesuet algoritm dlya pogloz'ayuz'ih chastits, to on u menya
est' kak v paskalevskom, tak i v assemblernom variante.
Privozhu timing dlya protsessora INTEL(R) CELERON(TM)-MMX - 333

Vremya rascheta integral'nyh harakteristik i kompleksnyh amplitud
na 1 chastitse s parametrom difraktsii 'x' v N_fi - tochkah po uglu [0,180 deg.]


build-in ASM (0.046 + 0.0016*N_fi) * N_x/10000 sec.
Turbo PAS 6.0 (0.080 + 0.0092*N_fi) * N_x/10000 sec.

gde N_x - chislo chlenov ryada ( N_x poryadka x + 4*x**1/3)

}

Const
nfimax = 2001; { chislo tochek po uglu v intervale 0 - 180 gradusov
konstanta v Turbo 6.0 mozhet prinimat' znacheniya
ot 2 do 4074 }

nHfimax = (nfimax+1) shr 1; { chislo tochek po uglu v intervale 0 - 90
gradusov. Vvidu simmetrii Pi i Tau
funktsii dostatochno poloviny intervala}

Type

{ Tipy dlya uglovyh funktsiy, poskol'ku problem s uglovymi funktsiyami
kak pravilo net zdes' ispol'zuetsya 64-bitovoe predstavlenie
chisel s plavayuz'ey tochkoy }

AngleType = array[0..nfimax-1] of double;
{Tip dlya indikatrisy }

SAngleType = array[0..nHfimax-1] of double;
{ dlya vspomogatel'nyh uglovyh funktsiy, zadavaemyh pri cos(teta) > =0 }

PTType = array[0..(nHfimax shl 1) - 1] of double;
{ dlya fuktsiy Pi i Tau zapisannyh v vide massiva massiva
(P[0],T[0],P[1],T[1],...P[i],T[i],...,P[n],T[n]) ,
gde indeks i sootvetstvuet nomeru ugla }

FourType = array[1..4] of double;
{ vspomogatel'nyy tip dlya optimizatsii summirovaniya ryadov
sootvetstvuet massivu Re(i1),Im(i1),Re(i2),Im(i2) }

EOType = array[0..nHfimax-1] of FourType;
{ massiv tipa FourType dlya chetnyh i nechetnyh uglovyh funktsiy }


VAR
s1r,s1i,
s2r,s2i : ^AngleType; { ssylki na deystvitel'nyy i mnimye chasti
amplitud i1 , i2 kak funktsii ugla }

i,nfi,j,i_x : integer;
n : longint;
i1,i2,i3,i4,
i0,inten,
DegPol,PhFun,
qg,qb,t,
qs,qe,y,la,
x,m1,m2 : double;
ft : text; { text file }

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

{ Protsedura raschityvaet sleduyuz'ie znacheniya uglovyh funktsiy
Pi,Tau dlya uglov, kosinusy kotoryh ravny mu[i].
nhfi - chislo uglov v intervale [0..90]
n - poryadok funktsiy (ili indeks ryada),
mu - massiv kosinusov
p0 - na vhode znacheniya Pi-funktsii ((n-2)-go poryadka)
na vyhode (n-1)-go
pt - na vhode predyduz'ie znacheniya Pi i Tau -funktsiy ((n-1)-go poryadka)
na vyhode n-go poryadka.
}

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);
{
Raschityvaet sleduyuz'ee znachenie funktsii Rikkati-Besselya.
V glavnoy protsedure ispol'zuetsya kak dlya Psi = x*Jn,
tak i dlya Xi = - x*Yn ( Kci = Psi -i*Xi )
twonm1 = 2*(n-1)
xinv = 1/x
brf1 - na vhode Bess_Fun(n-2)
na vyhode Bess_Fun(n-1)
brf2 - na vhode Bess_Fun(n-1)
na vyhode Bess_Fun(n)
Poskol'ku pri n>x funktsiya Psi nachinaet rashodit'sya, to
ispol'zuetsya tip extended (80 bit)
}
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);
{
Vychislyaet sleduyuz'ee znachenie funktsii logarifmicheskoy proizvodnoy
D_n(m*x) = d_ln(Psi(m*x))/d_(m*x).
Esli Im(m)= 0, to nachinaya s opredelennogo znacheniya
indeksa n pryamaya rekursiya rashoditsya. Znachenie indeksa pri kotorom
nachinaetsya rashodimost' opredelyaetsya kak Im(m), x , tak i Re(m).
Parametry protsedury:
rxr - Re(1/(m*x)
rxi - Im(1/(m*x)
dnr - Re(D_n)
dni - Im(D_n)
}
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);
{
Eta protsedura - osnovnaya chast' moego varianta vychisleniya
kompleksnyh amplitud (i1,i2). Razbienie na chetnuyu i
nechetnuyu chast' uglovyh funktsiy daet dvuhkratnyy vyigrysh pri
simmetrichnoy setke po uglu ( Teta(nfi-i) = 180 - Teta(i) )
po sravneniyu s drugimi Mi kodami. Zametim, chto
osnovnoe vremya raboty glavnoy Mi-protsedury BigP2 zatrachivaetsya
na rabotu etoy protsedury i protsedury NEXTPT8.

nh - chislo tochek ot 0 do 90°,
ar - Re(a_n)
ai - Im(a_n)
br - Re(b_n)
bi - Im(b_n)
b - priznak chetnosti
Pt - Pi,Tau funktsii n-go poryadka P[0],T[0],..P[i],T[i] dlya raznyh uglov
Ev - Chetnaya chast' ((i1[0],i2[0]),.......(i1[i],i2[i])
Od - Nechetnaya chast' ..................................
i1,i2 - kompleksny, i - sootvetsvuet uglu.

V ASM-variante 'b' menyalos' putem invertirovaniya bita, a dlya
Paskalya luchshe perepisat' protseduru tak:

Procedure Name(......; var b: Boolean ; ....);
VAR
...
Begin
....
if b then
..........
else
...........;
b:= not(b)
End; }

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);

{ Vychislyaet koeffitsienty a_n ili b_n , v zavisimosti ot br,bi

dlya a_n br= Re( D_n(m*x)/m + n/x), bi= Im( D_n(m*x)/m + n/x);
dlya b_n br= Re( m*D_n(m*x) + n/x), bi= Im( m*D_n(m*x) + n/x)
psi0 - psi(n-1)
psi1 - psi(n)
x0 - xi(n-1)
x1 - xi(n)
cr - Re( a_n ili b_n)
ci - Im( a_n ili b_n)
}
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;

{-------------------------------}
Procedure BigP2(m1,m2,x:extended;nfi:integer; var qe,qs,qg,qb:double;
var b1,b2,b3,b4 :AngleType); { x > 0.01, m2<<m1, m1>1 }
{ Raschet kompleksnyh amplitud i faktorov effektivnosti dlya odnorodnogo
shara s parametrom difraktsii "x".
m1 - Re(m)
m2 - Im(m), pogloz'eniyu sootvetstvuet polozhitel'noe znachenie m2.
nfi - chislo tochek ravnomernoy po uglu setki v intervale 0-180
b1 - Re(i1), b2 - Im(i1), b3 - Re(i2), b4 - Im(i2), - raschityvaemye
znacheniya amplitud v tochkah uglovoy setki.
.......................................

QG = 1/x/x * MOD( SUM( (2n+1) (-1)**n (AN-BN) ) )**2 }

Const
Eps = 1e-10; { pogreshnost' vychisleniya secheniya rasseyaniya. Mozhno umen'shit'.}
Var
Pe,Po : ^EOType; { chetnaya i nechetnaya chast' }
PiTau : PTType; { Pi,Tau - (Pi(0),Tau(0),..Pi(i),Tau(i) }
mua,PiPrev : SAngleType; { kosinusy i predyduz'ie znacheniya
funktsii Pi v intervale 0 - 90 }
hnr,hni,
nix,cf,
ani,anr,
bni,bnr,
ix,y,mi,mr,
imi,imr,
imxi,imxr,
xi1,xi0,
psi1,psi0,
d1i,d1r : extended;
{ 80 - bitovoe opisanie peremennyh uvelichivaet stabil'nost' algoritma i
diapazon ustoychivosti pryamyh rekursiy }

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);
{ writeln(memavail); }

{ Nachal'nye znacheniya uglovyh funktsiy }
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));
{ maksimal'noe chislo chlenov ryada. Vmesto 4 v raznyh programmah
vstrechayut'sya znacheniya 5 i 6. S uvelicheniem etogo koeffitsienta
tochnost' rascheta rasseyanie na ugly > 90 uvelichivaetsya. Odnako v svyazi
s rashozhdeniem algoritma pryamoy rekursii uvelichivat' etot koeffitsient
ne sleduet }


maxerror := EPS/nbrmax;
{ Etot operator ya nedavno izmenil !!
Pri bol'shih "x" , nbrmax primerno ravno "x", sledovatel'no
1/nbrmax - otsenivaet srednyuyu dolyu vklada n-go chlena v rasseyanie }

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;

{Vychislenie nachal'nyh znacheniy logarifmicheskoy proizvodnoy,
peremennye xi1,psi0,psi1 - zdes' ne imeyut otnosheniya k funktsiyam Besselya }
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;

{ Nachal'nye znacheniya funktsiy Rikkati-Besselya}
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);
{ V printsipe oprator repeat .... until mozhno zamenit' na tsikl
for L:= 1 to nbrmax ......,
poskol'ku:
- pri L>x ryady bystro shodyatsya;
- dlya bol'shih "x" n_max poryadka "x",
- a ot rashodimosti REPEAT UNTIL ne lechit. }

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);
m2:= 0;
repeat
write(' nfi[2..',nfimax,'] =');
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,'out1.txt');
{$I-}
append(ft); { Dobavit' zapic' v file }
{$I+}
if IOResult<>0 then
rewrite(ft);
writeln(ft);
writeln(ft,' x=',x:17,' Re(m)=',m1:17,' Im(m)=',m2:17);
writeln(ft,' Ext=',qe:17,' Sca=',qs:17);
writeln(ft,' gQsc=',qg:17,' Sback=',qb:17);
writeln(ft,' Angles Re(I1) Im(I1)',' Re(I2) Im(I2)');
for j:= 0 to nfi-1 do
begin
t:= j*180.0/(nfi-1);
writeln(ft,t:7:2,' ',s1r^[j]:16,' ',s1i^[j]:16,
' ',s2r^[j]:16,' ',s2i^[j]:16)
end;
close(ft);
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,' Ext=',qe:17,' Sca=',qs:17);
writeln(ft,' gQsc=',qg:17,' Sback=',qb:17);
writeln(ft,' Angles PhaseFun Intens Deg_Pol',
' Stoks1 Stoks2 Stoks3 Stoks4');

y:= 4/sqr(x)/qs;
if nfi > 1 then
for j:= 0 to nfi-1 do
begin
t:= j*180.0/(nfi-1);
i1:= (sqr(s1r^[j])+sqr(s1i^[j]));
i2:= (sqr(s2r^[j])+sqr(s2i^[j]));
i3:= (s1r^[j]*s2r^[j] + s1i^[j]*s2i^[j]);
i4:= (s1r^[j]*s2i^[j] - s2r^[j]*s1i^[j]);
i0:= i1+i2;
inten:= 0.5*i0;
DegPol:= (i2-i1)/i0;
PhFun:= Inten*y/(2*Pi);
writeln(ft,t:7:2,' ',Phfun:16,' ',Inten:16,
' ',degpol:16,
' ',i1*y:16,' ',i2*y:16,
' ',i3*y:16,' ',i4*y:16)
end;
close(ft);

End.