%Description: 
%
% This program generates the summer energy tables accompanying this
% submission. Summer energies (in Giga-Joules/m^2) are reported at each five
% degree latitude increment from 90S to 90N in an individual table.  Rows in
% each table correspond to kila-years before present (0 to 5000ky BP) and columns to
% different threshold values (0 to 600 W/m^2 at 25 W/m^2 increments).  Note
% that computation time, depending upon the machine, can take many hours.


clear;
lat=-90:5:90;   %latitude
tau=0:25:600;   %threshold values in W/m^2.
kyear=0:5000;   %thousands of years before present 
day=1:.1:365.2; %day of the year, where day 1 is January 1st.
[kyear_grid day_grid]=meshgrid(kyear,day);

for ct_lat=1:length(lat),
  lat(ct_lat),
  clear J; 
  
  %Calculate integrated summer energy
  Fsw=daily_insolation(kyear_grid,lat(ct_lat),day_grid);
  for ct_tau=1:length(tau),
    Fsw_temp=Fsw; 
    Fsw_temp(Fsw_temp<tau(ct_tau))=0;
    J(ct_tau,:)=sum(Fsw_temp,1)*(1/10)*(24*3600)/1e9; 
    %values are in (W/m^2) * (1 day / 10 increment) * (24*3600 seconds/day)
    %*(1 Giga-Joule / 1e9 Joules) gives output in Giga-Joules/m^2. 
  end;
  
  %add column indicating kila-year
  J=[kyear' J'];
  
  %write file
  if lat(ct_lat)<0,
    eval(['fid = fopen(''J_',num2str(abs(lat(ct_lat))),'South.txt'',''wt'');']);
  else,
    eval(['fid = fopen(''J_',num2str(abs(lat(ct_lat))),'North.txt'',''wt'');']);
  end;
    fprintf(fid,'Table of integrated summer insolation. \n');
    fprintf(fid,'The name of the file indicates the latitude. \n'); 
    fprintf(fid,'Rows correspond to time and columns to diurnal average insolation thresholds. \n'); 
    fprintf(fid,'Time (kilayears before present) is indicated in the first column and thresholds (W/m^2) in the first row. \n\n');
    fprintf(fid,'Peter Huybers, Harvard University, 2006.  (phuybers@fas.harvard.edu)\n\n\n');
    
    fprintf(fid,'%5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f %5.0f\n',[NaN tau]');
  fprintf(fid,'%5.0f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n',J');
  fclose(fid);
end;