%% Exercise 9 Part 1
%needs to be run after Exercise 8 Part 1
close all;
clc;
diff = [];
ModelCount = [];
h = [];
decay=@(tal,tpfr, x) -1/tal*(x-tpfr) + log(1-exp(-0.5/tal));
for i=1:44 %Loop through all the peaks
    y=[]; %Reset x and y on every iteration
    x=[];
    k = 0;
    l=0;
    N=0; %Find the value for N for each peak
    tpfr(i)=Time(TriggeredOnsetTime(i))-TrigTime(i);
    for j=1:(OTimes(i+1,1)-PTimes(i,1))*2 %Time between peak of the current and onset of the next, in half seconds
        if Counts(j+TriggeredPeakTime(i))~=0 %Only count nonzero points
            y=[y; (Counts(j+TriggeredPeakTime(i)))]; %Store log of fraction of total particles in CSTR2
           	x=[x; 1, Time(j+TriggeredPeakTime(i))-TrigTime(i)]; %Store negative of the time
            N = N + Counts(j+TriggeredPeakTime(i)); %calculates total count over one peak
        end
    end
    y(:,1) = y(:, 1)/N; %normalizes the actual counts
    for j =1:(OTimes(i,1)-TrigTime(i, 1))*2
        Tsb(j) = Counts(j+TriggeredTime(i)); %extracts counts in bins before onset of peak
    end
    Sbaseline = mean(Tsb); %calculates baseline value
    ModelCount = [];
    h = [];
    for j =1:(OTimes(i,1)-TrigTime(i, 1))*2 %calculates differences between model and actual for points before onset
        if Counts(j+TriggeredTime(i)) ~=0
            k = k+1;
            ModelCount(k) = Sbaseline/N;
            diff(k) = (Counts(j+TriggeredTime(i)) - ModelCount(k));
        end
    end

    %calculates difference between actual and model after onset
    for j=1:(OTimes(i+1,1)-PTimes(i,1))*2 %Time between peak of the current and onset of the next, in half seconds
        if Counts(j+TriggeredPeakTime(i))~=0 %Only count nonzero points
            k = k+1;
            l = l+1;
            ModelCount(k) = exp(decay(T(1,i), tpfr(i), l));
            diff(k) = (y(l) - ModelCount(k));
        end
    end
        h(i)= chi2gof(diff, 'Alpha', 0.05);
end