function varargout = Find_Ith(varargin)
global zeng zeng2 ZengIth
if nargin==0
action='initial';
else
action=varargin{1};
end
switch action
case 'initial'
%model setup-----------------------------------------------------
[Filename Path]=uigetfile([zeng2.cd zeng2.FILESEP 'model' zeng2.FILESEP '*.m']);
if strcmp('0',num2str(Filename))
return
end
zexst('open model file',Filename);
%/model setup-----------------------------------------------------
zeng2.options.vestim=1; % 1=calculated, 2=import
zexst('change options','vestim',1);
%zexst('DeleteFcn')
%------------------------------------
if length(zeng2.varlist)>1
zexst('err','Sorry.. This version support only one "axon" in the model file.');
return
else
for tempzeng=1:length(zeng2.varlist)
feval('global',zeng2.varlist{tempzeng})
end
end
%------------------------------------
zeng.dtmode=1;
zeng.dt=0.5; %us
zeng.tf=0.25; %ms
ZengIth.cal.Ierr=1; %
ZengIth.Istart= 500;
ZengIth.Ilimit= 200*ZengIth.Istart;
%-----------------------------------------
Fired_para(1)=0.1; %ms %[Detect_Start Detect_Threthold ]
eval(['Fired_para(2)=75+' zeng2.var{1}.name '.vini;']) %[Detect_Start Detect_Threthold]
%-----------------------------------------
if 0
ZengIth.cal.Dia=10;
else
% ZengIth.cal.Dia=1:0.5:20;
% ZengIth.cal.Dia=5:1:15;
ZengIth.cal.Dia=[14 15 21:25];
% ZengIth.cal.Dia=5:5:10;
end
ZengIth.cal.Dia_length=length(ZengIth.cal.Dia);
%-----------------------------------------
ZengIth.cal.X_number=5;
ZengIth.cal.X=zeros(ZengIth.cal.X_number,ZengIth.cal.Dia_length);
ZengIth.cal.X=zeros(ZengIth.cal.X_number,ZengIth.cal.Dia_length);
if 1 %systematic 0-L
ZengIth.cal.X(:,1)=linspace(0,1,ZengIth.cal.X_number)';
for i=2:ZengIth.cal.Dia_length
ZengIth.cal.X(:,i)=ZengIth.cal.X(:,1);
end
elseif 0 %systematic 0-L/2
ZengIth.cal.X(:,1)=linspace(0,0.5,ZengIth.cal.X_number)';
for i=2:ZengIth.cal.Dia_length
ZengIth.cal.X(:,i)=ZengIth.cal.X(:,1);
end
else %random
ZengIth.cal.X=rand(ZengIth.cal.X_number,ZengIth.cal.Dia_length)*2-1; %[-1 +1]
end
eval(['Linter_D=' zeng2.var{1}.name '.Linter_D;'])
for i=1:ZengIth.cal.Dia_length
ZengIth.cal.X(:,i)=ZengIth.cal.X(:,i)*Linter_D*ZengIth.cal.Dia(i);
end
%-----------------------------------------
%Don't sort(X), do sort(Z) only.
if 0 % random_z
temp=[50 200]; %range (um)
ZengIth.cal.Z=rand(ZengIth.cal.X_number,ZengIth.cal.Dia_length)*diff(temp)+temp(1);
ZengIth.cal.Z=sort(ZengIth.cal.Z);
else
ZengIth.cal.Z=[];
end
%-----------------------------------------
ZengIth.cal.Ith=zeros(ZengIth.cal.X_number,ZengIth.cal.Dia_length);
Iexstim_original=zeng.Iexstim;
temp=0;
Iexstim_n=length(zeng.Iexstim);
for k=1:Iexstim_n
temp=max(temp,abs(zeng.Iexstim(k).amp));
end
for k=1:Iexstim_n
zeng.Iexstim(k).amp=zeng.Iexstim(k).amp/temp;
end
Iexstim=zeng.Iexstim;
%-----------------------------------------
Total_axons=num2str(ZengIth.cal.Dia_length*ZengIth.cal.X_number);
Total_axons_i=0;
%-----------------------------------------
tic
%---------------------------
if 1 %for early termination
Early_Ter_check=1;
Early_Ter=figure(99);
text(0.25,0.5,'close me to stop')
else
Early_Ter_check=0;
end
for i=1:ZengIth.cal.Dia_length
%-----------------------------------------
if 0 %vary L
temp.L=400000; %um
eval([zeng2.var{1}.name '.dia =ZengIth.cal.Dia(i);'])
eval([zeng2.var{1}.name '.nseg =1+2*round(temp.L/2/ZengIth.cal.Dia(i)/' zeng2.var{1}.name '.Linter_D)']);
eval([zeng2.var{1}.name '.L =' zeng2.var{1}.name '.nseg*ZengIth.cal.Dia(i)*' zeng2.var{1}.name '.Linter_D']);
else %donot alter L
eval([zeng2.var{1}.name '.dia =ZengIth.cal.Dia(i);'])
end
%-----------------------------------------
for j=1:ZengIth.cal.X_number
Total_axons_i=Total_axons_i+1;
disp([num2str(Total_axons_i) '/' Total_axons])
Iminmaxmean=ZengIth.Istart;
Iminmax(1:2)=Iminmaxmean;
for k=1:Iexstim_n
zeng.Iexstim(k).amp =Iexstim(k).amp*Iminmaxmean;
zeng.Iexstim(k).xyz(1)=Iexstim(k).xyz(1)+ZengIth.cal.X(j,i);
if ~isempty(ZengIth.cal.Z)
zeng.Iexstim(k).xyz(3)=ZengIth.cal.Z(j,i);
end
end
zexst('Update GUI')
drawnow
[Fired]=zexst('calculate',Fired_para);
disp(['start with I = ' num2str([Iminmax]) ' - ' num2str(Fired)])
if Fired==-1
return
elseif Fired==0
while Fired==0
if Early_Ter_check
if ~ishandle(Early_Ter)
return
end
end
Iminmax(1)=Iminmax(2);
Iminmax(2)=Iminmax(2)*2;
for k=1:Iexstim_n
zeng.Iexstim(k).amp=Iminmax(2)*Iexstim(k).amp;
end
[Fired]=zexst('calculate',Fired_para);
disp([' I = ' num2str([mean(Iminmax)]) ' - ' num2str(Fired)])
if Iminmax(2)>ZengIth.Ilimit
Fired=1
end
end
elseif Fired==1
while Fired==1
if Early_Ter_check
if ~ishandle(Early_Ter)
return
end
end
Iminmax(2)=Iminmax(1);
Iminmax(1)=Iminmax(1)*0.75;
for k=1:Iexstim_n
zeng.Iexstim(k).amp=Iminmax(1)*Iexstim(k).amp;
end
[Fired]=zexst('calculate',Fired_para);
disp([' I = ' num2str([mean(Iminmax)]) ' - ' num2str(Fired)])
end
else
disp('err')
return
end
disp(['I found you at ' num2str(Iminmax)])
while 1
if Early_Ter_check
if ~ishandle(Early_Ter)
return
end
end
Iminmaxmean=mean(Iminmax);
if Iminmaxmean>ZengIth.Ilimit
ZengIth.cal.Ith(j,i)=ZengIth.Ilimit;
break
end
for k=1:Iexstim_n
zeng.Iexstim(k).amp=Iminmaxmean*Iexstim(k).amp;
end
[Fired]=zexst('calculate',Fired_para);
disp([' I = ' num2str([Iminmaxmean]) ' - ' num2str(Fired)])
if Fired==-1
return
elseif Fired==1
if (100*(Iminmaxmean-Iminmax(1))/Iminmaxmean)>ZengIth.cal.Ierr
Iminmax(2)=Iminmaxmean;
else
disp([' I = ' num2str([Iminmaxmean])])
ZengIth.cal.Ith(j,i)=Iminmaxmean;
break
end
elseif Fired==0
if (100*(Iminmax(2)-Iminmaxmean)/Iminmaxmean)>ZengIth.cal.Ierr
Iminmax(1)=Iminmaxmean;
else
disp([' I = ' num2str([Iminmax(2)])])
ZengIth.cal.Ith(j,i)=Iminmax(2);
break
end
end
end
end %for j=1:ZengIth.cal.X_number
Find_Ith('Save Threshold')
end% for i=1:ZengIth.cal.Dia_length
zeng.Iexstim=Iexstim_original;
Find_Ith('Save Threshold')
temp=toc;
figure
if size(ZengIth.cal.Ith,1)==1
plot(ZengIth.cal.Dia,ZengIth.cal.Ith);
else
errorbar(ZengIth.cal.Dia,mean(ZengIth.cal.Ith),std(ZengIth.cal.Ith),'b');
end
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
case 'Save Threshold'
if 1
ZengIth.Filename=[zeng2.file(1:(length(zeng2.file)-2)) '_Z-200' '_PW-' num2str(1000*zeng.Iexstim(1).dur) '_dt-' num2str(zeng.dt) '_D2-20.thr'];
ZengIth.Path=[zeng2.cd '\'];
else
[Filename Path]=uiputfile([thr.cd '\*.thr']);
end
if ~strcmp('0',num2str(ZengIth.Filename)) & ~strcmp('0',num2str(ZengIth.Path));
Have_dot=find(ZengIth.Filename=='.');
if isempty(Have_dot)
ZengIth.Filename=[ZengIth.Filename '.thr'];
end
dlmwrite([ZengIth.Path ZengIth.Filename],'');
FID=fopen([ZengIth.Path ZengIth.Filename],'w');
fprintf(FID,['I(mA) \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).amp) '\t' ]);
end
fprintf(FID,['\n PW(us) \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).dur*1000) '\t' ]);
end
fprintf(FID,['\n Delay(us) \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).delay*1000) '\t' ]);
end
fprintf(FID,['\n X \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).xyz(1)) '\t' ]);
end
fprintf(FID,['\n Y \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).xyz(2)) '\t' ]);
end
fprintf(FID,['\n Z \t']);
for k=1:length(zeng.Iexstim)
fprintf(FID,[num2str(zeng.Iexstim(k).xyz(3)) '\t' ]);
end
fprintf(FID,['\n Time step (us) \t']);
fprintf(FID,[num2str(zeng.dt)]);
fprintf(FID,['\n']);
for i=1:ZengIth.cal.Dia_length
fprintf(FID,[num2str(ZengIth.cal.Dia(i)) '\t' ]);
end
fprintf(FID,['\t']);
fprintf(FID,['\t' 'X']);
for i=1:ZengIth.cal.Dia_length
fprintf(FID,[ '\t' ]);
end
fprintf(FID,['\t' 'Z']);
fprintf(FID,['\n']);
for j=1:ZengIth.cal.X_number
for i=1:ZengIth.cal.Dia_length
fprintf(FID,[num2str(ZengIth.cal.Ith(j,i)) '\t' ]);
end
fprintf(FID,['\t']);
for i=1:ZengIth.cal.Dia_length
fprintf(FID,[num2str(ZengIth.cal.X(j,i)) '\t' ]);
end
fprintf(FID,['\t']);
if isempty(ZengIth.cal.Z)
else
for i=1:ZengIth.cal.Dia_length
fprintf(FID,[num2str(ZengIth.cal.Z(j,i)) '\t' ]);
end
end
fprintf(FID,['\n']);
end
Status=fclose(FID);
end
end
