/*
 *  test_hh_psc_alpha.sli
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

 /* BeginDocumentation
Name: testsuite::test_hh_psc_alpha - sli script for hh_psc_alpha model

Synopsis: (test_hh_psc_alpha) run -> 'rates.txt'


Description:

  test_hh_psc_alpha produces a rate-response (FI) curve of the Hogkin-Huxley
  neuron  in response to a range of different current (DC) stimulations.
  The result is stored in an extra file, rates.txt, and can be plotted,
  e.g. in matlab using
  x=load('hh_rates.txt'); plot(x(:,1),x(:,2))
  or
  plot "hh_rates.txt" with lines;
  in gnuplot
  
  Since a DC input affetcs only the neuron's channel dynamics, this routine
  does not yet check correctness of synaptic response.

Author: Schrader

SeeAlso: iaf_psc_exp, iaf_neuron, testsuite::test_iaf_i0, testsuite::test_iaf_i0_refractory, testsuite::test_iaf_dc
*/


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARAMETER SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  /simtime 1000 def % simulation time in ms - the longer the smoother the curve

  % amplitude range - in pA
  /dcfrom  0    def
  /dcstep  20   def
  /dcto    2000 def

  0.1 /h Set % simulation step size in mS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  /rateprint {
      f                           % file handle
      amp print (\t) print        % print amplidute with tab
      sd GetStatus /n_events get  % get number of spikes
      simtime 1000 div div        % convert to Hz
      print endl ;                % linefeed and pop
  } bind def

%%%%%%%%%%%%%%%%%%%%%%%%%%%% SIMULATION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % open datafile
  (hh_rates.txt) ofstream
  {
      /f Set
  }
  {
      /test_hh_psc_alpha /CouldNotOpenFile raiseerror
  } ifelse

  
   ResetKernel      % clear all existing network elements

   modeldict begin  
   synapsedict begin
   userdict begin

   /hh_psc_alpha Create   /neuron Set
   /spike_detector Create /sd     Set
      
   sd << /to_memory false >> SetStatus % spikes are not recorded but number of spikes (n_events)
      
   neuron sd 1.0 h Connect

   M_WARNING      setverbosity % hide multiple simulation calls
   
   % simulation loop
   dcfrom dcstep dcto {

       cvd /amp Set     % dc stimulation in pA
       neuron << /I_e amp >> SetStatus

       (Simulating with current I=)  amp cvs join  (pA) join =
	
	1000 Simulate                    % one second warm-up time for equilibrium state
	sd << /n_events 0 >> SetStatus   % then reset spike counts
	simtime Simulate                 % another simulation call to record neuron's 'clean' firing rate

	rateprint % print result to file

   } bind
   for

  end
  end
  end

  % close datafile
  f close

  (done.) =