Permalink
Switch branches/tags
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
77 lines (73 sloc) 2.06 KB
function [I,dph,th,thR,th0,ph0]=sincos(th,m1,m2,dom,pars)
% [I,dph,th,nthR,th0,ph0]=sincos(th,m1,m2,dom,pars)
%
% Calculates the longitudinal integral
% \int_{A}^{B}\sin(m1\phi)\cos(m2\phi)\.d\phi
% over a specified domain at some (Gauss-Legendre) integration point(s).
%
% INPUT:
%
% th Colatitude(s) where you want it evaluated
% m1,m2 Angular order of the cosine factors
% dom A matrix with the hatchings from PHICURVE, or
% 'patch' Spherical patch [default] with the following specs
% 'england' England, Scotland and Wales
% pars [th0,ph0,thR] for 'patch', with
% th0 Colatitude of the cap center, in radians
% ph0 Longitude of the cap center, in radians
% thR Radius of the cap, in radians
% N for 'england' with N the smoothness of the splining
%
% EXAMPLE:
%
% sincos('demo1');
%
% Last modified by fjsimons-at-alum.mit.edu, 04/20/2009
defval('th',[])
defval('m1',10);
defval('m2',4);
defval('dom','patch')
if isstr(dom)
switch dom
case 'patch'
defval('pars',[pi/4 pi/2 pi/9]);
th0=pars(1); ph0=pars(2); thR=pars(3);
defval('th',linspace(th0-thR,th0+thR,100));
phint=dphpatch(th,thR,th0,ph0);
case 'england'
defval('N',10)
phint=dphengland(th*180/pi,N);
phint=phint*pi/180;
otherwise
error('Specify valid domain')
end
else
phint=dom;
end
I=repmat(0,size(phint,1),1);
for index=1:size(phint,2)/2
A=phint(:,2*index-1); B=phint(:,2*index);
if m1~=m2 & m1~=0 & m2~=0 & m1~=-m2
I=I+1/2*(-cos((m1+m2)*B)*m1+...
cos((m1+m2)*B)*m2-...
cos((m1-m2)*B)*m1-...
cos((m1-m2)*B)*m2+...
cos((m1+m2)*A)*m1-...
cos((m1+m2)*A)*m2+...
cos((m1-m2)*A)*m1+...
cos((m1-m2)*A)*m2)/...
((m1-m2)*(m1+m2));
elseif m1~=m2 & m1==0
I=I+0;
elseif m1~=m2 & m2==0
I=I-(cos(m1*B)-cos(m1*A))/m1;
elseif (m1==m2 | (m1==-m2 & m1>0)) & m1~=0
m=m1;
I=I-1/2*(cos(m*B).^2-cos(m*A).^2)/m;
elseif (m1==m2 | (m1==-m2 & m2>0)) & m1~=0
m=m2;
I=I+1/2*(cos(m*B).^2-cos(m*A).^2)/m;
else
I=I+0;
end
end