function out = cm( task, opt)
% cm: a sample m-file for simulating conductance-based neuron models
%	It could be as a template for various dynamic simulations.
% Usage:
% >>cm init	
%	initialization
% >>cm run
%	run a simulation
% Use with: hh.m, ml.m
% August 1999 by K. Doya

% Declaration of global variables to be used repeatedly
% any variable not declared as "global" will be cleared after each call
global	model			% model name: 'hh' or 'ml'
global	ts ns vname		% time span, state dim, and variable names
global	Inp				% current input
global	T Xt It 		% time, state, and input tables
% If you want to check or change these variables from command line, type, e.g.
% >>global ts

% If you just type
% >>cm
% a demo script below will be run.
if nargin == 0
	cm init
	% the above is interpreted as "cm( 'init')" in matlab's hacky syntax.
	cm run
	return
end
% If a task is specified:
switch( task)
case 'init'				% Initialization
	if nargin >= 2, model = opt, end;	% change model
	if isempty(model), model = 'hh', end	% default model
	% default time span and initial state
	[ts,x] = feval( model,[],[],'init');
	% is equivalent to, e.g., hh([],[],'init');
	ns = size(x,1);		% state space dimension
	vname = feval( model,[],[],'name');	% variable names
	cm( 'figs');		% setup figures
	cm( 'pulse', 10);	% default stimulus
	T = 0;				% initialize time, state, and input trace
	Xt = x';
	It = 0;
case 'run'				% Run the simulation
	figure(1);			% real-time display of phase protrait
	% options for ode simulation: plot first and last variable
	odeopt = odeset('OutputFcn','odephas2','OutputSel',[1,ns]);
	% integration of 'stiff' differential equation
	% time steps are chosen automatically
	% Inp could be constant or (t,I) table
	if size(Inp,2)==2, ts=Inp([1,end],1); end
	[ T1, X1] = ode23s( model, ts, Xt(end,:)', odeopt, Inp);
	% concatenate with previous output
	T = [ T; T(end)+T1];	% time 
	Xt = [ Xt; X1];		% state
	if size(Inp,2)==2	% input
		It = [ It; interp1( Inp(:,1), Inp(:,2), T1)];	% table
	else
		It = [ It; repmat( Inp(1), size(T1))];		% const
	end
	figure(2); cm( 'wave');
	figure(1); cm( 'phase');
case 'figs'
	ss = get(0,'ScreenSize');
	dx = ss(3)/2; dy = ss(4)/2;		% figure offset
	sx = dx-24; sy = dy-48;			% figure size
	figure(1); set(1,'Position',[10,dy,sx,sy],'Name','Phase'); clf;
	grid('on'); xlabel( vname(1)); ylabel( vname(end));
	figure(2); set(2,'Position',[10+dx,dy,sx,sy],'Name','Wave'); clf;
case 'phase'		% Plot a phase portrait
	plot( Xt(:,1), Xt(:,end));		% plot x1 and xn
	grid('on');  xlabel( vname(1)); ylabel( vname(end));
case 'wave'			% Plot waveforms
	for i = 1:ns
		subplot(ns+1,1,i);	% plot i-th state variable
		plot( T, Xt(:,i));
		ylabel( vname(i,:));
	end
	subplot(ns+1,1,ns+1);	% plot current input
	plot( T, It(:));
	ylabel( 'I'); xlabel( 't (ms)');
case 'inp'			% Current input
	if ischar(opt), opt = eval(opt); end	% strig to value
	Inp = opt				% maybe constant or table
case 'pulse'			% Simple pulse
	if ischar(opt), opt = eval(opt); end	% strig to value
	t = [ 0; 0.2; 0.201; 0.25; 0.251; 1]*(ts(end)-ts(1)) + ts(1);
	i = [ 0; 0; 1; 1; 0; 0]*opt;
	Inp = [ t, i]
case 'save'			% Save all global variables, e.g.: >>cm save 'data1'
	save( opt);
case 'load'			% Relaod all global variables, including T, X
	load( opt);
otherwise
	error( [' invalid task: ', task]);
end

% end of cm.m