function [results] = calc_int_kernel(sm,params)
% MATLAB HELP:
%
% function [results] = calc_int_kernel(sm,params)
%
% compute the 1d travel-time kernel 
%
% inputs:
%  sm, the solar model
%  params, the parameters for the kernel
% outputs:
%  results, the kernel information
%
%  all of these structures are complicated, see codes.ps for details
%
% aaron birch, may 21, 2002




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

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

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

%
% if we are going 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
%
N=length(n);
delta=params.delta;


%
% do the computation 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;


%
% dont need these variables anymore
%
clear U  
clear V  
clear D

%
% do the time grid
%

time=linspace(params.T_MIN,params.T_MAX,params.T_STEPS);
F=amp;
SF=sqrt(F); 

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

%
% calculate the isolation filter
%
[f,df,ddf]=isolation_filter(w,F,legendreP(l),l,time);

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

%
%integrate ddf*f over time
%
P=integrate(time,ddf.*f);


%
% allocate space for the result
%
k=zeros(1,length(r_vector)); 


%
% loop over l
%
for (the_l = 1:max(l))
  
  %
  % find all the modes with this l
  %
  Ind=find(l==the_l);
  
  
  %
  % if we have modes at this l
  %
  if (length(Ind)>0 ) 
    the_l,
    
    %
    % l dependent amplitude
    %
    A = -1*(2*the_l+1)/4/pi/P*AsympLegendreP(the_l,delta); 
    
    
    %
    % and get the matrix G only using modes with given l
    %
    G=calcg(df,time,w(Ind));

    %
    % add and do the kernel, loop over radius
    %
    for (R_INDEX=1:length(r_vector))
      k(R_INDEX)=k(R_INDEX)+A*(SF(Ind)'.*UI(Ind,R_INDEX)')*G*(UI(Ind,R_INDEX).*SF(Ind));
    end; %% end loop over R_INDEX

  end;  %% end if we have modes at this l
end;  %% end loop over l
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%
% this is the kernel for the square of the sound speed
% multiply by two to get the kernel for the sound speed
%
k=k.*rho(index).*(c(index).^2)*R^3.*(  (r_vector).^2  );

%
% put the results in nice form
%
results.k=k;
results.c=c(index);
results.rho=rho(index);
results.f=f;
results.r=r_vector;
results.rho=rho(index);
results.n=n;
results.l=l;
results.w=w;