function I=sincosexp(m1,m2,mp)
% I=SINCOSEXP(m1,m2,mp)
%
% INPUT: 
%
% m1,m2,m3     Three spherical harmonic orders
% 
% Calculates the integral of 
%
% \int_{0}^{pi}\sin(m1\phi)\cos(m2\phi) exp(i*mp*phi)\,d\phi
%
% Silly me - this integral is no good: usually, we want 2*pi
%
% Last modified by fjsimons-at-alum.mit.edu, 08.08.2006

cospim1=cos(pi*m1); cospim2=cos(pi*m2);
sinpim1=sin(pi*m1); sinpim2=sin(pi*m2);
expimp=exp(mp*pi*i);

itop=-[-expimp*mp*m1^2*sinpim1*cospim2*i+expimp*m1^3*cospim1*cospim2+2* ...
	i*expimp*mp*m1*m2*cospim1*sinpim2-expimp*mp*m2^2*sinpim1*cospim2* ...
	i-expimp*m2^3*sinpim1*sinpim2+expimp*mp^3*sinpim1*cospim2*i- ...
	expimp*mp^2*m1*cospim1*cospim2+expimp*m1^2*m2*sinpim1*sinpim2- ...
	expimp*m1*m2^2*cospim1*cospim2+expimp*m2*mp^2*sinpim1*sinpim2+ ...
	mp^2*m1-m1^3+m1*m2^2];
ibot=[mp^4-2*mp^2*m1^2-2*mp^2*m2^2+m1^4-2*m1^2*m2^2+m2^4]; 

% Limiting case using L'Hopital
if abs(ibot)<eps
  disp('Using L''Hopital on m1')
  % Derivative wrt m1; could use m2 or mp, why not?
  itop=[expimp*mp*m2^2*cospim1*pi*cospim2*i+2*i*expimp*mp*m1*sinpim1* ...
	cospim2-3*expimp*m1^2*cospim1*cospim2+expimp*m1^3*sinpim1*pi*cospim2+ ...
	2*i*expimp*mp*m1*m2*sinpim1*pi*sinpim2+expimp*mp*m1^2*cospim1*pi* ...
	cospim2*i-expimp*mp^3*cospim1*pi*cospim2*i+expimp*m2^3*cospim1*pi* ...
	sinpim2-2*i*expimp*mp*m2*cospim1*sinpim2+expimp*mp^2*cospim1* ...
	cospim2-expimp*mp^2*m1*sinpim1*pi*cospim2-2*expimp*m1*m2*sinpim1* ...
	sinpim2-expimp*m1^2*m2*cospim1*pi*sinpim2+expimp*m2^2*cospim1*cospim2- ...
	expimp*m1*m2^2*sinpim1*pi*cospim2-expimp*m2*mp^2*cospim1*pi*sinpim2- ...
	mp^2+3*m1^2-m2^2];
  ibot=[-4*mp^2*m1+4*m1^3-4*m1*m2^2];
end

I=itop/ibot;