function [varargout] = zrun(varargin)
%the performance can be increased by precombining some operations.
global zeng zeng2 t
%varargout{1};
% NONE for normal calculation
% 0=not fired (For firing detection)
% 1=    fired (For firing detection)
%------------------------------------------------------------------------
%It is stupid, but it fixes bugs in Matlab 
%convert from ms -> ns
%dt cannot be smaller than 1 ns
if isempty(zeng2.setup.stimIin)
    %used when dt is variable
    stimIintime=[0 0];
    stimIintimerange=[0 0];
else
    stimIintime     =round(1000000*zeng2.setup.stimIin.time);
    stimIintimerange=round(1000000*zeng2.setup.stimIin.timerange);
end
if isempty(zeng2.setup.stimVe)
    %used when dt is variable
    stimVetime=[0 0];
    stimVetimerange=[0 0];
else
    stimVetime      =round(1000000*zeng2.setup.stimVe.time);
    stimVetimerange =round(1000000*zeng2.setup.stimVe.timerange);
end
%\It is stupid, but it fixes bugs in Matlab\
%------------------------------------------------------------------------
i_iterant=1; %used when dt is variable
i_sim=1;
n_sim=length(t);
if isempty(varargin)
    %for normal calculation: calculate to the end (tf)
    hWait=[];
    if zeng2.rundisp
        tic
        hWait_i_step=round(n_sim/20);
        ScreenSize=max(reshape(get(0,'screensize'),2,2)');
        hWait_position=[ScreenSize(1)/2-360/2 ScreenSize(2)/2-80/2 360 80];
        hWait=waitbar(i_sim/n_sim,['Please wait or Close this window to stop']);
        set(hWait,'unit','pixel','position',hWait_position)
    end
else
    %for firing detection: stop when fired or reach the tf.
    varargout{1}=0; %0=not fired, 1=fired
    hWait=-1;
    Detect_Start=varargin{1}(1);
    Detect_Threshold=varargin{1}(2);
end
%------------------------------------------------------------------------
if zeng.dtmode==1       
    %fix
    dt=zeng.dt/1000; %us->ms
else
    %var coarse
    dt=zeng.dtstart(zeng.dtmode)/1000; %us->ms
    dv=zeng.dv(zeng.dtmode,:);
end
%------------------------------------------------------------------------
Temperature_old=zeng.temperature;
dt_original=dt;
Gamma_dt=zeng2.setup.Gamma_dt;
Omega_dt=zeng2.setup.Omega_dt;
%---------------------------------------------
Stim_run=1;%for both Iin and Ve
Stim_PW_adjust=0;
stim_ini=zeros(zeng2.setup.nseg_total,1);
stimIin=stim_ini;
if isempty(zeng.Iinstim)
    StimIin_run=0;
    StimIinSP_run=0;
    StimIinWF_run=0;
else 
    StimIin_run=1;
    if max(zeng2.setup.stimIin.timerange)==0
        StimIinSP_run=0;     %for the square pulse
    else
        StimIinSP_run=1;     %for the square pulse
    end
    if isempty(zeng2.setup.stimIin.WaveformIndex)
        StimIinWF_run=0; %for the wave form
    else
        StimIinWF_run=1; %for the wave form
    end
end

stimVe=stim_ini;
if isempty(zeng.Iexstim)
    StimVe_run=0;
    StimVeSP_run=0;  
    StimVeWF_run=0;
else
    StimVe_run=1;
    if max(zeng2.setup.stimVe.timerange)==0
        StimVeSP_run=0;     %for the square pulse
    else
        StimVeSP_run=1;     %for the square pulse
    end
    if isempty(zeng2.setup.stimVe.WaveformIndex)
        StimVeWF_run=0; %for the wave form
    else
        StimVeWF_run=1; %for the wave form
    end
end
%---------------------------------------------
dummymodel_length=length(zeng2.dummymodel);
for s=1:dummymodel_length
    if ~(zeng2.dummymodel{s}.const_G & zeng2.dummymodel{s}.const_E)
        disp(['variable G and E are not supported yet'])
        return
    end
end
%---------------------------------------------
while t(i_sim)<zeng.tf
    for s=1:dummymodel_length
        temp_Gamma_dt=zeros(zeng2.dummymodel{s}.nodes_all_length,1);
        temp_Omega_dt=temp_Gamma_dt;
        for q=1:zeng2.dummymodel{s}.Inumber
            temp_G=zeng2.dummymodel{s}.G(:,q);
            temp_E=zeng2.dummymodel{s}.E(:,q);
            if ~isempty(zeng2.dummymodel{s}.I{q}.gate)  %detecting any current
         	    for r=1:length(zeng2.dummymodel{s}.I{q}.gate(:,1)) 
                    temp_G=temp_G.*(zeng2.dummygate{s}{zeng2.dummymodel{s}.I{q}.gate(r,1)}(:,i_sim).^zeng2.dummymodel{s}.I{q}.gate(r,2));
                end
            end
            temp_Gamma_dt=temp_Gamma_dt+temp_G;
            temp_Omega_dt=temp_Omega_dt+temp_G.*temp_E;
        end
        Gamma_dt(zeng2.dummymodel{s}.nodes_all)=temp_Gamma_dt.*zeng2.dummymodel{s}.i_C;
        Omega_dt(zeng2.dummymodel{s}.nodes_all)=temp_Omega_dt.*zeng2.dummymodel{s}.i_C;
    end
    t(i_sim+1)=t(i_sim)+dt;
    if t(i_sim+1)>zeng.tf
        t(i_sim+1)=zeng.tf;
        dt=t(i_sim+1)-t(i_sim);
    end
    if Stim_run | zeng.dtmode~=1 %dt=variable
        %When zeng.dtmode==2, t(i_sim) can go backward. The Stim_run check alone is useless.
        t_sim1=round(1000000*t(i_sim));
        t_sim2=round(1000000*t(i_sim+1));
        if Stim_PW_adjust & (zeng.dtmode==1)%fix
            %to solve the problem when dt was adjusted at the edges of PW
            dt=dt_original;
            t(i_sim+1)=t(i_sim)+dt;
            Stim_PW_adjust=0;
        end
        %Iin------------------------------------------------------------------------- 
        if StimIin_run & ~isempty(zeng2.setup.stimIin) %The 2nd part is for zeng.dtmode~=1 %dt=variable
            stimIin=stim_ini;
            if StimIinSP_run %Square Pulse Stimulation
                if t_sim1 >= stimIintimerange(2) % >=, rather than ==, is needed when dt is variable
                    if zeng.dtmode==1 %Don't combine this line with "if t_sim1 >= stimIintimerange(2)"
                        StimIinSP_run=0;
                        if StimIinWF_run==0 
                            StimIin_run=0;
                            if StimVe_run==0
                                Stim_run=0;
                            end
                        end
                    end
                elseif t_sim2 <= stimIintimerange(1)
                    
                elseif t_sim1 < stimIintimerange(1)
                    %set t+1 = the edge    
                    dt=min(zeng2.setup.stimIin.time(:,1))-t(i_sim);
                    Stim_PW_adjust=1;
                    t(i_sim+1)=t(i_sim)+dt;
                else %within the stimIint_range
                    time_stim1=find(stimIintime(:,1)<=t_sim1 & stimIintime(:,2)> t_sim1);
                    time_stim2=find(stimIintime(:,1)< t_sim2 & stimIintime(:,2)>=t_sim2);
                    if isempty(time_stim1)
                        %This can happen when dt=variable
                        %or when t(i_sim) is in between two pulses
                    else
                       for s=1:length(time_stim1)
                            stimIin(zeng2.setup.stimIin.node_amp{time_stim1(s)}(1))=zeng2.setup.stimIin.node_amp{time_stim1(s)}(2)*dt/zeng2.var{zeng2.setup.stimIin.var(time_stim1(s))}.setting.cm_area;
                       end  
                       if isempty(time_stim2)
                            dt=min(zeng2.setup.stimIin.time(time_stim1,2))-t(i_sim);
                   		    Stim_PW_adjust=1;
                			t(i_sim+1)=t(i_sim)+dt;
                       else
                            t_sim3=min(stimIintime(time_stim2,1));
	                        if t_sim3<t_sim2 & t_sim3>t_sim1
                                dt=min(zeng2.setup.stimIin.time(time_stim2,1))-t(i_sim);
                    			Stim_PW_adjust=1;
                    			t(i_sim+1)=t(i_sim)+dt;
                            end
                       end 
                    end
                end %if t_sim1 >= stimIintimerange(2)
            end %%if StimIinSP_run
            if StimIinWF_run
                temp_havesome=0;
                for i=1:length(zeng2.setup.stimIin.WaveformIndex)
                    temp_amp=feval(zeng.Iinstim(zeng2.setup.stimIin.WaveformIndex(i)).amp,round(t(i_sim+1)*1000000)/1000);
                    if ~isempty(temp_amp)
                        temp_havesome=1;
                        stimIin(zeng2.setup.stimIin.node_amp{zeng2.setup.stimIin.WaveformIndex(i)}(1))=zeng2.amp*temp_amp*dt/zeng2.var{zeng2.setup.stimIin.var(zeng2.setup.stimIin.WaveformIndex(i))}.setting.cm_area;
                    end
                end
                if temp_havesome==0 & zeng.dtmode==1
                    StimIinWF_run=0;
                    if StimIinSP_run==0;
                        StimIin_run=0;
                        if StimVe_run==0;
                            Stim_run=0;
                        end
                    end
                end %if temp_havesome==0
            end %if StimVeWF_run
        end%if StimIin_run
        %Ve--------------------------------------------------------------------------
        if StimVe_run & ~isempty(zeng2.setup.stimVe) %The 2nd part is for zeng.dtmode~=1 %dt=variable
            stimVe=stim_ini;
            if StimVeSP_run 
                if t_sim1 >= stimVetimerange(2)
                    if zeng.dtmode==1 %Don't combine this line with "if t_sim1 >= stimIintimerange(2)"
                        StimVeSP_run=0;
                        if StimVeWF_run==0
                            StimVe_run=0;
                            if StimIin_run==0
                                Stim_run=0;
                            end
                        end
                    end
                elseif t_sim2 <= stimVetimerange(1)
                elseif t_sim1 < stimVetimerange(1)
                    %set t+1 = the edge    
                    dt=min(zeng2.setup.stimVe.time(:,1))-t(i_sim);
                    Stim_PW_adjust=1;
                    t(i_sim+1)=t(i_sim)+dt;
                else %within the stimIext_range
                    time_stim1=find(stimVetime(:,1)<=t_sim1 & stimVetime(:,2)> t_sim1);
                    time_stim2=find(stimVetime(:,1)< t_sim2 & stimVetime(:,2)>=t_sim2);
                    if isempty(time_stim1)
                        %This can happen when dt=variable
                        %or when t(i_sim) is in between two pulses
                    else
                        stimVe=sum(zeng2.setup.stimVe.vestim(:,time_stim1),2);
                        if isempty(time_stim2)
                            dt=min(zeng2.setup.stimVe.time(time_stim1,2))-t(i_sim);
                   		    Stim_PW_adjust=1;
                		    t(i_sim+1)=t(i_sim)+dt;
                        else
                            t_sim3=min(stimVetime(time_stim2,1));
	                        if t_sim3<t_sim2 & t_sim3>t_sim1
                                dt=min(zeng2.setup.stimVe.time(time_stim2,1))-t(i_sim);
                    			Stim_PW_adjust=1;
                    			t(i_sim+1)=t(i_sim)+dt;
                            end
                        end 
                    end
                end %if t_sim1 == stimVetimerange(2)
            end %if StimVeSP_run
            if StimVeWF_run
                temp_havesome=0;
                for i=1:length(zeng2.setup.stimVe.WaveformIndex)
                    temp_amp=feval(zeng.Iexstim(zeng2.setup.stimVe.WaveformIndex(i)).amp,round(t(i_sim+1)*1000000)/1000);
                    if ~isempty(temp_amp)
                        temp_havesome=1;
                        stimVe=stimVe+temp_amp*zeng2.setup.stimVe.vestim(:,zeng2.setup.stimVe.WaveformIndex(i));
                    end
                end
                if temp_havesome==0 & zeng.dtmode==1
                    StimVeWF_run=0;
                    if StimVeSP_run==0;
                        StimVe_run=0;
                        if StimIin_run==0;
                            Stim_run=0;
                        end
                    end
                end %if temp_havesome==0
            end %if StimVeWF_run
        end%if StimVe_run
    end %if Stim_run
    %[i_sim t(i_sim) stimVe(28)]
    %----------------------------------------------------------------------------
    Vout=feval(zeng2.method,zeng2.dummyvm(:,i_sim),dt,stimVe,stimIin,zeng2.setup.lamda_dt_matrix,Gamma_dt,Omega_dt);
    %----------------------------------------------------------------------------
    if isempty(varargin)
        %normal calculation
        if  zeng2.rundisp
            if ishandle(hWait) 
                if mod(i_sim,hWait_i_step)==0
    		        waitbar(i_sim/n_sim)
                end
            else
	            break
            end
        end
    else
        %for fired-detecting
        if t(i_sim)>Detect_Start
            if 1  %detect every node
                if ~isempty(find(Vout>Detect_Threshold))
                    varargout{1}=1;
                    return
                end
            else  %detect only certain nodes
                temp=round(length(Vout)/2);
                if ~isempty(find(Vout([-5:1:5]+temp)>Detect_Threshold))
                    varargout{1}=1;
                    figure(99)
            %       plot(zeng2.dummyvm(round(size(zeng2.dummyvm,1)/2),:))
                    mesh(zeng2.dummyvm(:,1:(i_sim+1)))
                    drawnow
                    %pause(1)
                    return
                end
            end
        end

    end
    if zeng.dtmode==1 %fix
        %/common----------------------------------------------------------------------------
    	i_sim=i_sim+1;
    	zeng2.dummyvm(:,i_sim)=Vout;
        if ~isempty(zeng.Iexstim)
            zeng2.dummyve(:,i_sim)=stimVe;
        end
        
    	for s=1:dummymodel_length
            %Enable this part if you want to Adjust q10 in real time
            %if Temperature_old~=zeng.temperature
          	%	feval(zeng2.dummymodel{s}.name,'setting',zeng.temperature); 
            %   Temperature_old=zeng.temperature
            %end
       		for q=1:zeng2.dummymodel{s}.gatenumber
         		zeng2.dummygate{s}{q}(:,i_sim)=feval(zeng2.dummymodel{s}.name,zeng2.dummymodel{s}.gate{q},Vout([zeng2.dummymodel{s}.nodes_all]),zeng2.dummygate{s}{q}(:,i_sim-1),dt);
       		end
        end
       	%common/----------------------------------------------------------------------------
	else %if zeng.cal.dtval==1 %fix
        i_iterant=i_iterant+1;
        if i_iterant>10000
            disp('----------------------------------------------------')
            disp(['XXX--Variable time step cannot work with this configuration.'])
            disp(['XXX--please use fixed time step.'])
            disp('----------------------------------------------------')
            break
        end
      	%----------------------------------------------------------------------------
      	%check if dt is too big
        temp=max(abs(Vout-zeng2.dummyvm(:,i_sim)));
        if temp>dv(2)
            %lower the dt and recalculate.
            dt=max(dt/4,zeng.dtmin); %cannot be the same as when temp<dv(1)
            %--------------------------------------------------------------
            %"Two step" back".
            i_sim=max(i_sim-1,1);
            %--------------------------------------------------------------
        else
            %There is no need to the find a new value of Vprevious for 1st order Backward.
            %check if dt is too small
            if	temp<dv(1)
	    		dt=min(dt*2,zeng.dtmax); 
            end
   	    	%/common----------------------------------------------------------------------------
	        i_sim=i_sim+1;
    	    zeng2.dummyvm(:,i_sim)=Vout;
            if ~isempty(zeng.Iexstim)
                zeng2.dummyve(:,i_sim)=stimVe;
            end
    		for s=1:dummymodel_length
                %Enable this part if you want to Adjust q10 in real time
                %if Temperature_old~=zeng.temperature
       			%    feval(zeng2.dummymodel{s}.name,'setting',zeng.temperature)
                %   Temperature_old=zeng.temperature
                %end
        		for q=1:zeng2.dummymodel{s}.gatenumber
            		zeng2.dummygate{s}{q}(:,i_sim)=feval(zeng2.dummymodel{s}.name,zeng2.dummymodel{s}.gate{q},Vout([zeng2.dummymodel{s}.nodes_all]),zeng2.dummygate{s}{q}(:,i_sim-1),dt);
            	end
            end
     	    %common/----------------------------------------------------------------------------
            if i_sim==n_sim
                %add allocate.(double)
        		zeng2.dummyvm=[zeng2.dummyvm zeros(zeng2.setup.nseg_total,n_sim)];
                if ~isempty(zeng.Iexstim)
                    zeng2.dummyve=[zeng2.dummyve zeros(zeng2.setup.nseg_total,n_sim)];
                end
        		for s=1:dummymodel_length
        		    for q=1:zeng2.dummymodel{s}.gatenumber
            		    zeng2.dummygate{s}{q}=[zeng2.dummygate{s}{q} zeros(zeng2.dummymodel{s}.nodes_all_length,n_sim)];
            	    end
                 end
                 n_sim=2*n_sim;
            end
	   	end %temp>dv(2)
        %check if more n_sim is required.
    end %if zeng.cal.dtval==1 %fix
end %t(i)<zeng.cal.tf

if isempty(varargin)
    %normal calculation
    if zeng.dtmode==1 %fix
    else
        %remove junk------------------------------------------
        t=t(1:i_sim);
        zeng2.dummyvm=zeng2.dummyvm(:,1:i_sim);   
        if ~isempty(zeng.Iexstim)
            zeng2.dummyve=zeng2.dummyve(:,1:i_sim);
        end
        for s=1:dummymodel_length
            for q=1:zeng2.dummymodel{s}.gatenumber
                zeng2.dummygate{s}{q}=zeng2.dummygate{s}{q}(:,1:i_sim);
            end
        end
        %----------------------------------------------------
    end
    if ishandle(hWait)
        close(hWait)
    end
    
    if zeng2.rundisp
        temp=toc;
        if temp<=60
            disp(['finish run with ' num2str(temp) ' seconds' ])
        elseif temp<=3600
            temp_minute=fix(temp/60);
            temp_second=mod(temp,60);
            disp(['finish run with ' num2str(temp_minute) ':' num2str(temp_second) ' min:sec' ])
        else
            temp_hour=fix(temp/3600);
            temp=temp-temp_hour*3600;
            temp_minute=fix(temp/60);
            temp_second=mod(temp,60);
            disp(['finish run with ' num2str(temp_hour) ':' num2str(temp_minute) ':' num2str(temp_second) ' h:m:sec' ])
        end
    end
else %isempty(varargin)
    %Firing detection

    figure(99)
   %plot(zeng2.dummyvm(round(size(zeng2.dummyvm,1)/2),:))
    mesh(zeng2.dummyvm)
    drawnow
    %pause(1)
end