function [STA, t] = stim_sta(T,U,c1,c2, bounds,dt)
%STIM_STA Spike triggered average of whisker stimulus waveform and optical
%input

% Initialization
uu = unique(U); % unique units

t = bounds(1):dt:bounds(2);
% wave_n = 1:1:numel(c1);
STA = zeros(numel(uu),size(t,2),2);
% Get windows

for j = 1:numel(uu)
    
    % Spike times for unit j
    TT = T(U == uu(j));
    
    bounds = bounds(ones(size(TT)),:);
    TT = [TT, TT];
    bound_n = ceil((bounds + TT)/dt);
    K = 0;
    
    % For all spikes occuring in at times TT
    for i = 1:size(bound_n,1)
        
        % Current window
        n = bound_n(i,1):1:bound_n(i,2);
        
        try
            % Optical power
            STA(j,:,1) = STA(j,:,1) + c1(n);
            
            % Whisker stimulus waveform
            STA(j,:,2) = STA(j,:,2) + c2(n);
            
            K = K + 1;
            
        catch
            continue
        end
        
    end
    
    % Normalize by number of spikes
    STA(j,:,:) = STA(j,:,:)./K;
end


end

