diff --git a/buildAnalysisParamFile.m b/buildAnalysisParamFile.m
index 2f6080c..9e84f61 100644
--- a/buildAnalysisParamFile.m
+++ b/buildAnalysisParamFile.m
@@ -5,9 +5,9 @@
 
 
 %%%%%%%  USER PARAMETERS, EDIT ROUTINELY %%%%%%%%
-runNum = '001';
-dateSubject = '171004ALAN'; 
-machine = 'laptop';
+runNum = '005';
+dateSubject = '171009ALAN'; 
+machine = 'rig';
 
 switch machine
   case 'rig'
@@ -78,7 +78,7 @@
 analogInParams.analogInChannels = [138,139,140]; 
 analogInParams.channelNames = {'eyeX','eyeY','eyeD'};
 analogInParams.analogInChannelScaleBy = [5/32764 5/32764 5/32764]; %converts raw values to volts
-analogInParams.decimateFactorPass1 = 1; %note: product of the two decimate factors should be 30, if 1 khz samples desired
+analogInParams.decimateFactorPass1 = 1; 
 analogInParams.decimateFactorPass2 = 1;
 analogInParams.samPerMS = 1; %THIS IS AFTER DECIMATION, and applies to analogIn (should be raw rate/productOfDecimateFactors)
 % see http://www.mathworks.com/help/signal/examples/filter-design-gallery.html
@@ -197,11 +197,11 @@
 plotSwitch.prefImRaster = 0;
 plotSwitch.prefImRasterEvokedOverlay = 0;
 plotSwitch.prefImMultiChRasterEvokedOverlay = 0;
-plotSwitch.imageTuningSorted = 1;
+plotSwitch.imageTuningSorted = 0;
 plotSwitch.stimPrefBarPlot = 1;
 plotSwitch.stimPrefBarPlotEarly = 0;
 plotSwitch.stimPrefBarPlotLate = 0;
-plotSwitch.tuningCurves = 1;
+plotSwitch.tuningCurves = 0;
 plotSwitch.tuningCurvesEarly = 0;
 plotSwitch.tuningCurvesLate = 0;
 plotSwitch.RF = 1;
diff --git a/checkAnalysisParamFile.m b/checkAnalysisParamFile.m
new file mode 100644
index 0000000..2cf56c0
--- /dev/null
+++ b/checkAnalysisParamFile.m
@@ -0,0 +1,394 @@
+function [ output_args ] = checkAnalysisParamFile( analysisParamFilename )
+%UNTITLED3 Summary of this function goes here
+%   Detailed explanation goes here
+load(analysisParamFilename);
+logString = 'Assigned default values to: \n';
+
+% runNum
+assert(logical(exist('runNum','var')),'Invalid analysis parameter file: must specify a run number as runNum.');
+assert(ischar(runNum),'Invalid analysis parameter file: runNum must be a string'); 
+
+% dateSubject
+assert(logical(exist('dateSubject','var')),'Invalid analysis parameter file: must specify a date and subject as dateSubject.');
+assert(ischar(dateSubject),'Invalid analysis parameter file: dateSubject must be a string'); 
+
+% ephysVolume
+assert(logical(exist('ephysVolume','var')),'Invalid analysis parameter file: must specify the ephys data folder as ephysVolume.');
+assert(ischar(ephysVolume),'Invalid analysis parameter file: ephysVolume must be a string'); 
+
+% stimulusLogVolume
+assert(logical(exist('stimulusLogVolume','var')),'Invalid analysis parameter file: must specify the stimulus log folder as stimulusLogVolume.');
+assert(ischar(stimulusLogVolume),'Invalid analysis parameter file: stimulusLogVolume must be a string');
+
+% outputVolume
+assert(logical(exist('runNum','var')),'Invalid analysis parameter file: must specify a run number as runNum.');
+assert(ischar(runNum),'Invalid analysis parameter file: runNum must be a string'); 
+
+% stimParamsFilename
+assert(logical(exist('outputVolume','var')),'Invalid analysis parameter file: must specify an output folder as outputVolume.');
+assert(ischar(outputVolume),'Invalid analysis parameter file: outputVolume must be a string');
+
+% stimParamsFilename
+assert(logical(exist('outputVolume','var')),'Invalid analysis parameter file: must specify an output folder as outputVolume.');
+assert(ischar(outputVolume),'Invalid analysis parameter file: outputVolume must be a string');
+
+%%% SAVE SWITCHES %%%
+% saveFig
+if ~logical(exist('saveFig','var')) || ~ismember(saveFig,[0 1])
+  saveFig = 1; %#ok
+  logString = strcat(logString,'saveFig\n');
+end
+
+% closeFig
+if ~logical(exist('closeFig','var')) || ~ismember(closeFig,[0 1])
+  closeFig = 0; %#ok
+  logString = strcat(logString,'closeFig\n');
+end
+
+% exportFig
+if ~logical(exist('exportFig','var')) || ~ismember(exportFig,[0 1])
+  exportFig = 0; %#ok
+  logString = strcat(logString,'exportFig\n');
+end
+
+% saveFigData
+if ~logical(exist('saveFigData','var')) || ~ismember(saveFigData,[0 1])
+  saveFigData = 0; %#ok
+  logString = strcat(logString,'saveFigData\n');
+end
+
+% savePreprocessed
+if ~logical(exist('savePreprocessed','var')) || ~ismember(savePreprocessed,[0 1])
+  savePreprocessed = 0; %#ok
+  logString = strcat(logString,'savePreprocessed\n');
+end
+
+%%% VERBOSITY %%%
+if ~logical(exist('verbosity','var')) || ~ischar(verbostity) || ~ismember(verbosity,{'INFO','DEBUG','VERBOSE'})
+  verbosity = 'INFO'; %#ok
+  logString = strcat(logString,'verbosity\n');
+end
+
+%%% EPHYS PARAMS %%%
+if ~logical(exist('ephysParams','var')) || ~isstruct(ephysParams)
+  ephysParams.needLFP = 0;
+  ephysParams.needSpikes = 0;
+  logString = strcat(logString,'ephysParams\n');
+else
+  if ~isfield(ephysParams,'needLFP') || ~ismember(ephysParams.needLFP,[0 1])
+    ephysParams.needLFP = 0;
+    logString = strcat(logString,'ephysParams.needLFP\n');
+  end
+  if ~isfield(ephysParams,'needSpikes') || ~ismember(ephysParams.needSpikes,[0 1])
+    ephysParams.needSpikes = 0;
+    logString = strcat(logString,'ephysParams.needSpikes\n');
+  end
+  if ~isfield(ephysParams,'spikeChannels') || ~isnumeric(ephysParams.spikeChannels)
+    ephysParams.spikeChannels = [];
+    logString = strcat(logString,'ephysParams.spikeChannels\n');
+  end
+  if ~isfield(ephysParams,'lfpChannels') || ~isnumeric(ephysParams.lfpChannels)
+    ephysParams.lfpChannels = [];
+    logString = strcat(logString,'ephysParams.lfpChannels\n');
+  end
+  if ephysParams.needLFP && ephysParams.needSpikes
+    assert(length(ephysParams.spikeChannels) == length(ephysParams.lfpChannels),...
+      'Invalid analysis parameter file: spikeChannels and lfpChannels must be the same length if analyzing both');
+    assert(all(ephysParams.spikeChannels == ephysParams.lfpChannels),...
+      'Invalid analysis parameter file: spikeChannels and lfpChannels must be in the same order, if analyzing both');
+  end
+  % channelNames
+  if (ephysParams.needLFP || ephysParams.needSpikes)
+    numChannels = max(ephysParams.needLFP*length(ephysParams.lfpChannels),ephysParams.needSpikes*length(ephysParams.spikeChannels));
+    if isfield(ephysParams,'channelNames')
+      assert(length(channelNames) == numChannels,'Invalid analysis parameter file: channelNames length must match length of channels to analyze, or be zero.')
+    else
+      ephysParams.channelNames = cell(numChannels,1);
+      for channel_i = 1:numChannels
+        ephysParams.channelNames{channel_i} = sprintf('ch%d',channel_i);  
+      end
+      logString = strcat(logString,'ephysParams.channelNames\n');
+    end
+  end
+  % lfpChannelScaleBy
+  if ~isfield(ephysParams,'lfpChannelScaleBy') || ~isnumeric(ephysParams.lfpChannelScaleBy)
+    ephysParams.lfpChannelScaleBy = ones(size(ephysParams.lfpChannels));
+    logString = strcat(logString,'ephysParams.lfpChannelScaleBy (set to one)\n');
+  end
+  if length(ephysParams.lfpChannelScaleBy) == 1 && length(ephysParams.lfpChannels) > 1
+    ephysParams.lfpChannelScaleBy = ephysParams.lfpChannelScaleBy*ones(size(ephysParams.lfpChannels));
+    logString = strcat(logString,'ephysParams.lfpChannelScaleBy (applied single given value to all channels)\n');
+  end
+  assert(isfield(ephysParams,'cPtCal') && isnumeric(ephysParams.),'Invalid analysis parameter file: must specify conversion from spike time units to decimated LFP indices');
+  if ~isfield(ephysParams,'decimateFactorPass1') || ~isnumeric(ephysParams.decimateFactorPass1)
+    ephysParams.decimateFactorPass1 = 6;
+    logString = strcat(logString,'ephysParams.decimateFactorPass1\n');
+  end
+  if ~isfield(ephysParams,'decimateFactorPass2') || ~isnumeric(ephysParams.decimateFactorPass2)
+    ephysParams.decimateFactorPass2 = 5;
+    logString = strcat(logString,'ephysParams.decimateFactorPass2\n');
+  end
+  if ~isfield(ephysParams,'samPerMS') || ~isnumeric(ephysParams.samPerMS)
+    ephysParams.samPerMS = 1;
+    logString = strcat(logString,'ephysParams.samPerMS\n'); 
+  else
+    if ephysParams.samPerMS ~= 1
+      disp('Warning: lfp samples per ms ~= 1. In current implementation, this will lead to axes in samples incorrectly labeled as in ms');
+    end
+  end
+  if ~isfield(ephysParams,'unitsToUnsort') || ~iscell(ephysParams.unitsToUnsort) || ~length(ephysParams.unitsToUnsort) == length(ephysParams.spikeChannels)
+    ephysParams.unitsToUnsort = cell(length(ephysParams.spikeChannels));
+    logString = strcat(logString,'ephysParams.unitsToUnsort\n'); 
+  end
+  if ~isfield(ephysParams,'unitsToDiscard') || ~iscell(ephysParams.unitsToDiscard) || ~length(ephysParams.unitsToDiscard) == length(ephysParams.spikeChannels)
+    ephysParams.unitsToDiscard = cell(length(ephysParams.spikeChannels));
+    logString = strcat(logString,'ephysParams.unitsToDiscard\n'); 
+  end
+  if ~isfield(ephysParams,'spikeWaveformPca') || ~ismember(ephysParams.spikeWaveformPca,[0 1])
+    ephysParams.spikeWaveformPca = 0;
+    logString = strcat(logString,'ephysParams.spikeWaveformPca\n');
+  end
+  if ~isfield(ephysParams,'plotSpikeWaveforms') || ~ismember(ephysParams.plotSpikeWaveforms,[0 1])
+    ephysParams.needSpikes = 0;
+    logString = strcat(logString,'ephysParams.plotSpikeWaveforms\n');
+  end
+  if ~isfield(ephysParams,'shiftSpikeWaveforms') || ~ismember(ephysParams.shiftSpikeWaveforms,[0 1])
+    ephysParams.shiftSpikeWaveforms = 0;
+    logString = strcat(logString,'ephysParams.shiftSpikeWaveforms\n');
+  end
+  if ~isfield(ephysParams,'filter') %todo: add check for digitalFilter object type?
+    ephysParams.filter = '';
+    logString = strcat(logString,'ephysParams.filter\n');
+  end
+  if ~isfield(ephysParams,'plotFilterResult') || ~ismember(ephysParams.plotFilterResult,[0 1])
+    ephysParams.shiftSpikeWaveforms = 0;
+    logString = strcat(logString,'ephysParams.plotFilterResult\n');
+  end
+end
+
+%%% AnalogInParams %%%
+analogInParams.needAnalogIn = 0;
+analogInParams.analogInChannels = [138,139,140]; 
+analogInParams.channelNames = {'eyeX','eyeY','eyeD'};
+analogInParams.analogInChannelScaleBy = [5/32764 5/32764 5/32764]; %converts raw values to volts
+analogInParams.decimateFactorPass1 = 1; 
+analogInParams.decimateFactorPass2 = 1;
+analogInParams.samPerMS = 1; %THIS IS AFTER DECIMATION, and applies to analogIn (should be raw rate/productOfDecimateFactors)
+% see http://www.mathworks.com/help/signal/examples/filter-design-gallery.html
+butter200Hz_v1 = designfilt('lowpassiir', 'PassbandFrequency', 120, 'StopbandFrequency', 480, 'PassbandRipple', 1,...
+  'StopbandAttenuation', 60, 'SampleRate', 1000, 'MatchExactly', 'passband');  %this returns a 3rd order iir filter
+analogInParams.filters = {0,0,0};%{butter200Hz_v1;butter200Hz_v1;butter200Hz_v1}; %filter channel i if filters{i} is digital filter or 1x2 numeric array
+analogInParams.plotFilteredSignal = 1; %#ok
+
+photodiodeParams.needPhotodiode = 0;
+photodiodeParams.channels = [1;2]; %#ok
+
+% parameters preprocessLogFile, see function for details
+stimSyncParams.usePhotodiode = 0;        %#ok
+%
+eyeCalParams.needEyeCal = 0;
+eyeCalParams.method = 'zeroEachFixation';
+eyeCalParams.makePlots = 0;
+eyeCalParams.eyeXChannelInd = 1;
+eyeCalParams.eyeYChannelInd = 2;
+eyeCalParams.eyeDChannelInd = 3;
+eyeCalParams.gainX = 112;
+eyeCalParams.gainY = 107;
+eyeCalParams.flipX = 1;
+eyeCalParams.flipY = 1; 
+eyeCalParams.offsetX = -6.4;
+eyeCalParams.offsetY = -5.6; 
+eyeCalParams.minFixZeroTime = 1000; %#ok
+
+accelParams.needAccelCal = 0;
+accelParams.accelChannels = {[4;5;6]};
+accelParams.channelGains = {[1/.666 1/.666 1/.666]};
+accelParams.calMethods = {'hardcode'}; %other option is 'calFile'; calibration method
+accelParams.calFiles = {''}; %if method is 'calFile', an ns2 filename
+
+% parameters for excludeStimuli, see function for details
+excludeStimParams.fixPre = 100; %ms
+excludeStimParams.fixPost = 100; %ms
+excludeStimParams.flashPre = 0;  %ms
+excludeStimParams.flashPost = 0; %ms
+excludeStimParams.juicePre = 0; % optional, ms
+excludeStimParams.juicePost = 0; % optional, ms
+excludeStimParams.DEBUG = 0; % makes exclusion criterion plots if true
+% additional optional excludeStimParams: accel1, accel2, minStimDur (ms)
+
+psthParams.psthPre = 100; % if e.g. +200, then start psth 200ms before trial onset; 
+psthParams.psthImDur = 0;  % only need to set this for variable length stim runs; else, comes from log file
+psthParams.psthPost = 300;
+psthParams.smoothingWidth = 10;  %psth smoothing width, in ms
+psthParams.errorType = 1; %chronux convention: 1 is poisson, 2 is trialwise bootstrap, 3 is across trial std for binned spikes, bootstrap for spike times 
+psthParams.errorRangeZ = 1; %how many standard errors to show
+psthParams.bootstrapSamples = 100;
+psthParams.psthColormapFilename = 'cocode2.mat'; % a file with one variable, a colormap called 'map'
+
+
+% TW=3 with T=.2, then W = 15 Hz (5 tapers)
+% TW=1.5 with T=.1, then W = 15 Hz (2 tapers)
+% TW = 1.5 with T=.2, then W = 7.5 Hz (2 tapers)
+chronuxParams.tapers = [3 5]; %[3 5] is chronux default; 
+chronuxParams.pad = 1;
+chronuxParams.fs = 1;
+chronuxParams.trialave = 1;
+chronuxParams.err = [1 .05];  %note: first entry will be automatically switched to 2 if calcSwitch.useJacknife == 1
+chronuxParams.fpass = [0 .1]; 
+tfParams.movingWin = [200 5]; 
+tfParams.specgramRowAve = 0;
+
+correlParams.maxShift = []; % a number, or empty
+correlParams.matchTimeRanges = 1;
+correlParams.timeDifferenceBound = [0,200];
+correlParams.normalize = 1;
+correlParams.useJacknife = 0;
+correlParams.jacknifeDraws = 100;
+switch machine
+  case 'laptop'
+    correlParams.jacknifeParallelWorkers = 0;    
+  case 'hopper'
+    correlParams.jacknifeParallelWorkers = 0;   
+  case 'turing'
+    correlParams.jacknifeParallelWorkers = 20;    
+end
+spikeCorrelSmoothingWidth = 5; %ms
+filterPoints = -20*spikeCorrelSmoothingWidth:20*spikeCorrelSmoothingWidth;
+smoothingFilter = exp(-1*filterPoints.^2/(2*spikeCorrelSmoothingWidth^2));
+correlParams.smoothingFilter = smoothingFilter/sum(smoothingFilter); %#ok
+%
+lfpAlignParams.samPerMS = 1; % because this is after decimation
+lfpAlignParams.msPreAlign = psthParams.psthPre+tfParams.movingWin(1)/2; 
+lfpAlignParams.msPostAlign = psthParams.psthImDur+psthParams.psthPost+tfParams.movingWin(1)/2;
+%
+spikeAlignParams.preAlign = psthParams.psthPre+3*psthParams.smoothingWidth;
+spikeAlignParams.postAlign = psthParams.psthImDur+psthParams.psthPost+3*psthParams.smoothingWidth;   %#ok
+% for lfps, constrain first and (optional) last [n m] samples to 0 mean
+useDCSUB = 0;
+if useDCSUB
+  %lfpAlignParams.DCSUB_SAM = [lfpAlignParams.msPreAlign, lfpAlignParams.msPreAlign+10; 0, 0 ]; % 0-th order 
+  lfpAlignParams.DCSUB_SAM = [lfpAlignParams.msPreAlign, lfpAlignParams.msPreAlign+10;lfpAlignParams.msPreAlign, lfpAlignParams.msPreAlign+10 ]; % 1st order 
+else
+  lfpAlignParams.DCSUB_SAM = 0;   %#ok
+end
+% firing rate calculation epochs. Can provide either time (ms from stim onset),
+% or function handle, which will receive the minimum stimulus duration in the run as an input
+frEpochsCell = {{60, @(stimDur) stimDur+60};...
+                {60, 260}; ...
+                {260, @(stimDur) stimDur+60}}; %#ok
+
+
+% Boolean variables to specify which computations to perform; TODO: read
+% from config file, eventually with conditional on log file info
+calcCoherenceRFcpt = 0;  %#ok
+calcCoherenceRFcc = 0;   %#ok
+calcCoherenceRFptpt = 0; %#ok
+calcGrangerRF = 0;       %#ok 
+
+
+
+%%%% note: all analysisGroups cell arrays are nx1, NOT 1xn
+analysisGroups.selectivityIndex.groups = {{'face';'nonface'},{'face';'object'},{'face';'body'}};
+%
+analysisGroups.stimPrefBarPlot.groups = {{{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'};{'face';'object';'body'}}};
+analysisGroups.stimPrefBarPlot.colors  = {{{'b';'c';'y';'g';'m';'r';'k'};{'b';'g';'r'}}};
+analysisGroups.stimPrefBarPlot.names = {'fobPlus'};
+%
+analysisGroups.stimulusLabelGroups.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
+analysisGroups.stimulusLabelGroups.names = {'fobPlus'};
+analysisGroups.stimulusLabelGroups.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
+%
+analysisGroups.evokedPotentials.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
+analysisGroups.evokedPotentials.names = {'fobPlus'};
+analysisGroups.evokedPotentials.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
+%
+analysisGroups.analogInPotentials.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
+analysisGroups.analogInPotentials.channels = {[1; 2]};
+analysisGroups.analogInPotentials.names = {'eyePositions,fobPlus'};
+analysisGroups.analogInPotentials.units = {'degrees visual angle'};
+analysisGroups.analogInPotentials.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
+%
+analysisGroups.analogInDerivatives.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
+analysisGroups.analogInDerivatives.channels = {[1; 2]};
+analysisGroups.analogInDerivatives.names = {'eyeVelocity,fobPlus'};
+analysisGroups.analogInDerivatives.units = {'degrees visual angle/sec'};
+analysisGroups.analogInDerivatives.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
+%
+analysisGroups.colorPsthEvoked.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
+analysisGroups.colorPsthEvoked.names = {'fobPlus'};
+analysisGroups.colorPsthEvoked.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
+%
+analysisGroups.linePsthEvoked.groups = {{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'techno'}};
+analysisGroups.linePsthEvoked.names = {'fobPlus'};
+analysisGroups.linePsthEvoked.colors = {{'b';'c';'y';'g';'m';'r';'k'}};
+%
+analysisGroups.evokedPsthOnePane.groups = {{'face';'nonface'}};
+analysisGroups.evokedPsthOnePane.names = {'faceVnon'};
+%
+analysisGroups.tuningCurves.groups = {{'humanFaceL90','humanFaceL45','humanFaceFront','humanFaceR45','humanFaceR90'},...
+  {'monkeyFaceL90','monkeyFaceL45','monkeyFaceFront','monkeyFaceR45','monkeyFaceR90'}}; %can be images or categories
+analysisGroups.tuningCurves.paramValues = {[-90 -45 0 45 90], [-90 -45 0 45 90]};
+analysisGroups.tuningCurves.paramLabels = {'viewing angle (degrees)','viewing angle (degrees)'};
+analysisGroups.tuningCurves.names = {'Human face view','Monkey face view'};
+%
+analysisGroups.spectraByCategory.groups = {{'face';'nonface'}};  %todo: add spectra diff?
+analysisGroups.spectraByCategory.names = {'faceVnon'};
+analysisGroups.spectraByCategory.colors = {{'r';'b'}};
+%
+analysisGroups.tfSpectraByCategory.groups = {{'face'};{'nonface'}};%{'object'};{'body'}      %todo: add tf spectra diff?
+analysisGroups.tfSpectraByCategory.names = {'face','nonface'};%'nonface';'object';'body'
+%
+analysisGroups.lfpSingleTrialsByCategory.groups = {{'face';'nonface'}};
+analysisGroups.lfpSingleTrialsByCategory.names = {'faceVnon'};
+%
+analysisGroups.coherenceByCategory.groups = {{'face';'nonface'}}; %{'face';'object';'body'};{'humanFace';'monkeyFace';'place';'fruit';'humanBody';'monkeyBody';'hand';'techno'}
+analysisGroups.coherenceByCategory.colors = {{'r';'b'}}; %{'r';'g';'b'};{'b';'c';'y';'g';'m';'r';'k';'k'}
+analysisGroups.coherenceByCategory.names = {'faceVnon'}; %'fob';'slimCats'
+%
+analysisGroups.tfCouplingByCategory.groups = {{'face'};{'nonface'};{'object'};{'body'}};
+
+analysisGroups.byImage = {};      %#ok
+analysisGroupColors.byImage = {}; %#ok
+%%%%%
+
+calcSwitch.categoryPSTH = 1;
+calcSwitch.imagePSTH = 1;
+calcSwitch.faceSelectIndex = 1;
+calcSwitch.faceSelectIndexEarly = 1;
+calcSwitch.faceSelectIndexLate = 1;
+calcSwitch.inducedTrialMagnitudeCorrection = 0;
+calcSwitch.evokedSpectra = 0;
+calcSwitch.inducedSpectra = 0;
+calcSwitch.evokedImageTF = 0;
+calcSwitch.inducedImageTF = 0;
+calcSwitch.evokedCatTF = 0;
+calcSwitch.inducedCatTF = 0;
+calcSwitch.meanEvokedTF = 0;
+calcSwitch.trialMeanSpectra = 0;
+calcSwitch.coherenceByCategory = 0;
+calcSwitch.spikeTimes = 0;
+calcSwitch.useJacknife = 0;      
+
+if calcSwitch.useJacknife
+  chronuxParams.err(1) = 2; %#ok
+end
+
+%%% set paths and directories, EDIT RARELY %%%
+analogInFilename = sprintf('%s/%s/%s%s.ns2',ephysVolume,dateSubject,dateSubject,runNum);   %#ok
+lfpFilename = sprintf('%s/%s/%s%s.ns5',ephysVolume,dateSubject,dateSubject,runNum);        %#ok
+spikeFilename = sprintf('%s/%s/%s%s.nev',ephysVolume,dateSubject,dateSubject,runNum); %note that this file also contains blackrock digital in events
+taskFilename = sprintf('%s/%s/%s0%s.log',stimulusLogVolume,dateSubject,dateSubject,runNum); %information on stimuli and performance
+outDir = sprintf('%s/%s/%s/%s/',outputVolume,dateSubject,analysisLabel,runNum);
+analysisParamsFilename = strcat(outDir,analysisParamsFilenameStem);
+preprocessedDataFilename = strcat(outDir,preprocessedDataFilenameStem);                     %#ok
+%
+load('cocode2.mat');
+psthColormap = map;  %#ok
+%
+if ~exist(outDir,'dir')
+  mkdir(outDir);
+end
+save(analysisParamsFilename);
+
+end
+
diff --git a/dependencies/cart2polDeg.m b/dependencies/cart2polDeg.m
new file mode 100644
index 0000000..0a560e2
--- /dev/null
+++ b/dependencies/cart2polDeg.m
@@ -0,0 +1,7 @@
+function [ theta, r ] = cart2polDeg( x,y )
+%UNTITLED2 Summary of this function goes here
+%   Detailed explanation goes here
+[theta,r] = cart2pol(x,y);
+theta = rad2deg(theta);
+end
+
diff --git a/preprocessPhotodiodeStrobe.m b/preprocessPhotodiodeStrobe.m
index 907d6a4..02545cc 100644
--- a/preprocessPhotodiodeStrobe.m
+++ b/preprocessPhotodiodeStrobe.m
@@ -1,10 +1,385 @@
 function [ analogInData ] = preprocessPhotodiodeStrobe( analogInData, params )
-%preprocessPhotodiodeStrobe cleans a photodiode signal
-%   Detailed explanation goes here
+%preprocessPhotodiodeStrobe extracts frame times from a photodiode signal 
+%   Inputs:
+%   - analogInData: data array returned by preprocessAnalogIn (channels x samples)
+%   - params: struct with the fields
+%     - needPhotodiode: return immediately if 0 (type: logical)
+%     - analogInChannelInd: this is the index of the photodiode channel in the list of analog in
+%                           channels specified in the paramFile (type: int)
+%     - centerCornerOffset: this is the offset, in ms, of the photodiode
+%                           trigger from the screen-center frame time
+%                           (type: float)
+%     - levelCalibType: (type: string) options are:
+%       - hardcode
+%       -
+%       -
+%
+%
+%     - numLevels: (type: int)
+%     - hardcodeFromFile: (type: logical)
+%     - saveCalibFile: (type: logical)
+%     - inputCalibrationFile (required only if hardcodeFromFile): (type: string)
+%     - outputCalibrationFile (required only if saveCalibFile): (type:string)
 if ~params.needPhotodiode
   return
 end
-error('photodiode strobe preprocessing not yet implemented');
 
+diodeTrace = analogInData(params.analogInChannelInd,:); 
+diodeTrace = -1*(diodeTrace - mean(diodeTrace)); %this makes frames peaks
+centerCornerOffset = params.centerCornerOffset; 
+
+frameNumCeiling = ceil(length(diodeTrace)*params.framerate/params.sampleRate);
+% Preprocessing logic:
+% 1) Find peaks in the inverted diode trace. A peak occurs when the cathode ray passes the diode,
+%    even on black pixels. However, additional small-prominence peaks occur due to noise
+% 2) Sort peaks by their height: since both the the signal and its
+%    second derivative are high at true peaks, they should be the largest
+rawPeaks = findpeaks(diodeTrace);
+rawPeaksSorted = sort(diodeTrace(rawPeaks.loc),'descend');
+rawPeaksSorted = rawPeaksSorted(1:frameNumCeiling);
+
+% additional todo: add frame count and alternation checks for manual levels
+
+if strcmp(params.levelCalibType,'hardcode')  
+  if params.hardcodeFromFile
+    switch params.numLevels
+      case 1
+        load(params.inputCalibrationFile,'highThreshold');
+      case 2
+        load(params.inputCalibrationFile,'highThreshold','lowThreshold');
+      case 3
+        load(params.inputCalibrationFile,'highThreshold','midThreshold','lowThreshold');
+    end
+  else
+    switch params.numLevels
+      case 1
+        highThreshold = params.levelHigh;
+      case 2
+        highThreshold = params.levelHigh;
+        lowThreshold = params.levelLow;
+      case 3
+        highThreshold = params.levelHigh;
+        midThreshold = params.levelMid;
+        lowThreshold = params.levelLow;
+    end
+  end
+end
+              
+if any(strcmp(params.levelCalibType, {'auto','autoAndPlot','autoAndCheck'}))
+  % Auto-calibration logic: 
+  %  The largest jumps in the plot of sorted peak heights should occur between true peak types,  
+  %  or between the smallest peak type and the largest noise peak. So find and sort those jumps, 
+  %  then place level cut-offs at the mid-point of the largest ones. 
+
+  peakDiffs = diff(rawPeaksSorted);
+  [peakDiffsSorted,peakDiffsSortInds] = sort(peakDiffs);
+  if ~(strcmp(params.levelCalibType,'hardcodeAndPlot') || strcmp(params.levelCalibType,'hardcodeAndCheck'))
+    switch params.numLevels
+      case 1
+        highThresholdDiffInd = peakDifSortInds(1);
+        highThreshold = (rawPeaksSorted(highThresholdDiffInd + 1) - rawPeaksSorted(highThresholdDiffInd))/2;
+      case 2
+        highThresholdDiffInd = min(peakDiffSortInds(1:2)); %note: these end up on the high side of the threshold,i.e. the low index side
+        lowThresholdDiffInd = max(peakDiffSortInds(1:2));
+        highThreshold = (rawPeaksSorted(highThresholdDiffInd + 1) - rawPeaksSorted(highThresholdDiffInd))/2;
+        lowThreshold = (rawPeaksSorted(lowThresholdDiffInd + 1) - rawPeaksSorted(lowThresholdDiffInd))/2;
+      case 3
+        highThresholdDiffInd = min(peakDiffSortInds(1:3)); %note: these end up on the high side of the threshold,i.e. the low index side
+        lowThresholdDiffInd = max(peakDiffSortInds(1:3));
+        midThresholdDiffInd = setdiff(peakDiffSortInds(1:3),[highThresholdDiffInd, lowThresholdDiffInd]);
+        highThreshold = (rawPeaksSorted(highThresholdDiffInd + 1) - rawPeaksSorted(highThresholdDiffInd))/2;
+        midThreshold = (rawPeaksSorted(midThresholdDiffInd + 1) - rawPeaksSorted(midThresholdDiffInd))/2;
+        lowThreshold = (rawPeaksSorted(lowThresholdDiffInd + 1) - rawPeaksSorted(lowThresholdDiffInd))/2;
+    end
+  end
+end
+
+% if two levels, check that frames alternate appropriately
+if params.numLevels == 2
+  frameSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > lowThreshold);
+  highPeakSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > highThreshold);
+  lowPeakSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > lowThreshold && diodeTrace(rawPeaks.loc) < highThreshold);
+  frameCountDiff = length(highPeakSampleInds) - length(lowPeakSampleInds);
+  % with correct alternation, smallest interval between subsequent high peaks must tbe greater than longest frame; same for low peaks 
+  frameAlternationError = (min(min(diff(highPeakSampleInds)),min(diff(lowPeakSampleInds))) > max(diff(frameSampleInds)));
+end
+    
+%make plots, if requested, or if auto calib failed
+if any(strcmp(params.levelCalibType,{'autoAndPlot','autoAndCheck','hardcodeAndPlot','hardcodeAndCheck'})) || (params.numLevels == 2 && (abs(frameCountDiff) > 1 || frameAlternationError))
+  if params.numLevels == 2 && abs(frameCountDiff) > 1
+    fprintf('Frame count mismatch: Found %d bright-square frames and %d dark-square frames. Starting manual check.\n',...
+      length(highPeakSampleInds),length(lowPeakSampleInds));
+  end
+  if params.numLevels == 2 && frameAlternationError
+    fprintf('Found light-dark frame alternation error. Starting manual check.\n');
+  end
+  if any(strcmp(params.levelCalibType, {'hardcodeAndPlot','hardcodeAndCheck','hardcode'}))
+    calibSuffix = 'hard code';
+  else
+    calibSuffix = 'auto';
+  end
+
+  figure();
+  subplot(1,2,1);
+  framesToPlot = 10;  %magic number; should move to top
+  sampleRange = ceil(frameSampleInds(1) - 5*params.sampleRate/params.framerate):ceil(frameSampleInds(1) + framesToPlot*params.sampleRate/params.framerate)
+  plot(sampleRange,diodeTrace(sampleRange));
+  hold on
+  h = get(gca,'xlim');
+  plot([h(1) h(2)], [highThreshold highThreshold]);
+  forLegend121 = {'raw trace',sprintf('high threshold %s',calibSuffix)}; %note: 121 refers to the subplot
+  if params.numLevels == 3
+    plot([h(1) h(2)], [midThreshold midThreshold]);
+    forLegend121 = horzcat(forLegend121, {sprintf('mid threshold %s',calibSuffix)});
+  end
+  if params.numLevels > 1
+    plot([h(1) h(2)], [lowThreshold lowThreshold]);
+    forLegend121 = horzcat(forLegend121, {sprintf('low threshold %s',calibSuffix)});
+  end
+  plot(frameSampleInds(1:framesToPlot), diodeTrace(frameSampleInds(1:framesToPlot)));
+  forLegend121 = horzcat(forLegend121, {'frameTimes'});
+  legend(forLegend121);  
+  subplot(1,2,2);
+  plot(rawPeaksSorted);
+  hold on
+  h = get(gca,'xlim');
+  plot([h(1) h(2)], [highThreshold highThreshold]);
+  forLegend122 = {'peak levels',sprintf('high threshold %s',calibSuffix)}; %note: 122 refers to the subplot
+  if params.numLevels == 3
+    plot([h(1) h(2)], [midThreshold midThreshold]);
+    forLegend122 = horzcat(forLegend122, {sprintf('mid threshold %s',calibSuffix)});
+  end
+  if params.numLevels > 1
+    plot([h(1) h(2)], [lowThreshold lowThreshold]);
+    forLegend122 = horzcat(forLegend122, {sprintf('low threshold %s',calibSuffix)});
+  end
+  legend(forLegend122);
+  drawnow();
+
+  if any(strcmp(params.levelCalibType,{'autoAndCheck','hardcodeAndCheck'})) || params.numLevels == 2
+    tmp = input('Accept high threshold? Enter y or n: ','s');
+    if strcmp(tmp,'n')
+      tmp = input('Enter new high threshold, or ''quit'' to continue without photodiode sync:','s');
+      while isempty(str2num(tmp)) && ~strcmp(tmp,'quit')
+        tmp = input('Invalid input: enter new high threshold, or ''quit'' to continue without photodiode sync:','s');
+      end
+      if strcmp(tmp,'quit')
+        return %todo: need to handle this better
+      else
+        highThreshold = str2num(tmp);
+        subplot(1,2,1);
+        h = get(gca,'xlim');
+        plot([h(1) h(2)], [highThreshold highThreshold]);
+        forLegend121 = horzcat(forLegend121,{'high threshold manual'});
+        legend(forLegend121);
+        subplot(1,2,2);
+        h = get(gca,'xlim');
+        plot([h(1) h(2)], [highThreshold highThreshold]);
+        forLegend122 = horzcat(forLegend121,{'high threshold manual'});
+        legend(forLegend122);
+      end
+    end
+    if params.numLevels == 3
+      tmp = input('Accept mid threshold? Enter y or n: ','s');
+      if strcmp(tmp,'n')
+        tmp = input('Enter new mid threshold, or ''quit'' to continue without photodiode sync:','s');
+        while isempty(str2num(tmp)) && ~strcmp(tmp,'quit')
+          tmp = input('Invalid input: enter new mid threshold, or ''quit'' to continue without photodiode sync:','s');
+        end
+        if strcmp(tmp,'quit')
+          return %todo: need to handle this better
+        else
+          midThreshold = str2num(tmp);
+          subplot(1,2,1);
+          h = get(gca,'xlim');
+          plot([h(1) h(2)], [midThreshold midThreshold]);
+          forLegend121 = horzcat(forLegend121,{'mid threshold manual'});
+          legend(forLegend121);
+          subplot(1,2,2);
+          h = get(gca,'xlim');
+          plot([h(1) h(2)], [midThreshold midThreshold]);
+          forLegend122 = horzcat(forLegend121,{'mid threshold manual'});
+          legend(forLegend122);
+        end
+      end
+    end
+    if params.numLevels > 1
+      tmp = input('Accept low threshold? Enter y or n: ','s');
+      if strcmp(tmp,'n')
+        tmp = input('Enter new low threshold, or ''quit'' to continue without photodiode sync:','s');
+        while isempty(str2num(tmp)) && ~strcmp(tmp,'quit')
+          tmp = input('Invalid input: enter new low threshold, or ''quit'' to continue without photodiode sync:','s');
+        end
+        if strcmp(tmp,'quit')
+          return %todo: need to handle this better
+        else
+          lowThreshold = str2num(tmp);
+          subplot(1,2,1);
+          h = get(gca,'xlim');
+          plot([h(1) h(2)], [lowThreshold lowThreshold]);
+          forLegend121 = horzcat(forLegend121,{'low threshold manual'});
+          legend(forLegend121);
+          subplot(1,2,2);
+          h = get(gca,'xlim');
+          plot([h(1) h(2)], [lowThreshold lowThreshold]);
+          forLegend122 = horzcat(forLegend121,{'low threshold manual'});
+          legend(forLegend122);
+        end
+      end
+    end
+  end
+end
+
+% if two levels, check that frames alternate appropriately
+if params.numLevels == 2
+  frameSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > lowThreshold);
+  highPeakSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > highThreshold);
+  lowPeakSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > lowThreshold && diodeTrace(rawPeaks.loc) < highThreshold);
+  frameCountDiff = length(highPeakSampleInds) - length(lowPeakSampleInds);
+  % with correct alternation, smallest interval between subsequent high peaks must tbe greater than longest frame; same for low peaks 
+  frameAlternationError = (min(min(diff(highPeakSampleInds)),min(diff(lowPeakSampleInds))) > max(diff(frameSampleInds)));
+  notValid = (abs(frameCountDiff) > 1 || frameAlternationError);
+  if notValid
+    if abs(frameCountDiff) > 1
+      tmp = input(sprintf('Frame count mismatch: Found %d bright-square frames and %d dark-square frames. Press ENTER to start manual check, or q to quit.',...
+        length(highPeakSampleInds),length(lowPeakSampleInds)));
+    end
+    if frameAlternationError
+      tmp = input(sprintf('Found light-dark frame alternation error. Press ENTER to start manual check, or q to quit.'));
+    end
+    while ~strcmp('tmp',{'','q'})
+      tmp = imput('Invalid input. Press ENTER to re-start manual check, or q to quit.');
+    end
+    if isempty(tmp)
+      notValid = 1;
+    else
+      notValid = 0;
+    end
+  end
+else
+  notValid = 0;
+end
+
+if strcmp(params.levelCalibType,'manual') || notValid
+  notValid = 1;
+  while notValid
+    figure();
+    subplot(1,2,1);
+    framesToPlot = 10;  %magic number; should move to top
+    sampleRange = ceil(length(diodeTrace)/2 - (framesToPlot/2)*params.sampleRate/params.framerate):ceil(length(diodeTrace)/2 - (framesToPlot/2)*params.sampleRate/params.framerate)
+    plot(sampleRange,diodeTrace(sampleRange));
+    forLegend121 = {'raw trace'};
+    legend(forLegend121);
+    subplot(1,2,2);
+    plot(rawPeaksSorted);
+    forLegend122 = {'peak levels'};
+    legend(forLegend122);
+    drawnow();
+    tmp = input('Enter the high threshold, or ''quit'' to continue without photodiode sync:','s');
+    while isempty(str2num(tmp)) && ~strcmp(tmp,'quit')
+      tmp = input('Invalid input: enter the high threshold, or ''quit'' to continue without photodiode sync:','s');
+    end
+    if strcmp(tmp,'quit')
+      return %todo: need to handle this better
+    else
+      highThreshold = str2num(tmp);
+      subplot(1,2,1);
+      h = get(gca,'xlim');
+      plot([h(1) h(2)], [highThreshold highThreshold]);
+      forLegend121 = horzcat(forLegend121,{'high threshold manual'});
+      legend(forLegend121);
+      subplot(1,2,2);
+      h = get(gca,'xlim');
+      plot([h(1) h(2)], [highThreshold highThreshold]);
+      forLegend122 = horzcat(forLegend121,{'high threshold manual'});
+      legend(forLegend122);
+    end
+    if params.numLevels == 3
+      tmp = input('Enter the mid threshold, or ''quit'' to continue without photodiode sync:','s');
+      while isempty(str2num(tmp)) && ~strcmp(tmp,'quit')
+        tmp = input('Invalid input: enter the mid threshold, or ''quit'' to continue without photodiode sync:','s');
+      end
+      if strcmp(tmp,'quit')
+        return %todo: need to handle this better
+      else
+        midThreshold = str2num(tmp);
+        subplot(1,2,1);
+        h = get(gca,'xlim');
+        plot([h(1) h(2)], [midThreshold midThreshold]);
+        forLegend121 = horzcat(forLegend121,{'mid threshold manual'});
+        legend(forLegend121);
+        subplot(1,2,2);
+        h = get(gca,'xlim');
+        plot([h(1) h(2)], [midThreshold midThreshold]);
+        forLegend122 = horzcat(forLegend121,{'mid threshold manual'});
+        legend(forLegend122);
+      end
+    end
+    if params.numLevels > 1
+      tmp = input('Enter the low threshold, or ''quit'' to continue without photodiode sync:','s');
+      while isempty(str2num(tmp)) && ~strcmp(tmp,'quit')
+        tmp = input('Invalid input: enter the low threshold, or ''quit'' to continue without photodiode sync:','s');
+      end
+      if strcmp(tmp,'quit')
+        return %todo: need to handle this better
+      else
+        lowThreshold = str2num(tmp);
+        subplot(1,2,1);
+        h = get(gca,'xlim');
+        plot([h(1) h(2)], [lowThreshold lowThreshold]);
+        forLegend121 = horzcat(forLegend121,{'low threshold manual'});
+        legend(forLegend121);
+        subplot(1,2,2);
+        h = get(gca,'xlim');
+        plot([h(1) h(2)], [lowThreshold lowThreshold]);
+        forLegend122 = horzcat(forLegend121,{'low threshold manual'});
+        legend(forLegend122);
+      end
+    end
+    % if two levels, check that frames alternate appropriately
+    if params.numLevels == 2
+      frameSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > lowThreshold);
+      highPeakSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > highThreshold);
+      lowPeakSampleInds = rawPeaks.loc(diodeTrace(rawPeaks.loc) > lowThreshold && diodeTrace(rawPeaks.loc) < highThreshold);
+      frameCountDiff = length(highPeakSampleInds) - length(lowPeakSampleInds);
+      % with correct alternation, smallest interval between subsequent high peaks must tbe greater than longest frame; same for low peaks 
+      frameAlternationError = (min(min(diff(highPeakSampleInds)),min(diff(lowPeakSampleInds))) > max(diff(frameSampleInds)));
+      notValid = (abs(frameCountDiff) > 1 || frameAlternationError);
+      if notValid
+        if abs(frameCountDiff) > 1
+          tmp = input(sprintf('Frame count mismatch: Found %d bright-square frames and %d dark-square frames. Press ENTER to re-start manual check, or q to quit.',...
+            length(highPeakSampleInds),length(lowPeakSampleInds)));
+        end
+        if frameAlternationError
+          tmp = input(sprintf('Found light-dark frame alternation error. Press ENTER to re-start manual check, or q to quit.'));
+        end
+        while ~strcmp('tmp',{'','q'})
+          tmp = imput('Invalid input. Press ENTER to re-start manual check, or q to quit.');
+        end
+        if isempty(tmp)
+          continue
+        else
+          break
+        end
+      end
+    else
+      notValid = 0;
+    end
+  end
+end
+
+if params.saveCalibFig
+  %saveFigure(outDir,sprintf(,runNum), figData, saveFig, exportFig, saveFigData, sprintf('%s, Run %s',dateSubject,runNum));
+end
+if params.saveCalibFile
+  switch params.numLevels
+    case 1
+      save(params.outputCalibrationFilename, 'highThreshold');
+    case 2
+      save(params.outputCalibrationFilename, 'highThreshold', 'lowThreshold');
+    case 3
+      save(params.outputCalibrationFilename, 'highThreshold', 'midThreshold', 'lowThreshold');
+  end
+end  
 end
 
diff --git a/preprocessSpikes.m b/preprocessSpikes.m
index ed93b64..003635f 100644
--- a/preprocessSpikes.m
+++ b/preprocessSpikes.m
@@ -15,6 +15,11 @@
 NEV = openNEV(spikeFilename,'read','nosave','noparse'); %note: add param 'report' for verbose
 Output.VERBOSE('parsing serial IO packets');
 taskTriggers = NEV.Data.SerialDigitalIO;
+if ~ephysParams.needSpikes
+  spikesByChannel = {};
+  channelUnitNames = {};
+  return
+end
 
 %%%%% remove spike data from non-spike channels (e.g. reference electrodes), unsort low quality units, and remove noise units
 spikesByChannel = repmat(struct('times',[],'units',[],'waveforms',[]),length(params.spikeChannels),1);
diff --git a/runAnalyses.m b/runAnalyses.m
index 612a97d..f8f47c8 100644
--- a/runAnalyses.m
+++ b/runAnalyses.m
@@ -189,7 +189,7 @@
   end
 end
 
-if plotSwitch.imagePsth
+if isfield(plotSwitch,'imagePsth') && plotSwitch.imagePsth
   for channel_i = 1:length(channelNames)
     for unit_i = 1:length(channelUnitNames{channel_i})
       if length(channelUnitNames{channel_i}) == 2 && unit_i == 1 %if no isolated unit defined, plot just MUA, not also unsorted (since it's identical)
@@ -206,7 +206,7 @@
   end
 end
 
-if plotSwitch.categoryPsth
+if isfield(plotSwitch,'categoryPsth') && plotSwitch.categoryPsth
   for channel_i = 1:length(channelNames)
     for unit_i = 1:length(channelUnitNames{channel_i})
       if length(channelUnitNames{channel_i}) == 2 && unit_i == 1 %if no isolated unit defined, plot just MUA, not also unsorted (since it's identical)
@@ -298,7 +298,7 @@
         fprintf('%d) %s: %.2f +/- %.2f Hz\n',i,sortedImageLabels{end-i},imageSortedRates(end-i), imFrErrSorted(end-i));
       end
       % preferred images raster plot
-      if plotSwitch.prefImRaster
+      if isfield(plotSwitch,'prefImRaster') && plotSwitch.prefImRaster
         fh = figure();
         raster(spikesByImage(imageSortOrder(1:10)), sortedImageLabels(1:10), psthPre, psthPost, psthImDur, stimTiming.ISI, channel_i, unit_i, colors);
         title(sprintf('Preferred Images, %s %s',channelNames{channel_i},channelUnitNames{channel_i}{unit_i}));
@@ -308,7 +308,7 @@
         end
       end
       % preferred images raster-evoked overlay
-      if plotSwitch.prefImRasterEvokedOverlay
+      if isfield(plotSwitch,'prefImRasterEvokedOverlay') && plotSwitch.prefImRasterEvokedOverlay
         fh = figure();
         rasterEvoked(spikesByImage(imageSortOrder(1:10)), lfpByImage(imageSortOrder(1:10)), sortedImageLabels(1:10), psthPre, psthPost, psthImDur, stimTiming.ISI, lfpPaddedBy, channel_i, colors, 1)
         title(sprintf('Preferred Images, from top, %s %s',channelNames{channel_i},channelUnitNames{channel_i}{unit_i}));
@@ -318,7 +318,7 @@
         end
       end
       % preferred images raster-evoked overlay, with other channels
-      if plotSwitch.prefImMultiChRasterEvokedOverlay
+      if isfield(plotSwitch,'prefImMultiChRasterEvokedOverlay') && plotSwitch.prefImMultiChRasterEvokedOverlay
         fh = figure();
         rasterEvokedMultiCh(spikesByImage(imageSortOrder(1:10)), lfpByImage(imageSortOrder(1:10)), sortedImageLabels(1:10), psthPre, psthPost, psthImDur, stimTiming.ISI, lfpPaddedBy, 1:length(lfpChannels), channelNames, colors)
         title(sprintf('Preferred Images, from top, %s %s',channelNames{channel_i},channelUnitNames{channel_i}{unit_i}));
@@ -328,7 +328,7 @@
         end
       end
       % image preference barplot
-      if plotSwitch.imageTuningSorted
+      if isfield(plotSwitch,'imageTuningSorted') && plotSwitch.imageTuningSorted
         for group_i = 1:length(analysisGroups.stimulusLabelGroups.groups)
           fh = figure();
           groupName = analysisGroups.stimulusLabelGroups.names{group_i};
@@ -433,7 +433,8 @@
   end
   
   % category preference bar plot
-  calcSwitches = [plotSwitch.stimPrefBarPlot, plotSwitch.stimPrefBarPlotEarly, plotSwitch.stimPrefBarPlotLate];
+  calcSwitches = [isfield(plotSwitch,'stimPrefBarPlot') && plotSwitch.stimPrefBarPlot, isfield(plotSwitch,'stimPrefBarPlotEarly') && plotSwitch.stimPrefBarPlotEarly,...
+    isfield(plotSwitch,'stimPrefBarPlotLate') && plotSwitch.stimPrefBarPlotLate];
   for calc_i = 1:length(calcSwitches)
     if ~calcSwitches(calc_i)
       continue
@@ -501,7 +502,8 @@
   end
 
   % tuning curves
-  calcSwitches = [plotSwitch.tuningCurves, plotSwitch.tuningCurvesEarly, plotSwitch.tuningCurvesLate];  
+  calcSwitches = [isfield(plotSwitch,'tuningCurves') && plotSwitch.tuningCurves, isfield(plotSwitch,'tuningCurvesEarly') && plotSwitch.tuningCurvesEarly, ...
+    isfield(plotSwitch,'tuningCurvesLate') && plotSwitch.tuningCurvesLate];  
   for calc_i = 1:length(calcSwitches)
     if ~calcSwitches(calc_i)
       continue
@@ -571,7 +573,8 @@
   % these are the x and y values at which to interpolate RF values
   xi = linspace(min(rfGrid(:,1)),max(rfGrid(:,1)),200);
   yi = linspace(min(rfGrid(:,2)),max(rfGrid(:,2)),200);
-  calcSwitches = [plotSwitch.RF, plotSwitch.rfEarly, plotSwitch.rfLate];  
+  calcSwitches = [isfield(plotSwitch,'RF') && plotSwitch.RF, isfield(plotSwitch,'rfEarly') && plotSwitch.rfEarly,...
+    isfield(plotSwitch,'rfLate') && plotSwitch.rfLate];  
   for calc_i = 1:length(calcSwitches)
     if ~calcSwitches(calc_i)
       continue
@@ -618,11 +621,11 @@
           meanRF = meanRF + imageRF;
           display_map(rfGrid(:,1),rfGrid(:,2),imageRF,xi,yi,2.2857*gridsize,0,saveFig,sprintf('%s %s, %s RF',channelNames{channel_i},channelUnitNames{channel_i}{unit_i},pictureLabels{image_i}),...
             [outDir sprintf('RF_%s_%s_%s_Run%s.png',channelNames{channel_i},channelUnitNames{channel_i}{unit_i},pictureLabels{image_i},runNum)]);
-          if plotSwitch.latencyRF
+          if isfield(plotSwitch,'latencyRF') && plotSwitch.latencyRF
             display_map(rfGrid(:,1),rfGrid(:,2),spikeLatencyRF,xi,yi,2.2857*gridsize,0,saveFig,sprintf('%s %s, %s Latency RF',channelNames{channel_i},channelUnitNames{channel_i}{unit_i},pictureLabels{image_i}),...
               [outDir sprintf('LatencyRF_%s_%s_%s_Run%s.png',channelNames{channel_i},channelUnitNames{channel_i}{unit_i},pictureLabels{image_i},runNum)]);
           end
-          if plotSwitch.evokedPowerRF
+          if isfield(plotSwitch,'evokedPowerRF') && plotSwitch.evokedPowerRF
             display_map(rfGrid(:,1),rfGrid(:,2),evokedPowerRF,xi,yi,2.2857*gridsize,0,saveFig,sprintf('%s %s, %s Evoked Power RF',channelNames{channel_i},channelUnitNames{channel_i}{unit_i},pictureLabels{image_i}),...
               [outDir sprintf('EvokedPowerRF_%s_%s_%s_Run%s.png',channelNames{channel_i},channelUnitNames{channel_i}{unit_i},pictureLabels{image_i},runNum)]);
           end
@@ -718,7 +721,7 @@
 end
 
 %%%% evoked potential plots
-if plotSwitch.evokedPsthMuaMultiCh
+if isfield(plotSwitch,'evokedPsthMuaMultiCh') && plotSwitch.evokedPsthMuaMultiCh
   for group_i = 1:length(analysisGroups.evokedPsthOnePane.groups)
     group = analysisGroups.evokedPsthOnePane.groups{group_i};
     groupName = analysisGroups.evokedPsthOnePane.names{group_i};
@@ -785,7 +788,7 @@
 
 %evoked potentials
 for channel_i = 1:length(lfpChannels)
-  if plotSwitch.evokedByCategory
+  if isfield(plotSwitch,'evokedByCategory') && plotSwitch.evokedByCategory
     for group_i = 1:length(analysisGroups.evokedPotentials.groups)
       group = analysisGroups.evokedPotentials.groups{group_i};
       groupName = analysisGroups.evokedPotentials.names{group_i};
@@ -830,7 +833,7 @@
 end
 
 %evoked analog in potentials
-if plotSwitch.analogInByItem
+if isfield(plotSwitch,'analogInByItem') && plotSwitch.analogInByItem
   for group_i = 1:length(analysisGroups.analogInPotentials.groups)
     group = analysisGroups.analogInPotentials.groups{group_i};
     groupChannels = analysisGroups.analogInPotentials.channels{group_i};
@@ -878,7 +881,7 @@
 end
 
 %evoked analog in derivatives
-if plotSwitch.analogInDerivativesByItem
+if isfield(plotSwitch,'analogInDerivativesByItem') && plotSwitch.analogInDerivativesByItem
   for group_i = 1:length(analysisGroups.analogInPotentials.groups)
     group = analysisGroups.analogInDerivatives.groups{group_i};
     groupChannels = analysisGroups.analogInDerivatives.channels{group_i};
@@ -928,7 +931,7 @@
 
 % lfp - (color) psth subplot
 for channel_i = 1:length(lfpChannels)  
-  if plotSwitch.colorPsthEvoked
+  if isfield(plotSwitch,'colorPsthEvoked') && plotSwitch.colorPsthEvoked
     for group_i = 1:length(analysisGroups.colorPsthEvoked.groups)
       group = analysisGroups.colorPsthEvoked.groups{group_i};
       groupName = analysisGroups.colorPsthEvoked.names{group_i};
@@ -1005,7 +1008,7 @@
   else
     numSubplots = length(channelUnitNames{channel_i}) + 1;
   end
-  if plotSwitch.linePsthEvoked
+  if isfield(plotSwitch,'linePsthEvoked') && plotSwitch.linePsthEvoked
     for group_i = 1:length(analysisGroups.colorPsthEvoked.groups)
       group = analysisGroups.colorPsthEvoked.groups{group_i};
       groupName = analysisGroups.colorPsthEvoked.names{group_i};
@@ -1090,7 +1093,7 @@
 
 
 
-if plotSwitch.runSummary
+if isfield(plotSwitch,'runSummary') && plotSwitch.runSummary
   for channel_i = 1:length(channelNames)
     fh = figure();
     numSubplots = length(imSpikeCounts{channel_i})-1+4;
@@ -1179,7 +1182,7 @@
   end
 end
   
-if plotSwitch.runSummaryImMeanSub
+if isfield(plotSwitch,'runSummaryImMeanSub') && plotSwitch.runSummaryImMeanSub
   for channel_i = 1:length(channelNames)
     fh = figure();
     numSubplots = length(imSpikeCounts{channel_i})-1+4;
@@ -1268,7 +1271,7 @@
   end
 end
   
-if plotSwitch.runSummaryImMeanSubDiv
+if isfield(plotSwitch,'runSummaryImMeanSubDiv') && plotSwitch.runSummaryImMeanSubDiv
   for channel_i = 1:length(channelNames)
     fh = figure();
     numSubplots = length(imSpikeCounts{channel_i})-1+4;
@@ -1358,7 +1361,7 @@
 end
 
 %%% scatter plots %%%
-if plotSwitch.lfpPowerMuaScatter
+if isfield(plotSwitch,'lfpPowerMuaScatter') && plotSwitch.lfpPowerMuaScatter
   for epoch_i = 1:size(frEpochs,1)
     for channel_i = 1:length(channelNames)
       for unit_i = 1:length(channelUnitNames{channel_i})
@@ -1400,7 +1403,7 @@
   end
 end
   
-if plotSwitch.lfpPeakToPeakMuaScatter
+if isfield(plotSwitch,'lfpPeakToPeakMuaScatter') && plotSwitch.lfpPeakToPeakMuaScatter
   for epoch_i = 1:size(frEpochs,1)
     for channel_i = 1:length(channelNames)
       for unit_i = 1:length(channelUnitNames{channel_i})
@@ -1440,12 +1443,12 @@
   end
 end 
   
-if plotSwitch.lfpLatencyMuaLatency
+if isfield(plotSwitch,'lfpLatencyMuaLatency') && plotSwitch.lfpLatencyMuaLatency
   for channel_i = 1:length(channelNames)
   end
 end
   
-if plotSwitch.lfpPowerAcrossChannels && channel_i < length(lfpChannels)
+if isfield(plotSwitch,'lfpPowerAcrossChannels') && plotSwitch.lfpPowerAcrossChannels && channel_i < length(lfpChannels)
   for epoch_i = 1:size(frEpochs,1)
     for channel_i = 1:length(channelNames)
       for channel2_i = channel_i+1:length(lfpChannels)
@@ -1490,7 +1493,7 @@
 end
   
 % lfp p2p amplitude vs lfp p2p amplitude, across channels
-if plotSwitch.lfpPeakToPeakAcrossChannels && channel_i < length(lfpChannels)
+if isfield(plotSwitch,'lfpPeakToPeakAcrossChannels') && plotSwitch.lfpPeakToPeakAcrossChannels && channel_i < length(lfpChannels)
   for epoch_i = 1:size(frEpochs,1)
     for channel_i = 1:length(channelNames)
       for channel2_i = channel_i+1:length(lfpChannels)
@@ -1531,7 +1534,7 @@
 end
   
 % lfp latency vs lfp latency, across channels, defined as max dot product with mean
-if plotSwitch.lfpLatencyShiftAcrossChannels && channel_i < length(lfpChannels)
+if isfield(plotSwitch,'lfpLatencyShiftAcrossChannels') && plotSwitch.lfpLatencyShiftAcrossChannels && channel_i < length(lfpChannels)
   for channel_i = 1:length(channelNames)
     for channel2_i = channel_i+1:length(lfpChannels)
       for group_i = 1:length(analysisGroups.stimulusLabelGroups.groups)
@@ -1623,7 +1626,7 @@
     end
     % single trial evoked lfps; better as multi-channel subplot?
     % todo: programatize on category
-    if plotSwitch.singleTrialLfpByCategory
+    if isfield(plotSwitch,'singleTrialLfpByCategory') && plotSwitch.singleTrialLfpByCategory
       for group_i = 1:length(analysisGroups.lfpSingleTrialsByCategory.groups)
         group = analysisGroups.lfpSingleTrialsByCategory.groups{group_i};
         groupName = analysisGroups.lfpSingleTrialsByCategory.names{group_i};
@@ -1661,7 +1664,7 @@
     end
       
     
-    if plotSwitch.lfpSpectraByCategory
+    if isfield(plotSwitch,'lfpSpectraByCategory') && plotSwitch.lfpSpectraByCategory
       for group_i = 1:length(analysisGroups.spectraByCategory.groups)
         group = analysisGroups.spectraByCategory.groups{group_i};
         groupName = analysisGroups.spectraByCategory.names{group_i};
@@ -1878,7 +1881,7 @@
     end
      
 
-    if plotSwitch.spikeSpectraByCategory
+    if isfield(plotSwitch,'spikeSpectraByCategory') && plotSwitch.spikeSpectraByCategory
       for group_i = 1:length(analysisGroups.spectraByCategory.groups)
         group = analysisGroups.spectraByCategory.groups{group_i};
         groupName = analysisGroups.spectraByCategory.names{group_i};
@@ -2137,7 +2140,7 @@
       lfpByImageEvokedTmp = lfpByImage;
       lfpByImage = lfpByImageInduced;
     end
-    if plotSwitch.spikeSpectraTfByImage
+    if isfield(plotSwitch,'spikeSpectraTfByImage') && plotSwitch.spikeSpectraTfByImage
       for channel_i = 1:length(channelNames)
         for image_i = 1:length(pictureLabels)
           % todo: put in dB conversion and specgramrowave option
@@ -2210,7 +2213,7 @@
         end
       end
     end
-    if plotSwitch.lfpSpectraTfByImage
+    if isfield(plotSwitch,'lfpSpectraTfByImage') && plotSwitch.lfpSpectraTfByImage
       for channel_i = 1:length(channelNames)
         %todo: add Serr plot, significance plot, effect size plot, or similar
         %todo: add support for image calc groups, incl 'pref' wildcard
@@ -2266,7 +2269,7 @@
       lfpByCategory = lfpByCategoryInduced;
     end
     
-    if plotSwitch.tfSpectraByCategory %plots spectra individually by category/channel, and as nChannels x nCategories in group subplot
+    if isfield(plotSwitch,'tfSpectraByCategory') && plotSwitch.tfSpectraByCategory %plots spectra individually by category/channel, and as nChannels x nCategories in group subplot
       for group_i = 1:length(analysisGroups.tfSpectraByCategory.groups)
         group = analysisGroups.tfSpectraByCategory.groups{group_i};
         groupName = analysisGroups.tfSpectraByCategory.names{group_i};
@@ -2701,7 +2704,7 @@
                 close(fh);
               end
               
-              if plotSwitch.couplingPhasesUnwrapped
+              if isfield(plotSwitch,'couplingPhasesUnwrapped') && plotSwitch.couplingPhasesUnwrapped
                 %replot the phases with the implicit mod pi operation undone (so easy to see slope across frequencies
                 fh = figure();
                 lineProps.width = 3;
@@ -2723,7 +2726,7 @@
                 end
               end
               
-              if plotSwitch.couplingPhasesAsOffsets
+              if isfield(plotSwitch,'couplingPhasesAsOffsets') && plotSwitch.couplingPhasesAsOffsets
                 %replot the unwrapped phases as latency differences
                 fh = figure();
                 lineProps.width = 3;
@@ -2745,7 +2748,7 @@
                 end
               end
               
-             if plotSwitch.couplingPhasesPolar
+             if isfield(plotSwitch,'couplingPhasesPolar') && plotSwitch.couplingPhasesPolar
                 %replot the unwrapped phases as latency differences
                 fh = figure();
                 [unwrappedPhases,~] = unwrapPhases(phases);
@@ -2845,7 +2848,7 @@
                   close(fh);
                 end
                 
-                if plotSwitch.couplingPhasesUnwrapped
+                if isfield(plotSwitch,'couplingPhasesUnwrapped') && plotSwitch.couplingPhasesUnwrapped
                   %replot the phases with the implicit mod pi operation undone (so easy to see slope across frequencies
                   fh = figure();
                   lineProps.width = 3;
@@ -2865,7 +2868,7 @@
                   end
                 end
                 
-                if plotSwitch.couplingPhasesAsOffsets
+                if isfield(plotSwitch,'couplingPhasesAsOffsets') && plotSwitch.couplingPhasesAsOffsets
                   %replot the unwrapped phases as latency differences
                   fh = figure();
                   lineProps.width = 3;
@@ -2990,7 +2993,7 @@
                 end
 
 
-                if plotSwitch.tfErrs
+                if isfield(plotSwitch,'tfErrs') && plotSwitch.tfErrs
                   fh = figure(); 
                   subplot(2,3,1);
                   imagesc(t,f,C'); axis xy
@@ -3360,7 +3363,7 @@
               end
 
 
-              if plotSwitch.tfErrs
+              if isfield(plotSwitch,'tfErrs') && plotSwitch.tfErrs
                 fh = figure(); 
                 subplot(2,3,1);
                 imagesc(t,f,C'); axis xy
@@ -3796,7 +3799,7 @@
                 end
 
 
-                if plotSwitch.tfErrs
+                if isfield(plotSwitch,'tfErrs') && plotSwitch.tfErrs
                   fh = figure(); 
                   subplot(2,3,1);
                   imagesc(t,f,C'); axis xy