%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                  %
%  This routine plots the Alfven speed in the solar wind over the  %
%  following ranges:                                               %
%  Bo:  0-20 nT                                                    %
%  ni:  5-25 cm^-3                                                 %
%                                                                  %
%  Author:                                                         %
%  S. A. Ledvina, Univ. of California, Berkeley Nov. 27, 2023      %
%                                                                  %
%  Notes, Comments, Changes:                                       % 
%  S. A. Ledvina   11/27/23    Initial Implimentation              %
%                                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

set(0,'defaultaxesfontsize',12);

% Set consts.
mi = 1.67e-27           
mu0 = 4 * 3.1415926 * 1e-7


% Set range of IMF and SW density
b0 = 1.0:1.0:20.0;    % |B| in nT
n0 = 5.0:1.0:25.0;    %  n  is cm^-3


% Set up the grid
[b1,n1] = meshgrid(b0, n0);

b = b1*1.0e-9;    % convert to mks
ni = n1*1e6;

% Calculate the Alfven speed

mu0rho = mu0*ni*mi 
  for j = 1:20
      for i = 1:21
		va(i,j) = b(i,j) / sqrt(mu0rho(i,j)) / 1000.0 ; % m/s -> km/s
      end
  end

		va1 = 20*1.0e-9/sqrt(mu0*5*1.0e6*mi) /1000
		va2 = 20*1.0e-9/sqrt(mu0*25*1.0e6*mi) /1000
		va3 = 1*1.0e-9/sqrt(mu0*5*1.0e6*mi) /1000
		va4 = 1*1.0e-9/sqrt(mu0*25*1.0e6*mi) / 1000
		va5 = 10*1.0e-9/sqrt(mu0*14*1.0e6*mi) / 1000

% Plot the Alfven speed contours
  
  pcolor(n1,b1,va)
  colorbar
  title('Alven Speed in the Solar Wind')
  xlabel('Proton Density (cm^-^3)')
  ylabel('Magnetic Field Strength (nT)')

hold on

% Set and plot the Orbital speed and the Trans-Alfvenic envolope  
  lb = 0.8 * 35;
  ub = 1.2 * 35;
  llb = [lb,lb];
  uub = [ub,ub];
  vv = [35,35]
    contour(n1,b1,va,llb,'c--','linewidth',1.8)
    contour(n1,b1,va,vv ,'m'  ,'linewidth',1.8)
    contour(n1,b1,va,uub,'c--','linewidth',1.8)
    [c0,h0] = contour(n1,b1,va)
    clabel(c0,h0,'FontSize',20,'Color','c')


% Label the Alfven speed contours
    contour(n1,b1,va,[20,40,60,80,100,120,140,160],'w:','linewidth',1.8)
    [c,h] = contour(n1,b1,va)
    clabel(c,h,'FontSize',20,'Color','w')

% Label the Orbital speed and Trans-Alfvenic envolope contours
text(23,7.8,'35','FontSize',20,'Color','m','FontWeight','bold')
text(20,6,'28','FontSize',20,'Color','c','FontWeight','bold')    
text(20,9,'42','FontSize',20,'Color','c','FontWeight','bold')