########################################
Params file information

All scripts require a runname_params.txt file to exist in the data directory. This file contains all necessary information about the run to compute things correctly and produce normalized quantities. See example_params.txt for the expected structure. IMPORTANT NOTE the formatting of the name of the params file is important and must follow the naming format of your data files. For instance, g0 data use a - rather than _ between the base filename and elc/ion/field/etc. So, for these cases, the params file must be named runname-params.txt


##########################################
Table of choices for varid. i,j, k = x,y,z. spec = species as in filename
Note that case is ignored
-Particle distribution: dist_spec, specify dist0_spec or dist1_spec for PKPM
-EM fields: ei or bi
-Scalar potentials: phi, psi 
-Magnitude E or B: mage or magb
-Div E or B: dive or divb
-Density: n_spec
-Flow velocity: ui_spec
-Stream function: stream_spec
-Eparallel: epar
-Current: ji
-Species current: ji_spec
-Parallel current: jpar
-Species parallel current: jpar_spec
-J.E, Work: work
-Parallel work, Jpar.Epar: workpar
-J_spec.E, species work: work_spec
-Species parallel work, Jpar_spec.Epar: workpar_spec
-Component work, JiEi: worki
-Species component work, Ji_spec Ei: worki_spec
-Pressure: pij_spec NOTE: Vlasov data with only M2 data and all 5M data will always return trace pressure / dimV. For 5M data, p_spec is preferred, but pii_spec will return the same value.
-Parallel pressure: ppar_spec
-Perp pressure: pperp_spec
-Tr(P): trp_spec
-Temperature: temp_spec
-Parallel temp: temppar_spec
-Perp temp: tempperp_spec
-Tperp / Tpar: tempperppar_spec
-Tpar / Tperp: tempparperp_spec
-Pressure agyrotropy: agyro_spec
-Heat flux: qi_spec for M3i data or qijk_spec for M3ijk data. DO NOTE USE params["restFrame"] = 1 for M3i data if M2ij data does not exist!  
-Heat flux continued: for PKPM, qpar_spec or qperp_spec
-Fourth moment, for PKPM only: rparperp_spec or rperpperp_spec
-Beta: beta_spec
-Beta par or perp: betapar_spec, betaperp_spec
-Magnetic moment, p_perp_spec / B n_spec: mu_spec
-Gyroradius, sqrt(2 T_perp / m) / (|q| B / m): rho_spec
-Inertial length, c / omega_P: inertiallength_spec
-ExB drift component: exbDrifti
-Poynting vector: poyntingi or poynting for |ExB|
-Cross helicity, u_spec.b / .5(u_spec^2 + b^2): crossHelicity_spec
-Firehose instability threshold (pperp_spec - ppar_spec)/ppar_spec + 2./(betapar_spec) < 0 unstable: firehose_spec
-Mirror instability threshold (ppar_spec - pperp_spec)/ppar_spec + 1./(betaperp_spec) < 0 unstable: mirror_spec

--PiD related quantities, see Cassak PoP 29, 122306 (2022) (https://doi.org/10.1063/5.0125248)
-Pressure dilatation, -pTheta (Tr(p_spec) / 3)*div u_spec: ptheta_spec
-PiD, -Pi_jk_spec*D_jk_spec: pid_spec
-PiD normal component, equation 9a of Cassak: pidnormal_spec
-PiD shear component, equation 9b of Cassak: pidshear_spec


--Drifts: Add work to the name to calculate q n E dotted into work, e.g., curvatureDriftWork_elc
-Curvature drift component: curvatureDrifti_spec (F13 term 4 in Juno et al 2021)
-Curvature drift component: curvatureDriftv2i_spec (F13 only ppar part of term 4 in Juno 2021)
-Curvature drift component: curvatureDriftv0i_spec (Guiding center curvature drift)

-GradB drift: gradBDrifti_spec (F13 term 3 in Juno et al 2021)
-Magnetiziation drift: magDrifti_spec (F13 term 2 in Juno et al 2021)
-Agyrotropic drift: agyroDrifti_spec (F13 term 5 in Juno et al 2021)
-Diamagnetic drift: diamagDrifti_spec (F14 in Juno et al 2021)
-Betatron drift: betaDrifti_spec (F13 only pperp part of term 4 in Juno et al 2021)

-Dispersion relation data from linear solver: dispersion


###############################################
General use instructions
First, one must import gkData
from utils import gkData

See plotExample.py for basic usage example. The possible varid choices are defined at the top of this readme. Several other example input files are also included to demonstrate other use cases.

Once a gkData object has been read in (var in the plotExample.py script), you can access the following

var.data = The data array corresponding to the varid choice
var.coords = The coordinate array in code units, e.g., x, y, z
var.dx = The axis spacing array in code units
var.max = max(var.data)
var.min = min(var.data)
var.time = time in code units
var.mu = Species mass array, order defined by params file
var.n = Reference species density defined in params file
var.B0 = Reference B0 defined by params file
var.q = Species charge array
var.beta = species beta array
var.vA = species Alfven speed array
var.vt = species thermal speed array (vt = sqrt(2T_spec/m_spec))
var.omegaC = species gyrofrequency array
var.omegaP = species plasma frequency array
var.d = species inertial length array, c / omega_p
var.rho = species gyroradius array (rho = vt / omegaC)
var.debye = species debye array
var.mu0 = mu0
var.eps0 = eps0
var.c = speed of light

As an example, if you want your x axis in units of the inertial length of the first species in the params file, simply do var.coords[0] / var.d[0]

You can also access all of the params[*] variables, e.g.,
var.params["axesNorm"] = Array of axis normilization optionally defined by user input.

Plotting a gkData object directly can be done as plt.gkPlot(var) starting with the import command
from utils import gkPlot as plt

See plotExample.py for the various plotting related params[*] that gkPlot will use.
Note gkPlot takes optional arguments
show = 0 or 1 to display the figure
save = 0 or 1 to save the figure

var.integrate() behaves like the python integrate function on a gkData object, e.g.,
var.integrate() integrates over all axes, var.integrate(axis=0) integrates over the first axis, var.integrate(axis=(0,2)) integrates over axes 0 and 2. Note that axes integrated over are squeezed from the data array, and the corresponding coords components are removed. 

Simple operations can also be performed on gkData objects, e.g., for objects A, B, C
D = 2*A*B + C**2 will return a new gkData object D. 

###############################################
Field-particle correaltion (FPC) usage
from utils import FPC as FPC to import FPC diagnostics

[coords, FPC, t] = FPC.computeFPC(paramFile,ts,spec,params) 
returns the coords read in, the FPC, and the time (in code units) of the FPC
FPC = -0.5*q*v^2*dfdv*E

[coords, FPC, t] = FPC.computeFPCPKPM(paramFile,ts,spec,params) 
returns the coords read in, the kinetic pressure work term, and the time (in code units) of the FPC
FPC = 0.5*m_s vpar^2 (bb:grad u) d (vpar f) / d vpar

This works for arbitrary dimensionality


Options for the FPC
spec = species name as in params file

params["fieldAlign"] = 0 or 1 #Align FPC to the local magnetic field. Only use for 3V data.

params["driftAlign"] = 'curvatureDrift' #Rotate FPC with B and drift (see varid list for drifts). Only use for 3V data.

params["frameXform"] = [ux,uy,uz] #Transform frames, including electric field. This can only be used for single point FPC. Note that this Xform depends on time and so must be defined at each call for FPC.computeFPC at different times.

If placed in a time loop, one can do a time average easily using FPC.computeFPCAvg(fpc, nTau, avgDir). For example,

fpc  = []
for it in range(nt):
        print('Working on frame {0} of {1}'.format(it+1,nt))
        [coords, fpcTmp, t[it]] = FPC.computeFPC(paramFile,ts[it],spec,params)
        fpc.append(fpcTmp)
if nTau > 0:
	fpc = FPC.computeFPCAvg(fpc, nTau, avgDir)

Here
nTau = 3 #Frames over which to average. 0 or 1 does no averaging. Note centered ==> nTau must be odd and >= 3
avgDir = 0 #Backwards (-1), forward (1), or centered (0). Endpoints treated with telescoping windows.

For a complete example, see the included FPCExample_0xNv.py.


###################################################
Entropy usage
from utils import entropy

[coords, entropy, t] = entropy.getEntropy(paramFile,fileNum,spec,params,type=0)
returns the coords read in, the type of entropy density chosen, and the time (in code units)

Options for the FPC
spec = species name as in params file
type = the type of entropy 
0=conventional entropy -Int f ln f dv
1=relatitive entropy  - Int f ln (f / fMaxwelian) dv
2=position space entropy - n ln n 
3=velocity space entropy -Int f ln (f / n) dv
Note that type0 = type2 + type3

params["absVal"] = applies absolute value to the argument of the ln function. Sometimes necessary for non-positive values of f.