%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                              %
%  This routine plots the Sound speed in the solar wind for    %
%  Te: 5-20 eV                                                 %
%                                                              %
%  Author:                                                     %
%  S. A. Ledvina, Univ. of California, Berkeley Nov. 27, 2023  %
%                                                              %
%  Notes, Comments, Changes:                                   %
%  S. A. Ledvina   11/27/23    Initial Implimentation          %
%                                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set default font size
set(0,'defaultaxesfontsize',12);

% set constants
   gamma = 1.66
   mp = 1.67e-27
   kb = 1.38064852e-23    % {m^2 kg}/{s^2 K}
   z = 1

% select Te range   
   te1 = 1:1:20;            % eV
   te = te1 * 11604;        % K

% find sound speed   
   cs = sqrt( gamma*z*kb*te/mp) / 1000.0 ; 

% Calculate trans-sonic limits lb & ub and define the array for the
% orbital speed of Venus (os).
   tt = te1./te1;
   lb = 0.8 * cs;
   ub = 1.2 * cs;
   os = 1.0 * 35.0 * tt;

% Plot Te
   plot(te1,cs,'k','linewidth',2);
   hold on

% Plot the Venus Trans-sonic envolope and the orbital speed
   plot(te1,lb,'m--','linewidth',1.5);      
   plot(te1,os,'b'  ,'linewidth',2  );        % Venus orbital speed         
   plot(te1,ub,'r--','linewidth',1.5);

% Label the plot
   xlabel('Electron Temperature (eV)')
   ylabel('Sound Speed (km/s)')
   title ('Speed of Sound in Solar Wind')

% Plot legend:  
   legend({'C_s','M_s = 0.8','M_s = 1.0','M_s = 1.2'},'location',...
           'northwest','fontsize',10)

% Scale the plot size to make room for the labels  
   scalefactor = .99;              
   g = get(gca,'Position');
   g(1:2) = g(1:2) + (1-scalefactor)/2*g(3:4);
   g(3:4) = scalefactor*g(3:4);
   set(gca,'Position',g);