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
