% file: tetp_equilib_bmin
define tetp_equilib_bmin ( vshock, tau, bmin )
% Returns beta = Te/Tp;
% Inputs:  vshock in km/s; tau in s/cm^3; bmin NULL or value 0.0-1.0
{
   variable tetp_atshock, tau_start, tau_equil;

   % Return an approximation to T_e/T_p based on the shock velocity
   % and the current tau of the plasma.
   
   % Start with T_e/T_p at the shock from Ghavamian+ 2007;
   % use abs(vshock) in case it is negative.
   tetp_atshock = (400.0/max([400.0,abs(vshock)]))^2;
   % do not allow it to go below the m_e/m_p value (unlikely: v_s 17000 km/s)
   tetp_atshock = max([tetp_atshock, 1.0/1836.0]);

   % Also ...  clip it to estimate based on Adelsberg+ 2008:
   % Looks like ~ 0.06 from their plot - default if not given on input
   if (bmin == NULL) bmin = 0.06;
   tetp_atshock = max([tetp_atshock, bmin]);

   % Effect of tau in causing T_e/T_p to approach 1.0,
   % have it go log-log'ly from the T_e/T_p-at-shock value
   % getting to 1.0 at larger tau value:

   % Initial relationship used:
   %%return tetp_atshock^(1.0-(min([max([log10(tau),8.5]),12.0])-8.5)/3.5);
   
   % Have the tau_equil depend on vshock:
   %  500km/s-->10.2, 3500km/s-->12.8
   % use 9.9 at 400 km/s:
   tau_equil = 9.9 + log10( (max([400.0,abs(vshock)]) / 400.0)^3.1 );
   % keep a constant log(beta)/log(tau) slope of 0.44, so set the
   % tau_start appropriately:
   tau_start = tau_equil + log10(tetp_atshock)/0.44;
   return tetp_atshock^(1.0-
      ( min([max([log10(tau),tau_start]),tau_equil]) - tau_start ) / (tau_equil-tau_start) );

}

#iffalse

% Simple commands to plot this function a la Michael '02:
variable taus, betas, vshock, it;

label("Tau (s/cm^3)","beta = T_e / T_shock",
 "T_e/T_shock vs tau, for values of v_shock (km/s)");
xrange(0.9e7,1.1e14); xlog;
yrange(1.e-3,10.1); ylog;  
set_frame_line_width(3);
set_line_width(3);

taus = 10^[7.0:14.0:#101];
betas = 0.0*taus;
it=0;

vshock=5000.0;
_for it (0,length(taus)-1,1) {betas[it]=tetp_equilib(vshock, taus[it]);};
plot(taus,betas,1);        
color(1);
xylabel(2e7, 1.2*(vshock/400.0)^(-2.0), string(vshock));

vshock=2400.0;
_for it (0,length(taus)-1,1) {betas[it]=tetp_equilib(vshock, taus[it]);};
oplot(taus,betas,2);
color(2);
xylabel(2e7, 1.2*(vshock/400.0)^(-2.0), string(vshock));

vshock=1500.0;
_for it (0,length(taus)-1,1) {betas[it]=tetp_equilib(vshock, taus[it]);};
oplot(taus,betas,8);
color(8);
xylabel(2e7, 1.2*(vshock/400.0)^(-2.0), string(vshock));

vshock=800.0;
_for it (0,length(taus)-1,1) {betas[it]=tetp_equilib(vshock, taus[it]);};
oplot(taus,betas,4);
color(4);
xylabel(2e7, 1.2*(vshock/400.0)^(-2.0), string(vshock));

#endif