mirror of
https://github.com/JHUAPL/kaipy.git
synced 2026-01-09 14:28:02 -05:00
Merged development into ewinter-doc_fixes
This commit is contained in:
Eric Winter
@@ -91,7 +91,7 @@ Mode: CLI
|
||||
---------
|
||||
The Kaipy package also comes with a command line interface that allows you to run scripts to analyze MAGE model data. The CLI is a great way to automate the analysis of large datasets. The CLI is run from the terminal and has a variety of options to customize the analysis.
|
||||
|
||||
A complete list of the available scripts can be found at the `Scripts documentation <https://kaipy.readthedocs.io/en/latest/scripts.html>`_.
|
||||
A complete list of the available scripts can be found at the `Scripts documentation <https://kaipy-docs.readthedocs.io/en/latest/scripts.html>`_.
|
||||
|
||||
The quicklook directory has numerous scripts that can be used to generate plots and movies of the MAGE model output. For example the `msphpic.py` command makes a summary movie of the magnetosphere while the `mixpic.py` command makes a summary movie of the ionosphere.
|
||||
|
||||
|
||||
@@ -181,8 +181,55 @@ def genFatEgg(Ni=Ni0,Nj=Nj0,Rin=3.0,Rout=30.0,xtail=250,NumSph=5,TINY=1.0e-8,A=0
|
||||
tuple: A tuple containing the x-coordinates and y-coordinates of the grid points.
|
||||
|
||||
"""
|
||||
# Rest of the code...
|
||||
#Get ellipse parameters
|
||||
x0,a,b = Egglipses(Ni,Nj,Rin,Rout,xtail,NumSph)
|
||||
|
||||
xSun = np.max(x0+a) #Forward boundary
|
||||
xBack = np.abs(np.min(x0-a)) #Back boundary
|
||||
|
||||
#Eta is index (0,1) for theta mapping
|
||||
eta = np.linspace(0+TINY,1.0-TINY,Nj+1)
|
||||
L = np.pi
|
||||
xC = 0.5*np.pi #Any warping is symmetric
|
||||
|
||||
#theta = np.linspace(0+TINY,np.pi-TINY,Nj+1)
|
||||
xx = np.zeros((Ni+1,Nj+1))
|
||||
yy = np.zeros((Ni+1,Nj+1))
|
||||
tau0 = 0.35
|
||||
|
||||
di = Ni/8
|
||||
i1 = Ni-di
|
||||
|
||||
AScl = 3*A
|
||||
|
||||
for i in range(Ni+1):
|
||||
#Calculate theta profile for this shell
|
||||
xTi = np.abs(x0[i]-a[i])
|
||||
|
||||
theta = L*eta + A*(xC-L*eta)*(1-eta)*eta #Lines of constant phi
|
||||
if (xTi>xSun):
|
||||
#Outer region, keep tailward near-axis region packed tightly
|
||||
Ai = A - (A-AScl)*RampUp(xTi,xSun,xBack)
|
||||
thTail = L*eta + Ai*(xC-L*eta)*(1-eta)*eta
|
||||
#Only use back half (past half-way)
|
||||
theta[Nj//2:] = thTail[Nj//2:]
|
||||
|
||||
for j in range(Nj+1):
|
||||
th = theta[j]
|
||||
Ae = (np.cos(th)/a[i])**2.0 + (np.sin(th)/b[i])**2.0
|
||||
Be = -2*np.cos(th)*x0[i]/a[i]**2.0
|
||||
Ce = x0[i]**2.0/a[i]**2.0 - 1.0
|
||||
r = (-Be+np.sqrt(Be*Be-4*Ae*Ce))/(2*Ae)
|
||||
|
||||
xx[i,j] = r*np.cos(th)
|
||||
|
||||
taui = min(tau0*RampUp(1.0*i,i1,di),1.0)
|
||||
tauscl = taui*(xx[i,j])/a[i]
|
||||
tau = min(1.0+tauscl,1)
|
||||
yy[i,j] = r*np.sin(th)/np.sqrt(tau)
|
||||
|
||||
|
||||
return xx,yy
|
||||
|
||||
def RampUp(r, rC, lC):
|
||||
"""
|
||||
|
||||
@@ -2,11 +2,97 @@ import h5py
|
||||
import numpy as np
|
||||
import os, sys, subprocess
|
||||
from alive_progress import alive_bar, alive_it
|
||||
import datetime
|
||||
from astropy.time import Time
|
||||
#Import kaipy modules
|
||||
import kaipy.kdefs as kdefs
|
||||
import kaipy.kaiTools as ktools
|
||||
|
||||
|
||||
#------
|
||||
# Data Containers
|
||||
#------
|
||||
|
||||
class H5Info(object):
|
||||
"""
|
||||
Class to hold model data information from h5 output files
|
||||
|
||||
Args:
|
||||
h5fname (str): .h5 filename to scrape info from
|
||||
noSubsec (bool): Whether or not UTs should include subseconds (optional, default: True)
|
||||
|
||||
Attributes:
|
||||
fname (str): .h5 filename the rest of the info came from
|
||||
steps (List[int]): list of step numbers
|
||||
stepStrs (List[str]): list of step strings in Step#X h5.Group format
|
||||
Nt (int): Number of time steps
|
||||
times (List[float]): List of times corresponding to each step
|
||||
Units: Unchanged from file, likely seconds
|
||||
MJDs (List[float]): List of MJDs corresponding to each step
|
||||
UTs (List[float]): List of UTs corresponding to each step
|
||||
May or may not include subseconds based on constructor args
|
||||
"""
|
||||
|
||||
def __init__(self, h5fname: str, noSubsec:bool=True):
|
||||
# h5fname = h5fname.split('/')[-1]
|
||||
self.fname = h5fname
|
||||
self.Nt, self.steps = cntSteps(self.fname)
|
||||
self.stepStrs = ['Step#'+str(s) for s in self.steps]
|
||||
self.times = getTs(self.fname, self.steps, "time")
|
||||
self.MJDs = getTs(self.fname, self.steps, "MJD" )
|
||||
|
||||
f5 = h5py.File(h5fname)
|
||||
if noSubsec:
|
||||
ut_arr = np.zeros(self.Nt, dtype=datetime.datetime)
|
||||
for i,mjd in enumerate(self.MJDs):
|
||||
UTStr = Time(mjd,format='mjd').isot
|
||||
utTrim = UTStr.split('.')[0]
|
||||
UT = datetime.datetime.strptime(utTrim,'%Y-%m-%dT%H:%M:%S')
|
||||
ut_arr[i] = UT
|
||||
else:
|
||||
ut_arr = ktools.MJD2UT(self.MJDs)
|
||||
|
||||
self.UTs = ut_arr
|
||||
|
||||
f5.close()
|
||||
|
||||
def printStepInfo(self) -> None:
|
||||
"""
|
||||
Prints a summary of step info for the contained data
|
||||
|
||||
Args: None
|
||||
|
||||
Returns: None
|
||||
"""
|
||||
print("Step Info for {}:".format(self.fname))
|
||||
print(" Step start/end/stride : {} / {} / {}".format(self.steps[0], self.steps[-1], self.steps[1]-self.steps[0]))
|
||||
print(" Time start/end/stride : {:1.2f} / {:1.2f} / {:1.2f}".format(self.times[0], self.times[-1], self.times[1]-self.times[0]))
|
||||
print(" MJD start/end : {:1.8f} / {:1.8f}".format(self.MJDs[0], self.MJDs[-1]))
|
||||
print(" UT start/end : {} / {}".format(self.UTs[0], self.UTs[-1]))
|
||||
|
||||
|
||||
class TPInfo(H5Info):
|
||||
"""
|
||||
Extension of H5Info for test particle files
|
||||
|
||||
Attributes:
|
||||
Ntp (int): Number of test particles
|
||||
id2idxMap (Dict{int:}): Dict to help map TP id to its index
|
||||
"""
|
||||
|
||||
def __init__(self, fname:str, noSubsec:bool=True):
|
||||
super().__init__(fname, noSubsec)
|
||||
with h5py.File(self.fname) as tp5:
|
||||
ids = tp5[self.stepStrs[0]]['id']
|
||||
self.Ntp = ids.shape[0]
|
||||
self.id2idxMap = {}
|
||||
for i in range(self.Ntp):
|
||||
self.id2idxMap[int(ids[i])] = int(i)
|
||||
|
||||
#------
|
||||
# General functions
|
||||
#------
|
||||
|
||||
#Generate MPI-style name
|
||||
def genName(bStr, i, j, k, Ri, Rj, Rk, nRes=None):
|
||||
'''
|
||||
|
||||
@@ -62,6 +62,37 @@ def MJD2UT(mjd):
|
||||
else:
|
||||
return [datetime.datetime.strptime(UT[n], isotfmt) for n in range(len(UT))]
|
||||
|
||||
def utIdx(utList,ut):
|
||||
"""
|
||||
Finds the index of the closest datetime within an array of datetimes
|
||||
|
||||
Parameters:
|
||||
utList (List[datetime.datetime]): List of datetimes to find index within
|
||||
ut (datetime.datetime): datetime to find closest index to
|
||||
|
||||
Returns:
|
||||
Single integer value of the closest index within 'utList' to 'ut'
|
||||
"""
|
||||
return np.array([np.abs((utList[n]-ut).total_seconds()) for n in range(len(utList))]).argmin()
|
||||
|
||||
def pntIdx_2D(X2D, Y2D, pnt):
|
||||
"""
|
||||
Finds i, j, index within 'X2D' and 'Y2D' closest to 'pnt'
|
||||
Assumes cartesian coordinate system
|
||||
|
||||
Parameters:
|
||||
X2D (ndarray[float]): 2D array of X values
|
||||
Y2D (ndarray[float]): 2D array of Y values
|
||||
pnt (List[float]) : [x, y] point to find closest index to
|
||||
|
||||
Returns:
|
||||
Tuple of closest i, j index to point 'pnt'
|
||||
"""
|
||||
distSq = (X2D-pnt[0])**2 + (Y2D-pnt[1])**2 # Each source point's euclidian distance from 'pnt'
|
||||
aMinFlattened = distSq.argmin() # Index of min distance, if it was a flattened array
|
||||
i, j = np.unravel_index(aMinFlattened, X2D.shape) # Convert back into i,j location
|
||||
return i, j
|
||||
|
||||
def getRunInfo(fdir, ftag):
|
||||
"""
|
||||
Get information about the run.
|
||||
|
||||
@@ -696,6 +696,32 @@ def get_aspect(ax):
|
||||
|
||||
return disp_ratio #/ data_ratio
|
||||
|
||||
def setBndsByAspect(ax, bnds, axis='x', aspect=1):
|
||||
"""
|
||||
Given bounds of one axis, set the other bounds of Axes object such that aspect ratio and Axes box size is perserved
|
||||
|
||||
Parameters:
|
||||
ax: Axes object to set bounds for
|
||||
bnds (List[float]): [axMin, axMax] min/max bounds for fixed axis
|
||||
axis (str): ('x' or 'y') axis that fixed bounds are applied to
|
||||
aspect: x:y aspect ratio to preserve
|
||||
"""
|
||||
boxRatio = get_aspect(ax)
|
||||
da = bnds[1] - bnds[0]
|
||||
|
||||
if axis=='x':
|
||||
db = aspect*boxRatio*da
|
||||
bMin = -db/2
|
||||
bMax = db/2
|
||||
ax.set_xlim(bnds)
|
||||
ax.set_ylim([bMin, bMax])
|
||||
else:
|
||||
db = 1/aspect*boxRatio*da
|
||||
bMin = -db/2
|
||||
bMax = db/2
|
||||
ax.set_xlim([bMin, bMax])
|
||||
ax.set_ylim(bnds)
|
||||
|
||||
|
||||
# Methods specific for comparing gamhelio output to data from heliospheric
|
||||
# spacecraft.
|
||||
|
||||
@@ -114,7 +114,7 @@ def getRootVars(fname, gDims):
|
||||
doV = True
|
||||
if ("Step" in vID or kdefs.grpTimeCache in vID):
|
||||
doV = False
|
||||
if ((vID == "X") or (vID == "Y") or (vID == "Z")):
|
||||
if ((vID == "X") or (vID == "Y") or (vID == "Z") or (type(hf[k]) != h5py._hl.group.Group)): # Ignore root geom and root Groups
|
||||
doV = False
|
||||
if (doV):
|
||||
Nv = hf[k].shape
|
||||
|
||||
@@ -90,8 +90,8 @@ class AlamParams:
|
||||
# These can help to determine some things for higher-level lambda generators
|
||||
#tiote: Optional[float] = None # Ratio of ion temperature to electron temperature
|
||||
#ktMax: Optional[float] = None # Energy in eV
|
||||
emine: Optional[float]
|
||||
eminp: Optional[float]
|
||||
emaxe: Optional[float]
|
||||
emaxp: Optional[float]
|
||||
L_kt: Optional[float]
|
||||
emine: Optional[float] = None
|
||||
eminp: Optional[float] = None
|
||||
emaxe: Optional[float] = None
|
||||
emaxp: Optional[float] = None
|
||||
L_kt : Optional[float] = None
|
||||
|
||||
@@ -562,8 +562,64 @@ class remix:
|
||||
if not self.Initialized:
|
||||
sys.exit("Variables should be initialized for the specific hemisphere (call init_var) prior to efield calculation.")
|
||||
|
||||
# Rest of the code...
|
||||
|
||||
Psi = self.variables['potential']['data'] # note, these are numbers of cells. self.ion['X'].shape = Nr+1,Nt+1
|
||||
Nt,Np = Psi.shape
|
||||
|
||||
# Aliases to keep things short
|
||||
x = self.ion['X']
|
||||
y = self.ion['Y']
|
||||
|
||||
# note the change in naming convention from above
|
||||
# i.e., theta is now the polar angle
|
||||
# and phi is the azimuthal (what was theta)
|
||||
# TODO: make consistent throughout
|
||||
theta = np.arcsin(self.ion['R'])
|
||||
phi = self.ion['THETA']
|
||||
|
||||
# interpolate Psi to corners
|
||||
Psi_c = np.zeros(x.shape)
|
||||
Psi_c[1:-1,1:-1] = 0.25*(Psi[1:,1:]+Psi[:-1,1:]+Psi[1:,:-1]+Psi[:-1,:-1])
|
||||
|
||||
# fix up periodic
|
||||
Psi_c[1:-1,0] = 0.25*(Psi[1:,0]+Psi[:-1,0]+Psi[1:,-1]+Psi[:-1,-1])
|
||||
Psi_c[1:-1,-1] = Psi_c[1:-1,0]
|
||||
|
||||
# fix up pole
|
||||
Psi_pole = Psi[0,:].mean()
|
||||
Psi_c[0,1:-1] = 0.25*(2.*Psi_pole + Psi[0,:-1]+Psi[0,1:])
|
||||
Psi_c[0,0] = 0.25*(2.*Psi_pole + Psi[0,-1]+Psi[0,0])
|
||||
Psi_c[0,-1] = 0.25*(2.*Psi_pole + Psi[0,-1]+Psi[0,0])
|
||||
|
||||
# fix up low lat boundary
|
||||
# extrapolate linearly just like we did for the coordinates
|
||||
# (see genOutGrid in src/remix/mixio.F90)
|
||||
# note, neglecting the possibly non-uniform spacing (don't care)
|
||||
Psi_c[-1,:] = 2*Psi_c[-2,:]-Psi_c[-3,:]
|
||||
|
||||
# now, do the differencing
|
||||
# for each cell corner on the original grid, I have the coordinates and Psi_c
|
||||
# need to find the gradient at cell center
|
||||
# the result is the same size as Psi
|
||||
|
||||
# first etheta
|
||||
tmp = 0.5*(Psi_c[:,1:]+Psi_c[:,:-1]) # move to edge center
|
||||
dPsi = tmp[1:,:]-tmp[:-1,:]
|
||||
tmp = 0.5*(theta[:,1:]+theta[:,:-1])
|
||||
dtheta = tmp[1:,:]-tmp[:-1,:]
|
||||
etheta = dPsi/dtheta/ri # this is in V/m
|
||||
|
||||
# now ephi
|
||||
tmp = 0.5*(Psi_c[1:,:]+Psi_c[:-1,:]) # move to edge center
|
||||
dPsi = tmp[:,1:]-tmp[:,:-1]
|
||||
tmp = 0.5*(phi[1:,:]+phi[:-1,:])
|
||||
dphi = tmp[:,1:]-tmp[:,:-1]
|
||||
tc = 0.25*(theta[:-1,:-1]+theta[1:,:-1]+theta[:-1,1:]+theta[1:,1:]) # need this additionally
|
||||
ephi = dPsi/dphi/np.sin(tc)/ri # this is in V/m
|
||||
|
||||
if returnDeltas:
|
||||
return (-etheta,-ephi,dtheta,dphi) # E = -grad Psi
|
||||
else:
|
||||
return (-etheta,-ephi) # E = -grad Psi
|
||||
|
||||
def joule(self):
|
||||
"""
|
||||
|
||||
Reference in New Issue
Block a user