subroutine sphrd(lambda,xrot,xperp,rev,eps0,zeps)
C--------/---------/---------/---------/---------/---------/---------/--
C	>>>  lambda,xc,xb,eps0,zeps
C	<<<  sext(3)
C
C  xrot ... the half-length of the spheroid along the rotational z-axis
C  xperp ... the half-length of the spheroid along the perpendicular axis
C  rev ... equal-volume-sphere radius 	 
C--------/---------/---------/---------/---------/---------/---------/--
	implicit none
	integer NOUT
C ::: number of the output unit for cross sections and scattering matrix
      PARAMETER (NOUT=35)

	integer ij,ns,npol
	real*8 lambda,xrot,xperp,rev,eps0,pi,sext(5),xk,xmn,xmj,xx,xe,
     &      xlz,xlx,xdz,xdx,xvol,xlp,xdp,xcr,xfx,xapl,pf,qf
	complex*16 ci,cone,zalph(5),zeps

      DATA PI/3.141592653589793d0/
	data ci/(0.d0,1.d0)/,cone/(1.d0,0.d0)/

	xk=2.d0*pi*sqrt(eps0)/lambda
	xcr=pi*rev**2        !an effective geom. cross section
	xvol=xrot*xperp**2/3.d0   !V/(4.d0*pi)

	if (xrot.gt.xperp) then
	ns=1         !prolate spheroid
	else
	ns=2         !oblate spheroid
	end if

	xmj=max(xrot,xperp)   !major semiaxis
	xmn=min(xrot,xperp)   !minor semiaxis

	xe=(xmj**2-xmn**2)/xmj**2
	xe=sqrt(xe)            !eccentricity

	if (ns.eq.1) then ! prolate

	xlz=log((1.d0+xe)/(1.d0-xe))
	xlz=(1.d0-xe**2)*(-1.d0 + xlz/(2.d0*xe))/xe**2
*
	xdz=1.d0 + xlz*(1.d0+xe**2)/(1.d0-xe**2)
	xdz=3.d0*xdz/4.d0
	xdx=(3.d0*log((1.d0+xe)/(1.d0-xe))/(2.d0*xe) - xdz)/2.d0 

	else if (ns.eq.2) then ! oblate

	xlz=1.d0 - sqrt(1.d0-xe**2)*dasin(xe)/xe
	xlz=xlz/xe**2
*
	xdz=1.d0 + xlz*(1.d0-2.d0*xe**2)
	xdz=3.d0*xdz/4.d0
	xdx=(3.d0*sqrt(1.d0-xe**2)*dasin(xe)/xe - xdz)/2.d0 

	end if

	xlx=(1.d0-xlz)/2.d0

	do 20 npol=1,2

	if (npol.eq.1) then        !polarization along the rotation axis
        xlp=xlz
	  xdp=xdz
	  xfx=1.d0-(xk*xperp)**2/10.d0
	  xapl=xrot
	else if (npol.eq.2) then   !perpendicular polarization 
	  xlp=xlx
	  xdp=xdx
	  xfx=1.d0-(xk*xrot)**2/10.d0
	  xapl=xperp
	end if

* Static Rayleigh polarizability:
	zalph(1)=xvol*(zeps-cone)/(cone+xlp*(zeps-cone)) 
	
* MLWA polarizability
	zalph(2)=zalph(1)/(cone-zalph(1)*xk**2/xmj 
     & - ci*2.d0*xk**3*zalph(1)/3.d0)

* MLWA polarizability with a ddepol factor:
	zalph(3)=zalph(1)/(cone-xdp*zalph(1)*xk**2/xapl 
     &  - ci*2.d0*xk**3*zalph(1)/3.d0)

	pf=0.37d0
	qf=1.d0-pf

* MLWA polarizability with averaged ddepol factor:
	zalph(4)=zalph(1)/(cone-(qf*xdp/xrot +pf/xapl)*zalph(1)*xk**2
     &  - ci*2.d0*xk**3*zalph(1)/3.d0)

* MLWA polarizability with a ddepol factor + averaging:
	zalph(5)=xfx*zalph(1)/(cone-xdp*xfx*zalph(1)*xk**2/xrot 
     &  - ci*2.d0*xk**3*xfx*zalph(1)/3.d0)

	do ij=1,5
        sext(ij)=4.d0*pi*xk*imag(zalph(ij))/xcr    ! extinction cross section
* 
* xcr=pi*rev**2 - an effective geom. cross section
 	
      enddo

	if (npol.eq.1) write(NOUT+16,1107) lambda, sext 
	if (npol.eq.2) write(NOUT+17,1107) lambda, sext 

 20	continue

 1107 FORMAT (F8.2,5(5X,D12.6))
	return
	end