function out=hedgeSummerSample(w,cMu,cVar,N)
% hedgeSummerSample generates lots of random summers with daily temperature
% deviations and cloud covers, measures the performance of the bet-hedging
% model and the adaptive tracking model under these conditions, and
% collects a number of metrics about the summer. Function is fairly
% hardwired to the parameters needed to make Figure 5D,E
%=========================================================================
% Usage:
% out=hedgeSummerSample(w,cMu,cVar,N)
% Inputs:
%       * w = struct of weather data, e.g. wp in hedgeWeatherPack.mat
%       * cMu = mean light choice probability
%       * cVar = variance in light choice probability
%       * N = number of summers to simulate
% Output, out=matrix of model performance and summer characteristics

% toggle this to 1 or 0 to either randomly vary the length of the season or not
varySummerLengthBool=1;

% b and d parameters from '001' weather conditions
params=[0.04480243 0.012755];

% standard temp diff between shade and sun
shadeDiff=7;

%initialize output matrix
if varySummerLengthBool==1
    out={};
else
    out=[];
end

%for each summer to be simulated
for i=1:N
    disp(i); %comment this line in to monitor output pace
    
    if varySummerLengthBool==1 %if varying season length
        numDays=ceil(rand*258+107); %number of days in this season [107, 365]
        normsOrig=w.normal(w.flySeasonInterval);
        cloudOrig=w.cloud(w.flySeasonInterval);
        w.normal=interp1(1:length(normsOrig),normsOrig,linspace(1,length(normsOrig),numDays))'; %interpolate normals
        w.cloud=interp1(1:length(cloudOrig),cloudOrig,linspace(1,length(cloudOrig),numDays))'; %interpolate clouds
        w.flySeasonInterval=1:length(w.normal);
    end
    
    %run analytic model for adaptive tracking model, with new daily temp
    %deviations and cloud cover vector
    tempH=hedgeAnalytic(1,w,cMu,cVar,params(1),params(2),0,0,shadeDiff,1,'111',0);
    
    %capture the daily devs and cloud data from that run
    wPack.tempHist=tempH.tempHist;
    wPack.cloudCover=tempH.cloudCover;
    
    %run analytic model for bet-hedging model, in pre mode, passing in the
    %weather data from above run
    tempP=hedgeAnalytic(0,w,cMu,cVar,params(1),params(2),0,0,shadeDiff,1,'pre',wPack);
    
    %calculate the autocorrelation function of the daily temp deviations
    XC=xcorr(wPack.tempHist-w.normal(w.flySeasonInterval),'coeff');
    
    %final populations of each model
    Hpop=sum(tempH.pops(:,end));
    Ppop=sum(tempP.pops(:,end));
    
    %percentage difference between them
    Pperc=(Ppop-Hpop)/Hpop;
    
    %mean and std of the temperatures and cloud covers
    mTemp=mean(wPack.tempHist);
    sTemp=std(wPack.tempHist);
    mCloud=mean(wPack.cloudCover);
    sCloud=std(wPack.cloudCover);
    
    %standard deviation of the "sun  potential" = temp + shadeDiff*cloud
    %cover fraction
    sTShade=std(wPack.tempHist+shadeDiff*wPack.cloudCover);
    
    %calculate power in three bands of the autocorrelation function
    if varySummerLengthBool~=1
        xc1=sum(XC(215:215+7));
        xc2=sum(XC(215:215+30));
        xc3=sum(XC(215+7:215+30));
    end
    
    %calculate histogram of daily temperatures in 2 degree bins
    num0=sum(wPack.tempHist>-2)-sum(wPack.tempHist>0);
    num1=sum(wPack.tempHist>0)-sum(wPack.tempHist>2);
    num2=sum(wPack.tempHist>2)-sum(wPack.tempHist>4);
    num3=sum(wPack.tempHist>4)-sum(wPack.tempHist>6);
    num4=sum(wPack.tempHist>6)-sum(wPack.tempHist>8);
    num5=sum(wPack.tempHist>8)-sum(wPack.tempHist>10);
    num6=sum(wPack.tempHist>10)-sum(wPack.tempHist>12);
    num7=sum(wPack.tempHist>12)-sum(wPack.tempHist>14);
    num8=sum(wPack.tempHist>14)-sum(wPack.tempHist>16);
    num9=sum(wPack.tempHist>16)-sum(wPack.tempHist>18);
    num10=sum(wPack.tempHist>18)-sum(wPack.tempHist>20);
    num11=sum(wPack.tempHist>20)-sum(wPack.tempHist>22);
    num12=sum(wPack.tempHist>22)-sum(wPack.tempHist>24);
    num13=sum(wPack.tempHist>24)-sum(wPack.tempHist>26);
    num14=sum(wPack.tempHist>26)-sum(wPack.tempHist>28);
    num15=sum(wPack.tempHist>28)-sum(wPack.tempHist>30);
    num16=sum(wPack.tempHist>30)-sum(wPack.tempHist>32);
    num17=sum(wPack.tempHist>32)-sum(wPack.tempHist>34);
    
    % calculate correlation scores between:
    corl0=corr2(wPack.tempHist,w.normal(w.flySeasonInterval));          % the temp history and the standard profile
    corl1=corr2(wPack.tempHist,wPack.cloudCover);                       % the temp history and the cloud cover
    corl2=corr2(wPack.cloudCover,w.normal(w.flySeasonInterval));        % the cloud cover and the standard profile
    corl3=corr2(wPack.tempHist-w.normal(w.flySeasonInterval),w.normal(w.flySeasonInterval));    %the temp deviations and the standard profile
    
    % calculate binned differences in:
    if varySummerLengthBool~=1
        sum1=sum(wPack.tempHist(1:107))-sum(wPack.tempHist(108:end));                                      %  total heat in the first and second half of summer
        sum2=sum(wPack.tempHist(54:160))-sum(wPack.tempHist(160:end))-sum(wPack.tempHist(1:53));           % middle of the summer minus start and end
        sum3=sum(wPack.cloudCover(1:107))-sum(wPack.cloudCover(108:end));                                  % cloud cover in halves
        sum4=sum(wPack.cloudCover(54:160))-sum(wPack.cloudCover(160:end))-sum(wPack.cloudCover(1:53));     % cluod cover in middle vs end
    end
    
    % save all these calculated terms into a new vector
    temp=[Hpop Ppop Pperc mTemp sTemp mCloud sCloud ];%1-7
    temp=[temp sTShade];
    if varySummerLengthBool~=1; temp=[temp xc1 xc2 xc3]; end  %9-11
    temp=[temp num0 num1 num2 num3 num4 num5 num6 num7 num8 num9]; %12-21
    temp=[temp num10 num11 num12 num13 num14 num15 num16 num17]; %22-29
    temp=[temp corl0 corl1 corl2 corl3]; %30-33
   if varySummerLengthBool~=1;  temp=[temp sum1 sum2 sum3 sum4]; end;  %34-37
    
    %save this row of scores into the output matrix.
    if varySummerLengthBool==1
        out=[out;[temp NaN numDays NaN {wPack.tempHist'}]];
    else
        out=[out;[temp NaN XC' NaN periodogram(wPack.tempHist-w.normal(w.flySeasonInterval))' NaN wPack.tempHist']];
    end
end

% sound(0.2*sin(.5*(1:1000))) %comment this in to beep when done.