%function I=get_inso(ky,phi2); 
%  
%Updates insolation forcing using a linear
%combination of the adjacent Ky intervals.

function I=get_inso(ky,phi2); 

  %Convert years in ky before present
  ky=-ky/1000;   
  
  if nargin<2, 
    phi2=65; 
  end;
  
  if ky==round(ky); 
    I=load_inso(ky,phi2); 
  else,
    I1=load_inso(floor(ky),phi2); 
    I2=load_inso(ceil(ky),phi2);
    I=(ceil(ky)-ky)*I1+(ky-floor(ky))*I2; 
  end;
  
 
function I=load_inso(ky,phi2); 

  if exist(['Inso_',num2str(ky),'x',num2str(length(phi2)),'.mat'])==2;
    eval(['load Inso_',num2str(ky),'x',num2str(length(phi2)),'.mat I']);
  else
    display('Computing and saving insolation');
    I=inso2(ky,phi2,1:365);
    eval(['save Inso_',num2str(ky),'x',num2str(length(phi2)),'.mat I']);
  end;