Bstart = 0.01;
Bend = 100.0;
Astart =0.0001;
Aend = 0.9999;
Bprec = 200000; %2000
Aprec = 200; %200
SH = -1.0*ones(Aprec+1,Bprec+1);
%t_sh= -1.0*ones(Aprec+1,Bprec+1);
%flux= -1.0*ones(Aprec+1,Bprec+1);
%%throughput= -1.0*ones(Aprec+1,Bprec+1);
%throughput_per_energy= -1.0*ones(Aprec+1,Bprec+1);
%rho_integral= -1.0*ones(Aprec+1,Bprec+1);
%energy= -1.0*ones(Aprec+1,Bprec+1);
%throughput_per_energy_08= -1.0*ones(Aprec+1,Bprec+1);
%throughput_per_energy_06= -1.0*ones(Aprec+1,Bprec+1);
%throughput_per_energy_04= -1.0*ones(Aprec+1,Bprec+1);



%size(SH)
%k=waitforbuttonpress

A_vector = linspace(Astart,Aend,Aprec+1);
B_vector = linspace(Bstart,Bend,Bprec+1);

for B = Bstart:((Bend-Bstart)/Bprec):Bend

    for A  = Astart:((Aend-Astart)/Aprec):Aend
        i = round(1+(A-Astart)/(Aend-Astart)*Aprec);
        j = round(1+(B-Bstart)/(Bend-Bstart)*Bprec);
        %flux(i,j) = 1-A/2; 
       % k=waitforbuttonpress
        if (A*B)<1-sqrt(A*(2-A))
            SH(i,j) = 0.0; %cake shutdown only possible
            %t_sh(i,j) = (2*B)/((2-A)^2);

        elseif (A*B)<1
            SH(i,j) = 1.0; % filter shutdown first
            %t_sh(i,j) = ((A-1)^2 + 2*(A*B-1) +2*abs(1-A*B)*sqrt(abs(A*(2-A))))/(A^2*B*(2-A)^2/2);

        else
            SH(i,j) = 1.0; % filter shutdown first
            %t_sh(i,j) = ((A-1)^2 + 2*(A*B-1) -2*abs(1-A*B)*sqrt(abs(A*(2-A))))/(A^2*B*(2-A)^2/2);

        end
        %throughput(i,j) = flux(i,j)*t_sh(i,j);
        %t_sh(i,j)
        %t = linspace(0.0,t_sh(i,j)-0.00000000001,1000);
        %A
        %F1 = flux(i,j)*(1+A*B*flux(i,j)*t).^2+ B*(1-A*B)^2*flux(i,j)*t + (B-A*B*B)*(1+A*B*flux(i,j)*t);
        %F2 = 1+A*B*flux(i,j)*t;
        %F3 = B*(1-A*B)^2*(2*flux(i,j)*t+2*A*B*flux(i,j)*flux(i,j)*t.^2);
        %F4 = B^2*(1-A*B)^2*(1+flux(i,j)*t).^2;
        %F5 = (F3-2*F1.*F2 + sign(1-A*B)* realsqrt(F3.^2 + 4*F4.*F2.^2 -4*F1.*F2.*F3))./(2*F2.^2);
        %rho = 1/A*(1-realsqrt(1+2*A*F5));
        %rho_integral(i,j) = trapz(t,rho);
        %energy(i,j)=rho_integral(i,j)*flux(i,j);
        %throughput_per_energy(i,j) = throughput(i,j)/(energy(i,j)^(1.0)); % change the exponent from 1.0 down to 0.5
        %throughput_per_energy_08(i,j) = throughput(i,j)/(energy(i,j)^(0.8)); % change the exponent from 1.0 down to 0.5
        %throughput_per_energy_06(i,j) = throughput(i,j)/(energy(i,j)^(0.6)); % change the exponent from 1.0 down to 0.5
        %throughput_per_energy_04(i,j) = throughput(i,j)/(energy(i,j)^(0.4)); % change the exponent from 1.0 down to 0.5


        
        %rho_integral(i,j)
        %k=waitforbuttonpress
    end
end
%size((Astart:((Aend-Astart)/Aprec):Aend)(:))
%size((Bstart:((Bend-Bstart)/Bprec):Bend)(:))
%size(A_vector);
%size(B_vector);
%size(SH);
%size(SH');
gc_vector = 1./B_vector
h = pcolor(A_vector,gc_vector,SH');
set(h, 'EdgeColor', 'none');
colormap summer;
hold on
gf = linspace(0.0,0.99, 1000);
%plot(Anew, 1./Anew,'--','Color','r');
plot(gf, gf./(1-sqrt(gf.*(2-gf))),'Color','k','Linewidth',3);


%figure(2)
%h2 = pcolor(A_vector,gc_vector, throughput');
%colormap jet;
%%colorbar
%caxis([0 6]);
%set(h2, 'EdgeColor', 'none');
%hold on
%plot(Anew, 1./Anew,'--','Color','g');
%plot(Anew, (1-sqrt(Anew.*(2-Anew)))./Anew,'Color','k');


%pointsize = 10;
%scatter(Astart:((Aend-Astart)/Aprec):Aend, Bstart:((Bend-Bstart)/Bprec):Bend, pointsize, SH);

