% Jon Newman
% Turn OU process into analog stimulation waveform for ChR2 excitation

%% Exp. Parameters
% We will make dt 1 millisecond so we get ff = 1000 Hz;
mea = 'gimli-19305';
rep = 100;
IappW = 100; %msec
F = 2; %Hz
start_time = 0; % sec
dt = 0.001; % discretisation in time, needs to be adjusted to sampling frequency (1/ff)
stim_time = 10; % X seconds of stimulation
tauI = 5; % correlation time of the OU process
u_noise  = 2.5; % sets mean control voltage in the center of the input range
s_noise = 0.625; % ensures we will not staturate the light ofte
wakeUpTime = 0; %seconds
fid = [ datestr(now,29) '_' mea '_OU-stim_tau'  num2str(tauI) 'sec-u-' num2str(u_noise) '-s-' num2str(s_noise) '_00'];

%% Derive OU snippet
N=ceil(stim_time/dt); % number of noise points
rset=randn(N,1);
t_noise = (0:dt:stim_time-dt);
e_noise=zeros(1,N);

% make unit std deviation noise
A_e = exp(-dt/tauI);
sigma=1; % std of the noise
for s = 1:(N-1)
    e_noise(s+1) = e_noise(s)*A_e + sigma*sqrt(1-A_e^2)*rset(s);
end
% scale and shift
e_noise = s_noise*e_noise + u_noise;

% storage
T = [0 start_time];
C = [0 0];
A = [-0.5 -0.5];

Td = [0];
Sd = [0];

%% Make repeated OU/current injection processes
for i = 1:rep
%     if mod(i,2)
        T = [T T(end) + t_noise];
        A = [A e_noise];
        C = [C zeros(1,length(e_noise))];
        
        b_trig_on = dec2bin(1,32);
        b_trig_off = dec2bin(9,32);
        
        % Trigger a new stimulation cycle on bit 0
        Td = [Td T(end) + t_noise(1)];
        Sd = [Sd 1];
        Td = [Td T(end) + t_noise(1) + 0.0001];
        Sd = [Sd 0];
        
        
%     else
%         b_space = dec2bin(0,31);
%         b_IappOn = ['1' b_space];
%         b_IappOff = ['0' b_space];
%         int_on = bin2dec(b_IappOn);
%         int_off = bin2dec(b_IappOff);
%         t1 = 1/F:1/F:stim_time;
%         t0 = t1 + IappW/1000;
%         [td ind] = sort([t1 t0]);
%         
%         s1 = int_on * ones(size(t1));
%         s0 = int_off * ones(size(t1));
%         s = [s1 s0];
%         s = s(ind);
%         
%         % Trigger a new stimulation cycle on bit 0
%         Td = [Td T(end) + t_noise(1)];
%         Sd = [Sd 1];
%         Td = [Td T(end) + t_noise(1) + 0.0001];
%         Sd = [Sd 0];
%         
%         % Mark current injections
%         Td = [Td T(end)+td];
%         Sd = [Sd s];
%         T = [T T(end) + t_noise];
%         A = [A e_noise];
%         C = [C zeros(1,length(e_noise))];
%     end
    
    
    
end

A(end) = -0.5;

%% Create a muliplier function to taper the onset of the OU over the first 60 seconds
T1 = T(T<=60);
T2 = T(T>60);
KK = [(1.0/60)*(T1-start_time) ones(size(T2))];
A = A.*KK;


%% Plot
close all
figure
subplot(211)
plot(T,A)
xlabel('time (sec)')
ylabel('V_{CTL} (volts)')
subplot(212)
stairs(Td,Sd)
xlabel('time (sec)')
ylabel('Digital port state (uint32)')

%% Make the NeuroRighter open loop stimulation files

makedigfile(fid,Td',Sd');
makeauxfile(fid,T',C',A');





