function stdp3
%used to induce spike-timing-dependent long-term synaptic plasticity  -- presyn cell hold in Vclamp, postsyn cell switched to Iclamp
%stdp3 is the file used to output pairs, pre-post or post-pre, or triplets, either as pre-post-pre or as post-pre-post

hdfml('closeall')
clear all
%____________________
%CHECK GAIN FOR THE CURRENT-CLAMPED CHANNEL
output_channel_to_presynaptic=2 %tells which electrode is being treated as on the presynaptic side in MATLAB channel numbering (output channel 1 is hardware channel 0 etc.)

stimulation_pattern=1          %1=pre-post, 2=post-pre, 3=pre-post-pre, 4=post-pre-post
pulse_duration_presyn=2.2;       %in milliseconds
pulse_duration_postsyn=2;    %in milliseconds

timing_offset=10                %in ms, how much time btwn stims
OVSF_presyn=5                  %output voltage scaling factor to voltage-clamped (presynaptic) neuron
OVSF_postsyn=5                 %SET at 5 for 2nA current injection to postsynaptic neuron
mcc_output_gain_presyn=5       %set to match the output gain setting selected on multiclamp commander
mcc_output_gain_postsyn=5
recording_status=1;
number_of_cycles=60           %60 for standard induction protocol
cycle_duration=1;              %in seconds
sampling_rate=10000;           %in kHz

%____________________
samples_per_ms=sampling_rate/1000;
sample_offset=timing_offset*samples_per_ms;samples_presyn=pulse_duration_presyn*samples_per_ms;samples_postsyn=pulse_duration_postsyn*samples_per_ms;
samples_per_cycle=cycle_duration*sampling_rate;
cycle_midpoint=samples_per_cycle/2;
output_channel_to_postsynaptic=1/(output_channel_to_presynaptic/2); %or use =bitcmp(output_channel_to_presynaptic,1);
picoamp_converter=2000/mcc_output_gain_presyn;
millivolt_converter=1000/mcc_output_gain_postsyn;
arbitrary_postsyn_scaler=10;        %used in plotting presynaptic membrane current and postsynaptic membrane potential on the same plot

stdpstim_presyn=zeros(samples_per_cycle,1); 
stdpstim_postsyn=zeros(samples_per_cycle,1);
if stimulation_pattern==1
    stdpstim_presyn(cycle_midpoint-(sample_offset/2)+1:cycle_midpoint-(sample_offset/2)+samples_presyn)=1;
    stdpstim_postsyn(cycle_midpoint+(sample_offset/2)+1:cycle_midpoint+(sample_offset/2)+samples_postsyn)=1;
elseif stimulation_pattern==2
    stdpstim_presyn(cycle_midpoint+(sample_offset/2)+1:cycle_midpoint+(sample_offset/2)+samples_presyn)=1;
    stdpstim_postsyn(cycle_midpoint-(sample_offset/2)+1:cycle_midpoint-(sample_offset/2)+samples_postsyn)=1;
elseif stimulation_pattern==3
    stdpstim_presyn(cycle_midpoint-sample_offset+1:cycle_midpoint-sample_offset+samples_presyn)=1;
    stdpstim_presyn(cycle_midpoint+sample_offset+1:cycle_midpoint+sample_offset+samples_presyn)=1;
    stdpstim_postsyn(cycle_midpoint+1:cycle_midpoint+samples_postsyn)=1;
elseif stimulation_pattern==4
    stdpstim_presyn(cycle_midpoint+1:cycle_midpoint+samples_presyn)=1;
    stdpstim_postsyn(cycle_midpoint-sample_offset+1:cycle_midpoint-sample_offset+samples_postsyn)=1;
    stdpstim_postsyn(cycle_midpoint+sample_offset+1:cycle_midpoint+sample_offset+samples_postsyn)=1;
elseif stimulation_pattern==5
    stdpstim_presyn(cycle_midpoint+1:cycle_midpoint+samples_presyn)=1;
    stdpstim_postsyn(cycle_midpoint+sample_offset+1:cycle_midpoint+sample_offset+samples_postsyn)=1;
    stdpstim_postsyn(cycle_midpoint+2*sample_offset+1:cycle_midpoint+2*sample_offset+samples_postsyn)=1;
    stdpstim_postsyn(cycle_midpoint+3*sample_offset+1:cycle_midpoint+3*sample_offset+samples_postsyn)=1;
    stdpstim_postsyn(cycle_midpoint+4*sample_offset+1:cycle_midpoint+4*sample_offset+samples_postsyn)=1;
end
stdpstim_presyn=stdpstim_presyn*OVSF_presyn;
stdpstim_postsyn=stdpstim_postsyn*OVSF_postsyn;
stdp_output=[stdpstim_presyn stdpstim_postsyn];

if output_channel_to_presynaptic==2
       stdp_output=fliplr(stdp_output);
end

cyclesize=size(stdp_output)
stdp_aatput=repmat(stdp_output,number_of_cycles,1);
aatputsize=size(stdp_aatput)

figure(15)
set(15,'Position',[-49 -145 672 504])
plot(stdp_output)
axis([cycle_midpoint-(2*sample_offset) cycle_midpoint+(2*sample_offset) -10 10]) 
% figure(21)
% set(21,'Position',[641 -144 672 504])
% plot(stdp_aatput)

%cleanup the workspace before getting underway
synch=[];
delete(synch)
clear synch
stdp_output_obj=[];
delete(stdp_output_obj)
clear stdp_output_obj
stdp_input_obj=[];
delete(stdp_input_obj)
clear stdp_input_obj

%building the output object
stdp_output_obj=analogoutput('nidaq');
addchannel(stdp_output_obj,0);            
addchannel(stdp_output_obj,1);            
%set(stdp_output_obj,'TransferMode','SingleDMA')
set(stdp_output_obj,'TransferMode','Interrupts')
%set(stdp_output_obj,'TriggerType','Immediate')
set(stdp_output_obj,'TriggerType','Manual')
set(stdp_output_obj,'SampleRate',sampling_rate)
%actualout=get(stdp_output_obj,'SampleRate')           %these two lines are for debugging/checks on stdp_output_obj
%get(stdp_output_obj)


%building the input object
stdp_input_obj=analoginput('nidaq');
addchannel(stdp_input_obj,0);
addchannel(stdp_input_obj,1); 
addchannel(stdp_input_obj,2);               %use to collect temperature data
set(stdp_input_obj,'SampleRate',sampling_rate)       
set(stdp_input_obj,'SamplesPerTrigger',samples_per_cycle*number_of_cycles)
set(stdp_input_obj,'ChannelSkewMode','Minimum')
set(stdp_input_obj,'TransferMode','SingleDMA') 
%set(stdp_input_obj,'TransferMode','Interrupts')
%set(stdp_input_obj,'TriggerType','Immediate')
set(stdp_input_obj,'TriggerType','Manual')
if recording_status                         %this is part of matlab DAQ logging, the files created for which will be used for the .hdf patch log
    set(stdp_input_obj,'LoggingMode','Disk&Memory')
else
    set(stdp_input_obj,'LoggingMode','Memory')
end

%building a vector of input and output objects to synch them with a trigger
synch=[stdp_output_obj stdp_input_obj];

putdata (stdp_output_obj,stdp_aatput)
if recording_status
    patchlog_name=create_patchlog_name;
    stdplog_name=strcat('stdp_',patchlog_name) 
    sd_id_stdp=hdfsd('start',stdplog_name,'create');        %this sd file will be written to by update_patchlog.m
    statupdate_end_sd_id_stdp=hdfsd('end',sd_id_stdp);                    %will be opened update_patchlog when it is called below
    trial = 1;
    time_since_prev_trial=stimulation_pattern;                     %for use by generate_log_filename  -- beginning aug 1 2004 will input stimulation_pattern, instead of 0, here
    log_filename=generate_log_filename(time_since_prev_trial)
    set(stdp_input_obj,'LogFileName',log_filename)
end
start(synch)
tic
trigger(synch)

i = 1;
while stdp_input_obj.SamplesAcquired < stdp_input_obj.SamplesPerTrigger
    while stdp_input_obj.SamplesAcquired < samples_per_cycle*i                   %wait out the first cycle
    end
    datums=peekdata(stdp_input_obj,samples_per_cycle);
    presyn_peek=datums(:,output_channel_to_presynaptic)*picoamp_converter;       %converting to picoamps
    postsyn_peek=datums(:,output_channel_to_postsynaptic)*millivolt_converter;   %converting to millivolts
    figure(16)
    clf
    plot(presyn_peek,'r')
    hold on
    plot(postsyn_peek*arbitrary_postsyn_scaler,'b')
    xlim([.4*samples_per_cycle .55*samples_per_cycle])
    ylim([-8000 1000])
    xlabel('sample')
    ylabel(strcat('presynaptic mem current (pA) and postsynaptic mem potential (mV) *',num2str(arbitrary_postsyn_scaler)))
    %grid on
    %set(16,'doublebuffer','on')
    title([sprintf('Peekdata calls (presynaptic red, postsynaptic black): '), num2str(i)])
    legend(stdplog_name)
    drawnow
    i = i + 1;
end
% while strcmp(stdp_input_obj.Running,'On')
% end
toc
%get(synch,'Running')

%saving the induction process (a single 'trial', containing 'number_of_cycles' cyclic repetitions)
if recording_status
    update_patchlog(trial,time_since_prev_trial,log_filename,OVSF_presyn,OVSF_postsyn,pulse_duration_presyn,pulse_duration_postsyn,sampling_rate,...
    timing_offset,output_channel_to_presynaptic,stdplog_name)
end

%cleanup
synch=[];
delete(synch)
clear synch
stdp_output_obj=[];
delete(stdp_output_obj)
clear stdp_output_obj
stdp_input_obj=[];
delete(stdp_input_obj)
clear stdp_input_obj
hdfml('closeall')
clear