; the procedure merlyn calculates amplitude scattering matrix 
; elements and efficiencies for extinction, total scattering
; and backscattering for a given size parameter and relative
; refractive index

pro merlyn,x,refrel,qext,qsca,qback,qabs,qratio
d=dcomplexarr(100000)
dx=x
y=x*refrel

; series terminated after nstop terms

xstop=x+4d0*x^(1d0/3d0)+2d0
nstop=long(xstop)
ymod=double(abs(y))
nmx=long(max([xstop,ymod]))+15

; logarithmic derivative d(j) calculated by downward recurrence
; beginning with initial value 0.0 + i*0.0 at j=nmx

d(nmx)=dcomplex(0d0,0d0)
nn=nmx-1
for n=1,nn do begin
  rn=nmx-n+1
  d(nmx-n)=(rn/y)-(1d0/(d(nmx-n+1)+rn/y))
  endfor

; riccati-bessel functions with real argument x calculated by
; upward recurrence

psi0=cos(dx)
psi1=sin(dx)
chi0=-sin(x)
chi1=cos(x)
apsi0=psi0
apsi1=psi1
xi0=dcomplex(apsi0,-chi0)
xi1=dcomplex(apsi1,-chi1)
qsca=0d0
qext=0d0
qabs=0d0
qback=0d0
back=dcomplex(0d0,0d0)
n=1
midpoint: dn=n
rn=n
fn=(2d0*rn+1d0)/(rn*(rn+1d0))
psi=(2d0*dn-1d0)*psi1/dx-psi0
apsi=psi
chi=(2d0*rn-1d0)*chi1/x-chi0
xi=dcomplex(apsi,-chi)
an=(d(n)/refrel+rn/x)*apsi-apsi1
an=an/((d(n)/refrel+rn/x)*xi-xi1)
bn=(refrel*d(n)+rn/x)*apsi-apsi1
bn=bn/((refrel*d(n)+rn/x)*xi-xi1)
qsca=qsca+(2d0*rn+1d0)*(double(abs(an))*double(abs(an))$
            +double(abs(bn))*double(abs(bn)))
qext=qext+(2d0*rn+1d0)*double(an+bn)
back=back+(2d0*rn+1d0)*((-1d0)^rn)*(an-bn)
psi0=psi1
psi1=psi
apsi1=psi1
chi0=chi1
chi1=chi
xi1=dcomplex(apsi1,-chi1)
n=n+1
rn=n
if n lt nstop then goto,midpoint
qsca=(2./(x*x))*qsca
qext=(2./(x*x))*qext
qback=(1./(x*x))*double(abs(back))*double(abs(back))
qabs=qext-qsca
qratio=qabs/qext
return
end
;