# HG changeset patch
# User Ray Meddis <rmeddis@essex.ac.uk>
# Date 1310129307 -3600
# Node ID b51bf546ca3f6dead1d0ac2191e6c250e02943cb
# Parent  02aa9826efe07052d80d42a2a1a0c4eedd78a47b
physiologyProb

diff -r 02aa9826efe0 -r b51bf546ca3f Help and reference data/MAP 1_14 DEVELOPMENT log.doc
Binary file Help and reference data/MAP 1_14 DEVELOPMENT log.doc has changed
diff -r 02aa9826efe0 -r b51bf546ca3f Help and reference data/MAP1_14 quick reference.doc
Binary file Help and reference data/MAP1_14 quick reference.doc has changed
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/expGUI_MT.fig
Binary file multithreshold 1.46/expGUI_MT.fig has changed
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/expGUI_MT.m
--- a/multithreshold 1.46/expGUI_MT.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/expGUI_MT.m	Fri Jul 08 13:48:27 2011 +0100
@@ -1,9 +1,9 @@
 % expGUI_MT = 'experimenter GUI for multiThreshold
 % allows the experimenter to design experiments.
 % The *running* of experiments is left to subjGUI.m
-% 
+%
 % There are three kinds of experiments known as 'earOptions':
-% 1. Measurements using real listeners 
+% 1. Measurements using real listeners
 %  'left', 'right',  'diotic', 'dichoticLeft', 'dichoticRight'
 % 2. Measurements using the MAP model as the subject
 %   'MAPmodelListen',  'MAPmodelMultiCh', 'MAPmodelSingleCh'
@@ -88,7 +88,7 @@
 setLocationOfGUIs(handles)
 
 function setLocationOfGUIs(handles)
-global checkForPreviousGUI  % holds screen positioning across repeated calls 
+global checkForPreviousGUI  % holds screen positioning across repeated calls
 scrnsize=get(0,'screensize');
 checkForPreviousGUI=[];
 % if isstruct(checkForPreviousGUI)...
@@ -140,6 +140,7 @@
 % addpath (['..' filesep 'modules'], ['..' filesep 'utilities'], ...
 %     ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
 %     ['..' filesep 'testPrograms'])
+addpath (['..' filesep 'testPrograms'])
 
 % specify all variables that  need to be set on the GUI
 variableNames={'stimulusDelay','maskerDuration','maskerLevel',...
@@ -314,7 +315,7 @@
 paradigm=paradigmNames{chosenOption};
 experiment.paradigm=paradigm;
 
-        %Paradigm: read in all relevant parameters
+%Paradigm: read in all relevant parameters
 % a file must exist with this name 'paradigm_<paradigm>'
 % 'handles' are only occasionally used
 addpath ('paradigms')
@@ -347,7 +348,7 @@
 eval (cmd);
 
 % establish popup menus on the basis of the paradigm file
-set(handles.popupmenuRandomize,'value', betweenRuns.randomizeSequence)	
+set(handles.popupmenuRandomize,'value', betweenRuns.randomizeSequence)
 set(handles.popupmenuPhase,'string', stimulusParameters.maskerPhase)
 if stimulusParameters.includeCue
     set(handles.popupmenuCueNoCue,'value', 1)
@@ -383,14 +384,14 @@
 
 % ------------------------------------------------------ aResetPopupMenus
 function aResetPopupMenus(handles)
-global   stimulusParameters betweenRuns variableNames 
+global   stimulusParameters betweenRuns variableNames
 global targetTypes maskerTypes experiment backgroundTypes
 
 switch experiment.threshEstMethod
     case {'MaxLikelihood','oneIntervalUpDown'}
         set(handles.editstopCriteriaBox, 'string', ...
             num2str(experiment.singleIntervalMaxTrials))
-        
+
     case {'2I2AFC++','2I2AFC+++'}
         set(handles.editstopCriteriaBox, 'string', ...
             num2str(experiment.stopCriteria2IFC))
@@ -410,7 +411,7 @@
 for i=1:length(variableNames)
     if strcmp(variableNames{i},betweenRuns.variableName1)
         variableParameter1ID=i;
-    end   
+    end
     if strcmp(variableNames{i},betweenRuns.variableName2)
         variableParameter2ID=i;
     end
@@ -453,7 +454,7 @@
         set(handles.editStatsModel, 'visible', 'on')
         set(handles.textStatsModel, 'visible', 'on')
         set(handles.pushbuttonStop, 'visible', 'on')
-showModelPushButtons(handles, 0)
+        showModelPushButtons(handles, 0)
         set(handles.editCatchTrialRate, 'visible', 'off')
         set(handles.textCatchTrials, 'visible', 'off')
         set(handles.editcalibrationdB, 'visible', 'off')
@@ -462,25 +463,25 @@
         set(handles.textCue, 'visible', 'off')
         set(handles.editMusicLevel,'visible', 'off')
         set(handles.textMusicLevel,'visible', 'off')
-        
+
     case {'MAPmodel',  'MAPmodelListen', 'MAPmodelMultiCh', 'MAPmodelSingleCh'}
         set(handles.popupmenuCueNoCue, 'visible', 'off')
         set(handles.editStatsModel, 'visible', 'off')
         set(handles.textStatsModel, 'visible', 'off')
         set(handles.pushbuttonStop, 'visible', 'on')
-showModelPushButtons(handles, 1)
+        showModelPushButtons(handles, 1)
         set(handles.editcalibrationdB, 'visible', 'off')
         set(handles.textcalibration, 'visible', 'off')
         set(handles.textCue, 'visible', 'off')
         set(handles.editMusicLevel,'visible', 'off')
         set(handles.textMusicLevel,'visible', 'off')
-        
+
     otherwise
         % i.e. using real subjects (left, right, diotic, dichotic)
         set(handles.editStatsModel, 'visible', 'off')
         set(handles.textStatsModel, 'visible', 'off')
         set(handles.pushbuttonStop, 'visible', 'off')
-showModelPushButtons(handles, 0)
+        showModelPushButtons(handles, 0)
         set(handles.editCatchTrialRate, 'visible', 'on')
         set(handles.textCatchTrials, 'visible', 'on')
         set(handles.editcalibrationdB, 'visible', 'on')
@@ -497,7 +498,7 @@
         set(handles.textCue,'visible', 'on')
         set(handles.editstopCriteriaBox, 'string', ...
             num2str(experiment.singleIntervalMaxTrials))
-        
+
         if stimulusParameters.includeCue==0
             set(handles.editcueTestDifference,'visible', 'off')
             set(handles.textcueTestDifference,'visible', 'off')
@@ -505,7 +506,7 @@
             set(handles.editcueTestDifference,'visible', 'on')
             set(handles.textcueTestDifference,'visible', 'on')
         end
-        
+
     case {'2I2AFC++','2I2AFC+++'}
         set(handles.editCatchTrialRate, 'visible', 'off')
         set(handles.textCatchTrials, 'visible', 'off')
@@ -609,7 +610,7 @@
 
 % ------------------------------------------------ pushbuttonRun_Callback
 function pushbuttonRun_Callback(hObject, eventdata, handles)
-global checkForPreviousGUI % holds screen positioning across repeated calls 
+global checkForPreviousGUI % holds screen positioning across repeated calls
 global experiment betweenRuns paradigmNames errormsg
 checkForPreviousGUI.GUIposition=get(handles.figure1,'position');
 experiment.singleShot=0;
@@ -621,8 +622,8 @@
         set(handles.edittargetDuration,'string', num2str(0.25))
         set(handles.editstopCriteriaBox,'string','10') % nTrials
         run (handles)
-        
-        if ~isempty(errormsg)
+
+        if strcmp(errormsg,'manually stopped')
             disp(errormsg)
             optionNo=strmatch('profile',paradigmNames);
             set(handles.popupmenuParadigm,'value',optionNo);
@@ -630,14 +631,14 @@
             aParadigmSelection(handles)
             return
         end
-        
+
         global resultsTable
         longTone=resultsTable(2:end,2:end);
-        
+
         set(handles.edittargetDuration,'string', num2str(0.016))
         set(handles.editstopCriteriaBox,'string','10') % nTrials
         run (handles)
-        if ~isempty(errormsg)
+        if strcmp(errormsg,'manually stopped')
             disp(errormsg)
             optionNo=strmatch('profile',paradigmNames);
             set(handles.popupmenuParadigm,'value',optionNo);
@@ -646,7 +647,7 @@
             aParadigmSelection(handles)
             return
         end
-        
+
         shortTone=resultsTable(2:end,2:end);
 
         % use these threshold for TMC
@@ -656,10 +657,10 @@
         aParadigmSelection(handles)
         set(handles.edittargetLevel,'string', thresholds16ms+10);
         set(handles.editstopCriteriaBox,'string','10')  % nTrials
-        pause(.1) 
+        pause(.1)
         run (handles)
-        
-        if ~isempty(errormsg)
+
+        if strcmp(errormsg,'manually stopped')
             disp(errormsg)
             optionNo=strmatch('profile',paradigmNames);
             set(handles.popupmenuParadigm,'value',optionNo);
@@ -672,17 +673,29 @@
         TMC=resultsTable(2:end,2:end);
         gaps=resultsTable(2:end,1);
         BFs=resultsTable(1, 2:end);
-        
+
         % use these threshold for IFMC
         optionNo=strmatch('IFMC',paradigmNames);
         set(handles.popupmenuParadigm,'value',optionNo);
         aParadigmSelection(handles)
         set(handles.edittargetLevel,'string', thresholds16ms+10);
         set(handles.editstopCriteriaBox,'string','10')  % nTrials
-        pause(.1) 
+        pause(.1)
         run (handles)
 
-        if ~isempty(errormsg)
+
+        IFMCs=resultsTable(2:end,2:end);
+        offBFs=resultsTable(2:end,1);
+
+        % reset original paradigm
+        optionNo=strmatch('profile',paradigmNames);
+        set(handles.popupmenuParadigm,'value',optionNo);
+        aParadigmSelection(handles)
+
+        save profile longTone shortTone gaps BFs TMC offBFs IFMCs
+        plotProfile(longTone,shortTone,gaps,BFs,TMC,offBFs,IFMCs)
+
+        if strcmp(errormsg,'manually stopped')
             disp(errormsg)
             optionNo=strmatch('profile',paradigmNames);
             set(handles.popupmenuParadigm,'value',optionNo);
@@ -690,19 +703,7 @@
             experiment.stop=-0;
             aParadigmSelection(handles)
             return
-        end       
-
-        IFMCs=resultsTable(2:end,2:end);
-        offBFs=resultsTable(2:end,1);
-        
-        % reset original paradigm
-        optionNo=strmatch('profile',paradigmNames);
-        set(handles.popupmenuParadigm,'value',optionNo);
-        aParadigmSelection(handles)
-
-        save profile longTone shortTone gaps BFs TMC offBFs IFMCs
-        plotProfile(longTone,shortTone,gaps,BFs,TMC,offBFs,IFMCs)
-        
+        end
     otherwise
         run (handles)
         experiment.stop=0;
@@ -713,7 +714,7 @@
 tic
 expGUIhandles=handles;
 set(handles.pushbuttonStop, 'backgroundColor', [.941 .941 .941])
-set(handles.editparamChanges,'visible','off')
+% set(handles.editparamChanges,'visible','off')
 
 % message box white (removes any previous error message)
 set(handles.textMSG,...
@@ -742,7 +743,7 @@
 
 % --- Executes on button press in pushbuttonSingleShot.
 function pushbuttonSingleShot_Callback(hObject, eventdata, handles)
-global experiment 
+global experiment
 experiment.singleShot=1;
 
 % special test for spontaneous activity
@@ -764,7 +765,7 @@
 % ------------------------------------------aReadAndCheckParameterBoxes
 function errorMsg=aReadAndCheckParameterBoxes(handles)
 global experiment  stimulusParameters betweenRuns  statsModel
-global variableNames  LevittControl
+global variableNames  LevittControl paramChanges
 % When the program sets the parameters all should be well
 % But when the user changes them...
 
@@ -780,8 +781,8 @@
 
 switch experiment.ear
     case { 'MAPmodel', 'MAPmodelMultiCh', ...
-            'MAPmodelSingleCh', 'MAPmodelListen'}         
-        % MAPmodel writes forced parameter settings to the screen 
+            'MAPmodelSingleCh', 'MAPmodelListen'}
+        % MAPmodel writes forced parameter settings to the screen
         %  so that they can be read from there
         set(handles.popupmenuRandomize,'value',2)       % fixed sequence
         set(handles.editstimulusDelay,'string','0.01')  % no stimulus delay
@@ -840,13 +841,13 @@
 end
 
 % calibration
-%  this is used to *reduce* the output signal from what it otherwise 
+%  this is used to *reduce* the output signal from what it otherwise
 %  would be
 % signal values are between 1 - 2^23
 %  these are interpreted as microPascals between -29 dB and 128 dB SPL
-% calibrationdB adjusts these values to compensate for equipment 
+% calibrationdB adjusts these values to compensate for equipment
 %  characteristics
-%  this will change the range. e.g. a 7 dB calibration will yield 
+%  this will change the range. e.g. a 7 dB calibration will yield
 %   a range of -36 to 121 dB SPL
 % Calibration is not used when modelling. Values are treated as dB SPL
 stimulusParameters.calibrationdB=...
@@ -946,7 +947,7 @@
         % start value for step until reduced
         LevittControl.startLevelStep= stimulusParameters.WRVsteps(1);
         % reduced step size
-        LevittControl.steadyLevittStep= stimulusParameters.WRVsteps(2); 
+        LevittControl.steadyLevittStep= stimulusParameters.WRVsteps(2);
         LevittControl.TurnsToSmallSteps= 2;
         LevittControl.useLastNturns= 2*experiment.peaksUsed;
         LevittControl.minReversals= ...
@@ -1006,7 +1007,7 @@
 
 % Check WRVstartValues for length and compatibility with randomization
 % only one start value supplied so all start values are the same
-if length(stimulusParameters.WRVstartValues)==1     
+if length(stimulusParameters.WRVstartValues)==1
     stimulusParameters.WRVstartValues= ...
         repmat(stimulusParameters.WRVstartValues, 1, ...
         length(betweenRuns.variableList1)...
@@ -1051,7 +1052,7 @@
                 stimulusParameters.maskerDuration)
             addToMsg(...
                 'Warning: masker and target duration not the same.',1,1)
-        end        
+        end
         if ~isequal(stimulusParameters.maskerLevel, ...
                 stimulusParameters.targetLevel)
             addToMsg(['Warning: masker and target level different.'...
@@ -1059,6 +1060,10 @@
         end
 end
 
+% identify model parameter changes if any
+paramChanges=get(handles.editparamChanges,'string');
+eval(paramChanges);
+
 % -------------------------------------------- aSetSampleRate
 function aSetSampleRate(sampleRate, handles)
 global  stimulusParameters
@@ -1106,24 +1111,24 @@
     case {'statsModelLogistic', 'statsModelRareEvent'}
         set(handles.editStatsModel, 'visible', 'off')
         set(handles.textStatsModel, 'visible', 'off')
-        
+
         % default psychometric bin size and logistic slopes
         set(handles.popupmenuRandomize,'value',2)   % fixed sequence
         set(handles.editName,'string', 'statsModel')
         %         experiment.headphonesUsed=0;
-        
+
     case {'MAPmodelListen', 'MAPmodelMultiCh', 'MAPmodelSingleCh'}
         stimulusParameters.includeCue=0;						 % no cue
         set(handles.popupmenuCueNoCue,'value', 2)
-        
+
         set(handles.editCatchTrialRate,'string','0 0  2 ');%no catch trials
         set(handles.editName,'string', 'Normal')			% force name
         experiment.name=get(handles.editName,'string');	% read name back
-        set(handles.editcalibrationdB,'string','0')	
-        
+        set(handles.editcalibrationdB,'string','0')
+
         set(handles.popupmenuRandomize,'value',2)       % fixed sequence
         set(handles.editstimulusDelay,'string','0')
-        
+
         set(handles.editSaveData,'string', '0')
         set(handles.editSubjectFont,'string', '10');
         experiment.MacGThreshold=0; % num MacG spikes to exceed threshold
@@ -1178,7 +1183,7 @@
             case 'training'
                 experiment.possLogSlopes=0.5;
         end
-        
+
     case 'oneIntervalUpDown'
         experiment.functionEstMethod='logisticLS';
         set(handles.textstopCriteria,'string', 'stop criteria \ maxTrials')
@@ -1189,7 +1194,7 @@
             otherwise
                 experiment.allowCatchTrials= 1;
         end
-        
+
     case {'2I2AFC++',  '2I2AFC+++'}
         LevittControl.rule='++'; %  e.g. '++' or '+++'
         experiment.singleIntervalMaxTrials=experiment.stopCriteria2IFC;
@@ -1209,7 +1214,7 @@
 % NB responsibility for this is now transferred to the paradigm file
 switch experiment.threshEstMethod
     % only one value required for level change
-    case {'2I2AFC++', '2A2AIFC+++'}		
+    case {'2I2AFC++', '2A2AIFC+++'}
         stimulusParameters.subjectText=...
             'did the tone occur in window 1 or 2?';
     case {'MaxLikelihood',  'oneIntervalUpDown'};
@@ -1399,41 +1404,49 @@
 drawnow
 
 function pushbuttonOME_Callback(hObject, eventdata, handles)
-global experiment
+global experiment paramChanges
 aReadAndCheckParameterBoxes(handles);
-testOME(experiment.name);
+testOME(experiment.name, paramChanges);
 
 function pushbuttonBM_Callback(hObject, eventdata, handles)
-global  stimulusParameters experiment
+global  stimulusParameters experiment paramChanges
 aReadAndCheckParameterBoxes(handles);
 relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];
-
+AN_spikesOrProbability='probability';
 testBM(stimulusParameters.targetFrequency, ...
-    experiment.name,relativeFrequencies);
+    experiment.name,relativeFrequencies, AN_spikesOrProbability, ...
+    paramChanges);
 
 function pushbuttonAN_Callback(hObject, eventdata, handles)
-global stimulusParameters
+global stimulusParameters experiment paramChanges
 aReadAndCheckParameterBoxes(handles);
 % now carry out tests
 showPSTHs=0;
 targetFrequency=stimulusParameters.targetFrequency(1);
 BFlist=targetFrequency;
 
-testAN(targetFrequency,BFlist);
+testAN(targetFrequency,BFlist,-10:10:90,experiment.name, paramChanges);
 
 function pushbuttonPhLk_Callback(hObject, eventdata, handles)
+global experiment
 aReadAndCheckParameterBoxes(handles);
-testPhaseLocking
+testPhaseLocking(experiment.name)
 
 function pushbuttonSYN_Callback(hObject, eventdata, handles)
+global stimulusParameters experiment paramChanges
 aReadAndCheckParameterBoxes(handles);
-testSynapse
+% now carry out tests
+showPSTHs=0;
+targetFrequency=stimulusParameters.targetFrequency(1);
+BFlist=targetFrequency;
+testSynapse(BFlist,experiment.name, paramChanges)
 
 function pushbuttonFM_Callback(hObject, eventdata, handles)
-global stimulusParameters experiment
+global stimulusParameters experiment paramChanges
 aReadAndCheckParameterBoxes(handles);
 showPSTHs=1;
-testFM(stimulusParameters.targetFrequency(1),experiment.name, showPSTHs)
+testFM(stimulusParameters.targetFrequency(1),experiment.name, ...
+    showPSTHs, paramChanges)
 
 function popupmenuPhase_Callback(hObject, eventdata, handles)
 global stimulusParameters
@@ -1447,10 +1460,9 @@
 function pushbuttonParams_Callback(hObject, eventdata, handles)
 global experiment stimulusParameters
 aReadAndCheckParameterBoxes(handles);
-% print model parameters using the 'name' box (e.g. CTa -> MAPparamsCTa)
 showParams=1; BFlist=-1;
 paramChanges=get(handles.editparamChanges,'string');
-eval(paramChanges)
+eval(paramChanges);
 
 paramFunctionName=['method=MAPparams' experiment.name ...
     '(BFlist, stimulusParameters.sampleRate, showParams,paramChanges);'];
@@ -1459,10 +1471,10 @@
 
 % --- Executes on button press in pushbuttonRP.
 function pushbuttonRP_Callback(hObject, eventdata, handles)
-global experiment stimulusParameters
+global experiment stimulusParameters paramChanges
 aReadAndCheckParameterBoxes(handles);
 % now carry out test
-testRP(stimulusParameters.targetFrequency,experiment.name)
+testRP(stimulusParameters.targetFrequency,experiment.name, paramChanges)
 
 % function handles % ??
 
@@ -1755,7 +1767,7 @@
 end
 
 
-    
+
 
 
 function editparamChanges_Callback(hObject, eventdata, handles)
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/nextStimulus.m
--- a/multithreshold 1.46/nextStimulus.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/nextStimulus.m	Fri Jul 08 13:48:27 2011 +0100
@@ -7,13 +7,13 @@
 errormsg='';
 
 % interrupt by 'stop' button
-if experiment.stop
-    disp('******** experiment manually stopped  *****************')
-    experiment.status= 'waitingForStart';
-    errormsg='manually stopped';
-    addToMsg(errormsg,1)
-    return
-end
+% if experiment.stop
+%     disp('******** experiment manually stopped  *****************')
+%     experiment.status= 'waitingForStart';
+%     errormsg='manually stopped';
+%     addToMsg(errormsg,1)
+%     return
+% end
 
 % -----------------------------------------choose catch trials at random
 % catch trials are for subject threshold measurements only
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/old files/MAPmodel.m
--- a/multithreshold 1.46/old files/MAPmodel.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/old files/MAPmodel.m	Fri Jul 08 13:48:27 2011 +0100
@@ -4,7 +4,7 @@
 global outerMiddleEarParams DRNLParams AN_IHCsynapseParams
 
 savePath=path;
-addpath('..\MAP')
+addpath(['..' filesep 'MAP'], ['..' filesep 'utilities'])
 modelResponse=[];
 MacGregorResponse=[];
 
@@ -25,12 +25,12 @@
         MAPparamsName, AN_spikesOrProbability);
     
 if showPlotsAndDetails
-    options.showModelParameters=0;
+    options.printModelParameters=0;
     options.showModelOutput=1;
     options.printFiringRates=1;
     options.showACF=0;
     options.showEfferent=1;
-    UTIL_showMAP(options, paramChanges)
+    UTIL_showMAP(options)
 end
 
 % No response,  probably caused by hitting 'stop' button
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/paradigms/paradigm_IFMC.m
--- a/multithreshold 1.46/paradigms/paradigm_IFMC.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/paradigms/paradigm_IFMC.m	Fri Jul 08 13:48:27 2011 +0100
@@ -4,7 +4,7 @@
 paradigm_training(handles) % default
 
 stimulusParameters.WRVname='maskerLevel';
-stimulusParameters.WRVstartValues=-10;
+stimulusParameters.WRVstartValues=50;
 stimulusParameters.WRVsteps= [-10 -2];
 stimulusParameters.WRVlimits=[-30 110];
 
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/paradigms/paradigm_TMC.m
--- a/multithreshold 1.46/paradigms/paradigm_TMC.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/paradigms/paradigm_TMC.m	Fri Jul 08 13:48:27 2011 +0100
@@ -4,7 +4,7 @@
 paradigm_training(handles) % default
 
 stimulusParameters.WRVname='maskerLevel';
-stimulusParameters.WRVstartValues=-10;
+stimulusParameters.WRVstartValues=50;
 stimulusParameters.WRVsteps= [-10 -4];
 stimulusParameters.WRVlimits=[-30 110];
 
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/paradigms/paradigm_profile.m
--- a/multithreshold 1.46/paradigms/paradigm_profile.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/paradigms/paradigm_profile.m	Fri Jul 08 13:48:27 2011 +0100
@@ -5,7 +5,6 @@
 
 betweenRuns.variableName1='targetFrequency';
 betweenRuns.variableList1=[250 500 1000 2000 4000 8000 ];
-betweenRuns.variableList1=str2num(get(handles.edittargetFrequency,'string'));
 betweenRuns.variableName2='targetDuration';
 betweenRuns.variableList2=0.016;
 betweenRuns.randomizeSequence=1; % 'random sequence'
@@ -15,7 +14,7 @@
 
 stimulusParameters.targetType='tone';
 stimulusParameters.targetPhase='sin';
-stimulusParameters.targetFrequency=1000;
+stimulusParameters.targetFrequency=[250 500 1000 2000 4000 8000 ];
 stimulusParameters.targetDuration=0.016;
 stimulusParameters.targetLevel=stimulusParameters.WRVstartValues(1);
 
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/plotProfile.m
--- a/multithreshold 1.46/plotProfile.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/plotProfile.m	Fri Jul 08 13:48:27 2011 +0100
@@ -5,63 +5,118 @@
     load profile
 end
 
-normLongTone=[	11.4	1.55	-13.5	-6.35	-6.4	7.45];
-normShortTone=[	23.85	18.9	9.85	10.6	9.55	21.9];
+% comparison data (e.g. participants)
+% rows are BFs
 
-normGaps=0.01:0.01:0.09;
-normTMC=	[
-28.5	35.0	49.3	70.1	80.5	85.5    NaN      NaN    NaN;			
-31.1	44.3	48.4	59.5	56.4	76.7	70.2	82.4	76.3;
-33.4	38.4	48.8	55.8	64.5	78.7	84.2	88.3	90.3;
-25.4	37.0	49.2	49.7	58.2	69.6	87.7	95.8	93.0;
-18.2	23.5	27.4	41.5	64.3	82.1	86.7	91.2	NaN;
-32.5	35.8	43.5	52.1	69.1	78.6	86.6	86.0	NaN;
-    ];
-normTMC=normTMC';
+% -------------------------------------------JSan
+compareBFs=[250	500	1000	2000	3000
+];
+compareLongTone=  [67.46485	56.95655	65.01985	61.46655	73.33265
+];
+compareShortTone=[	72.3185	63.2818	69.0373	65.2853	76
+];
 
-normIFMC=[
-50	42	34	35	34	33	37;
-58	51	38	33	28	41	49;
-57	41	27	20	28	37	66;
-61	49	27	20	34	68	79;
-67	45	27	22	46	74	87;
-62	62	43	22	47	56	83;
+compareGaps=[0.02 0.05 0.08];
+compareTMC=	[
+    84	69	77	75	93
+88	73	81	79	97
+95	79	85	83	98
 ];
-normIFMC=normIFMC';
+
+compareMaskerFreqs=[0.7  0.9 1 1.1  1.3 ];
+compareIFMCs=[83.1698	77.3165	79.8474	82.9074	82.3294
+80.9667	73.6653	80.9446	80.7005	79.0022
+82.0135	71.2284	78.7345	74.3342	84.126
+79.3348	70.3347	78.5769	79.7274	91.9849
+79.2308	76.3164	83.6881	86.2105	NaN
+  ];
+
+% % -------------------------------------------JE
+% BFs=[250 500 1000 2000 4000 8000];
+% compareLongTone=  [32 30	31	40	54 NaN];
+% compareShortTone=[	49	50	47	56	63	NaN];
+% 
+% compareGaps=0.01:0.01:0.09;
+% compareTMC=	[
+%     69	83	82	NaN	NaN	NaN	NaN	NaN	NaN
+%     61	68	79	88	93	NaN	NaN	NaN	NaN
+%     63	69	79	84	92	NaN	NaN	NaN	NaN
+%     67	71	75	80	82	84	88	93	NaN
+%     82	82	86	86	NaN	83	88	90	75
+%     NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
+%     ];
+% compareTMC=compareTMC';
+% 
+% compareMaskerFreqs=[0.5	0.7	0.9	1	1.1	1.3	1.6];
+% compareIFMCs=[
+%     64	60	58	58	57	60	64
+%     65	63	58	55	54	59	69
+%     68	64	60	59	62	73	79
+%     76	75	71	67	68	71	77
+%     79	71	68	69	73	75	77
+%     76	73	75	75	76	80	NaN
+%     ];
+% compareIFMCs=compareIFMCs';
+
+% -------------------------------------------CMR
+% CMR 
+% BFs=[250 500 1000 2000 4000 8000];
+% compareLongTone=[	11.4	1.55	-13.5	-6.35	-6.4	7.45];
+% compareShortTone=[	23.85	18.9	9.85	10.6	9.55	21.9];
+% 
+% compareGaps=0.01:0.01:0.09;
+% compareTMC=	[
+%     28.5	35.0	49.3	70.1	80.5	85.5    NaN      NaN    NaN;
+%     31.1	44.3	48.4	59.5	56.4	76.7	70.2	82.4	76.3;
+%     33.4	38.4	48.8	55.8	64.5	78.7	84.2	88.3	90.3;
+%     25.4	37.0	49.2	49.7	58.2	69.6	87.7	95.8	93.0;
+%     18.2	23.5	27.4	41.5	64.3	82.1	86.7	91.2	NaN;
+%     32.5	35.8	43.5	52.1	69.1	78.6	86.6	86.0	NaN;
+%     ];
+% compareTMC=compareTMC';
+% 
+% compareMaskerFreqs=[0.5	0.7	0.9	1	1.1	1.3	1.6];
+% compareIFMCs=[
+%     50	42	34	35	34	33	37;
+%     58	51	38	33	28	41	49;
+%     57	41	27	20	28	37	66;
+%     61	49	27	20	34	68	79;
+%     67	45	27	22	46	74	87;
+%     62	62	43	22	47	56	83;
+%     ];
+% compareIFMCs=compareIFMCs';
 
 % absolute thresholds
 figure(90), clf
 subplot(2,1,2)
 semilogx(BFs,longTone,'ko-','lineWidth',2); hold on
 semilogx(BFs,shortTone,'bo-','lineWidth',2); hold on
-semilogx(BFs,normLongTone,'ko:'); hold on
-semilogx(BFs,normShortTone,'bo:'); hold on
+semilogx(compareBFs,compareLongTone,'ko:'); hold on
+semilogx(compareBFs,compareShortTone,'bo:'); hold on
 ylim([0 100])
 
 % TMC
 for BFno=1:length(BFs)
     subplot(2,6,BFno)
     plot(gaps,TMC(:,BFno)-longTone(BFno),'r','lineWidth',3), hold on
-    plot(normGaps,normTMC(:,BFno)-longTone(BFno),'k:')
-    ylim([-10 90])
+    plot(gaps,TMC(:,BFno),'b','lineWidth',3), hold on
+    ylim([-10 110])
     xlim([0.01 0.1])
-    grid on    
+    grid on
     if BFno==1
         ylabel('masker dB SL')
         xlabel('gap')
-        text(0.02,80,' TMC','backgroundColor','w')
+%         text(0.02,80,' TMC','backgroundColor','w')
     end
-        title([num2str(BFs(BFno)) ' Hz'])
-        set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'})
+    title([num2str(BFs(BFno)) ' Hz'])
+    set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'})
 end
 
 % IFMCs
 for BFno=1:length(BFs)
-    BF=BFs(BFno);
     freq=offBFs'*BFs(BFno);
     subplot(2,1,2)
     semilogx(freq,IFMCs(:,BFno),'r','lineWidth',3), hold on
-    semilogx(freq,normIFMC(:,BFno),'k:')
     ylim([0 100])
     xlim([100 12000])
     grid on
@@ -69,3 +124,29 @@
 xlabel('frequency (Hz)')
 ylabel('masker dB / probe dB')
 set(gca,'XTick',BFs)
+
+for BFno=1:length(compareBFs)
+    subplot(2,6,BFno)
+    plot(compareGaps,compareTMC(:,BFno)-longTone(BFno),'k:')
+    plot(compareGaps,compareTMC(:,BFno),'k:')
+    ylim([-10 110])
+    xlim([0.01 0.1])
+    grid on
+    if BFno==1
+        ylabel('masker dB SL')
+        xlabel('gap')
+%         text(0.02,80,' TMC','backgroundColor','w')
+    end
+    title([num2str(BFs(BFno)) ' Hz'])
+    set(gca,'XTick',[ 0.1],'xTickLabel', { '0.1'})
+end
+
+% IFMCs
+for BFno=1:length(compareBFs)
+    compareFreq=compareMaskerFreqs'*BFs(BFno);
+    subplot(2,1,2)
+    semilogx(compareFreq,compareIFMCs(:,BFno),'k:')
+    ylim([0 100])
+    xlim([100 12000])
+    grid on
+end
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/printReport.m
--- a/multithreshold 1.46/printReport.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/printReport.m	Fri Jul 08 13:48:27 2011 +0100
@@ -6,7 +6,7 @@
 
 global experiment stimulusParameters betweenRuns withinRuns statsModel   audio
 global LevittControl expGUIhandles
-
+global paramChanges
 
 global inputStimulusParams OMEParams DRNLParams
 global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
@@ -24,7 +24,9 @@
     printReportGuide.showTracks=experiment.printTracks;
     printReportGuide.fileName=[];
     if experiment.saveData
-        saveFileName=['savedData/' experiment.name '_' experiment.date '_' experiment.paradigm];
+        saveFileName=...
+            ['savedData/' experiment.name '_' ...
+            experiment.date '_' experiment.paradigm];
     else
         % save this data (just in case)
         saveFileName=['savedData/mostRecentResults'];
@@ -84,6 +86,7 @@
 msg=printTabTable(resultsTable,  headers);
 addToMsg(msg,0)
 fprintf('\n')
+disp(paramChanges)
 
 % sort tracks into the same order
 betweenRuns.levelTracks=betweenRuns.levelTracks(idx1);
@@ -181,6 +184,7 @@
 end
 
 fprintf('\nparadigm:\t%s\n ', experiment.paradigm)
+disp(paramChanges)
 
 % ------------------------------------------------------- sortTablesForPrinting
 function table= sortTablesForPrinting(idx1,idx2, var1values,var2values, x)
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/profile.mat
Binary file multithreshold 1.46/profile.mat has changed
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/savedData/mostRecentResults.mat
Binary file multithreshold 1.46/savedData/mostRecentResults.mat has changed
diff -r 02aa9826efe0 -r b51bf546ca3f multithreshold 1.46/subjGUI_MT.m
--- a/multithreshold 1.46/subjGUI_MT.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/multithreshold 1.46/subjGUI_MT.m	Fri Jul 08 13:48:27 2011 +0100
@@ -1349,11 +1349,11 @@
     %  without waiting for button press
     startNewRun(handles)
 
-    % show sample Rate on GUI; it must be set in MAPparams
+    % show sample Rate on GUI; it must be set in MAPparams ##??
     set(expGUIhandles.textsampleRate,'string',...
         num2str(stimulusParameters.sampleRate))
 
-    if experiment.singleShot
+    if experiment.singleShot % ##??
         AN_IHCsynapseParams.plotSynapseContents=1;
     else
         AN_IHCsynapseParams.plotSynapseContents=0;
@@ -1424,7 +1424,9 @@
 function [modelResponse, MacGregorResponse]=MAPmodel
 
 global experiment stimulusParameters audio withinRuns
-global outerMiddleEarParams DRNLParams AN_IHCsynapseParams
+% global outerMiddleEarParams DRNLParams AN_IHCsynapseParams
+global ICoutput ANdt dt savedBFlist ANprobRateOutput expGUIhandles
+global paramChanges
 
 savePath=path;
 addpath(['..' filesep 'MAP'], ['..' filesep 'utilities'])
@@ -1444,12 +1446,12 @@
 % ---------------------------------------------- run Model
 MAPparamsName=experiment.name;
 showPlotsAndDetails=experiment.MAPplot;
+
+% important buried constant ##??
 AN_spikesOrProbability='spikes';
 AN_spikesOrProbability='probability';
 
 % [response, method]=MAPsequenceSeg(audio, method, 1:8);
-global ICoutput ANdt dt savedBFlist ANprobRateOutput expGUIhandles
-global stimulusParameters experiment
 
 if sum(strcmp(experiment.ear,{'MAPmodelMultiCh', 'MAPmodelListen'}))
     % use BFlist specified in MAPparams file
@@ -1458,7 +1460,8 @@
     BFlist=stimulusParameters.targetFrequency;
 end
 paramChanges=get(expGUIhandles.editparamChanges,'string');
-eval(paramChanges)
+% convert from string to a cell array
+eval(paramChanges);
 
 MAP1_14(audio, stimulusParameters.sampleRate, BFlist,...
     MAPparamsName, AN_spikesOrProbability, paramChanges);
diff -r 02aa9826efe0 -r b51bf546ca3f parameterStore/MAPparamsEndo.m
--- a/parameterStore/MAPparamsEndo.m	Fri Jul 01 12:59:47 2011 +0100
+++ b/parameterStore/MAPparamsEndo.m	Fri Jul 08 13:48:27 2011 +0100
@@ -129,7 +129,6 @@
 %  #5 IHC_RP
 IHC_cilia_RPParams.Cab=	4e-012;         % IHC capacitance (F)
 % IHC_cilia_RPParams.Cab=	1e-012;         % IHC capacitance (F)
-IHC_cilia_RPParams.Et=	0.100;          % endocochlear potential (V)
 IHC_cilia_RPParams.Et=	0.05;          % endocochlear potential (V)
 
 IHC_cilia_RPParams.Gk=	2e-008;         % 1e-8 potassium conductance (S)
diff -r 02aa9826efe0 -r b51bf546ca3f parameterStore/MAPparamsJE.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/parameterStore/MAPparamsJE.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,302 @@
+function method=MAPparamsJE ...
+    (BFlist, sampleRate, showParams, paramChanges)
+% MAPparams<> establishes a complete set of MAP parameters
+% Parameter file names must be of the form <MAPparams> <name>
+%
+% input arguments
+%  BFlist     (optional) specifies the desired list of channel BFs
+%    otherwise defaults set below
+%  sampleRate (optional), default is 50000.
+%  showParams (optional) =1 prints out the complete set of parameters
+% output argument
+%  method passes a miscelleny of values
+
+global inputStimulusParams OMEParams DRNLParams IHC_cilia_RPParams
+global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
+global MacGregorParams MacGregorMultiParams  filteredSACFParams
+global experiment % used by calls from multiThreshold only
+
+
+currentFile=mfilename;                      % i.e. the name of this mfile
+method.parameterSource=currentFile(10:end); % for the record
+
+efferentDelay=0.010;
+method.segmentDuration=efferentDelay;
+
+if nargin<3, showParams=0; end
+if nargin<2, sampleRate=50000; end
+if nargin<1 || BFlist(1)<0 % if BFlist= -1, set BFlist to default
+    lowestBF=250; 	highestBF= 8000; 	numChannels=21;
+    % 21 chs (250-8k)includes BFs at 250 500 1000 2000 4000 8000
+    BFlist=round(logspace(log10(lowestBF),log10(highestBF),numChannels));
+end
+% BFlist=1000;
+
+% preserve for backward campatibility
+method.nonlinCF=BFlist; 
+method.dt=1/sampleRate; 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% set  model parameters
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%  #1 inputStimulus
+inputStimulusParams=[];
+inputStimulusParams.sampleRate= sampleRate; 
+
+%%  #2 outerMiddleEar
+OMEParams=[];  % clear the structure first
+% outer ear resonances band pass filter  [gain lp order  hp]
+OMEParams.externalResonanceFilters=      [ 10 1 1000 4000];
+                                       
+% highpass stapes filter  
+%  Huber gives 2e-9 m at 80 dB and 1 kHz (2e-13 at 0 dB SPL)
+OMEParams.OMEstapesLPcutoff= 1000;
+OMEParams.stapesScalar=	     45e-9;
+
+% Acoustic reflex: maximum attenuation should be around 25 dB Price (1966)
+% i.e. a minimum ratio of 0.056.
+% 'spikes' model: AR based on brainstem spiking activity (LSR)
+OMEParams.rateToAttenuationFactor=0.006;   % * N(all ICspikes)
+% OMEParams.rateToAttenuationFactor=0;   % * N(all ICspikes)
+
+% 'probability model': Ar based on AN firing probabilities (LSR)
+OMEParams.rateToAttenuationFactorProb=0.01;% * N(all ANrates)
+% OMEParams.rateToAttenuationFactorProb=0;% * N(all ANrates)
+
+% asymptote should be around 100-200 ms
+OMEParams.ARtau=.05; % AR smoothing function
+% delay must be longer than the segment length
+OMEParams.ARdelay=efferentDelay;  %Moss gives 8.5 ms latency
+OMEParams.ARrateThreshold=0;
+
+%%  #3 DRNL
+DRNLParams=[];  % clear the structure first
+DRNLParams.BFlist=BFlist;
+
+% DRNL nonlinear path
+DRNLParams.a=5e4;     % DRNL.a=0 means no OHCs (no nonlinear path)
+DRNLParams.a=2e4;     % DRNL.a=0 means no OHCs (no nonlinear path)
+
+DRNLParams.b=8e-6;    % *compression threshold raised compression
+% DRNLParams.b=1;    % b=1 means no compression
+
+DRNLParams.c=0.2;     % compression exponent
+% nonlinear filters
+DRNLParams.nonlinCFs=BFlist;
+DRNLParams.nonlinOrder=	3;           % order of nonlinear gammatone filters
+p=0.2895;   q=170;  % human  (% p=0.14;   q=366;  % cat)
+DRNLParams.nlBWs=  p * BFlist + q;
+DRNLParams.p=p;   DRNLParams.q=q;   % save p and q for printing only
+
+% DRNL linear path:
+DRNLParams.g=1000;     % linear path gain factor
+% linCF is not necessarily the same as nonlinCF
+minLinCF=153.13; coeffLinCF=0.7341;   % linCF>nonlinBF for BF < 1 kHz
+DRNLParams.linCFs=minLinCF+coeffLinCF*BFlist;
+DRNLParams.linOrder=	3;           % order of linear gammatone filters
+minLinBW=100; coeffLinBW=0.6531;
+DRNLParams.linBWs=minLinBW + coeffLinBW*BFlist; % bandwidths of linear  filters
+
+% DRNL MOC efferents
+DRNLParams.MOCdelay = efferentDelay;            % must be < segment length!
+
+% 'spikes' model: MOC based on brainstem spiking activity (HSR)
+DRNLParams.rateToAttenuationFactor = .01;  % strength of MOC
+%      DRNLParams.rateToAttenuationFactor = 0;  % strength of MOC
+% 'probability' model: MOC based on AN spiking activity (HSR)
+DRNLParams.rateToAttenuationFactorProb = .0055;  % strength of MOC
+% DRNLParams.rateToAttenuationFactorProb = .0;  % strength of MOC
+DRNLParams.MOCrateThresholdProb =70;                % spikes/s probability only
+
+DRNLParams.MOCtau =.1;                         % smoothing for MOC
+
+
+%% #4 IHC_cilia_RPParams
+
+IHC_cilia_RPParams.tc=	0.0003;   % 0.0003 filter time simulates viscocity
+% IHC_cilia_RPParams.tc=	0.0005;   % 0.0003 filter time simulates viscocity
+IHC_cilia_RPParams.C=	0.03;      % 0.1 scalar (C_cilia ) 
+IHC_cilia_RPParams.u0=	5e-9;       
+IHC_cilia_RPParams.s0=	30e-9;
+IHC_cilia_RPParams.u1=	1e-9;
+IHC_cilia_RPParams.s1=	1e-9;
+
+IHC_cilia_RPParams.Gmax= 6e-9;    % 2.5e-9 maximum conductance (Siemens)
+IHC_cilia_RPParams.Ga=	1e-9;  % 4.3e-9 fixed apical membrane conductance
+IHC_cilia_RPParams.Ga=	.8e-9;  % 4.3e-9 fixed apical membrane conductance
+
+%  #5 IHC_RP
+IHC_cilia_RPParams.Cab=	4e-012;         % IHC capacitance (F)
+% IHC_cilia_RPParams.Cab=	1e-012;         % IHC capacitance (F)
+IHC_cilia_RPParams.Et=	0.05;          % endocochlear potential (V)
+
+IHC_cilia_RPParams.Gk=	2e-008;         % 1e-8 potassium conductance (S)
+IHC_cilia_RPParams.Ek=	-0.08;          % -0.084 K equilibrium potential
+IHC_cilia_RPParams.Rpc=	0.04;           % combined resistances
+
+
+%%  #5 IHCpreSynapse
+IHCpreSynapseParams=[];
+IHCpreSynapseParams.GmaxCa=	14e-9;% maximum calcium conductance
+IHCpreSynapseParams.GmaxCa=	12e-9;% maximum calcium conductance
+IHCpreSynapseParams.ECa=	0.066;  % calcium equilibrium potential
+IHCpreSynapseParams.beta=	400;	% determine Ca channel opening
+IHCpreSynapseParams.gamma=	100;	% determine Ca channel opening
+IHCpreSynapseParams.tauM=	0.00005;	% membrane time constant ?0.1ms
+IHCpreSynapseParams.power=	3;
+% reminder: changing z has a strong effect on HF thresholds (like Et)
+IHCpreSynapseParams.z=	    2e42;   % scalar Ca -> vesicle release rate
+
+LSRtauCa=35e-6;            HSRtauCa=85e-6;            % seconds
+% LSRtauCa=35e-6;            HSRtauCa=70e-6;            % seconds
+IHCpreSynapseParams.tauCa= [LSRtauCa HSRtauCa]; %LSR and HSR fiber
+
+%%  #6 AN_IHCsynapse
+% c=kym/(y(l+r)+kl)	(spontaneous rate)
+% c=(approx)  ym/l  (saturated rate)
+AN_IHCsynapseParams=[];             % clear the structure first
+AN_IHCsynapseParams.M=	12;         % maximum vesicles at synapse
+AN_IHCsynapseParams.y=	4;          % depleted vesicle replacement rate
+AN_IHCsynapseParams.y=	6;          % depleted vesicle replacement rate
+
+AN_IHCsynapseParams.x=	30;         % replenishment from re-uptake store
+AN_IHCsynapseParams.x=	60;         % replenishment from re-uptake store
+
+% reduce l to increase saturated rate
+AN_IHCsynapseParams.l=	100; % *loss rate of vesicles from the cleft
+AN_IHCsynapseParams.l=	250; % *loss rate of vesicles from the cleft
+
+AN_IHCsynapseParams.r=	500; % *reuptake rate from cleft into cell
+% AN_IHCsynapseParams.r=	300; % *reuptake rate from cleft into cell
+
+AN_IHCsynapseParams.refractory_period=	0.00075;
+% number of AN fibers at each BF (used only for spike generation)
+AN_IHCsynapseParams.numFibers=	100; 
+AN_IHCsynapseParams.TWdelay=0.004;  % ?delay before stimulus first spike
+
+AN_IHCsynapseParams.ANspeedUpFactor=5; % longer epochs for computing spikes.
+
+%%  #7 MacGregorMulti (first order brainstem neurons)
+MacGregorMultiParams=[];
+MacGregorMultiType='chopper'; % MacGregorMultiType='primary-like'; %choose
+switch MacGregorMultiType
+    case 'primary-like'
+        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF
+        MacGregorMultiParams.type = 'primary-like cell';
+        MacGregorMultiParams.fibersPerNeuron=4;   % N input fibers
+        MacGregorMultiParams.dendriteLPfreq=200;  % dendritic filter
+        MacGregorMultiParams.currentPerSpike=0.11e-6; % (A) per spike
+        MacGregorMultiParams.Cap=4.55e-9;   % cell capacitance (Siemens)
+        MacGregorMultiParams.tauM=5e-4;     % membrane time constant (s)
+        MacGregorMultiParams.Ek=-0.01;      % K+ eq. potential (V)
+        MacGregorMultiParams.dGkSpike=3.64e-5; % K+ cond.shift on spike,S
+        MacGregorMultiParams.tauGk=	0.0012; % K+ conductance tau (s)
+        MacGregorMultiParams.Th0=	0.01;   % equilibrium threshold (V)
+        MacGregorMultiParams.c=	0.01;       % threshold shift on spike, (V)
+        MacGregorMultiParams.tauTh=	0.015;  % variable threshold tau
+        MacGregorMultiParams.Er=-0.06;      % resting potential (V)
+        MacGregorMultiParams.Eb=0.06;       % spike height (V)
+
+    case 'chopper'
+        MacGregorMultiParams.nNeuronsPerBF=	10;   % N neurons per BF
+        MacGregorMultiParams.type = 'chopper cell';
+        MacGregorMultiParams.fibersPerNeuron=10;  % N input fibers
+%         MacGregorMultiParams.fibersPerNeuron=6;  % N input fibers
+
+        MacGregorMultiParams.dendriteLPfreq=50;   % dendritic filter
+        MacGregorMultiParams.currentPerSpike=35e-9; % *per spike
+%         MacGregorMultiParams.currentPerSpike=30e-9; % *per spike
+        
+        MacGregorMultiParams.Cap=1.67e-8; % ??cell capacitance (Siemens)
+        MacGregorMultiParams.tauM=0.002;  % membrane time constant (s)
+        MacGregorMultiParams.Ek=-0.01;    % K+ eq. potential (V)
+        MacGregorMultiParams.dGkSpike=1.33e-4; % K+ cond.shift on spike,S
+        MacGregorMultiParams.tauGk=	0.0005;% K+ conductance tau (s)
+        MacGregorMultiParams.Th0=	0.01; % equilibrium threshold (V)
+        MacGregorMultiParams.c=	0;        % threshold shift on spike, (V)
+        MacGregorMultiParams.tauTh=	0.02; % variable threshold tau
+        MacGregorMultiParams.Er=-0.06;    % resting potential (V)
+        MacGregorMultiParams.Eb=0.06;     % spike height (V)
+        MacGregorMultiParams.PSTHbinWidth=	1e-4;
+end
+
+%%  #8 MacGregor (second-order neuron). Only one per channel
+MacGregorParams=[];                 % clear the structure first
+MacGregorParams.type = 'chopper cell';
+MacGregorParams.fibersPerNeuron=10; % N input fibers
+MacGregorParams.dendriteLPfreq=100; % dendritic filter
+MacGregorParams.currentPerSpike=120e-9;% *(A) per spike
+MacGregorParams.currentPerSpike=40e-9;% *(A) per spike
+
+MacGregorParams.Cap=16.7e-9;        % cell capacitance (Siemens)
+MacGregorParams.tauM=0.002;         % membrane time constant (s)
+MacGregorParams.Ek=-0.01;           % K+ eq. potential (V)
+MacGregorParams.dGkSpike=1.33e-4;   % K+ cond.shift on spike,S
+MacGregorParams.tauGk=	0.0005;     % K+ conductance tau (s)
+MacGregorParams.Th0=	0.01;       % equilibrium threshold (V)
+MacGregorParams.c=	0;              % threshold shift on spike, (V)
+MacGregorParams.tauTh=	0.02;       % variable threshold tau
+MacGregorParams.Er=-0.06;           % resting potential (V)
+MacGregorParams.Eb=0.06;            % spike height (V)
+MacGregorParams.debugging=0;        % (special)
+% wideband accepts input from all channels (of same fiber type)
+% use wideband to create inhibitory units
+MacGregorParams.wideband=0;         % special for wideband units
+% MacGregorParams.saveAllData=0;
+
+%%  #9 filteredSACF
+minPitch=	300; maxPitch=	3000; numPitches=60;    % specify lags
+pitches=100*log10(logspace(minPitch/100, maxPitch/100, numPitches));
+filteredSACFParams.lags=1./pitches;     % autocorrelation lags vector
+filteredSACFParams.acfTau=	.003;       % time constant of running ACF
+filteredSACFParams.lambda=	0.12;       % slower filter to smooth ACF
+filteredSACFParams.plotFilteredSACF=1;  % 0 plots unfiltered ACFs
+filteredSACFParams.plotACFs=0;          % special plot (see code)
+%  filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags
+filteredSACFParams.lagsProcedure=  'useAllLags';
+% filteredSACFParams.lagsProcedure=  'useBernsteinLagWeights';
+% filteredSACFParams.lagsProcedure=  'omitShortLags';
+filteredSACFParams.criterionForOmittingLags=3;
+
+% checks
+if AN_IHCsynapseParams.numFibers<MacGregorMultiParams.fibersPerNeuron
+    error('MacGregorMulti: too few input fibers for input to MacG unit')
+end
+
+
+%% now accept last minute parameter changes required by the calling program
+% paramChanges
+if nargin>3 && ~isempty(paramChanges)
+    nChanges=length(paramChanges);
+    for idx=1:nChanges
+        eval(paramChanges{idx})
+    end
+end
+
+
+%% write all parameters to the command window
+% showParams is currently set at the top of htis function
+if showParams
+    fprintf('\n %%%%%%%%\n')
+    fprintf('\n%s\n', method.parameterSource)
+    fprintf('\n')
+    nm=UTIL_paramsList(whos);
+    for i=1:length(nm)
+        %         eval(['UTIL_showStruct(' nm{i} ', ''' nm{i} ''')'])
+        if ~strcmp(nm(i), 'method')
+            eval(['UTIL_showStructureSummary(' nm{i} ', ''' nm{i} ''', 10)'])
+        end
+    end
+
+    % highlight parameter changes made locally
+    if nargin>3 && ~isempty(paramChanges)
+        fprintf('\n Local parameter changes:\n')
+        for i=1:length(paramChanges)
+            disp(paramChanges{i})
+        end
+    end
+end
+
+% for backward compatibility
+experiment.comparisonData=[];
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testAN.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testAN.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,333 @@
+function vectorStrength=testAN(targetFrequency,BFlist, levels, ...
+    paramsName,paramChanges)
+% testIHC used either for IHC I/O function ...
+%  or receptive field (doReceptiveFields=1)
+
+global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams
+global AN_IHCsynapseParams
+
+    global ANoutput ANdt CNoutput ICoutput ICmembraneOutput tauCas
+    global ARattenuation MOCattenuation
+
+dbstop if error
+restorePath=path;
+
+addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
+    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
+    ['..' filesep 'testPrograms'])
+
+if nargin<5
+    paramChanges=[];
+end
+
+if nargin<4
+    paramsName='Normal';
+end
+
+if nargin<3
+    levels=-10:10:80;
+    % levels=80:10:90;
+end
+
+nLevels=length(levels);
+
+toneDuration=.2;
+rampDuration=0.002;
+silenceDuration=.02;
+localPSTHbinwidth=0.001;
+
+% Use only the first frequency in the GUI targetFrequency box to defineBF
+% targetFrequency=stimulusParameters.targetFrequency(1);
+% BFlist=targetFrequency;
+
+AN_HSRonset=zeros(nLevels,1);
+AN_HSRsaturated=zeros(nLevels,1);
+AN_LSRonset=zeros(nLevels,1);
+AN_LSRsaturated=zeros(nLevels,1);
+CNLSRrate=zeros(nLevels,1);
+CNHSRsaturated=zeros(nLevels,1);
+ICHSRsaturated=zeros(nLevels,1);
+ICLSRsaturated=zeros(nLevels,1);
+vectorStrength=zeros(nLevels,1);
+
+AR=zeros(nLevels,1);
+MOC=zeros(nLevels,1);
+
+% ANoutput=zeros(200,200);
+
+figure(15), clf
+set(gcf,'position',[980   356   401   321])
+figure(5), clf
+set(gcf,'position', [980 34 400 295])
+drawnow
+
+%% guarantee that the sample rate is at least 10 times the frequency
+sampleRate=50000;
+while sampleRate< 10* targetFrequency
+    sampleRate=sampleRate+10000;
+end
+
+%% adjust sample rate so that the pure tone stimulus has an integer
+%% numver of epochs in a period
+dt=1/sampleRate;
+period=1/targetFrequency;
+ANspeedUpFactor=5;  %anticipating MAP (needs clearing up)
+ANdt=dt*ANspeedUpFactor;
+ANdt=period/round(period/ANdt);
+dt=ANdt/ANspeedUpFactor;
+sampleRate=1/dt;
+epochsPerPeriod=sampleRate*period;
+
+%% main computational loop (vary level)    
+levelNo=0;
+for leveldB=levels
+    levelNo=levelNo+1;
+
+    fprintf('%4.0f\t', leveldB)
+    amp=28e-6*10^(leveldB/20);
+
+    time=dt:dt:toneDuration;
+    rampTime=dt:dt:rampDuration;
+    ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+        ones(1,length(time)-length(rampTime))];
+    ramp=ramp.*fliplr(ramp);
+
+    silence=zeros(1,round(silenceDuration/dt));
+
+    % create signal (leveldB/ targetFrequency)
+    inputSignal=amp*sin(2*pi*targetFrequency'*time);
+    inputSignal= ramp.*inputSignal;
+    inputSignal=[silence inputSignal];
+
+    %% run the model
+    AN_spikesOrProbability='spikes';
+    showPlotsAndDetails=0;
+
+
+    MAP1_14(inputSignal, 1/dt, BFlist, ...
+        paramsName, AN_spikesOrProbability, paramChanges);
+
+    nTaus=length(tauCas);
+
+    %LSR (same as HSR if no LSR fibers present)
+    [nANFibers nTimePoints]=size(ANoutput);
+
+    numLSRfibers=nANFibers/nTaus;
+    numHSRfibers=numLSRfibers;
+
+    LSRspikes=ANoutput(1:numLSRfibers,:);
+    PSTH=UTIL_PSTHmaker(LSRspikes, ANdt, localPSTHbinwidth);
+    PSTHLSR=mean(PSTH,1)/localPSTHbinwidth;  % across fibers rates
+    PSTHtime=localPSTHbinwidth:localPSTHbinwidth:...
+        localPSTHbinwidth*length(PSTH);
+    AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window
+    AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end));
+
+    % HSR
+    HSRspikes= ANoutput(end- numHSRfibers+1:end, :);
+    PSTH=UTIL_PSTHmaker(HSRspikes, ANdt, localPSTHbinwidth);
+    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
+    AN_HSRonset(levelNo)= max(PSTH);
+    AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end));
+
+    figure(5), subplot(2,2,2)
+    hold off, bar(PSTHtime,PSTH, 'b')
+    hold on,  bar(PSTHtime,PSTHLSR,'r')
+    ylim([0 1000])
+    xlim([0 length(PSTH)*localPSTHbinwidth])
+    set(gcf,'name',[num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);
+
+    % AN - CV
+    %  CV is computed 5 times. Use the middle one (3) as most typical
+    cvANHSR=  UTIL_CV(HSRspikes, ANdt);
+
+    % AN - vector strength
+    PSTH=sum(HSRspikes);
+    [PH, binTimes]=UTIL_periodHistogram...
+        (PSTH, ANdt, targetFrequency);
+    VS=UTIL_vectorStrength(PH);
+    vectorStrength(levelNo)=VS;
+    disp(['sat rate= ' num2str(AN_HSRsaturated(levelNo)) ...
+        ';   phase-locking VS = ' num2str(VS)])
+    title(['AN HSR: CV=' num2str(cvANHSR(3),'%5.2f') ...
+        'VS=' num2str(VS,'%5.2f')])
+
+    % CN - first-order neurons
+
+    % CN LSR
+    [nCNneurons c]=size(CNoutput);
+    nLSRneurons=round(nCNneurons/nTaus);
+    CNLSRspikes=CNoutput(1:nLSRneurons,:);
+    PSTH=UTIL_PSTHmaker(CNLSRspikes, ANdt, localPSTHbinwidth);
+    PSTH=sum(PSTH)/nLSRneurons;
+    CNLSRrate(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;
+
+    %CN HSR
+    MacGregorMultiHSRspikes=...
+        CNoutput(end-nLSRneurons:end,:);
+    PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth);
+    PSTH=sum(PSTH)/nLSRneurons;
+    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
+
+    CNHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));
+
+    figure(5), subplot(2,2,3)
+    bar(PSTHtime,PSTH)
+    ylim([0 1000])
+    xlim([0 length(PSTH)*localPSTHbinwidth])
+    cvMMHSR= UTIL_CV(MacGregorMultiHSRspikes, ANdt);
+    title(['CN    CV= ' num2str(cvMMHSR(3),'%5.2f')])
+
+    % IC LSR
+    [nICneurons c]=size(ICoutput);
+    nLSRneurons=round(nICneurons/nTaus);
+    ICLSRspikes=ICoutput(1:nLSRneurons,:);
+    PSTH=UTIL_PSTHmaker(ICLSRspikes, ANdt, localPSTHbinwidth);
+    ICLSRsaturated(levelNo)=mean(PSTH(round(length(PSTH)/2):end))/localPSTHbinwidth;
+
+    %IC HSR
+    MacGregorMultiHSRspikes=...
+        ICoutput(end-nLSRneurons:end,:);
+    PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, ANdt, localPSTHbinwidth);
+    PSTH=sum(PSTH)/nLSRneurons;
+    PSTH=mean(PSTH,1)/localPSTHbinwidth; % sum across fibers (HSR only)
+
+    ICHSRsaturated(levelNo)=mean(PSTH(length(PSTH)/2:end));
+
+    AR(levelNo)=min(ARattenuation);
+    MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));
+
+    time=dt:dt:dt*size(ICmembraneOutput,2);
+    figure(5), subplot(2,2,4)
+    plot(time,ICmembraneOutput(2, 1:end),'k')
+    ylim([-0.07 0])
+    xlim([0 max(time)])
+    title(['IC  ' num2str(leveldB,'%4.0f') 'dB'])
+    drawnow
+
+    figure(5), subplot(2,2,1)
+    plot(20*log10(MOC), 'k'),
+    title(' MOC'), ylabel('dB attenuation')
+    ylim([-30 0])
+
+
+end % level
+figure(5), subplot(2,2,1)
+plot(levels,20*log10(MOC), 'k'),
+title(' MOC'), ylabel('dB attenuation')
+ylim([-30 0])
+xlim([0 max(levels)])
+
+fprintf('\n')
+toneDuration=2;
+rampDuration=0.004;
+silenceDuration=.02;
+nRepeats=200;   % no. of AN fibers
+fprintf('toneDuration  %6.3f\n', toneDuration)
+fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)
+fprintf('levels')
+fprintf('%6.2f\t', levels)
+fprintf('\n')
+
+
+% ---------------------------------------------------- display parameters
+
+figure(15), clf
+nRows=2; nCols=2;
+
+% AN rate - level ONSET functions
+subplot(nRows,nCols,1)
+plot(levels,AN_LSRonset,'ro'), hold on
+plot(levels,AN_HSRonset,'ko'), hold off
+ylim([0 1000]),  xlim([min(levels) max(levels)])
+ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];
+title( ttl)
+xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on
+text(0, 800, 'AN onset', 'fontsize', 16)
+
+% AN rate - level ADAPTED function
+subplot(nRows,nCols,2)
+plot(levels,AN_LSRsaturated, 'ro'), hold on
+plot(levels,AN_HSRsaturated, 'ko'), hold off
+ylim([0 400])
+set(gca,'ytick',0:50:300)
+xlim([min(levels) max(levels)])
+set(gca,'xtick',[levels(1):20:levels(end)])
+%     grid on
+ttl=[   'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')...
+    '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];
+title( ttl)
+xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
+text(0, 340, 'AN adapted', 'fontsize', 16), grid on
+
+% CN rate - level ADAPTED function
+subplot(nRows,nCols,3)
+plot(levels,CNLSRrate, 'ro'), hold on
+plot(levels,CNHSRsaturated, 'ko'), hold off
+ylim([0 400])
+set(gca,'ytick',0:50:300)
+xlim([min(levels) max(levels)])
+set(gca,'xtick',[levels(1):20:levels(end)])
+%     grid on
+ttl=[   'spont=' num2str(mean(CNHSRsaturated(1,:)),'%4.0f') '  sat=' ...
+    num2str(mean(CNHSRsaturated(end,1)),'%4.0f')];
+title( ttl)
+xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
+text(0, 350, 'CN', 'fontsize', 16), grid on
+
+% IC rate - level ADAPTED function
+subplot(nRows,nCols,4)
+plot(levels,ICLSRsaturated, 'ro'), hold on
+plot(levels,ICHSRsaturated, 'ko'), hold off
+ylim([0 400])
+set(gca,'ytick',0:50:300)
+xlim([min(levels) max(levels)])
+set(gca,'xtick',[levels(1):20:levels(end)]), grid on
+
+
+ttl=['spont=' num2str(mean(ICHSRsaturated(1,:)),'%4.0f') ...
+    '  sat=' num2str(mean(ICHSRsaturated(end,1)),'%4.0f')];
+title( ttl)
+xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
+text(0, 350, 'IC', 'fontsize', 16)
+set(gcf,'name',' AN CN IC rate/level')
+
+peakVectorStrength=max(vectorStrength);
+
+fprintf('\n')
+disp('levels vectorStrength')
+fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
+fprintf('\n')
+fprintf('Phase locking, max vector strength=\t %6.4f\n\n',...
+    max(vectorStrength))
+
+allData=[ levels'  AN_HSRonset AN_HSRsaturated...
+    AN_LSRonset AN_LSRsaturated ...
+    CNHSRsaturated CNLSRrate...
+    ICHSRsaturated ICLSRsaturated];
+fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');
+UTIL_printTabTable(round(allData))
+fprintf('VS (phase locking)= \t%6.4f\n\n',...
+    max(vectorStrength))
+
+UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')
+UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
+UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
+
+fprintf('\n')
+disp('levels vectorStrength')
+fprintf('%3.0f \t %6.4f \n', [levels; vectorStrength'])
+fprintf('\n')
+fprintf('Phase locking, max vector strength= \t%6.4f\n\n',...
+    max(vectorStrength))
+
+allData=[ levels'  AN_HSRonset AN_HSRsaturated...
+    AN_LSRonset AN_LSRsaturated ...
+    CNHSRsaturated CNLSRrate...
+    ICHSRsaturated ICLSRsaturated];
+fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');
+UTIL_printTabTable(round(allData))
+fprintf('VS (phase locking)= \t%6.4f\n\n',...
+    max(vectorStrength))
+
+path(restorePath)
\ No newline at end of file
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testANprob.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testANprob.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,199 @@
+function vectorStrength=testANprob(targetFrequency,BFlist, levels, ...
+    paramsName, paramChanges)
+% testIHC used either for IHC I/O function ...
+%  or receptive field (doReceptiveFields=1)
+
+global IHC_VResp_VivoParams  IHC_cilia_RPParams IHCpreSynapseParams
+global AN_IHCsynapseParams
+
+global ANprobRateOutput dt tauCas
+global ARattenuation MOCattenuation
+
+AN_spikesOrProbability='probability';
+
+dbstop if error
+addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
+    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
+    ['..' filesep 'testPrograms'])
+
+if nargin<5
+    paramChanges=[];
+end
+
+if nargin<4
+    paramsName='Normal';
+end
+
+if nargin<3
+    levels=-10:10:80;
+end
+
+nLevels=length(levels);
+
+toneDuration=.2;
+rampDuration=0.002;
+silenceDuration=.02;
+localPSTHbinwidth=0.001;
+
+% Use only the first frequency in the GUI targetFrequency box to defineBF
+% targetFrequency=stimulusParameters.targetFrequency(1);
+% BFlist=targetFrequency;
+
+AN_HSRonset=zeros(nLevels,1);
+AN_HSRsaturated=zeros(nLevels,1);
+AN_LSRonset=zeros(nLevels,1);
+AN_LSRsaturated=zeros(nLevels,1);
+
+AR=zeros(nLevels,1);
+MOC=zeros(nLevels,1);
+
+figure(15), clf
+set(gcf,'position',[607    17   368   321])
+drawnow
+
+%% guarantee that the sample rate is at least 10 times the frequency
+sampleRate=50000;
+while sampleRate< 10* targetFrequency
+    sampleRate=sampleRate+10000;
+end
+
+%% adjust sample rate so that the pure tone stimulus has an integer
+%% numver of epochs in a period
+dt=1/sampleRate;
+period=1/targetFrequency;
+
+%% main computational loop (vary level)
+levelNo=0;
+for leveldB=levels
+    levelNo=levelNo+1;
+
+    fprintf('%4.0f\t', leveldB)
+    amp=28e-6*10^(leveldB/20);
+
+    time=dt:dt:toneDuration;
+    rampTime=dt:dt:rampDuration;
+    ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+        ones(1,length(time)-length(rampTime))];
+    ramp=ramp.*fliplr(ramp);
+
+    silence=zeros(1,round(silenceDuration/dt));
+
+    % create signal (leveldB/ targetFrequency)
+    inputSignal=amp*sin(2*pi*targetFrequency'*time);
+    inputSignal= ramp.*inputSignal;
+    inputSignal=[silence inputSignal];
+
+    %% run the model
+    showPlotsAndDetails=0;
+
+
+    MAP1_14(inputSignal, 1/dt, BFlist, ...
+        paramsName, AN_spikesOrProbability, paramChanges);
+
+    nTaus=length(tauCas);
+
+    %LSR (same as HSR if no LSR fibers present)
+    [nANFibers nTimePoints]=size(ANprobRateOutput);
+
+    numLSRfibers=1;
+    numHSRfibers=numLSRfibers;
+
+    LSRspikes=ANprobRateOutput(1:numLSRfibers,:);
+    PSTH=UTIL_PSTHmaker(LSRspikes, dt, localPSTHbinwidth);
+    PSTHLSR=PSTH/(localPSTHbinwidth/dt);  % across fibers rates
+    PSTHtime=localPSTHbinwidth:localPSTHbinwidth:...
+        localPSTHbinwidth*length(PSTH);
+    AN_LSRonset(levelNo)= max(PSTHLSR); % peak in 5 ms window
+    AN_LSRsaturated(levelNo)= mean(PSTHLSR(round(length(PSTH)/2):end));
+
+    % HSR
+    HSRspikes= ANprobRateOutput(end- numHSRfibers+1:end, :);
+    PSTH=UTIL_PSTHmaker(HSRspikes, dt, localPSTHbinwidth);
+    PSTH=PSTH/(localPSTHbinwidth/dt); % sum across fibers (HSR only)
+    AN_HSRonset(levelNo)= max(PSTH);
+    AN_HSRsaturated(levelNo)= mean(PSTH(round(length(PSTH)/2): end));
+
+    figure(15), subplot(2,2,4)
+    hold off, bar(PSTHtime,PSTH, 'b')
+    hold on,  bar(PSTHtime,PSTHLSR,'r')
+    ylim([0 1000])
+    xlim([0 length(PSTH)*localPSTHbinwidth])
+    set(gcf,'name',[num2str(BFlist), ' Hz: ' num2str(leveldB) ' dB']);
+
+    AR(levelNo)=min(ARattenuation);
+    MOC(levelNo)=min(MOCattenuation(length(MOCattenuation)/2:end));
+
+
+    figure(15), subplot(2,2,3)
+    plot(20*log10(MOC), 'k'),
+    title(' MOC'), ylabel('dB attenuation')
+    ylim([-30 0])
+
+
+end % level
+figure(15), subplot(2,2,3)
+plot(levels,20*log10(MOC), 'k'),
+title(' MOC'), ylabel('dB attenuation')
+ylim([-30 0])
+xlim([0 max(levels)])
+
+fprintf('\n')
+toneDuration=2;
+rampDuration=0.004;
+silenceDuration=.02;
+nRepeats=200;   % no. of AN fibers
+fprintf('toneDuration  %6.3f\n', toneDuration)
+fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)
+fprintf('levels')
+fprintf('%6.2f\t', levels)
+fprintf('\n')
+
+
+% ---------------------------------------------------- display parameters
+
+
+nRows=2; nCols=2;
+
+% AN rate - level ONSET functions
+subplot(nRows,nCols,1)
+plot(levels,AN_LSRonset,'ro'), hold on
+plot(levels,AN_HSRonset,'ko'), hold off
+ylim([0 1000]),  xlim([min(levels) max(levels)])
+ttl=['tauCa= ' num2str(IHCpreSynapseParams.tauCa)];
+title( ttl)
+xlabel('level dB SPL'), ylabel('peak rate (sp/s)'), grid on
+text(0, 800, 'AN onset', 'fontsize', 16)
+
+% AN rate - level ADAPTED function
+subplot(nRows,nCols,2)
+plot(levels,AN_LSRsaturated, 'ro'), hold on
+plot(levels,AN_HSRsaturated, 'ko'), hold off
+ylim([0 400])
+set(gca,'ytick',0:50:300)
+xlim([min(levels) max(levels)])
+set(gca,'xtick',[levels(1):20:levels(end)])
+%     grid on
+ttl=[   'spont=' num2str(mean(AN_HSRsaturated(1,:)),'%4.0f')...
+    '  sat=' num2str(mean(AN_HSRsaturated(end,1)),'%4.0f')];
+title( ttl)
+xlabel('level dB SPL'), ylabel ('adapted rate (sp/s)')
+text(0, 340, 'AN adapted', 'fontsize', 16), grid on
+
+allData=[ levels'  AN_HSRonset AN_HSRsaturated...
+    AN_LSRonset AN_LSRsaturated ];
+fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');
+UTIL_printTabTable(round(allData))
+
+
+UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')
+UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
+UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
+
+fprintf('\n')
+disp('levels vectorStrength')
+
+allData=[ levels'  AN_HSRonset AN_HSRsaturated...
+    AN_LSRonset AN_LSRsaturated ];
+fprintf('\n levels \tANHSR Onset \tANHSR adapted\tANLSR Onset \tANLSR adapted\tCNHSR\tCNLSR\tICHSR  \tICLSR \n');
+UTIL_printTabTable(round(allData))
+
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testBM.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testBM.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,164 @@
+function testBM (BMlocations, paramsName,...
+    relativeFrequencies, AN_spikesOrProbability, paramChanges)
+% testBM generates input output functions for DRNL model for any number
+% of locations.
+% Computations are bast on a single channel model (channelBFs=BF)
+% peak displacement (peakAmp) is measured.
+%  if DRNLParams.useMOC is chosen, the full model is run (slow)
+%  otherwise only DRNL is computed.
+% Tuning curves are generated based on a range of frequencies reletove to
+% the BF of the location.
+%
+
+global    DRNLParams
+
+if nargin<5
+    paramChanges=[];
+end
+
+if nargin<4
+    AN_spikesOrProbability='probability';
+end
+
+savePath=path;
+addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
+
+levels=-10:10:90;   nLevels=length(levels);
+% levels= 50;   nLevels=length(levels);
+
+% refBMdisplacement is the displacement of the BM at threshold
+% 1 nm disp at  threshold (9 kHz, Ruggero)
+% ? adjust for frequency
+refBMdisplacement= 1e-8; % adjusted for 10 nm at 1 kHz 
+
+toneDuration=.200;
+rampDuration=0.01;
+silenceDuration=0.01;
+
+sampleRate=30000;
+
+dbstop if error
+figure(3), clf
+set(gcf,'position',[280   350   327   326])
+set(gcf,'name','DRNL - BM')
+pause(0.1)
+
+finalSummary=[];
+nBFs=length(BMlocations);
+BFno=0; plotCount=0;
+for BF=BMlocations
+    BFno=BFno+1;
+    plotCount=plotCount+nBFs;
+    stimulusFrequencies=BF* relativeFrequencies;
+    nFrequencies=length(stimulusFrequencies);
+
+    peakAmpBM=zeros(nLevels,nFrequencies);
+    peakAmpBMdB=NaN(nLevels,nFrequencies);
+    peakEfferent=NaN(nLevels,nFrequencies);
+    peakAREfferent=NaN(nLevels,nFrequencies);
+
+
+    levelNo=0;
+    for leveldB=levels
+        disp(['level= ' num2str(leveldB)])
+        levelNo=levelNo+1;
+
+        freqNo=0;
+        for frequency=stimulusFrequencies
+            freqNo=freqNo+1;
+
+            % Generate stimuli
+            globalStimParams.FS=sampleRate;
+            globalStimParams.overallDuration=...
+                toneDuration+silenceDuration;  % s
+            stim.phases='sin';
+            stim.type='tone';
+            stim.toneDuration=toneDuration;
+            stim.frequencies=frequency;
+            stim.amplitudesdB=leveldB;
+            stim.beginSilence=silenceDuration;
+            stim.rampOnDur=rampDuration;
+            % no offset ramp
+            stim.rampOffDur=rampDuration;
+            doPlot=0;
+            inputSignal=stimulusCreate(globalStimParams, stim, doPlot);
+            inputSignal=inputSignal(:,1)';
+
+            %% run the model
+            MAPparamsName=paramsName;
+
+            global DRNLoutput MOCattenuation ARattenuation
+            MAP1_14(inputSignal, sampleRate, BF, ...
+                MAPparamsName, AN_spikesOrProbability, paramChanges);
+
+            DRNLresponse=DRNLoutput;
+            peakAmp=max(max(...
+                DRNLresponse(:, end-round(length(DRNLresponse)/2):end)));
+            peakAmpBM(levelNo,freqNo)=peakAmp;
+            peakAmpBMdB(levelNo,freqNo)=...
+                20*log10(peakAmp/refBMdisplacement);
+            peakEfferent(levelNo,freqNo)=min(min(MOCattenuation));
+            peakAREfferent(levelNo,freqNo)=min(min(ARattenuation));
+
+        end  % tone frequency
+    end  % level
+
+    %% analyses results and plot
+
+    % BM I/O plot (top panel)
+    figure(3)
+    subplot(3,nBFs,BFno), cla
+    plot(levels,peakAmpBMdB, 'linewidth',2)
+    hold on, plot(levels, repmat(refBMdisplacement,1,length(levels)))
+    hold off
+    title(['BF=' num2str(BF,'%5.0f') ' - ' paramsName])
+    xlabel('level')
+    %     set(gca,'xtick',levels),  grid on
+    if length(levels)>1,xlim([min(levels) max(levels)]), end
+    ylabel(['dB re:' num2str(refBMdisplacement,'%6.1e') 'm'])
+    ylim([-20 50])
+    set(gca,'ytick',[-10 0 10 20 40])
+    grid on
+    %     legend({num2str(stimulusFrequencies')}, 'location', 'EastOutside')
+    UTIL_printTabTable([levels' peakAmpBMdB], ...
+        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')
+    finalSummary=[finalSummary peakAmpBMdB];
+
+    % Tuning curve
+    if length(relativeFrequencies)>2
+        figure(3), subplot(3,nBFs, 2*nBFs+BFno)
+        %         contour(stimulusFrequencies,levels,peakAmpBM,...
+        %             [refBMdisplacement refBMdisplacement],'r')
+        contour(stimulusFrequencies,levels,peakAmpBM,...
+            refBMdisplacement.*[1 5 10 50 100])
+        title(['tuning curve at ' num2str(refBMdisplacement) 'm']);
+        ylabel('level (dB) at reference')
+        xlim([100 10000])
+        grid on
+        hold on
+        set(gca,'xscale','log')
+    end
+
+
+    % MOC contribution
+    figure(3)
+    subplot(3,nBFs,nBFs+BFno), cla
+    plot(levels,20*log10(peakEfferent), 'linewidth',2)
+    ylabel('MOC (dB attenuation)'), xlabel('level')
+    title(['peak MOC: model= ' AN_spikesOrProbability])
+    grid on
+    if length(levels)>1, xlim([min(levels) max(levels)]), end
+
+    % AR contribution
+    hold on
+    plot(levels,20*log10(peakAREfferent), 'r')
+    hold off
+
+end % best frequency
+
+UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10)
+
+    UTIL_printTabTable([levels' finalSummary], ...
+        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')
+
+path(savePath);
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testFM.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testFM.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,255 @@
+function testFM(BFlist,paramsName,showPSTHs,paramChanges)
+%   specify whether you want AN 'probability' or 'spikes'
+%       spikes is more realistic but takes longer
+%       refractory effect is included only for spikes
+%
+
+% specify the AN ANresponse bin width (normally 1 ms).
+%      This can influence the measure of the AN onset rate based on the
+%      largest bin in the ANresponse
+%
+% Demonstration is based on Harris and Dallos
+
+global inputStimulusParams outerMiddleEarParams DRNLParams 
+global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
+
+dbstop if error
+
+restorePath=path;
+addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
+    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
+    ['..' filesep 'testPrograms'])
+
+if nargin<3
+    paramChanges=[];
+end
+
+% masker and probe levels are relative to this threshold
+thresholdAtCF=10; % dB SPL
+thresholdAtCF=5; % dB SPL
+
+if nargin<1, showPSTHs=1;end
+
+sampleRate=50000;
+
+% fetch BF from GUI: use only the first target frequency
+maskerFrequency=BFlist;
+maskerDuration=.1;
+
+targetFrequency=maskerFrequency;
+probeLeveldB=20+thresholdAtCF;	% H&D use 20 dB SL/ TMC uses 10 dB SL
+probeDuration=0.008; % HD use 15 ms probe (fig 3).
+probeDuration=0.015; % HD use 15 ms probe (fig 3).
+
+rampDuration=.001;  % HD use 1 ms linear ramp
+initialSilenceDuration=0.02;
+finalSilenceDuration=0.05;      % useful for showing recovery
+
+maskerLevels=[-80   10 20 30 40 60 ] + thresholdAtCF;
+% maskerLevels=[-80   40 60 ] + thresholdAtCF;
+
+dt=1/sampleRate;
+
+figure(7), clf
+set(gcf,'position',[613    36   360   247])
+set(gcf,'name', ['forward masking: thresholdAtCF=' num2str(thresholdAtCF)])
+
+if showPSTHs
+    figure(8), clf
+    set(gcf,'name', 'Harris and Dallos simulation')
+    set(gcf,'position',[980    36   380   249])
+end
+
+% Plot Harris and Dallos result for comparison
+gapDurations=[0.001	0.002	0.005	0.01	0.02	0.05	0.1	0.3];
+HDmaskerLevels=[+10	+20	+30	+40	+60];
+HDresponse=[
+    1 1 1 1 1 1 1 1;
+    0.8  	0.82	0.81	0.83	0.87	0.95	1	    1;
+    0.48	0.5	    0.54	0.58	0.7	    0.85	0.95	1;
+    0.3	    0.31	0.35	0.4	    0.5	    0.68	0.82	0.94;
+    0.2 	0.27	0.27	0.29	0.42	0.64	0.75	0.92;
+    0.15	0.17	0.18	0.23	0.3	     0.5	0.6	    0.82];
+
+figure(7)
+semilogx(gapDurations,HDresponse,'o'), hold on
+legend(strvcat(num2str(maskerLevels')),'location','southeast')
+title([ 'masker dB: ' num2str(HDmaskerLevels)])
+
+%% Run the trials
+maxProbeResponse=0;
+levelNo=0;
+resultsMatrix=zeros(length(maskerLevels), length(gapDurations));
+summaryFiringRates=[];
+
+% initial silence
+time=dt: dt: initialSilenceDuration;
+initialSilence=zeros(1,length(time));
+
+% probe
+time=dt: dt: probeDuration;
+amp=28e-6*10^(probeLeveldB/20);
+probe=amp*sin(2*pi.*targetFrequency*time);
+% ramps
+% catch rampTime error
+if rampDuration>0.5*probeDuration, rampDuration=probeDuration/2; end
+rampTime=dt:dt:rampDuration;
+% raised cosine ramp
+ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+    ones(1,length(time)-length(rampTime))];
+%  onset ramp
+probe=probe.*ramp;
+%  offset ramp makes complete ramp for probe
+ramp=fliplr(ramp);
+% apply ramp mask to probe. Probe does not change below
+probe=probe.*ramp;
+
+% final silence
+time=dt: dt: finalSilenceDuration;
+finalSilence=zeros(1,length(time));
+
+PSTHplotCount=0;
+longestSignalDuration=initialSilenceDuration + maskerDuration +...
+    max(gapDurations) + probeDuration + finalSilenceDuration ;
+for maskerLeveldB=maskerLevels
+    levelNo=levelNo+1;
+    allDurations=[];
+    allFiringRates=[];
+
+    %masker
+    time=dt: dt: maskerDuration;
+    masker=sin(2*pi.*maskerFrequency*time);
+    % masker ramps
+    if rampDuration>0.5*maskerDuration
+        % catch ramp duration error
+        rampDuration=maskerDuration/2;
+    end
+    % NB masker ramp (not probe ramp)
+    rampTime=dt:dt:rampDuration;
+    % raised cosine ramp
+    ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi))...
+        ones(1,length(time)-length(rampTime))];
+    %  onset ramp
+    masker=masker.*ramp;
+    %  offset ramp
+    ramp=fliplr(ramp);
+    % apply ramp
+    masker=masker.*ramp;
+    amp=28e-6*10^(maskerLeveldB/20);
+    maskerPa=amp*masker;
+
+    for gapDuration=gapDurations
+        time=dt: dt: gapDuration;
+        gap=zeros(1,length(time));
+
+        inputSignal=...
+            [initialSilence maskerPa gap probe finalSilence];
+
+        % **********************************  run MAP model
+        
+        global  ANprobRateOutput  tauCas  ...
+
+    showPlotsAndDetails=0;
+
+AN_spikesOrProbability='probability';
+
+MAP1_14(inputSignal, 1/dt, targetFrequency, ...
+    paramsName, AN_spikesOrProbability, paramChanges);
+ 
+    [nFibers c]=size(ANprobRateOutput);
+    nLSRfibers=nFibers/length(tauCas);
+            ANresponse=ANprobRateOutput(end-nLSRfibers:end,:);
+            ANresponse=sum(ANresponse)/nLSRfibers;
+    
+        % analyse results
+        probeStart=initialSilenceDuration+maskerDuration+gapDuration;
+        PSTHbinWidth=dt;
+        responseDelay=0.005;
+        probeTimes=probeStart+responseDelay:...
+            PSTHbinWidth:probeStart+probeDuration+responseDelay;
+        probeIDX=round(probeTimes/PSTHbinWidth);
+        probePSTH=ANresponse(probeIDX);
+        firingRate=mean(probePSTH);
+        % NB this only works if you start with the lowest level masker
+        maxProbeResponse=max([maxProbeResponse firingRate]);
+        allDurations=[allDurations gapDuration];
+        allFiringRates=[allFiringRates firingRate];
+        
+        %% show PSTHs
+        if showPSTHs
+            nLevels=length(maskerLevels);
+            nDurations=length(gapDurations);
+            figure(8)
+            PSTHbinWidth=1e-3;
+            PSTH=UTIL_PSTHmaker(ANresponse, dt, PSTHbinWidth);
+            PSTH=PSTH*dt/PSTHbinWidth;
+            PSTHplotCount=PSTHplotCount+1;
+            subplot(nLevels,nDurations,PSTHplotCount)
+            probeTime=PSTHbinWidth:PSTHbinWidth:...
+                PSTHbinWidth*length(PSTH);
+            bar(probeTime, PSTH)
+            if strcmp(AN_spikesOrProbability, 'spikes')
+                ylim([0 500])
+            else
+                ylim([0 500])
+            end
+            xlim([0 longestSignalDuration])
+            grid on
+
+            if PSTHplotCount< (nLevels-1) * nDurations+1
+                set(gca,'xticklabel',[])
+            end
+
+            if ~isequal(mod(PSTHplotCount,nDurations),1)
+                set(gca,'yticklabel',[])
+            else
+                ylabel([num2str(maskerLevels...
+                    (round(PSTHplotCount/nDurations) +1))])
+            end
+
+            if PSTHplotCount<=nDurations
+                title([num2str(1000*gapDurations(PSTHplotCount)) 'ms'])
+            end
+        end % showPSTHs
+
+    end     % gapDurations duration
+    summaryFiringRates=[summaryFiringRates allFiringRates'];
+
+    figure(7), hold on
+    semilogx(allDurations, allFiringRates/maxProbeResponse)
+    ylim([0 1]), hold on
+    ylim([0 inf]), xlim([min(gapDurations) max(gapDurations)])
+    xlim([0.001 1])
+    pause(0.1) % to allow for CTRL/C interrupts
+    resultsMatrix(levelNo,:)=allFiringRates;
+end          % maskerLevel
+
+disp('delay/ rates')
+disp(num2str(round( [1000*allDurations' summaryFiringRates])))
+
+% replot with best adjustment
+% figure(34), clf% use for separate plot
+figure(7), clf
+peakProbe=max(max(resultsMatrix));
+resultsMatrix=resultsMatrix/peakProbe;
+semilogx(gapDurations,HDresponse,'o'), hold on
+title(['FrMsk: probe ' num2str(probeLeveldB)...
+    'dB SL: peakProbe=' num2str(peakProbe,'%5.0f') ' sp/s'])
+xlabel('gap duration (s)'), ylabel ('probe response')
+semilogx(allDurations, resultsMatrix'), ylim([0 1])
+ylim([0 inf]),
+xlim([0.001 1])
+legend(strvcat(num2str(maskerLevels'-thresholdAtCF)), -1)
+
+% ------------------------------------------------- display parameters
+disp(['parameter file was: ' paramsName])
+fprintf('\n')
+% UTIL_showStruct(inputStimulusParams, 'inputStimulusParams')
+% UTIL_showStruct(outerMiddleEarParams, 'outerMiddleEarParams')
+% UTIL_showStruct(DRNLParams, 'DRNLParams')
+% UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams')
+UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
+UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
+
+
+path(restorePath);
\ No newline at end of file
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testOME.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testOME.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,90 @@
+function testOME(paramsName, paramChanges)
+
+savePath=path;
+addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
+
+if nargin<2
+    paramChanges=[];
+end
+
+sampleRate=50000;
+
+dt=1/sampleRate;
+leveldBSPL=80; 		% dB SPL as used by Huber (may trigger AR)
+amp=10^(leveldBSPL/20)*28e-6;
+duration=.05;
+time=dt: dt: duration;
+
+%% Comparison data (human)
+% These data are taken directly from Huber 2001 (Fig. 4)
+HuberFrequencies=[600	  800	 1000	 2000	 3000	4000 6000 8000];
+HuberDisplacementAt80dBSPL=[1.5E-9	1.5E-09	1.5E-09	1.0E-09	7.0E-10	...
+    3.0E-10	2.0E-10	1.0E-10]; % m;
+% HuberVelocityAt80dBSPL= 2*pi*HuberFrequencies.*HuberDisplacementAt80dBSPL;
+
+figure(2), clf, subplot(2,1,1)
+set(2,'position',[5   349   268   327])
+semilogx(HuberFrequencies, 20*log10(HuberDisplacementAt80dBSPL/1e-10),...
+    'ko', 'MarkerFaceColor','k', 'Marker','o', 'markerSize',6)
+hold on
+
+%% Generate test stimulus .................................................................
+
+% independent test using discrete frequencies
+peakResponses=[];
+peakTMpressure=[];
+frequencies=[200 400 HuberFrequencies 10000];
+for toneFrequency=frequencies
+    inputSignal=amp*sin(2*pi*toneFrequency*time);
+
+    showPlotsAndDetails=0;
+    AN_spikesOrProbability='probability';
+    % switch off AR & MOC (Huber's patients were deaf)
+    paramChanges{1}='OMEParams.rateToAttenuationFactorProb=0;';
+    paramChanges{2}='DRNLParams.rateToAttenuationFactorProb = 0;';
+
+    global OMEoutput  OMEextEarPressure TMoutput ARattenuation
+    % BF is irrelevant
+    MAP1_14(inputSignal, sampleRate, -1, ...
+        paramsName, AN_spikesOrProbability, paramChanges);
+
+    peakDisplacement=max(OMEoutput(floor(end/2):end));
+    peakResponses=[peakResponses peakDisplacement];
+
+    peakTMpressure=[peakTMpressure max(OMEextEarPressure)];
+    disp([' AR attenuation (dB):   ' num2str(20*log10(min(ARattenuation)))])
+end
+
+%% Report
+disp('frequency displacement(m)')
+% disp(num2str([frequencies' peakResponses']))
+fprintf('%6.0f \t%10.3e\n',[frequencies' peakResponses']')
+
+% stapes peak displacement
+figure(2), subplot(2,1,1), hold on
+semilogx(frequencies, 20*log10(peakResponses/1e-10), 'r', 'linewidth',4)
+set(gca,'xScale','log')
+% ylim([1e-11 1e-8])
+xlim([100 10000]), ylim([0 30])
+grid on
+title(['stapes at ' num2str(leveldBSPL) ' (NB deaf)'])
+ylabel('disp: dB re 1e-10m')
+xlabel('stimulus frequency (Hz)')
+legend({'Huber et al','model'},'location','southWest')
+set(gcf,'name','OME')
+
+% external ear resonance
+figure(2), subplot(2,1,2),hold off
+semilogx(frequencies, 20*log10(peakTMpressure/28e-6)-leveldBSPL,...
+    'k', 'linewidth',2)
+xlim([100 10000]) %, ylim([-10 30])
+grid on
+title(['External ear resonances' ])
+ylabel('dB')
+xlabel('frequency')
+set(gcf,'name','OME: external resonances')
+% ---------------------------------------------------------- display parameters
+disp(['parameter file was: ' paramsName])
+fprintf('\n')
+
+path(savePath);
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testPeriphery.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testPeriphery.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,161 @@
+function testPeriphery (BMlocations, paramsName)
+% testBM generates input output functions for DRNL model for any number
+% of locations.
+% Computations are bast on a single channel model (channelBFs=BF)
+% peak displacement (peakAmp) is measured.
+%  if DRNLParams.useMOC is chosen, the full model is run (slow)
+%  otherwise only DRNL is computed.
+% Tuning curves are generated based on a range of frequencies reletove to
+% the BF of the location.
+%
+
+global    DRNLParams
+
+savePath=path;
+
+addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
+
+levels=-10:10:90;   nLevels=length(levels);
+% levels= 50;   nLevels=length(levels);
+
+% relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];
+relativeFrequencies=1;
+
+% refBMdisplacement is the displacement of the BM at threshold
+% 1 nm disp at  threshold (9 kHz, Ruggero)
+refBMdisplacement= 1e-8;            % adjusted for 10 nm at 1 kHz
+
+toneDuration=.200;
+rampDuration=0.01;
+silenceDuration=0.01;
+
+sampleRate=30000;
+
+dbstop if error
+figure(3), clf
+% set(gcf,'position',[276    33   331   645])
+set(gcf,'name','DRNL - BM')
+
+finalSummary=[];
+nBFs=length(BMlocations);
+BFno=0; plotCount=0;
+for BF=BMlocations
+    BFno=BFno+1;
+    plotCount=plotCount+nBFs;
+    stimulusFrequencies=BF* relativeFrequencies;
+    nFrequencies=length(stimulusFrequencies);
+
+    peakAmpBM=zeros(nLevels,nFrequencies);
+    peakAmpBMdB=NaN(nLevels,nFrequencies);
+    peakEfferent=NaN(nLevels,nFrequencies);
+    peakAREfferent=NaN(nLevels,nFrequencies);
+
+
+    levelNo=0;
+    for leveldB=levels
+        disp(['level= ' num2str(leveldB)])
+        levelNo=levelNo+1;
+
+        freqNo=0;
+        for frequency=stimulusFrequencies
+            freqNo=freqNo+1;
+
+            % Generate stimuli
+            globalStimParams.FS=sampleRate;
+            globalStimParams.overallDuration=...
+                toneDuration+silenceDuration;  % s
+            stim.type='tone';
+            stim.phases='sin';
+            stim.toneDuration=toneDuration;
+            stim.frequencies=frequency;
+            stim.amplitudesdB=leveldB;
+            stim.beginSilence=silenceDuration;
+            stim.rampOnDur=rampDuration;
+            stim.rampOffDur=rampDuration;
+            doPlot=0;
+            inputSignal=stimulusCreate(globalStimParams, stim, doPlot);
+            inputSignal=inputSignal(:,1)';
+
+            %% run the model
+            MAPparamsName=paramsName;
+            AN_spikesOrProbability='probability';
+            % spikes are slow but can be used to study MOC using IC units
+            % AN_spikesOrProbability='spikes';
+
+            global DRNLoutput MOCattenuation ARattenuation IHCoutput
+            MAP1_14(inputSignal, sampleRate, BF, ...
+                MAPparamsName, AN_spikesOrProbability);
+
+            DRNLresponse=IHCoutput;
+            peakAmp=max(max(...
+                DRNLresponse(:, end-round(length(DRNLresponse)/2):end)));
+            peakAmpBM(levelNo,freqNo)=peakAmp;
+            if peakAmp>0
+            peakAmpBMdB(levelNo,freqNo)=...
+                20*log10(peakAmp/refBMdisplacement);
+            else
+                peakAmpBMdB(levelNo,freqNo)=peakAmp;
+            end
+            peakEfferent(levelNo,freqNo)=min(min(MOCattenuation));
+            peakAREfferent(levelNo,freqNo)=min(min(ARattenuation));
+
+        end  % tone frequency
+    end  % level
+
+    %% analyses results and plot
+
+    % BM I/O plot (top panel)
+    figure(3)
+    subplot(3,nBFs,BFno), cla
+    plot(levels,peakAmpBMdB, 'linewidth',2)
+    hold on, plot(levels, repmat(refBMdisplacement,1,length(levels)))
+    hold off
+    title(['BF=' num2str(BF,'%5.0f') ' - ' paramsName])
+    xlabel('level')
+    %     set(gca,'xtick',levels),  grid on
+    if length(levels)>1,xlim([min(levels) max(levels)]), end
+    ylabel(['dB re:' num2str(refBMdisplacement,'%6.1e') 'm'])
+    ylim([-20 50])
+    set(gca,'ytick',[-10 0 10 20 40])
+    %     legend({num2str(stimulusFrequencies')}, 'location', 'EastOutside')
+    UTIL_printTabTable([levels' peakAmpBMdB], ...
+        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')
+    finalSummary=[finalSummary peakAmpBMdB];
+
+    % Tuning curve
+    if length(relativeFrequencies)>2
+        figure(3), subplot(3,nBFs, nBFs+BFno)
+        %         contour(stimulusFrequencies,levels,peakAmpBM,...
+        %             [refBMdisplacement refBMdisplacement],'r')
+        contour(stimulusFrequencies,levels,peakAmpBM,...
+            refBMdisplacement.*[1 5 10 50 100])
+        title(['tuning curve at ' num2str(refBMdisplacement) 'm']);
+        ylabel('level (dB) at reference')
+        xlim([100 10000])
+        hold on
+        set(gca,'xscale','log')
+    end
+
+
+    % MOC contribution
+    figure(3)
+    subplot(3,nBFs,2*nBFs+BFno), cla
+    plot(levels,20*log10(peakEfferent), 'linewidth',2)
+    ylabel('MOC (dB attenuation)'), xlabel('level')
+    title(['peak MOC: model= ' AN_spikesOrProbability])
+    grid on
+    if length(levels)>1, xlim([min(levels) max(levels)]), end
+
+    % AR contribution
+    hold on
+    plot(levels,20*log10(peakAREfferent), 'r')
+    hold off
+
+end % best frequency
+
+UTIL_showStructureSummary(DRNLParams, 'DRNLParams', 10)
+
+    UTIL_printTabTable([levels' finalSummary], ...
+        num2str([0 stimulusFrequencies]','%5.0f'), '%5.0f')
+
+path(savePath);
\ No newline at end of file
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testPhaseLocking.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testPhaseLocking.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,51 @@
+function testPhaseLocking(paramsName, paramChanges)
+
+if nargin<2
+    paramChanges=[];
+end
+
+testFrequencies=[250 500 1000 2000 4000 8000];
+levels=50:10:80;
+figure(14), clf
+set(gcf,'position', [980    36   383   321])
+set(gcf,'name', 'phase locking')
+allStrengths=zeros(length(testFrequencies), length(levels));
+peakVectorStrength=zeros(1,length(testFrequencies));
+freqCount=0;
+for targetFrequency=testFrequencies;
+    %single test
+    freqCount=freqCount+1;
+    vectorStrength=...
+        testAN(targetFrequency,targetFrequency, levels,...
+        paramsName, paramChanges);
+    allStrengths(freqCount,:)=vectorStrength';
+    peakVectorStrength(freqCount)=max(vectorStrength');
+end
+%% plot results
+figure(14)
+subplot(2,1,2)
+plot(levels,allStrengths)
+xlabel('levels')
+ylabel('vector strength')
+legend (num2str(testFrequencies'),'location','eastOutside')
+
+subplot(2,1,1)
+semilogx(testFrequencies,peakVectorStrength)
+grid on
+title ('peak vector strength')
+xlabel('frequency')
+ylabel('vector strength')
+
+johnson=[250	0.79
+500	0.82
+1000	0.8
+2000	0.7
+4000	0.25
+5500	0.05
+];
+hold on
+plot(johnson(:,1),johnson(:,2),'o')
+legend({'model','Johnson 80'},'location','eastOutside')
+hold off
+
+
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testPhysiology.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testPhysiology.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,21 @@
+function testPhysiology(BF,paramsName, paramChanges)
+
+restorePath=path;
+addpath (['..' filesep 'MAP'])
+
+if nargin<3
+    paramChanges=[];
+end
+
+testOME(paramsName,paramChanges)
+relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];
+testBM (BF, paramsName,relativeFrequencies,'spikes', paramChanges)
+testRP(BF,paramsName,paramChanges)
+testSynapse(BF,paramsName,paramChanges)
+testFM(BF,paramsName,1,paramChanges)
+testPhaseLocking(paramsName,paramChanges)
+testAN(BF,BF, -10:10:80,paramsName,paramChanges);
+figure(14)
+figure(15)
+
+path(restorePath)
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testPhysiologyProb.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testPhysiologyProb.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,19 @@
+function testPhysiologyProb(BF,paramsName, paramChanges)
+
+restorePath=path;
+addpath (['..' filesep 'MAP'])
+
+if nargin<3
+    paramChanges=[];
+end
+
+testOME(paramsName, paramChanges)
+relativeFrequencies=[0.25    .5   .75  1  1.25 1.5    2];
+testBM (BF, paramsName,relativeFrequencies,'probability', paramChanges)
+testRP(BF,paramsName, paramChanges)
+testSynapse(BF,paramsName, paramChanges)
+testANprob(BF,BF, -10:10:80,paramsName, paramChanges);
+
+figure(4)
+
+path(restorePath)
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testRF.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testRF.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,295 @@
+function testRF
+% testIHC used either for IHC I/O function or receptive field (doReceptiveFields=1)
+
+global experiment method stimulusParameters expGUIhandles
+global inputStimulusParams  IHC_ciliaParams
+global IHC_VResp_VivoParams IHCpreSynapseParams  AN_IHCsynapseParams
+dbstop if error
+% set(expGUIhandles.pushbuttonStop, 'backgroundColor', [.941 .941 .941])
+
+addpath (['..' filesep 'MAP'], ['..' filesep 'utilities'], ...
+    ['..' filesep 'parameterStore'],  ['..' filesep 'wavFileStore'],...
+    ['..' filesep 'testPrograms'])
+
+targetFrequency=stimulusParameters.targetFrequency(1);
+
+sampleRate=50000;
+doReceptiveFields=1;
+
+toneDuration=.05;
+rampDuration=0.004;
+silenceDuration=.02;
+
+nRepeats=100;   % no. of AN fibers
+
+plotGraphsForIHC=1;
+% number of MacGregor units is set in the parameter file.
+
+if doReceptiveFields
+    % show all receptive field
+    frequencies=targetFrequency*    [  0.5         0.7         0.9     1       1.1         1.3         1.6];
+    levels=0:20:80; nLevels=length(levels);
+    figure(14), clf
+    figure(15), clf
+else
+    % show only I/O function at BF
+    frequencies=targetFrequency;
+    levels=-20:10:90;
+    %         levels=10:.25:13;
+    %         levels=-20:1:-15
+    nLevels=length(levels);
+%     figure(13), clf,
+%     set (gcf, 'name', ['IHC/AN  input/output' num2str(AN_IHCsynapseParams.numFibers) ' repeats'])
+%     drawnow
+end
+nFrequencies=length(frequencies);
+
+IHC_RP_peak=zeros(nLevels,nFrequencies);
+IHC_RP_min=zeros(nLevels,nFrequencies);
+IHC_RP_dc=zeros(nLevels,nFrequencies);
+AN_HSRonset=zeros(nLevels,nFrequencies);
+AN_HSRsaturated=zeros(nLevels,nFrequencies);
+AN_LSRonset=zeros(nLevels,nFrequencies);
+AN_LSRsaturated=zeros(nLevels,nFrequencies);
+CNLSRsaturated=zeros(nLevels,nFrequencies);
+CNHSRsaturated=zeros(nLevels,nFrequencies);
+ICHSRsaturated=zeros(nLevels,nFrequencies);
+ICLSRsaturated=zeros(nLevels,nFrequencies);
+
+
+levelNo=0; PSTHplotCount=0;
+for leveldB=levels
+    fprintf('%4.0f\t', leveldB)
+    levelNo=levelNo+1;
+    amp=28e-6*10^(leveldB/20);
+    
+    freqNo=0;
+    for frequency=frequencies
+
+        paramFunctionName=['method=MAPparams' experiment.name ...
+            '(' num2str(targetFrequency) ');' ];
+        eval(paramFunctionName);  % read parameters afresh each pass
+
+        dt=method.dt;
+        time=dt:dt:toneDuration;        
+        rampTime=dt:dt:rampDuration;
+        ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))];
+        ramp=ramp.*fliplr(ramp);
+        
+        silence=zeros(1,round(silenceDuration/dt));
+        
+        toneStartptr=length(silence)+1;
+        toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
+        toneEndptr=toneStartptr+round(toneDuration/dt) -1;
+        
+        % create signal (leveldB/ frequency)
+        freqNo=freqNo+1;
+        inputSignal=amp*sin(2*pi*frequency'*time);
+        inputSignal= ramp.*inputSignal;
+        inputSignal=[silence inputSignal silence];
+        
+        if doReceptiveFields  % receptive field
+            method.plotGraphs=	0; % plot only PSTHs
+        else
+            method.plotGraphs=	plotGraphsForIHC; % show progress
+        end
+        
+        targetChannelNo=1;
+        
+        % force parameters
+         % the number of AN fibers at each BF
+        AN_IHCsynapseParams.numFibers=	nRepeats;
+        AN_IHCsynapseParams. mode= 'spikes';
+        AN_IHCsynapseParams.plotSynapseContents=0;
+        AN_IHCsynapseParams.PSTHbinWidth=.001;
+        
+        method.DRNLSave=1;
+        method.IHC_cilia_RPSave=1;
+        method.PSTHbinWidth=1e-3; % useful 1-ms default for all PSTHs
+        method.AN_IHCsynapseSave=1;
+        method.MacGregorMultiSave=1;
+        method.MacGregorSave=1;
+        method.dt=dt;
+              
+        moduleSequence=[1:8];       
+
+        global ANdt ARAttenuation TMoutput OMEoutput ...
+    DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
+    IHCoutput ANprobRateOutput ANoutput savePavailable tauCas  ...
+    CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates MOCattenuation
+
+AN_spikesOrProbability='spikes';
+AN_spikesOrProbability='probability';
+MAPparamsName='Normal';
+
+MAP1_14(inputSignal, 1/dt, targetFrequency, ...
+    MAPparamsName, AN_spikesOrProbability);
+        
+        % RP
+        IHC_RPData=IHC_cilia_output;
+        IHC_RPData=IHCoutput(targetChannelNo,:);
+        IHC_RP_peak(levelNo,freqNo)=max(IHC_RPData(toneStartptr:toneEndptr));
+        IHC_RP_min(levelNo,freqNo)=min(IHC_RPData(toneStartptr:toneEndptr));
+        IHC_RP_dc(levelNo,freqNo)=mean(IHC_RPData(toneStartptr:toneEndptr));
+        
+        % AN next
+        AN_IHCsynapseAllData=ANoutput;
+        method.PSTHbinWidth=0.001;
+        
+        nTaus=length(tauCas);
+        numANfibers=size(ANoutput,1);
+        numLSRfibers=numANfibers/nTaus;
+        
+        %LSR (same as HSR if no LSR fibers present)
+        channelPtr1=(targetChannelNo-1)*numANfibers+1;
+        channelPtr2=channelPtr1+numANfibers-1;
+        LSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);
+        method.dt=method.AN_IHCsynapsedt;
+        PSTH=UTIL_PSTHmaker(LSRspikes, method);
+        PSTH=sum(PSTH,1); % sum across fibers (HSR only)
+        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
+        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
+        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
+        AN_LSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);
+        AN_LSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);
+        
+        % HSR
+        channelPtr1=numLSRfibers+(targetChannelNo-1)*method.numANfibers+1;
+        channelPtr2=channelPtr1+method.numANfibers-1;
+        HSRspikes=AN_IHCsynapseAllData(channelPtr1:channelPtr2,:);
+        method.dt=method.AN_IHCsynapsedt;
+        PSTH=UTIL_PSTHmaker(HSRspikes, method);
+        PSTH=sum(PSTH,1); % sum across fibers (HSR only)
+        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
+        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
+        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
+        AN_HSRonset(levelNo,freqNo)=max(max(PSTH))/(method.PSTHbinWidth*method.numANfibers);
+        AN_HSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr))/(method.numANfibers*toneDuration/2);
+        [cvANHSR, cvTimes, allTimeStamps, allISIs]=  UTIL_CV(HSRspikes, method.AN_IHCsynapsedt);
+        
+        PSTHplotCount=PSTHplotCount+1;
+        if doReceptiveFields  % receptive field for HSR only
+            figure(14), set(gcf,'name','AN')
+            plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
+            ylim([0 method.numANfibers])
+            xlabel(['CV= ' num2str(max(cvANHSR),'%4.2f')],'fontsize',8)
+        end % doReceptiveFields
+        
+        % CN
+        MacGregorMultiAllData=method.MacGregorMultiData;
+        numLSRfibers=method.McGMultinNeuronsPerBF*length(method.nonlinCF)* (nTaus-1);
+        
+        %LSR (same as HSR if no LSR fibers present)
+        channelPtr1=(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1;
+        channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1;
+        MacGregorMultiLSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:);
+        method.dt=method.MacGregorMultidt;
+        PSTH=UTIL_PSTHmaker(MacGregorMultiLSRspikes, method);
+        PSTH=sum(PSTH,1); % sum across fibers (HSR only)
+        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
+        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
+        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
+        CNLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
+        CNLSRsaturated(levelNo,freqNo)=CNLSRsaturated(levelNo,freqNo)...
+            /((toneDuration/2)*method.McGMultinNeuronsPerBF);
+        
+        %HSR
+        channelPtr1=numLSRfibers+(targetChannelNo-1)*method.McGMultinNeuronsPerBF+1;
+        channelPtr2=channelPtr1+method.McGMultinNeuronsPerBF-1;
+        MacGregorMultiHSRspikes=MacGregorMultiAllData(channelPtr1:channelPtr2,:);
+        method.dt=method.MacGregorMultidt;
+        PSTH=UTIL_PSTHmaker(MacGregorMultiHSRspikes, method);
+        PSTH=sum(PSTH,1); % sum across fibers (HSR only)
+        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
+        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
+        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
+        CNHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
+        CNHSRsaturated(levelNo,freqNo)=CNHSRsaturated(levelNo,freqNo)...
+            /((toneDuration/2)*method.McGMultinNeuronsPerBF);
+        [cvMMHSR, cvTimes, allTimeStamps, allISIs]=  UTIL_CV(MacGregorMultiHSRspikes, method.MacGregorMultidt);
+        
+        if doReceptiveFields  % receptive field
+            figure(15), set(gcf,'name','CN HSR input')
+            plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
+            ylim([0 method.McGMultinNeuronsPerBF])
+            xlabel(['CV= ' num2str(max(cvMMHSR),'%4.2f')],'fontsize',8)
+        end
+        
+        MacGregorAllData=method.MacGregorData;
+        numLSRfibers=length(method.nonlinCF)* (nTaus-1);
+        
+        %LSR (same as HSR if no LSR fibers present)
+        channelPtr1=targetChannelNo;
+        MacGregorLSR=MacGregorAllData(channelPtr1,:);
+        method.dt=method.MacGregordt;
+        PSTH=UTIL_PSTHmaker(MacGregorLSR, method);
+        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
+        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
+        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
+        ICLSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
+        ICLSRsaturated(levelNo,freqNo)=ICLSRsaturated(levelNo,freqNo)/(toneDuration/2);
+        
+        %LSR (same as HSR if no LSR fibers present)
+        channelPtr1=numLSRfibers+targetChannelNo;
+        MacGregorHSR=MacGregorAllData(channelPtr1,:);
+        method.dt=method.MacGregordt;
+        PSTH=UTIL_PSTHmaker(MacGregorHSR, method);
+        PSTHStartptr=round(silenceDuration/method.PSTHbinWidth)+1;
+        PSTHMidptr=PSTHStartptr+round(toneDuration/(2*method.PSTHbinWidth)) -1;
+        PSTHEndptr=PSTHStartptr+round(toneDuration/method.PSTHbinWidth) -1;
+        ICHSRsaturated(levelNo,freqNo)=sum(PSTH(PSTHMidptr:PSTHEndptr));
+        ICHSRsaturated(levelNo,freqNo)=ICHSRsaturated(levelNo,freqNo)/(toneDuration/2);
+        [cvICHSR, cvTimes, allTimeStamps, allISIs]=  UTIL_CV(MacGregorHSR, method.MacGregordt);
+        
+%         if doReceptiveFields  % receptive field
+%             figure(16), set(gcf,'name','IC HSR input')
+%             plotReceptiveFields(method, PSTH, PSTHplotCount, levels, frequencies)
+%             ylim([0 method.McGMultinNeuronsPerBF])
+%             xlabel(['CV= ' num2str(max(cvICHSR),'%4.2f')],'fontsize',8)
+%         end
+    end % frequency
+end % level
+fprintf('\n')
+toneDuration=2;
+rampDuration=0.004;
+silenceDuration=.02;
+nRepeats=200;   % no. of AN fibers
+fprintf('toneDuration  %6.3f\n', toneDuration)
+fprintf(' %6.0f  AN fibers (repeats)\n', nRepeats)
+fprintf('levels')
+fprintf('%6.2f\t', levels)
+fprintf('\n')
+
+
+% ---------------------------------------------------------- display parameters
+disp(['parameter file was: ' experiment.name])
+fprintf('\n')
+UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_cilia_RPParams')
+UTIL_showStruct(IHCpreSynapseParams, 'IHCpreSynapseParams')
+UTIL_showStruct(AN_IHCsynapseParams, 'AN_IHCsynapseParams')
+
+
+
+function plotReceptiveFields(method, PSTH, PSTHplotCount,  levels, frequencies)
+
+% show PSTH for each level/frequency combination
+nLevels=length(levels);
+nFrequencies=length(frequencies);
+
+PSTHtime=method.PSTHbinWidth:method.PSTHbinWidth:method.PSTHbinWidth*length(PSTH);
+subplot(nLevels,nFrequencies,PSTHplotCount)
+bar(PSTHtime, PSTH)
+xlim([0 max(PSTHtime)])
+% write axis labels only at left and bottom
+if PSTHplotCount< (nLevels-1) * nFrequencies+1
+    set(gca,'xticklabel',[])
+end
+if ~isequal(mod(PSTHplotCount,nFrequencies),1)
+    set(gca,'yticklabel',[])
+else
+    ylabel([num2str(levels(round(PSTHplotCount/nFrequencies) +1)) ' dB'])
+end
+% add titles only on top row
+if PSTHplotCount<=nFrequencies
+    title([num2str(frequencies(PSTHplotCount)) ' Hz'])
+end
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testRP.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testRP.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,205 @@
+function testRP(BFs,MAPparamsName,paramChanges)
+% testIHC used for IHC I/O function
+
+global experiment method inputStimulusParams
+global stimulusParameters IHC_VResp_VivoParams IHC_cilia_RPParams
+savePath=path;
+addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
+dbstop if error
+
+figure(4), clf,
+set (gcf, 'name', ['IHC'])
+set(gcf,'position',[613   354   360   322])
+drawColors='rgbkmcy';
+drawnow
+
+if nargin<3
+    paramChanges=[]; 
+end
+
+levels=-20:10:100;
+nLevels=length(levels);
+toneDuration=.05;
+silenceDuration=.01;
+sampleRate=50000;
+dt=1/sampleRate;
+
+allIHC_RP_peak=[];
+allIHC_RP_dc=[];
+
+for BFno=1:length(BFs)
+    BF=BFs(BFno);
+    targetFrequency=BF;
+    % OR
+    %Patuzzi and Sellick test (see ELP & AEM, 2006)
+    % targetFrequency=100;
+
+    IHC_RP_peak=zeros(nLevels,1);
+    IHC_RP_min=zeros(nLevels,1);
+    IHC_RP_dc=zeros(nLevels,1);
+
+    time=dt:dt:toneDuration;
+
+    rampDuration=0.004;
+    rampTime=dt:dt:rampDuration;
+    ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+        ones(1,length(time)-length(rampTime))];
+    ramp=ramp.*fliplr(ramp);
+
+    silence=zeros(1,round(silenceDuration/dt));
+
+    toneStartptr=length(silence)+1;
+    toneMidptr=toneStartptr+round(toneDuration/(2*dt)) -1;
+    toneEndptr=toneStartptr+round(toneDuration/dt) -1;
+
+    levelNo=0;
+    for leveldB=levels
+        levelNo=levelNo+1;
+        % replicate at all levels
+        amp=28e-6*10^(leveldB/20);
+
+        %% create signal (leveldB/ frequency)
+        inputSignal=amp*sin(2*pi*targetFrequency'*time);
+        inputSignal= ramp.*inputSignal;
+        inputSignal=[silence inputSignal silence];
+        inputStimulusParams.sampleRate=1/dt;
+%         global IHC_ciliaParams
+
+        %% disable efferent for fast processing
+        method.DRNLSave=1;
+        method.IHC_cilia_RPSave=1;
+        method.IHCpreSynapseSave=1;
+        method.IHC_cilia_RPSave=1;
+        method.segmentDuration=-1;
+        moduleSequence=1:4;
+
+        %% run the model
+        global  DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
+            IHCoutput
+        AN_spikesOrProbability='probability';
+
+        MAP1_14(inputSignal, sampleRate, BF, ...
+            MAPparamsName, AN_spikesOrProbability, paramChanges);
+
+        % DRNL
+        DRNLoutput=DRNLoutput;
+        DRNL_peak(levelNo,1)=max(DRNLoutput(toneMidptr:toneEndptr));
+        DRNL_min(levelNo,1)=min(DRNLoutput(toneMidptr:toneEndptr));
+        DRNL_dc(levelNo,1)=mean(DRNLoutput(toneMidptr:toneEndptr));
+
+        % cilia
+        IHC_ciliaData=IHC_cilia_output;
+        IHC_ciliaData=IHC_ciliaData;
+        IHC_cilia_peak(levelNo,1)=...
+            max(IHC_ciliaData(toneMidptr:toneEndptr));
+        IHC_cilia_min(levelNo,1)=...
+            min(IHC_ciliaData(toneMidptr:toneEndptr));
+        IHC_cilia_dc(levelNo,1)=...
+            mean(IHC_ciliaData(toneMidptr:toneEndptr));
+
+        % RP
+        IHC_RPData=IHCoutput;
+        IHC_RPData=IHC_RPData;
+        IHC_RP_peak(levelNo,1)=...
+            max(IHC_RPData(toneMidptr:toneEndptr));
+        IHC_RP_min(levelNo,1)=...
+            min(IHC_RPData(toneMidptr:toneEndptr));
+        IHC_RP_dc(levelNo,1)=...
+            mean(IHC_RPData(toneMidptr:toneEndptr));
+
+    end % level
+
+
+    disp(['parameter file was: ' MAPparamsName])
+    fprintf('\n')
+
+    %% plot DRNL
+    subplot(2,2,1)
+%     referenceDisp= 9e-9*1000/targetFrequency;
+%     plot(levels,20*log10(DRNL_peak/referenceDisp), drawColors(BFno), ...
+%         'linewidth',2), hold on
+    referenceDisp=10e-9;
+    plot(levels,20*log10(DRNL_peak/referenceDisp), drawColors(BFno), ...
+        'linewidth',2), hold on
+    title([' DRNL peak:  ' num2str(BFs) ' Hz'])
+    ylabel ('log10DRNL(m)'), xlabel('dB SPL')
+    xlim([min(levels) max(levels)]), ylim([-10 50])
+    grid on
+
+    %% plot cilia displacement
+    figure(4)
+    subplot(2,2,2)
+    restingIHC_cilia=IHCrestingCiliaCond;
+    plot(levels, IHC_cilia_peak,'k', 'linewidth',2), hold on
+    plot(levels, IHC_cilia_min,'r', 'linewidth',2)
+    hold on,
+    plot([min(levels) max(levels)], ...
+        [restingIHC_cilia restingIHC_cilia], 'g')
+    title(' IHC apical cond.')
+    ylabel ('IHCcilia(conductance)'), xlabel('dB SPL')
+    xlim([min(levels) max(levels)])
+    grid on
+
+    %% plot receptor potentials
+    figure(4)
+    subplot(2,2,3)
+    % RP I/O function min and max
+    restingRP=IHC_RP_peak(1);
+    toPlot= [fliplr(IHC_RP_min(:,1)') IHC_RP_peak(:,1)'];
+    microPa=   28e-6*10.^(levels/20);
+    microPa=[-fliplr(microPa) microPa];
+    plot(microPa,toPlot, drawColors(BFno), 'linewidth',2)
+    % ylim([0 300])
+
+    %% Dallos and Harris data
+    dallosx=[-0.9	-0.1	-0.001	0.001	0.01	0.9];
+    dallosy=[-8	-7.8	-6.5	11	16.5	22]/1000 + restingRP;
+    hold on, plot(dallosx,dallosy, 'o')
+    plot([-1 1], [restingRP restingRP], 'r')
+    title(' Dallos(86) data at 800 Hz')
+    ylabel ('receptor potential(V)'), xlabel('Pa')
+    ylim([-0.08 -0.02]), xlim([-1 1])
+    grid on
+
+    %% RP I/O function min and max
+    figure(4)
+    subplot(2,2,4)
+    restingRP=IHC_RP_peak(1);
+    peakRP=max(IHC_RP_peak);
+    plot(levels, IHC_RP_peak,drawColors(BFno), 'linewidth',2)
+    hold on
+    plot(levels, IHC_RP_dc, [drawColors(BFno) ':'], 'linewidth',2)
+    hold on,
+    plot([min(levels) max(levels)], [restingRP restingRP], 'r')
+    xlim([min(levels) max(levels)])
+    % animal data
+    sndLevel=[5	15	25	35	45	55	65	75];
+RPanimal=restingRP+[0.5	2	4.6	5.8	6.4	7.2	8	10.2]/1000;
+% could be misleading when restingRP changes
+RPanimal=-0.060+[0.5	2	4.6	5.8	6.4	7.2	8	10.2]/1000;
+hold on, plot(sndLevel,RPanimal,'o')
+
+    grid on
+    title(['Et= ' num2str(IHC_cilia_RPParams.Et) ':  RP data 7 kHz Patuzzi'])
+    ylabel ('RP(V)'), xlabel('dB SPL')
+    ylim([-0.08 -0.04])
+    allIHC_RP_peak=[allIHC_RP_peak IHC_RP_peak];
+    allIHC_RP_dc=[allIHC_RP_dc IHC_RP_dc];
+
+    fprintf('level\t peak\t DC\n')
+    UTIL_printTabTable([levels' IHC_RP_peak IHC_RP_dc])
+    % disp(['restingIHC_cilia= ' num2str(restingIHC_cilia)])
+    fprintf('peakRP= \t%6.3f', peakRP)
+    fprintf('\nrestingRP= \t%6.3f', restingRP)
+    fprintf('\ndifference= \t%6.3f\n', (peakRP-restingRP))
+    drawnow
+end
+% UTIL_showStruct(IHC_VResp_VivoParams, 'IHC_VResp_VivoParams')
+UTIL_showStruct(IHC_cilia_RPParams, 'IHC_cilia_RPParams')
+    fprintf('level\t peak\n')
+    UTIL_printTabTable([levels' allIHC_RP_peak])
+    fprintf('level\t DC\n')
+    UTIL_printTabTable([levels' allIHC_RP_dc])
+
+path(savePath);
+disp(paramChanges)
diff -r 02aa9826efe0 -r b51bf546ca3f testPrograms/testSynapse.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/testPrograms/testSynapse.m	Fri Jul 08 13:48:27 2011 +0100
@@ -0,0 +1,87 @@
+function testSynapse(BFlist,paramsName, paramChanges)
+% testSynapse tracks the quantity of available transmitter vesicles
+%  the computations are single channel using the first frequency
+%  in the targetFrequency box of the expGUI.
+% For, speed this function uses only probability and HSR fibers.
+% This cannot be changed because of the way AN_IHCsynapse is coded.This).
+
+global experiment  method IHCpreSynapseParams
+global  AN_IHCsynapseParams  stimulusParameters
+savePath=path;
+addpath (['..' filesep 'utilities'],['..' filesep 'MAP'])
+
+if nargin<3
+    paramChanges=[];
+end
+
+figure(6),clf
+plotColors='rbgkrbgkrbgkrbgkrbgkrbgk';
+set(gcf,'position',[5    32   264   243])
+
+sampleRate=5e4; dt=1/sampleRate;
+
+maskerLevels=-0:10:100;
+
+targetFrequency=BFlist;
+
+silenceDuration=0.015;
+maskerDuration=0.1;
+recoveryDuration=0.15;
+rampDuration=0.004;
+
+maskerTime=dt:dt:maskerDuration;
+
+rampTime=dt:dt:rampDuration;
+ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
+    ones(1,length(maskerTime)-length(rampTime))];
+ramp=ramp.*fliplr(ramp);
+
+initialSilence=zeros(1,round(silenceDuration/dt));
+recoverySilence=zeros(1,round(recoveryDuration/dt));
+
+signal=sin(2*pi*targetFrequency'*maskerTime);
+signal= ramp.*signal;
+signal=[initialSilence signal  recoverySilence];
+
+levelCount=0;
+qtMatrix=[];
+for leveldB=maskerLevels
+    levelCount=levelCount+1;
+
+    amp=28e-6*10^(leveldB/20);
+    inputSignal=amp*signal;
+
+    AN_spikesOrProbability='probability';
+    showPlotsAndDetails=0;
+
+    global savePavailable
+   
+        MAP1_14(inputSignal, 1/dt, targetFrequency, ...
+            paramsName, AN_spikesOrProbability, paramChanges);
+
+    % ignore LSR channels (if any) at the top of the matrix
+    qt=savePavailable(end, :);
+
+    synapsedt=dt;
+    time=synapsedt:synapsedt:synapsedt*length(qt);
+
+    figure(6)
+    qtMatrix=[qtMatrix; qt];
+    plot(time,qt,  plotColors(levelCount))
+    hold on
+    xlim([0 max(time)])
+    ylim([0 AN_IHCsynapseParams.M])
+end
+
+set(gcf,'name','pre-synaptic available transmitter')
+title(['q - available vesicles:' num2str(BFlist) ' Hz'])
+legend(strvcat(num2str(maskerLevels')),'location','southeast')
+legend boxoff
+grid on
+
+figure(88), [c,H]=contour(time, maskerLevels,qtMatrix); clabel(c, H);
+set(gcf,'position',[ 276    31   328   246])
+xlabel('time'), ylabel('maskerLevels')
+title('contour plot of available transmitter')
+
+path(savePath);