function [results]=calc_kernel(sm,params);
% MATLAB HELP:
%
% function [results]=calc_kernel(sm,params);
%
% inputs:
%  sm, the solar model
%  params, the parameters for the kernel we want
% outputs
%  results, a structure containing the kernel and the grids, etc.
%
% these input/outputs are pretty complicated, see codes.ps or
% codes.pdf for more info
%
% aaron birch, may 21, 2002


%
% unpack what we need from the solar model structure
%
R=sm.R;
r=sm.r;
c=sm.c;
rho=sm.rho;
n=sm.n;
l=sm.l;
w=sm.w;
U=sm.U;
V=sm.V;
clear sm

%
% set up the grid in angle
%
phi_vector=linspace(-0.5*params.delta,1.5*params.delta,params.NPHI);

%
% set up the grid in radius
%
min_index=min(find(r>params.MINR));
index=round(linspace(min_index,length(r),params.NR));

%
% pick out the modes that we are going to compute with
%
[n,l,w,amp,I]=select_modes(n,l,w,params);

%
% if we wanted to skip some modes
%
SKIP=params.SKIP;
n=n(1:SKIP:length(n));
l=l(1:SKIP:length(l));
w=w(1:SKIP:length(w));
amp=amp(1:SKIP:length(amp));
I=I(1:SKIP:length(I));
U=U(I,:);
V=V(I,:);


%
% define for later convenience
%
N=length(n);
delta=params.delta;

%
% setup the grid in radius
% 
r_vector=r(index);

%
% do the computation of the divergence of the eigenfunctions
% 
UI=zeros(N,length(index));
D=zeros(size(U));
for (i=1:N)
	D(i,:)=gradient(U(i,:),R*r)+(r*R).^(-1).*(2*U(i,:)-sqrt(l(i)*(l(i)+1))*V(i,:));
	UI(i,:)=D(i,index);
end;


%
% don't need these variable any more
%
clear U  
clear V  
clear D

%
% do the time grid
%
time=linspace(params.T_MIN,params.T_MAX,params.T_STEPS);

%
% the filter
%
F=amp;
SF=sqrt(F);  

%
% calculate the legendre polynomials at cos(delta)
%
legendreP=calc_legendre(max(l),cos(delta));

%
% and then calculate the zero order signal at delta
%
[f,df,ddf]=isolation_filter(w,F,legendreP(l),l,time);

%
% if we wanted to only compute the isolation filter
%
if (params.CALC_F_ONLY)
  return
end;




%%%%%%%% introduce some angles %%%%%%%%%%

%
% for a cut containing the source and receiver
%
if (params.CUT==0)
for (PHI_INDEX=1:length(phi_vector) )
  if (phi_vector(PHI_INDEX)>0 & phi_vector(PHI_INDEX)<pi )
     phi_ps(PHI_INDEX)=phi_vector(PHI_INDEX);
  else
     phi_ps(PHI_INDEX)=-1*phi_vector(PHI_INDEX);
  end

  if (phi_vector(PHI_INDEX)>delta-pi & phi_vector(PHI_INDEX)<delta )
     phi_pr(PHI_INDEX)=delta-phi_vector(PHI_INDEX);
  else
     phi_pr(PHI_INDEX)=phi_vector(PHI_INDEX)-delta;
  end
end;
end;
%
%
%
if (params.CUT==1)
for (PHI_INDEX=1:length(phi_vector) )
 phi_ps(PHI_INDEX)=acos(sin(pi/2+(delta/2-phi_vector(PHI_INDEX)))*cos(delta/2));
 phi_pr(PHI_INDEX)=phi_ps(PHI_INDEX);
end;
end;




%
% and get the coupling matrix G
%
G=calcg(df,time,w);

%
% calculate P
%
P=sum( ddf.*f)*(time(2)-time(1)); %% integrate ddf*f over time

%
% allocate space for the results
%
k=zeros(length(phi_vector),length(r_vector));

%
% calculate all the legendre polynomials that we need
%
legendre_ps = calc_legendre(max(l),cos(phi_ps));
legendre_pr = calc_legendre(max(l),cos(phi_pr));



%
% loop over phi
%
for (PHI_INDEX=1:length(phi_vector) )  

  tic,

  GEO_PS=(2*l+1)/(4*pi).*legendre_ps(l,PHI_INDEX);
  GEO_PR=(2*l+1)/(4*pi).*legendre_pr(l,PHI_INDEX);

  %
  % loop over r
  %
  for (R_INDEX=1:length(r_vector))
    k(PHI_INDEX,R_INDEX)=(GEO_PS'.*SF'.*UI(:,R_INDEX)')*G*(UI(:,R_INDEX).*SF.*GEO_PR);
  end; 
  


  toc,
end; %% end loop over phi


%
% put in the factors that depend only on radius
%

%
% this is the kernel for the square of the sound speed
%
size(rho),size(c),size(k(1,:)),
for (PHI_INDEX=1:length(phi_vector))
	k(PHI_INDEX,:)=-1*rho(index).*(c(index).^2).*k(PHI_INDEX,:)/P;
end;

%
% stash all the results that we want into the results structure 
%
results.c=c(index);
results.rho=rho(index);
results.k=k;
results.n=n;
results.l=l;
results.w=w;
results.r=r(index);
results.phi=phi_vector;
results.CUT=params.CUT;
results.f=f;

%
% ALL DONE!
%