program ztbmie;
type complex=record r,i:real; end;

(************************************************************************)
(*                              Disclaimer                              *)
(* The authors, the publisher and/or their employees do not accept any  *)
(* responsibility or legal obligation for damage, harm, injury or loss  *)
(* caused by or due to the use of the program code presented.           *)
(************************************************************************)

const nmax=200; nmaxplusone=201;
   minreal=1.0E-38; maxreal=0.99E+38; precision=1.0E-06;

var i,n,maxn:integer;
   x,m,g,Qsca,theta:real;
   pie,tau,psi,psiprime,psiy,psiprimey:array [1..nmax] of real;
   zeta,zetaprime,zetay,zetaprimey:array [1..nmaxplusone] of complex;
   a,b:array [0..nmax] of complex;
   i1,i2:array [0..720] of real;

(************************************************************************)
(*                 ARITHMATIC FUNCTIONS AND PROCEDURES                  *)
(************************************************************************)

function pwrrs(x,y:real):real; (* x raised to the power y *)
   var scrint:integer;         (* defined for x>0 or (x<0 and y=integer) *)

   function is_int(y:real):boolean;
      begin is_int:=(y-trunc(y)=0.0); end;

   begin
   if ( (x<0)and(not is_int(y)) )
     then begin pwrrs:=0.0; writeln('error in function pwrrs'); end;
   if x=0 then begin pwrrs:=0.0; end;
   if ( (x<0) and (is_int(y)) )
      then begin
      if odd(round(y)) then pwrrs:=-pwrrs(-x,y)
                       else pwrrs:= pwrrs(-x,y);
      end;
   if x>0 then pwrrs:=exp(y*ln(x));
   end; (* pwrrs *)

function cabs(z:complex):real;               (* cabs:=|z| *)
   begin  cabs:=sqrt(sqr(z.r)+sqr(z.i)); end;

function creal(z:complex):real;               (* creal:=z.r *)
   begin  creal:=z.r; end;

procedure cadd(z1,z2:complex;var sum:complex); (* sum:=z1+z2 *)
   begin sum.r:=z1.r+z2.r; sum.i:=z1.i+z2.i; end;

procedure cmult(z1,z2:complex;var prod:complex); (* prod:=z1*z2 *)
   begin prod.r:=z1.r*z2.r-z1.i*z2.i;
         prod.i:=z1.r*z2.i+z1.i*z2.r; end;

procedure cdiv(z1,z2:complex;var ratio:complex); (* ratio:=z1/z2 *)
   begin
   if ((z2.r=0.0) and (z2.i=0.0))
      then begin
      ratio.r:=maxreal; ratio.i:=maxreal;
      write('error in cdiv');
      end
      else begin
      ratio.r:=(z1.r*z2.r+z1.i*z2.i)/(z2.r*z2.r+z2.i*z2.i);
      ratio.i:=(z1.i*z2.r-z1.r*z2.i)/(z2.r*z2.r+z2.i*z2.i);
      end;
   end; (* cdiv *)

(************************************************************************)
(*                        AUXILIARY PROCEDURES                          *)
(************************************************************************)

procedure calc_psi_zeta(x:real; maxn:integer);
   var n:integer;
   begin
   (* psi[n](x) = x*j[n](x)  HLST p123 , HBMF 10.1.11 *)
   psi[1]:=sin(x)/x-cos(x);
   psi[2]:=(3*pwrrs(x,-2)-1)*sin(x)-3*cos(x)/x;
   if maxn>2 then for n:=3 to maxn do
      psi[n]:=(2*n-1)*psi[n-1]/x-psi[n-2];              (* HBMF 10.1.19 *)

   (* psi'[n](x) *)
   psiprime[1]:=(1-pwrrs(x,-2))*sin(x)+cos(x)/x;        (* HBMF 10.1.21 *)
   if maxn>1 then for n:=2 to maxn do
         psiprime[n]:=psi[n-1]-n*psi[n]/x;

   (* zeta[n](x)=x*h(2)[n](x) complex HLST p123 *)
   zeta[1].r:=sin(x)/x-cos(x);                          (* HBMF 10.1.17 *)
   zeta[1].i:=sin(x)+cos(x)/x;
   zeta[2].r:=(-1+3*pwrrs(x,-2))*sin(x) - 3/x*cos(x);
   zeta[2].i:=3/x*sin(x) + (-1+3*pwrrs(x,-2))*cos(x);
   for n:=3 to maxn+1 do
      begin
      zeta[n].r:=(2*n-1)*zeta[n-1].r/x - zeta[n-2].r;
      zeta[n].i:=(2*n-1)*zeta[n-1].i/x - zeta[n-2].i;
      end;

   (* zeta'[n](x) complex *)
   for n:=1 to maxn do
      begin
      zetaprime[n].r:=(n+1)*zeta[n].r/x - zeta[n+1].r;
      zetaprime[n].i:=(n+1)*zeta[n].i/x - zeta[n+1].i;
      end;
   end; (* calc_psi_zeta *)

procedure mie_var(n:integer);
   (* calculates psi(x),psi'(x),zeta(x),zeta'(x) *)
   begin  (* and psi(y),psi'(y),zeta(y),zeta'(y) *)
   calc_psi_zeta(x*m,n);
   psiy[n]      :=psi[n];
   psiprimey[n] :=psiprime[n];
   zetay[n]     :=zeta[n];
   zetaprimey[n]:=zetaprime[n];
   calc_psi_zeta(x,n);
   end; (* mie_var *)

procedure calc_pie_tau(theta:real);
   var n:integer;
   begin
   (* pie[n](cos(theta))=dp[n]/dcos(theta) HLST p124 and HBMF p342 *)
   pie[1]:=-1;
   pie[2]:=-3*cos(theta);
   if maxn>2 then for n:=3 to maxn do
      pie[n]:=(((2*n-1)*cos(theta)*pie[n-1]-n*pie[n-2]))/(n-1);

   (* tau[n]=dP1[n]/d(theta)=dP1[n]/dcos(theta)*sin(theta) HLST p124 *)
   tau[1]:=-cos(theta);
   if maxn>1 then for n:=2 to maxn do
      tau[n]:=(n*cos(theta)*pie[n]-(n+1)*pie[n-1]);     (* HBMF 8.5.4 *)
   end; (* calc_pie_tau *)

(************************************************************************)
(*                      PROCEDURES FOR MAIN PROGRAM                     *)
(************************************************************************)

procedure nulling;
   begin
   for i:= 0 to 720 do begin i1[i]:=0.0; i2[i]:=0.0; end;
   for i:=1 to nmax do
      begin
      psiy[i]:=0.0;      psi[i]:=0.0;
      psiprimey[i]:=0.0; psiprime[i]:=0.0;
      end;
   for i:=0 to nmaxplusone do
      begin
      zetay[i].r:=0.0;       zeta[i].r:=0.0;
      zetay[i].i:=0.0;       zeta[i].i:=0.0;
      zetaprimey[i].r:=0.0;  zetaprime[i].r:=0.0;
      zetaprimey[i].i:=0.0;  zetaprime[i].i:=0.0;
      end;
   end;

procedure read_x_m;
   var d,lambda,n1,n2:real;
   begin
   n1:=0.0;
   write('x = ') ;readln(x);
   if x=0
      then begin
      write('lambda (nm) = '); readln(lambda);
      write('d      (nm) = '); readln(d);
      write('rfer.index surrounding = '); readln(n1);
      x:=2*pi*n1*d/(2*lambda);
      end;
   write('m = '); readln(m);
   if m=0
      then begin
      
 
      write('rfer.index    particle = '); readln(n2);
      if n1=0.0
         then begin
         write('rfer.index surrounding = '); readln(n1);
         end;
      m:=n2/n1;
      end;
   end; (* read_x_m *)

procedure mie_an_bn;    (* calculates a[n],b[n] HLST p123 *)
   var maxabsa,maxabsb:real; (* maxima of a[n] and b[n]   *)
     num,den:complex;        (* numerator and denominator *)

   procedure calc_an;
      begin
      num.r:=psiprimey[n]*psi[n]-m*psiy[n]*psiprime[n];
      num.i:=0;
      den.r:=psiprimey[n]*zeta[n].r-m*psiy[n]*zetaprime[n].r;
      den.i:=psiprimey[n]*zeta[n].i-m*psiy[n]*zetaprime[n].i;
      cdiv(num,den,a[n]);
      end;

   procedure calc_bn;
      begin
      num.r:=m*psiprimey[n]*psi[n]-psiy[n]*psiprime[n];
      num.i:=0;
      den.r:=m*psiprimey[n]*zeta[n].r-psiy[n]*zetaprime[n].r;
      den.i:=m*psiprimey[n]*zeta[n].i-psiy[n]*zetaprime[n].i;
      cdiv(num,den,b[n]);
      end;

   begin
   maxabsa:=minreal; maxabsb:=minreal;
   n:=0;
   repeat
      inc(n);
      mie_var(n);
      calc_an;
      if cabs(a[n])>maxabsa then maxabsa:=cabs(a[n]);
      calc_bn;
      if cabs(b[n])>maxabsb then maxabsb:=cabs(b[n])
      until (n>1)and(cabs(a[n])/maxabsa<precision) and
                    (cabs(b[n])/maxabsb<precision);
   maxn:=n;
   inc(n);
   mie_var(n);
   calc_an;  (* a[n+1] and b[n+1] needed to calculate g *)
   calc_bn;
   end; (* mie_an_bn *)


procedure mie_sigma(theta:real);
   (* calculates i1 en i2  HLST p35 by HLST p125 *)
   var s1,s2:complex;
   begin
   calc_pie_tau(theta);
   s1.r:=0.0  ;s1.i:=0.0;
   s2.r:=0.0  ;s2.i:=0.0;
   for n:=1 to maxn do
      begin
      s1.r:=s1.r + (2*n+1) * (a[n].r*pie[n]+b[n].r*tau[n]) / (n*(n+1));
      s1.i:=s1.i + (2*n+1) * (a[n].i*pie[n]+b[n].i*tau[n]) / (n*(n+1));
      s2.r:=s2.r + (2*n+1) * (b[n].r*pie[n]+a[n].r*tau[n]) / (n*(n+1));
      s2.i:=s2.i + (2*n+1) * (b[n].i*pie[n]+a[n].i*tau[n]) / (n*(n+1));
      end;
   i1[i]:=sqr(s1.r) + sqr(s1.i);
   i2[i]:=sqr(s2.r) + sqr(s2.i);
   end; (* mie_sigma *)

procedure calc_qsca; (* HLST p128 *)
   begin
   qsca:=0.0;
   for n:=maxn downto 1 do
      qsca:=qsca+2*(2*n+1)*(sqr(cabs(a[n]))+sqr(cabs(b[n])))/sqr(x);
   end;

procedure calc_g; (* HLST p128 *)
   var scr1,scr2,scr3:complex;
   begin
   g:=0.0;
   for n:=maxn downto 1 do
      
 
begin
      scr1.r:=a[n+1].r ;scr1.i:=-a[n+1].i ;cmult(a[n],scr1,scr1);
      scr2.r:=b[n+1].r ;scr2.i:=-b[n+1].i ;cmult(b[n],scr2,scr2);
      scr3.r:=b[n].r   ;scr3.i:=-b[n].i   ;cmult(a[n],scr3,scr3);
      cadd(scr1,scr2,scr1);
      g:=g+n*(n+2)*creal(scr1)/(n+1) + (2*n+1)*creal(scr3)/(n*(n+1));
      end;
   g:=4*g/(sqr(x)*qsca);
   end; (* calc_g *)


BEGIN (* main *)
nulling;
read_x_m;
mie_an_bn;
for i:=0 to 720 do mie_sigma(i*pi/720);
calc_qsca;
calc_g;
(* at this stage i1[0..720] and i2[0..720] contain the phase function *)
(* qsca is the relative scattering crosssection, g the anisotropy factor *)

END.  (* main *)