kappaTestScript

kappaTestScript.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 12 23:46:13 2017

@author: Pawel M. Kozlowski
"""


import numpy as np
import astropy.units as u
import scipy.integrate as spint
import PlasmaPy.plasmapy.physics.quantum as qm
from plasmapy.physics.parameters import (thermal_speed,
                                         kappa_thermal_speed)
from plasmapy.physics.distribution import (Maxwellian_velocity_3D,
                                           Maxwellian_speed_3D,
                                           kappa_velocity_1D,
                                           kappa_velocity_3D)

#%% testing kappa
temp = 1 * u.eV
kappa = 4
particle = 'e'
vTh = (kappa_thermal_speed(temp,
                           kappa,
                           particle=particle)).si.value
vels = np.arange(-10*vTh, 10*vTh, 0.1*vTh) * (u.m/u.s)
fs = kappa_velocity_1D(vels,
                       temp,
                       kappa,
                       particle=particle,
                       V_drift=0,
                       vTh=np.nan,
                       units="units")

#plt.plot(vels, fs)
#plt.savefig('kappa1D.png')

data = np.ones((len(vels),2))
data[:,0] = vels
data[:,1] = fs
np.savetxt('kappa1D.txt', data)

#%% testing value of RMS speed numerically
def rmsSpeed(T, kappa, particle):
    vTh = kappa_thermal_speed(T,
                              kappa=kappa,
                              particle=particle).si.value
    def rmsDistFunc(vx, vy, vz):
        """
        Closure for setting up triple integral of
        distribution function times squared velocity
        """
        vSq = vx ** 2 + vy ** 2 + vz **2
        func = kappa_velocity_3D(vx,
                                 vy,
                                 vz,
                                 T,
                                 kappa,
                                 particle=particle,
                                 Vx_drift=0,
                                 Vy_drift=0,
                                 Vz_drift=0,
                                 vTh=vTh,
                                 units="unitless")
        return vSq * func
    # setting up integration from -10*vTh to 10*vTh, which is close to Inf
    infApprox = (30 * vTh)
    # integrating, this should be close to 1
    integ = spint.tplquad(rmsDistFunc,
                          -infApprox,
                          infApprox,
                          lambda z: -infApprox,
                          lambda z: infApprox,
                          lambda z, y: -infApprox,
                          lambda z, y: infApprox,
                          args=(),
                          epsabs=1e-4,
                          epsrel=1e-4,
                          )
    return integ[0] ** 0.5
T = 1.0 * u.eV
kappa = 4
particle = 'H'

vRMSInteg = rmsSpeed(T, kappa, particle)
vRMSDirect = kappa_thermal_speed(T,
                                 kappa=kappa,
                                 particle=particle,
                                 method="rms")
vRMSMaxwell = thermal_speed(T,
                            particle=particle,
                            method="rms")
print(f"Numerically obtained vRMS: {vRMSInteg}")
print(f"Functionally obtained vRMS: {vRMSDirect}")
print(f"Maxwellian vRMS: {vRMSMaxwell}")


#%% testing value of mean magnitude speed numerically
def meanSpeed(T, kappa, particle):
    vTh = kappa_thermal_speed(T,
                              kappa=kappa,
                              particle=particle).si.value
    def meanDistFunc(vx, vy, vz):
        """
        Closure for setting up triple integral of
        distribution function times squared velocity
        """
        vMag = np.sqrt(vx ** 2 + vy ** 2 + vz ** 2)
#        vMag = vx + vy + vz
        func = kappa_velocity_3D(vx,
                                 vy,
                                 vz,
                                 T,
                                 kappa,
                                 particle=particle,
                                 Vx_drift=0,
                                 Vy_drift=0,
                                 Vz_drift=0,
                                 vTh=vTh,
                                 units="unitless")
        return vMag * func
    # setting up integration from -10*vTh to 10*vTh, which is close to Inf
    infApprox = (30 * vTh)
    # integrating, this should be close to 1
    integ = spint.tplquad(meanDistFunc,
                          -infApprox,
                          infApprox,
                          lambda z: -infApprox,
                          lambda z: infApprox,
                          lambda z, y: -infApprox,
                          lambda z, y: infApprox,
                          args=(),
                          epsabs=1e-4,
                          epsrel=1e-4,
                          )
    return integ[0]
T = 1.0 * u.eV
kappa = 4
particle = 'H'

vmeanInteg = meanSpeed(T, kappa, particle)
vmeanDirect = kappa_thermal_speed(T,
                                  kappa=kappa,
                                  particle=particle,
                                  method="mean_magnitude")
vmeanMaxwell = thermal_speed(T,
                             particle=particle,
                             method="mean_magnitude")
print(f"Numerically obtained vmean: {vmeanInteg}")
print(f"Functionally obtained vmean: {vmeanDirect}")
print(f"Maxwellian vmean: {vmeanMaxwell}")