In this problem we are given an input file which contains the solar wind protons velocity space distributions obtained from a particle simulation code for an observer located at 1.0 AU distance from the Sun. Simulated observer collects only the solar wind protons with directional velocity range from 0 to 8 *10^5 m/s and velocity resolution of 5*10^3 m/s in all three directions. The magnetic field is constant +5 nT in a with a magnetic vector along Y axis, and the solar wind flows along X axis.

clc;clear all;close all;

Constants

mp   = 1.6726E-27;      % Mass of a proton [kg]
kb   = 1.3806503e-23;   % Boltzmann constant [J/kg]

Load the file with wind protons velocity space distributions into MATLAB The velocity space distribution values are stored in ‘f’

load data
dv = 5e3;               %Velocity resolution in all directions [m/s]
% Velocity space [m/s]
v_x = -8E5:dv:8E5;
v_y = -8E5:dv:8E5;
v_z = -8E5:dv:8E5;

Solar wind protons number density,Zero moment

ns = sum(f(:))*(dv^3);
fprintf('Solar wind protons number density %4.2d m^-3 \n\n',ns);

Solar wind protons number density 3.00e+006 m^-3

Solar wind protons fluxes in all directions,First moment

Fi_x = 0; Fi_y = 0; Fi_z = 0;
len = length(v_x);
 for i = 1:len
       for j = 1:len
            for k = 1:len
                  Fi_x = Fi_x + v_x(i)*f(i,j,k);
                  Fi_y = Fi_y + v_y(j)*f(i,j,k);
                  Fi_z = Fi_z + v_z(k)*f(i,j,k);
            end
         end
    end
Fi_x = Fi_x*dv^3; Fi_y = Fi_y*dv^3; Fi_z = Fi_z*dv^3;

fprintf('Solar wind protons flux in X direction %4.2d m-2/s\n',Fi_x);
fprintf('Solar wind protons flux in Y direction %4.2d m-2/s\n',Fi_y);
fprintf('Solar wind protons flux in Z direction %4.2d m-2/s\n\n',Fi_z);
Fi = [Fi_x,Fi_y,Fi_z];

Solar wind protons flux in X direction -1.09e+012 m-2/s
Solar wind protons flux in Y direction 7.58e-007 m-2/s
Solar wind protons flux in Z direction -6.91e+004 m-2/s

Solar wind bulk flow velocity

Urt = Fi/ns/1000;       %[km/s]
fprintf('Solar wind bulk flow velocity U [%4.2f %4.2f %4.2f] km/s\n',Urt);

Solar wind bulk flow velocity U [-365.00 0.00 -0.00] km/s

Distribution function plot

Plot the distribution function in the vx – vy plane and compare it with distribution function in the vx – vz plane. The center of the distribution contours should clearly show the bulk velocity in each plane.

v_x = v_x/1000; v_y = v_y/1000; v_z = v_z/1000;
% Maximum of velocity space distribution function
[f_max,ind] = max(f(:));
size_f = [321 321 321];
% Location of maximum
[f_mi,f_mj,f_mk] = ind2sub(size_f,ind);

Distribution function in the vx – vy plane

figure(1);
hold on;grid on;
contour(v_x, v_y, f(:,:,f_mk) );
scatter(Urt(2),Urt(1),'*r');
% Matlab will auto adjust the axis and ellipse will appear as circle
% We adjust the axis manually to see the contour shape
axis([-120,120,-480,-240]);
title('Distribution function in the vx - vy plane');
ylabel('v_x [km/s]');
xlabel('v_y [km/s]');
hold off;

Distribution function in the vx – vz plane

figure(2);
hold on;grid on;
contour(v_x, v_z, squeeze(f(:,f_mj,:)) );
axis([-120,120,-480,-240]);
scatter(Urt(3),Urt(1),'*r');
title('Distribution function in the vx - vz plane');
ylabel('v_x [km/s]');
xlabel('v_z [km/s]');
hold off;

The type of the distribution function is Bi-Maxwellian, because the shape in vx – vy is ellipse.

Pressure tensor and scalar pressure

C = V – U and called random velocity

Ux = Urt(1);Uy = Urt(2);Uz = Urt(3);
% To save calculation time we use scalar variables instead of matrix
P11 = 0;P22 = 0; P33 = 0;
P12 =0; P23 = 0; P13 = 0;
for i = 1:len
       for j = 1:len
            for k = 1:len
                  c1 = v_x(i) - Ux;
                  c2 = v_y(j) - Uy;
                  c3 = v_z(k) - Uz;
                  P11 = P11 + c1^2*f(i,j,k);
                  P22 = P22 + c2^2*f(i,j,k);
                  P33 = P33 + c3^2*f(i,j,k);
                  P12 = P12 + c1*c2*f(i,j,k);
                  P23 = P23 + c2*c3*f(i,j,k);
                  P13 = P13 + c1*c3*f(i,j,k);
             end
        end
 end
P = [P11,P12,P13;
     P12,P22,P23;
     P13,P23,P33];
 P = P*mp*dv^3*1e6      % [kg*m^3*m^2/s^2]

P =

  1.0e-011 *

    0.3106   -0.0000    0.0000
   -0.0000    0.9319    0.0000
    0.0000    0.0000    0.3106

Scalar pressure

The scalar pressure ps is defined as one third of the trace of Pij

ps = 1/3*(P(1,1) + P(2,2) + P(3,3));
fprintf('Scalar pressure %4.2d nPa\n',ps*1e9)

Scalar pressure 5.18e-003 nPa

Solar wind thermal velocity and the solar wind kinetic temperature

Kinetic temperature of species can be given in a form of a 3 × 3 matrix

Ts = P/(ns*kb)      % K

Ts =

  1.0e+005 *

    0.7500   -0.0000    0.0000
   -0.0000    2.2500    0.0000
    0.0000    0.0000    0.7500

Thermal Velocity

Vth = ([Ts(1,1),Ts(2,2),Ts(3,3)].*(2*kb/mp)).^0.5;
fprintf('Thermal Velocity Vth [ %4.2f %4.2f %4.2f ] km/s\n',Vth/1000)
V_sp = norm(Vth) ;

Thermal Velocity Vth [ 35.19 60.95 35.19 ] km/s