function result = mie_esquare2(m, x, n)

% Test of term-wise equality of absorption efficiency
% n is term number, x size parameter,
% m refractive index of sphere, Sept. 13, 2002

nj=100*round(2+x+4*x.^(1/3))+100;
nu =(n+0.5); 
m1=real(m); m2=imag(m);e2=imag(m.*m);
abcd=mie_ab(m,x);
an=abcd(1,n);bn=abcd(2,n);
an2=abs(an).^2;
bn2=abs(bn).^2;

abcd=mie_cd(m,x);
cn=abcd(1,n);dn=abcd(2,n);
cn2=abs(cn).^2;
dn2=abs(dn).^2;
dx=x/nj;
es=[];
for j=1:nj,
    xj=dx.*j;
    z=m.*xj;
    sqz= sqrt(0.5*pi./z);
    bz = besselj(nu, z).*sqz;      % This is j_n(z)
    bz2=(abs(bz)).^2;
    b1z = besselj(nu-1, z).*sqz;   % This is j_n-1(z)
    az = b1z-n.*bz./z;
    az2=(abs(az)).^2;
    z2=(abs(z)).^2;
    n1 =n*(n+1);
    n2 =2*(2*n+1);
    mn=real(bz2.*n2);
    nn1=az2;
    nn2=bz2.*n1./z2;
    nn=n2.*real(nn1+nn2);
    es=[es 0.25*(cn2*mn+dn2*nn)];
end;
xxj=[0:dx:x];
een=[es(1) es];

plot(xxj,een);
title(sprintf('Sq. Ampl. Field in Sphere of Mode n=%g, m=%g+%gi, x=%g',n,m1,m2,x))
xlabel('r k')

x2=x.*x;
nj1=nj+1;
en1=0.5*een(nj1).*x2;     % End-Term correction in integral
enx=een*(xxj.*xxj)'-en1;    % Trapezoidal radial integration
inte=dx.*enx;
Qabs=4.*e2.*inte./x2;
Qab=(n2/x2)*(real(an)-an2+real(bn)-bn2);
result=[Qabs;Qab];