%  FIG.4B, FIG.4C

len1=90;
trails=3; 


v_p3e=linspace(0,3.01,len1)*1.0e-3;
v_sr=0.75e-3; 
v_es=1.799e-3;
delay2=0.06;
v_p3xi=0.101e-3;

% v_sr=1*1.0e-3;
v_p1xi=3.01e-4;
open1=0;
open2=0;
KK=1;

open3=0;
open4=1;
AA=1;




% matlabpool local 3


Pmax1=zeros(len1,trails);
Pmax2=zeros(len1,trails);
Pmin1=zeros(len1,trails);
Pmin2=zeros(len1,trails);
FR=zeros(len1,trails);
State=zeros(len1,trails);
FD=zeros(len1,trails);

for i=1:len1
    tic
    vv_p3e=v_p3e(i); 
    pmax1=[];
    pmax2=[];
    pmin1=[];
    pmin2=[];
    state=[];
    fd=[];
    parfor j=1:trails
        [p1,p2,p3,p4,p5,p6]=BGCT_subfun2(delay2,v_sr,v_p1xi,open1,open2,KK,v_p3xi,open3,open4,AA,v_es,vv_p3e);
        pmax1=[pmax1,p1];
        pmax2=[pmax2,p2];
        pmin1=[pmin1,p3];
        pmin2=[pmin2,p4];
        state=[state,p5];
        fd=[fd,p6];
    end
    
    Pmax1(i,:)=pmax1;
    Pmax2(i,:)=pmax2;
    Pmin1(i,:)=pmin1;
    Pmin2(i,:)=pmin2;
    State(i,:)=state;
    FD(i,:)=fd;
    toc
end

% matlabpool close
plot(linspace(0.0,3.01,len1),FD,'.')
% v_p3e=v_p3e*1000;
% plot(v_p3e,Pmax1,'.',v_p3e,Pmax2,'.',v_p3e,Pmin2,'.',v_p3e,Pmin1,'.')