function simulation(network,para,iinp,trial_num,hd_txt,OUTPUT_OPT)

%%DO NOT EDIT THIS FILE, UNLESS YOU ARE ABOSLUTE SURE WHAT YOU ARE DOING.
%%This is the template file use to generate the ODE simulation .m file
%%"simulation.m"
%%Written by Ning Jiang, Institute of Biomedical Engineering, Univesity of New
%%Brunswick, NB, Canada, E3B 5A3
%%Email: ning.jiang@unb.ca
%%Created Date: Nov 19, 2006
%%Rev. 1.1, Nov 20, 2006 --adding para to the saved variables

global inp;
duration=para.duration*1000; 
neuron_num=size(iinp,2);
if OUTPUT_OPT
    out_v=zeros(1,neuron_num);
    out_t=0;
end
initial=zeros(1,3*neuron_num);
initial(1:3:3*neuron_num-2)=rand(1,neuron_num)*5-60;
initial(2:3:3*neuron_num-1)=rand(1,neuron_num);
initial(3:3:3*neuron_num)=rand(1,neuron_num)*0.5;
options=odeset('InitialStep',para.stepsize);
firing_t=cell(1,neuron_num);
v=initial(1:3:3*neuron_num-2);
t_course=0;
t0=0; 
for n=1:size(iinp,1)
    inp=iinp(n,:);
    [t,y]=ode45(@hhn,[t0,t0+para.stepsize],initial,options);
    initial=y(size(y,1),:);
    v=[v;y((2:size(y,1)),(1:3:neuron_num*3-2))];
    t_course=[t_course;t(2:length(t))];
    if rem(n,1000/para.stepsize) == 0
        if hd_txt
            set(hd_txt,'string',['trial ',num2str(trial_num),':  ',num2str(duration/1000-n*para.stepsize/1000),' sec left']);
            pause(0.01);
        else
            disp(['trial ',num2str(trial_num),':  ',num2str(duration/1000-n*para.stepsize/1000),' sec left']);
        end
        for nn=1:neuron_num
            [ft,abnormal]=get_interval(v(:,nn),t_course,1);
            if ft ~=0
                firing_t{nn}=[firing_t{nn},ft];
            end
        end
        if OUTPUT_OPT
            out_v=[out_v;v(2:size(v,1),:)];
            out_t=[out_t;t_course(2:length(t_course))];
        end
        v=v(size(v,1),:);
        t_course=t_course(length(t_course));
    end
    t0=t(length(t));
end
if OUTPUT_OPT
    out_t=out_t(2:length(out_t));
    out_v=out_v(2:size(out_v,1),:);
    save([para.resultfile,'_',num2str(trial_num)],'firing_t','out_t','out_v','iinp','network','para');
else
    save([para.resultfile,'_',num2str(trial_num)],'firing_t','iinp','network','para');
end

function dydt=hhn(t,y)
global inp;