Merged in development (pull request #26)

Development
2026-01-09 15:58:04 -05:00 · 2025-12-10 16:24:41 +00:00
parent 4e98d39444 21c2150138
commit 01b89c8a4a
18 changed files with 1419 additions and 96 deletions
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -0,0 +1,64 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it using the metadata below."
+title: "Kaipy: Python package for analysis and visualization of MAGE and space weather modeling data"
+version: "1.1.2"
+date-released: 2025-07-03
+doi: 10.5281/zenodo.15801041  
+repository-code: "https://github.com/JHUAPL/kaipy"
+license: "BSD 3-Clause"
+
+authors:
+  - family-names: Wiltberger
+    given-names: Michael
+    orcid: https://orcid.org/0000-0002-4844-3148
+    affiliation: "NSF National Center for Atmospheric Research"
+  - family-names: Bao
+    given-names: Shanshan
+    orcid: https://orcid.org/0000-0002-5209-3988
+    affiliation: "Rice University"
+  - family-names: Garretson
+    given-names: Jeffery
+    orcid: https://orcid.org/0000-0003-3805-9860
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: McCubbin
+    given-names: Andrew
+    orcid: https://orcid.org/0000-0002-6222-3627
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: Merkin
+    given-names: Slava
+    orcid: https://orcid.org/0000-0003-4344-5424
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: Michael
+    given-names: Adam
+    orcid: https://orcid.org/0000-0003-2227-1242
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: Pham
+    given-names: Kevin
+    orcid: https://orcid.org/0000-0001-5031-5519
+    affiliation: "NSF National Center for Atmospheric Research"
+  - family-names: Provornikova
+    given-names: Elena
+    orcid: https://orcid.org/0000-0001-8875-7478
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: Rao
+    given-names: Nikhil
+    affiliation: "NSF National Center for Atmospheric Research"
+    orcid: https://orcid.org/0000-0003-2639-9892
+  - family-names: Sciola
+    given-names: Anthony
+    orcid: https://orcid.org/0000-0002-9752-9618
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: Sorathia
+    given-names: Kareem
+    orcid: https://orcid.org/0000-0002-6011-5470
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+  - family-names: Winter
+    given-names: Eric
+    orcid: https://orcid.org/0000-0001-5226-2107
+    affiliation: "Johns Hopkins University Applied Physics Laboratory"
+
+keywords:
+  - "space weather"
+  - "MAGE"
+  - "kaipy"
+  - "geospace modeling"
--- a/kaipy/gamera/gampp.py
+++ b/kaipy/gamera/gampp.py
@@ -275,7 +275,7 @@ class GameraPipe(object):

 			for data,vID in coords:
 				datasets = []
-				with alive_bar(NrX,title=f"{titStr}/{vID}".ljust(kdefs.barLab),length=kdefs.barLen,bar=kdefs.barDef) as bar, \
+				with alive_bar(NrX,title=f"{titStr}/{vID}".ljust(kdefs.barLab),length=kdefs.barLen,bar=kdefs.barDef2) as bar, \
 					ProcessPoolExecutor(max_workers=self.nWorkers) as executor:
 					futures = [executor.submit(kh5.PullVarLoc, fIn, vID, loc=loc) for fIn, loc in files]
 					for future in as_completed(futures):
@@ -328,7 +328,7 @@ class GameraPipe(object):
 					else:
 						titStr = None
 					NrX = max(self.Nr,1)
-					with alive_bar(NrX,title=titStr,length=kdefs.barLen) as bar:
+					with alive_bar(NrX,title=titStr,length=kdefs.barLen,bar=kdefs.barDef2) as bar:
 						for (i,j,k) in itertools.product(range(self.Ri),range(self.Rj),range(self.Rk)):
 							iS = i *self.dNi
 							jS = j *self.dNj
@@ -454,7 +454,7 @@ class GameraPipe(object):
 			else:
 				titStr = ''
 			NrX = max(self.Nr,1)
-			with alive_bar(NrX,title=titStr.ljust(kdefs.barLab),length=kdefs.barLen,disable=not doVerb) as bar:
+			with alive_bar(NrX,title=titStr.ljust(kdefs.barLab),length=kdefs.barLen,bar=kdefs.barDef,disable=not doVerb) as bar:
 				for (i,j,k) in itertools.product(range(self.Ri),range(self.Rj),range(self.Rk)):

 					iS = i*self.dNi
--- a/kaipy/gamera/msphViz.py
+++ b/kaipy/gamera/msphViz.py
@@ -74,10 +74,10 @@ def GetSizeBds(args):

 	return xyBds

-#Plot absolute error in the requested, or given, equatorial field
-def PlotEqErrAbs(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=True, doDeco=True, vMin=1e-9, vMax=1e-4, doLog=True, doVerb=True):
+#Plot absolute error in the requested, or given, equatorial or meridional field
+def PlotErrAbs(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=True, doDeco=True, vMin=1e-9, vMax=1e-4, doLog=True, doVerb=True, doEq=True):
 	"""
-	PlotEqErrAbs function plots the absolute error between two gsph objects.
+	PlotErrAbs function plots the absolute error between two gsph objects.

 	Args:
 		gsphP (gsph): The gsph object representing the predicted values.
@@ -93,6 +93,7 @@ def PlotEqErrAbs(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=T
 		vMax (float, optional): The maximum value for the colorbar. (default: 1e-4)
 		doLog (bool, optional): Whether to use logarithmic scale for the colorbar. (default: True)
 		doVerb (bool, optional): Whether to print verbose output. (default: True)
+		doEq (bool, optional): Whether to plot equatorial or meridional (default: True)

 	Returns:
 		dataAbs (ndarray): The absolute error data.
@@ -118,8 +119,8 @@ def PlotEqErrAbs(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=T
 		Ax.clear()
 	dataAbs = None
 	for fn in fieldNames:
-		dataP = gsphP.EggSlice(fn, nStp, doEq=True, doVerb=doVerb)
-		dataO = gsphO.EggSlice(fn, nStp, doEq=True, doVerb=doVerb)
+		dataP = gsphP.EggSlice(fn, nStp, doEq=doEq, doVerb=doVerb)
+		dataO = gsphO.EggSlice(fn, nStp, doEq=doEq, doVerb=doVerb)
 		if dataAbs is None:
 			dataAbs = np.square(dataO - dataP)
 		else:
@@ -132,13 +133,16 @@ def PlotEqErrAbs(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=T
 	if doDeco:
 		kv.addEarth2D(ax=Ax)
 		Ax.set_xlabel('SM-X [Re]')
+		if doEq:
 			Ax.set_ylabel('SM-Y [Re]')
+		else:
+			Ax.set_ylabel('SM-Z [Re]')
 	return dataAbs

-#Plot relative error in the requested, or given, equatorial field
-def PlotEqErrRel(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=True, doDeco=True, vMin=1e-16, vMax=1, doLog=True, doVerb=True):
+#Plot relative error in the requested, or given, equatorial or meridional field
+def PlotErrRel(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=True, doDeco=True, vMin=1e-16, vMax=1, doLog=True, doVerb=True, doEq=True):
 	"""
-	PlotEqErrRel function plots the relative error between two gsph objects.
+	PlotErrRel function plots the relative error between two gsph objects.

 	Args:
 		gsphP (gsph): The gsph object representing the predicted values.
@@ -154,6 +158,7 @@ def PlotEqErrRel(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=T
 		vMax (float, optional): The maximum value for the colorbar. (default: 1)
 		doLog (bool, optional): Whether to use logarithmic scale for the colorbar. (default: True)
 		doVerb (bool, optional): Whether to print verbose output. (default: True)
+		doEq (bool, optional): Whether to plot equatorial or meridional (default: True)

 	Returns:
 		dataRel (ndarray): The relative error data.
@@ -180,8 +185,8 @@ def PlotEqErrRel(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=T
 	dataAbs = None
 	dataBase = None
 	for fn in fieldNames:
-		dataP = gsphP.EggSlice(fn, nStp, doEq=True, doVerb=doVerb)
-		dataO = gsphO.EggSlice(fn, nStp, doEq=True, doVerb=doVerb)
+		dataP = gsphP.EggSlice(fn, nStp, doEq=doEq, doVerb=doVerb)
+		dataO = gsphO.EggSlice(fn, nStp, doEq=doEq, doVerb=doVerb)
 		if dataAbs is None:
 			dataBase = np.square(dataP)
 			dataAbs = np.square(dataO - dataP)
@@ -200,7 +205,10 @@ def PlotEqErrRel(gsphP, gsphO, nStp, xyBds, Ax, fieldNames, AxCB=None, doClear=T
 	if doDeco:
 		kv.addEarth2D(ax=Ax)
 		Ax.set_xlabel('SM-X [Re]')
+		if doEq:
 			Ax.set_ylabel('SM-Y [Re]')
+		else:
+			Ax.set_ylabel('SM-Z [Re]')
 	return dataRel

 #Plot absolute error along the requested logical axis
--- a/kaipy/kaiViz.py
+++ b/kaipy/kaiViz.py
@@ -15,6 +15,7 @@ from matplotlib.colors import Normalize
 from matplotlib.colors import SymLogNorm
 from matplotlib.patches import Wedge
 from matplotlib import ticker
+import shutil

 # Kaipy modules
 from kaipy.kdefs import *
@@ -288,6 +289,13 @@ def savePic(fOut, dpiQ=300, doTrim=True, bLenX=20, bLenY=None, doClose=False, do
        else:
            plt.close(saveFigure)

+#Checks to see if image magick is in the current path
+def isMagick():
+    path = shutil.which('magick')
+    if (path is None):
+        return False
+    else:
+        return True

 #Use imagemagick to trim whitespace off figure
 #doEven: Guarantee even number of pixels in X/Y
@@ -309,7 +317,9 @@ def trimFig(fName, bLenX=20, bLenY=None, doEven=True):
    if bLenY is None:
        bLenY = bLenX

-    ComS = 'convert -trim -border %dx%d -bordercolor "#FFFFFF" ' % (bLenX, bLenY) + fName + ' ' + fName
+    if (not isMagick()):
+        return
+    ComS = 'magick ' + fName + ' -trim -border %dx%d -bordercolor "#FFFFFF" '%(bLenX,bLenY) + ' ' + fName
    os.system(ComS)

    if doEven:
@@ -355,7 +365,8 @@ def ShaveX(fName):
    Returns:
        fName file is cropped by one pixel on the right side.
    """
-    ComS = 'convert -crop -1+0 +repage ' + fName + ' ' + fName
+
+    ComS = 'magick %s -crop -1+0 +repage %s'%(fName,fName)
    os.system(ComS)

 def ShaveY(fName):
@@ -368,7 +379,8 @@ def ShaveY(fName):
    Returns:
        fName file is cropped by one pixel on the top.
    """
-    ComS = 'convert -crop +0-1 +repage ' + fName + ' ' + fName
+
+    ComS = 'magick %s -crop +0-1 +repage %s'%(fName,fName)
    os.system(ComS)

 #---------------------------------
--- a/kaipy/kaixml.py
+++ b/kaipy/kaixml.py
@@ -0,0 +1,31 @@
+# Standard modules
+import os
+import json
+import datetime
+
+# Third-party modules
+import xml.etree.ElementTree as ET
+
+def getXmlElement(E, name, caseInsensitive=True):
+    if E is None:
+        return None
+    for element in E.iter():
+        if caseInsensitive:
+            if element.tag.casefold() == name.casefold():
+                return element
+        else:
+            if element.tag == name:
+                return element
+    return None
+
+def getXmlAttribute(E, name, caseInsensitive=True):
+    if E is None:
+        return None
+    for attr_name,attr_value in E.attrib.items():
+        if caseInsensitive:
+            if attr_name.casefold() == name.casefold():
+                return attr_value
+        else:
+            if attr_name == name:
+                return attr_value
+    return None
--- a/kaipy/kdefs.py
+++ b/kaipy/kdefs.py
@@ -111,7 +111,8 @@ vc_cgs = 2.99792458e10 # [cm/s] speed of light
 barLen = 30
 barLab = 30 
 #barDef = 'fish'
-barDef = alive_progress.animations.bars.bar_factory(tip="><('>", chars='∙',background='')
+barDef  = alive_progress.animations.bars.bar_factory(tip="><('>", chars="∙",background="")
+barDef2 = alive_progress.animations.bars.bar_factory('∙○⦿●',tip=" ><>")

 #------
 # I/O
--- a/kaipy/remix/remix.py
+++ b/kaipy/remix/remix.py
@@ -7,6 +7,8 @@ import h5py
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
 import numpy as np
+import kaipy.kaiViz as kv
+import cmasher as cmr

 # Kaipy modules
 from kaipy.kdefs import RionE, REarth
@@ -19,6 +21,7 @@ mu0o4pi = 1.e-7          # mu0=4pi*10^-7 => mu0/4pi=10^-7
 facMax = 1.5
 facCM = cm.RdBu_r
 flxCM = cm.inferno
+eflxCM = plt.get_cmap('cmr.seaweed')

 class remix:
 	"""
@@ -51,7 +54,7 @@ class remix:
 		calcFaceAreas(self, x, y)
 			Calculates the area of each face in a quad mesh.
 	
-		plot(self, varname, ncontours=16, addlabels={}, gs=None, doInset=False, doCB=True, doCBVert=True, doGTYPE=False, doPP=False)
+		plot(self, varname, ncontours=16, addlabels={}, gs=None, doInset=False, doCB=True, doCBVert=True, doGTYPE=False, nPP=0, doMNDF=False)
 			Plots the specified variable.
 	
 		efield(self, returnDeltas=False, ri=Ri*1e3)
@@ -115,7 +118,10 @@ class remix:
 			'Deflux': {'min': 0, 'max': 10.},
 			'Penergy': {'min': 0, 'max': 240},
 			'Pflux': {'min': 0, 'max': 1.e7},
-			'Peflux': {'min': 0, 'max': 1.}
+			'Peflux': {'min': 0, 'max': 1.},
+			'beta': {'min':0.0, 'max':0.5},
+			'deltae': {'min':0, 'max':20},
+			'amtype': {'min':0, 'max':3}
 		}

 	def get_data(self, h5file, step):
@@ -171,6 +177,9 @@ class remix:
 			self.variables['sigmah']['data'] = self.ion['Hall conductance ' + h]
 			self.variables['energy']['data'] = self.ion['Average energy ' + h]
 			self.variables['flux']['data'] = self.ion['Number flux ' + h]
+			self.variables['beta']['data'] = self.ion['RCM beta '+h]
+			self.variables['deltae']['data'] = self.ion['Mono potential drop '+h]
+			self.variables['amtype']['data'] = self.ion['Auroral model type '+h]
 			if 'RCM grid type ' + h in self.ion.keys():
 				self.variables['gtype']['data'] = self.ion['RCM grid type ' + h]
 			if 'RCM plasmasphere density ' + h in self.ion.keys():
@@ -201,6 +210,9 @@ class remix:
 			self.variables['sigmah']['data'] = self.ion['Hall conductance ' + h][:, ::-1]
 			self.variables['energy']['data'] = self.ion['Average energy ' + h][:, ::-1]
 			self.variables['flux']['data'] = self.ion['Number flux ' + h][:, ::-1]
+			self.variables['beta']['data'] = self.ion['RCM beta '+h][:, ::-1]
+			self.variables['deltae']['data'] = self.ion['Mono potential drop '+h][:, ::-1]
+			self.variables['amtype']['data'] = self.ion['Auroral model type '+h][:, ::-1]
 			if 'RCM grid type ' + h in self.ion.keys():
 				self.variables['gtype']['data'] = self.ion['RCM grid type ' + h][:, ::-1]
 			if 'RCM plasmasphere density ' + h in self.ion.keys():
@@ -305,7 +317,7 @@ class remix:


 	# TODO: define and consolidate allowed variable names
-	def plot(self, varname, ncontours=16, addlabels={}, gs=None, doInset=False, doCB=True, doCBVert=True, doGTYPE=False, doPP=False):
+	def plot(self, varname, ncontours=16, addlabels={}, gs=None, doInset=False, doCB=True, doCBVert=True, doGTYPE=False, nPP=0, doMNDF=False):
 		"""
 		Plot the specified variable on a polar grid.

@@ -318,7 +330,8 @@ class remix:
 			doCB (bool): Whether to show the colorbar (default: True).
 			doCBVert (bool): Whether to show the colorbar vertically (default: True).
 			doGTYPE (bool): Whether to overplot grid type contours (default: False).
-			doPP (bool): Whether to overplot polar cap boundary contours (default: False).
+			nPP (real): plasmaspheric number density contour (default: 0 and not shown).
+			doMNDF (bool): Whether to plot mono/diffuse using gree/blue colorbar (default: False).

 		Returns:
 			ax (AxesSubplot): The AxesSubplot object containing the plot.
@@ -412,7 +425,13 @@ class remix:
 				LW = 0.75
 				alpha = 1
 				tOff = np.pi / 2.
-			ax.contour(tc + tOff, rc, tmp, levels=con_level, colors=con_color, linewidths=LW, alpha=alpha)
+			contour = ax.contour(tc + tOff, rc, tmp, levels=con_level, colors=con_color, linewidths=LW, alpha=alpha)
+			if con_name=='npsp':
+				cl_con = 'red'
+				ax.clabel(contour, inline=False, fontsize=8, fmt=lambda val: f"Npsph={val:.1f}", colors=cl_con, use_clabeltext=True)
+			if con_name=='gtype':
+				cl_con = 'purple'
+				ax.clabel(contour, inline=False, fontsize=8, fmt=lambda val: f"GIMAG={val:.2f}", colors=cl_con, use_clabeltext=True)

 		if not self.Initialized:
 			sys.exit("Variables should be initialized for the specific hemisphere (call init_var) prior to plotting.")
@@ -431,6 +450,10 @@ class remix:
 					'energy'    : r'Energy [keV]',
 					'flux'      : r'Flux [1/cm$^2$s]',
 					'eflux'     : r'Energy flux [erg/cm$^2$s]',
+					'gtype'     : r'IMAG grid type',
+					'beta'      : r'$\beta$',
+					'deltae'    : r'$\Delta$E [kV]',
+					'amtype'    : r'Auroral model type',
 					'ephi'      : r'$E_\phi$ [mV/m]',
 					'etheta'    : r'$E_\theta$ [mV/m]',
 					'efield'    : r'|E| [mV/m]',
@@ -461,7 +484,7 @@ class remix:
 			upper = self.variables[varname]['data'].max()

 		# define number format string for min/max labels
-		if varname !='flux' and varname !='Pflux': 
+		if varname not in ['flux','Dflux','Mflux','Pflux']:
 			if varname == 'jped':
 				format_str = '%.2f'
 			else:
@@ -472,7 +495,10 @@ class remix:
 		if (varname == 'potential') or (varname == 'current'):
 			cmap=facCM
 		elif varname in ['flux','energy','eflux','joule','Mflux','Menergy','Meflux','Dflux','Denergy','Deflux','Pflux','Penergy','Peflux']:
+			if not doMNDF:
 				cmap=flxCM
+			else:
+				cmap=eflxCM
 			latlblclr = 'white'
 		elif (varname == 'velocity'): 
 			cmap=cm.YlOrRd
@@ -480,6 +506,12 @@ class remix:
 			cmap=None #cm.RdBu_r#None # cm.hsv
 		elif (varname == 'magpert'): 
 			cmap = cm.YlOrRd
+		elif (varname == 'beta'): 
+			cmap=cm.inferno
+		elif (varname == 'gtype'): 
+			cmap=cm.magma_r
+		elif (varname == 'deltae'): 
+			cmap=cm.viridis
 		else:
 			cmap=None # default is used
 		
@@ -518,7 +550,7 @@ class remix:
 		if varname == 'joule':
 			self.variables[varname]['data'] = self.joule()*1.e3  # convert to mW/m^2

-		variable = self.variables[varname]['data']
+#		variable = self.variables[varname]['data']

 		fig = plt.gcf()
 		
@@ -528,7 +560,39 @@ class remix:
 		else:
 			ax=fig.add_subplot(polar=True) 

+#		p=ax.pcolormesh(theta+tOff,r,variable,cmap=cmap,vmin=lower,vmax=upper)
+
+		if not doMNDF:
+			# traditional colorbar limits and maps
+			variable = self.variables[varname]['data']
 			p=ax.pcolormesh(theta+tOff,r,variable,cmap=cmap,vmin=lower,vmax=upper)
+		else:
+			# for mono/diffuse, use different colorbar limits and maps
+			efxAmax = 10 # limit in the Mono-diffuse asymmetric colorbar.
+			numAmax = 1e9
+			engAmax = 20
+			idiff = self.variables['amtype']['data']<=2
+			if varname in ['eflux','Meflux','Deflux','thmeflux']:
+				variable = self.variables[varname]['data']
+				variable[idiff] = -abs(variable[idiff])
+				# use linear scale if below 0.1 mW/m^2
+				# set color bar limit at efxAmax mW/m^2
+				vQ = kv.genNorm(efxAmax,doSymLog=False,linP=0.01*efxAmax)
+				p=ax.pcolormesh(theta+tOff,r,variable,cmap=cmap,norm=vQ)
+			elif varname in ['flux','Mflux','Dflux','thmflux']:
+				variable = self.variables[varname]['data']
+				variable[idiff] = -abs(variable[idiff])
+				# use linear scale if below 1e7 #/cm^2/s
+				# set color bar limit at numAmax #/cm^2/s
+				vQ = kv.genNorm(numAmax,doSymLog=False,linP=0.01*numAmax)
+				p=ax.pcolormesh(theta+tOff,r,variable,cmap=cmap,norm=vQ)
+			elif varname in ['energy','Menergy','Denergy']:
+				variable = self.variables[varname]['data']
+				variable[idiff] = -abs(variable[idiff])
+				# use linear scale if below 0.1 keV
+				# set color bar limit at engAmax keV
+				vQ = kv.genNorm(engAmax,doSymLog=False,linP=0.01*engAmax)
+				p=ax.pcolormesh(theta+tOff,r,variable,cmap=cmap,norm=vQ)

 		if (not doInset):
 			if (doCB):
@@ -538,8 +602,21 @@ class remix:
 					cb=plt.colorbar(p,ax=ax,pad=0.1,shrink=0.85,orientation="horizontal")
 				cb.set_label(cblabels[varname])

+				if not doMNDF:
 					ax.text(-75.*np.pi/180.,1.2*r.max(),('min: '+format_str+'\nmax: ' +format_str) % 
 					  (variable.min() ,variable.max()))
+				else:
+					if varname in ['flux','Mflux','Dflux']:
+						cticks = cb.get_ticks()
+						cticklab = [str('{:6.1e}'.format(abs(tick))) for tick in cticks]
+						cb.set_ticklabels(cticklab)
+					if varname in ['eflux','Meflux','Deflux','energy','Menergy','Denergy']:
+						cticks = cb.get_ticks()
+						cticklab = [str('{:.1f}'.format(abs(tick))) for tick in cticks]
+						cb.set_ticklabels(cticklab)
+				if varname in ['gtype']:
+					cb.set_ticks([0,1])
+					cb.set_ticklabels(['MHD','IMAG'])
 			
 		lines, labels = plt.thetagrids((0.,90.,180.,270.),hour_labels)
 		lines, labels = plt.rgrids(circles,lbls,fontsize=8,color=latlblclr)
@@ -550,8 +627,18 @@ class remix:
 		# ax.axis([0,2*np.pi,0,r.max()],'tight')
 		ax.axis([0,2*np.pi,0,r.max()])
 		if (not doInset):
-			ax.text(-75.*np.pi/180.,1.2*r.max(),('min: '+format_str+'\nmax: ' +format_str) % 
-				  (variable.min() ,variable.max()))
+			if doMNDF and varname in ['eflux','Meflux','Deflux','flux','Mflux','Dflux','energy','Menergy','Denergy']:
+				mono=variable.copy()
+				diff=variable.copy()
+				mono[ idiff]=0.0
+				diff[~idiff]=0.0
+#				print('mono min/max:',mono.min(),mono.max())
+#				print('diff min/max:',diff.min(),diff.max())
+				ax.text(-77.*np.pi/180.,1.1*r.max(),('Diff max: '+format_str) % (-diff.min()))
+				ax.text( 75.*np.pi/180.,1.1*r.max(),('Mono max: '+format_str) % ( mono.max()))
+			else:
+				ax.text(-75.*np.pi/180.,1.2*r.max(),('min: '+format_str+'\nmax: ' +format_str) % (variable.min() ,variable.max()))
+
 			if varname in ['eflux','Meflux','Deflux','Peflux']:
 				ax.text(75.*np.pi/180.,1.3*r.max(),('GW: '+format_str) % power)
 			if (varname == 'current'):
@@ -561,9 +648,11 @@ class remix:
 		if varname=='current': 
 			potential_overplot(doInset)
 		if doGTYPE and varname=='eflux': 
-			boundary_overplot('gtype',[0.01,0.99],'green',doInset)
-		if doPP and varname=='eflux':
-			boundary_overplot('npsp',[10],'cyan',doInset)
+			cl_con = 'purple'
+			boundary_overplot('gtype',[0.5],cl_con,doInset)
+		if nPP>0.0 and varname=='eflux':
+			cl_con = 'red'
+			boundary_overplot('npsp',[nPP],cl_con,doInset)

 		return ax
 	# mpl.rcParams['contour.negative_linestyle'] = 'solid'
--- a/kaipy/scripts/postproc/slimh5.py
+++ b/kaipy/scripts/postproc/slimh5.py
@@ -56,9 +56,10 @@ def createfile(iH5,fOut):
        sQ = str(Q)
        #Don't include stuff that starts with "Step"
        if "Step" not in sQ and cacheName not in sQ:
+            if isinstance(iH5[sQ], h5py.Group):
+                oH5.copy(iH5[sQ], oH5)
+            else:
                oH5.create_dataset(sQ,data=iH5[sQ])
-        if cacheName in sQ:
-            oH5.create_group(sQ)
    return oH5


@@ -172,6 +173,7 @@ def main():
        #Add the cache after steps, select the same steps for the cache that are contained in the
        #Ns:Ne:Nsk start,end,stride
        if cacheName in iH5.keys():
+            oH5.create_group(cacheName)
            for Q in iH5[cacheName].keys():
                sQ = str(Q)
               
--- a/kaipy/scripts/quicklook/gamerrVid.py
+++ b/kaipy/scripts/quicklook/gamerrVid.py
@@ -50,7 +50,7 @@ def makeMovie(frame_dir,movie_name):
 	subprocess.run(cmd, check=True)

 # python allows changes by reference to the errTimes,errListRel, errListAbs lists
-def makeImage(i,gsph1,gsph2,tOut,doVerb,xyBds,fnList,oDir,errTimes,errListRel,errListAbs,cv,dataCounter, vO, figSz, noMPI, noLog, fieldNames):
+def makeImage(i,gsph1,gsph2,tOut,doVerb,doEq,logAx,xyBds,fnList,oDir,errTimes,errListRel,errListAbs,cv,dataCounter, vO, figSz, noMPI, noLog, fieldNames):
 	if doVerb:
 		print("Making image %d"%(i))
 	#Convert time (in seconds) to Step #
@@ -77,23 +77,45 @@ def makeImage(i,gsph1,gsph2,tOut,doVerb,xyBds,fnList,oDir,errTimes,errListRel,er
 	AxB2.clear()
 	
 	#plot upper left msph error
-	mviz.PlotEqErrRel(gsph1,gsph2,nStp,xyBds,AxTL,fnList,AxCB=AxCT,doVerb=doVerb)
+	mviz.PlotErrRel(gsph1,gsph2,nStp,xyBds,AxTL,fnList,AxCB=AxCT,doVerb=doVerb,doEq=doEq)
+	if doEq:
 		AxTL.set_title("Equatorial Slice of Relative Error")
+	else:
+		AxTL.set_title("Meridional Slice of Relative Error")
 	
-	#plot upper right k-axis error
-	mviz.PlotLogicalErrRel(gsph1,gsph2,nStp,AxTR,fnList,2,doVerb=doVerb)
+	#plot upper right cumulative logical error
+	mviz.PlotLogicalErrRel(gsph1,gsph2,nStp,AxTR,fnList,logAx,doVerb=doVerb)
+	if logAx == 0:
+		AxTR.set_title("Per-Cell Relative Error along I-Axis")
+	elif logAx == 1:
+		AxTR.set_title("Per-Cell Relative Error along J-Axis")
+	else:
 		AxTR.set_title("Per-Cell Relative Error along K-Axis")
 	if (not noMPI):
 		#plot I-MPI decomp on logical plot
-		if(gsph2.Ri > 1):
+		if(gsph2.Ri > 1 and logAx != 0):
 			for im in range(gsph2.Ri):
 				i0 = im*gsph2.dNi
+				if logAx == 1:
+					AxTR.plot([i0, i0],[0, gsph2.Nk],"deepskyblue",linewidth=0.25,alpha=0.5)
+				else:
 					AxTR.plot([i0, i0],[0, gsph2.Nj],"deepskyblue",linewidth=0.25,alpha=0.5)
 		#plot J-MPI decomp on logical plot
-		if (gsph2.Rj>1):
+		if (gsph2.Rj>1 and logAx != 1):
 			for jm in range(1,gsph2.Rj):
 				j0 = jm*gsph2.dNj
+				if logAx == 0:
+					AxTR.plot([j0, j0],[0, gsph2.Nk],"deepskyblue",linewidth=0.25,alpha=0.5)
+				else:
 					AxTR.plot([0, gsph2.Ni],[j0, j0],"deepskyblue",linewidth=0.25,alpha=0.5)
+		#plot K-MPI decomp on logical plot
+		if (gsph2.Rk>1 and logAx != 2):
+			for km in range(1,gsph2.Rk):
+				k0 = km*gsph2.dNk
+				if logAx == 0:
+					AxTR.plot([0, gsph2.Nj],[k0, k0],"deepskyblue",linewidth=0.25,alpha=0.5)
+				else:
+					AxTR.plot([0, gsph2.Ni],[k0, k0],"deepskyblue",linewidth=0.25,alpha=0.5)
 	
 	#plot bottom line plot
 	etval = tOut[i]/60.0
@@ -168,15 +190,17 @@ def create_command_line_parser():
 	fdir2 = os.getcwd()
 	ftag2 = "msphere"
 	oDir = "vid2D"
-	ts = 0    #[min]
-	te = 200  #[min]
+	ts = 0.0	  #[min]
+	te = 200.0	#[min]
 	dt = 0.0 #[sec] 0 default means every timestep
+	logAx = 2 # axis to accumulate logical error along
 	Nth = 1 #Number of threads
 	noMPI = False # Don't add MPI tiling
 	noLog = False
 	fieldNames = "Bx, By, Bz"
 	doVerb = False
 	skipMovie = False
+	merid = False

 	MainS = """Creates simple multi-panel figure for Gamera magnetosphere run
 	Left Panel - Residual vertical magnetic field
@@ -189,14 +213,16 @@ def create_command_line_parser():
 	parser.add_argument('-d2',type=str,metavar="directory",default=fdir2,help="Directory to read second dataset from (default: %(default)s)")
 	parser.add_argument('-id2',type=str,metavar="runid",default=ftag2,help="RunID of second dataset (default: %(default)s)")
 	parser.add_argument('-o',type=str,metavar="directory",default=oDir,help="Subdirectory to write to (default: %(default)s)")
-	parser.add_argument('-ts' ,type=int,metavar="tStart",default=ts,help="Starting time [min] (default: %(default)s)")
-	parser.add_argument('-te' ,type=int,metavar="tEnd"  ,default=te,help="Ending time   [min] (default: %(default)s)")
-	parser.add_argument('-dt' ,type=int,metavar="dt"    ,default=dt,help="Cadence       [sec] (default: %(default)s)")
+	parser.add_argument('-ts' ,type=float,metavar="tStart",default=ts,help="Starting time [min] (default: %(default)s)")
+	parser.add_argument('-te' ,type=float,metavar="tEnd"	,default=te,help="Ending time	[min] (default: %(default)s)")
+	parser.add_argument('-dt' ,type=float,metavar="dt"	,default=dt,help="Cadence		[sec] (default: %(default)s)")
+	parser.add_argument('-logAx',type=int,metavar="logAx",default=logAx,help="Index of the axis to accumulate along in the upper-right plot (default: %(default)s)")
 	parser.add_argument('-Nth' ,type=int,metavar="Nth",default=Nth,help="Number of threads to use (default: %(default)s)")
 	parser.add_argument('-f',type=str,metavar="fieldnames",default=fieldNames,help="Comma-separated fields to plot (default: %(default)s)")
 	parser.add_argument('-linear',action='store_true', default=noLog,help="Plot linear line plot instead of logarithmic (default: %(default)s)")
 	parser.add_argument('-v',action='store_true', default=doVerb,help="Do verbose output (default: %(default)s)")
 	parser.add_argument('-skipMovie',action='store_true', default=skipMovie,help="Skip automatic movie generation afterwards (default: %(default)s)")
+	parser.add_argument('-merid',action='store_true', default=merid,help="Plot meridional instead of equatorial slice (default: %(default)s)")
 	#parser.add_argument('-nompi', action='store_true', default=noMPI,help="Don't show MPI boundaries (default: %(default)s)")


@@ -209,15 +235,17 @@ def main():
 	fdir2 = os.getcwd()
 	ftag2 = "msphere"
 	oDir = "vid2D"
-	ts = 0    #[min]
-	te = 200  #[min]
+	ts = 0.0	  #[min]
+	te = 200.0	#[min]
 	dt = 0.0 #[sec] 0 default means every timestep
+	logAx = 2
 	Nth = 1 #Number of threads
 	noMPI = False # Don't add MPI tiling
 	noLog = False
 	fieldNames = "Bx, By, Bz"
 	doVerb = False
 	skipMovie = False
+	doEq = True

 	parser = create_command_line_parser()
 	mviz.AddSizeArgs(parser)
@@ -231,11 +259,13 @@ def main():
 	ts	= args.ts
 	te	= args.te
 	dt	= args.dt
+	logAx = args.logAx
 	oSub = args.o
 	Nth = args.Nth
 	fieldNames = args.f
 	noLog = args.linear
 	doVerb = args.v
+	doEq = not args.merid
 	#noMPI = args.noMPI
 	
 	fnList = [item.strip() for item in fieldNames.split(',')]
@@ -275,8 +305,8 @@ def main():
 	Nt = len(tOut)
 	vO = np.arange(0,Nt)

-	print("Writing %d outputs between minutes %d and %d"%(Nt,ts,te))
-	print("Using %d threads"%(Nth))
+	print(f"Writing {Nt} outputs between minutes {ts} and {te}")
+	print(f"Using {Nth} threads")
 	
 	errTimes = []
 	errListRel = []
@@ -284,7 +314,7 @@ def main():

 	#Loop over sub-range
 	titstr = "Comparing '%s' to '%s'"%(fdir1,fdir2)
-	with alive_bar(Nt,title=titstr.ljust(kdefs.barLab),length=kdefs.barLen,disable=doVerb) as bar:
+	with alive_bar(Nt,title=titstr.ljust(kdefs.barLab),length=kdefs.barLen,bar=kdefs.barDef,disable=doVerb) as bar:
 		#with concurrent.futures.ThreadPoolExecutor(max_workers=Nth) as executor:
 		with concurrent.futures.ProcessPoolExecutor(max_workers=Nth) as executor:
 			m = multiprocessing.Manager()
@@ -293,8 +323,8 @@ def main():
 			met = m.list(errTimes)
 			melr = m.list(errListRel)
 			mela = m.list(errListAbs)
-			#imageFutures = {executor.submit(makeImage,i,gsph1,gsph2,tOut,doVerb,xyBds,fnList,oDir,errTimes,errListRel,errListAbs,cv): i for i in range(0,Nt)}
-			imageFutures = {executor.submit(makeImage,i,gsph1,gsph2,tOut,doVerb,xyBds,fnList,oDir,met,melr,mela,cv,dataCounter,vO,figSz,noMPI,noLog,fieldNames): i for i in range(0,Nt)}
+			#imageFutures = {executor.submit(makeImage,i,gsph1,gsph2,tOut,doVerb,doEq,logAx,xyBds,fnList,oDir,errTimes,errListRel,errListAbs,cv): i for i in range(0,Nt)}
+			imageFutures = {executor.submit(makeImage,i,gsph1,gsph2,tOut,doVerb,doEq,logAx,xyBds,fnList,oDir,met,melr,mela,cv,dataCounter,vO,figSz,noMPI,noLog,fieldNames): i for i in range(0,Nt)}
 			for future in concurrent.futures.as_completed(imageFutures):
 				try:
 					retVal = future.result()
--- a/kaipy/scripts/quicklook/gamerrpic.py
+++ b/kaipy/scripts/quicklook/gamerrpic.py
@@ -102,8 +102,8 @@ def main():
 	AxL.clear()
 	AxR.clear()

-	mviz.PlotEqErrRel(gsph1,gsph2,nStp,xyBds,AxL,fnList,AxCB=AxCL)
-	mviz.PlotEqErrAbs(gsph1,gsph2,nStp,xyBds,AxR,fnList,AxCB=AxCR)
+	mviz.PlotErrRel(gsph1,gsph2,nStp,xyBds,AxL,fnList,AxCB=AxCL)
+	mviz.PlotErrAbs(gsph1,gsph2,nStp,xyBds,AxR,fnList,AxCB=AxCR)

 	gsph1.AddTime(nStp,AxL,xy=[0.025,0.89],fs="x-large")

--- a/kaipy/scripts/quicklook/mageDash.py
+++ b/kaipy/scripts/quicklook/mageDash.py
--- a/kaipy/scripts/quicklook/mixpic.py
+++ b/kaipy/scripts/quicklook/mixpic.py
@@ -122,8 +122,12 @@ def create_command_line_parser():
        help="Show RCM grid type in the eflx plot (default: %(default)s)"
    )
    parser.add_argument(
-        '-PP', action='store_true', default=False,
-        help="Show plasmapause (10/cc) in the eflx/nflx plot (default: %(default)s)"
+        '--nPP', type=float, metavar="nPP", default=0,
+        help="Plasmasphere density contour to show (default: %(default)/cc not shown)"
+    )
+    parser.add_argument(
+        '-MNDF', action='store_true', default=False,
+        help="Show mono-diffuse eflx/nflx in linear/log green/blue colormap (default: %(default)s)"
    )
    parser.add_argument(
        '-vid', action='store_true', default=False,
@@ -156,7 +160,8 @@ def makePlot(i, remixFile, nStp, args, varDict):
    spacecraft = args.spacecraft
    verbose = args.verbose
    do_GTYPE = args.GTYPE
-    do_PP = args.PP
+    nPP = args.nPP
+    do_MNDF = args.MNDF
    do_vid = args.vid
    do_overwrite = args.overwrite
    do_hash = not args.nohash
@@ -233,12 +238,12 @@ def makePlot(i, remixFile, nStp, args, varDict):
        axs[0] = ion.plot('current', gs=gs[0, 0])
        axs[1] = ion.plot('sigmap', gs=gs[0, 1])
        axs[2] = ion.plot('sigmah', gs=gs[0, 2])
-        axs[3] = ion.plot('joule', gs=gs[1, 0])
-        axs[4] = ion.plot('energy', gs=gs[1, 1])
        if do_nflux:
-            axs[5] = ion.plot('flux', gs=gs[1, 2],doGTYPE=do_GTYPE,doPP=do_PP)
+            axs[3] = ion.plot('flux', gs=gs[1, 0],doGTYPE=do_GTYPE,nPP=nPP,doMNDF=do_MNDF)
        else:
-            axs[5] = ion.plot('eflux', gs=gs[1, 2],doGTYPE=do_GTYPE,doPP=do_PP)
+            axs[3] = ion.plot('joule', gs=gs[1, 0])
+        axs[4] = ion.plot('energy', gs=gs[1, 1],doGTYPE=do_GTYPE,nPP=nPP,doMNDF=do_MNDF)
+        axs[5] = ion.plot('eflux', gs=gs[1, 2],doGTYPE=do_GTYPE,nPP=nPP,doMNDF=do_MNDF)

        # If requested, plot the magnetic footprints for the specified
        # spacecraft.
@@ -314,7 +319,8 @@ def main():
    spacecraft = args.spacecraft
    verbose = args.verbose
    do_GTYPE = args.GTYPE
-    do_PP = args.PP
+    nPP = args.nPP
+    do_MNDF = args.MNDF
    do_vid = args.vid
    do_overwrite = args.overwrite
    do_hash = not args.nohash
--- a/kaipy/solarWind/OMNI.py
+++ b/kaipy/solarWind/OMNI.py
@@ -5,8 +5,13 @@ import re
 # Third-party modules
 import numpy
 from cdasws import CdasWs
-from pyspedas import kyoto
 from pytplot import get_data
+try:
+    # Try the newer import path first (for recent pyspedas versions)
+    from pyspedas.projects import kyoto
+except ImportError:
+    # Fall back to the older import path
+    from pyspedas import kyoto

 # Kaipy modules
 import kaipy.transform
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ include-package-data = true

 [project]
 name = "kaipy"
-version = "1.1.3"
+version = "1.1.4"
 authors = [
    {name = "Eric Winter", email = "eric.winter@jhuapl.edu"},
 ]
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,11 +4,16 @@ cartopy
 cdasws
 cmasher
 configparser
+dash
+dash_bootstrap_components
+flask_caching
 h5py
 jupyterlab
 matplotlib
 pandas
 progressbar
+pyqt5
+pyqtwebengine
 pyspedas
 pytest
 slack_sdk
--- a/setup.py
+++ b/setup.py
@@ -2,7 +2,7 @@ from setuptools import setup, find_packages

 setup(
    name='kaipy',
-    version='1.1.3',
+    version='1.1.4',
    description='Python software for CGS MAGE and other Kaiju models',
    author='Kaiju team',
    author_email='wiltbemj@ucar.edu',
@@ -25,11 +25,16 @@ setup(
        'cdasws',
        'cmasher',
        'configparser',
+        'dash',
+        'dash_bootstrap_components',
+        'flask_caching',
        'h5py',
        'jupyterlab',
        'matplotlib',
        'pandas',
        'progressbar',
+        'pyqt5',
+        'pyqtwebengine',
        'pyspedas',
        'pytest',
        'slack_sdk',
@@ -91,6 +96,7 @@ setup(
            'gamsphVid=kaipy.scripts.quicklook.gamsphVid:main',
            'heliomovie=kaipy.scripts.quicklook.heliomovie:main',
            'heliopic=kaipy.scripts.quicklook.heliopic:main',
+            'mageDash=kaipy.scripts.quicklook.mageDash:main',
            'mixpic=kaipy.scripts.quicklook.mixpic:main',
            'msphpic=kaipy.scripts.quicklook.msphpic:main',
            'raijupic=kaipy.scripts.quicklook.raijupic:main',
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -47,6 +47,12 @@ def write_mix_h5(file_path,mjds=None):
            grp.create_dataset("Average energy SOUTH", data=np.zeros((Nlat, Nlon)))
            grp.create_dataset("Number flux NORTH", data=np.zeros((Nlat, Nlon)))
            grp.create_dataset("Number flux SOUTH", data=np.zeros((Nlat, Nlon)))
+            grp.create_dataset("RCM beta NORTH", data=np.zeros((Nlat, Nlon)))
+            grp.create_dataset("RCM beta SOUTH", data=np.zeros((Nlat, Nlon)))
+            grp.create_dataset("Mono potential drop NORTH", data=np.zeros((Nlat, Nlon)))
+            grp.create_dataset("Mono potential drop SOUTH", data=np.zeros((Nlat, Nlon)))
+            grp.create_dataset("Auroral model type NORTH", data=np.zeros((Nlat, Nlon)))
+            grp.create_dataset("Auroral model type SOUTH", data=np.zeros((Nlat, Nlon)))

    return file_path

--- a/tests/gamera/test_mspViz.py
+++ b/tests/gamera/test_mspViz.py
@@ -7,7 +7,7 @@ import argparse
 from argparse import RawTextHelpFormatter

 import kaipy.remix.remix as remix
-from kaipy.gamera.msphViz import AddSizeArgs, GetSizeBds, PlotEqErrAbs, PlotEqErrRel, PlotLogicalErrAbs, PlotLogicalErrRel, CalcTotalErrAbs, CalcTotalErrRel, PlotEqB, PlotMerid, PlotJyXZ, PlotEqEphi, PlotMPI, AddIonBoxes, plotPlane, plotXY, plotXZ
+from kaipy.gamera.msphViz import AddSizeArgs, GetSizeBds, PlotErrAbs, PlotErrRel, PlotLogicalErrAbs, PlotLogicalErrRel, CalcTotalErrAbs, CalcTotalErrRel, PlotEqB, PlotMerid, PlotJyXZ, PlotEqEphi, PlotMPI, AddIonBoxes, plotPlane, plotXY, plotXZ


@pytest.fixture
@@ -26,22 +26,22 @@ def test_GetSizeBds():
    result = GetSizeBds(args)
    assert result == [-100.0, 20.0, -60.0, 60.0]

-def test_PlotEqErrAbs(gamera_pipe):
+def test_PlotErrAbs(gamera_pipe):
    fig, ax = plt.subplots()
    gsphP = gamera_pipe
    gsphO = gamera_pipe
    fieldNames = ['Bz']
    xyBds = [-100, 100, -100, 100]
-    result = PlotEqErrAbs(gsphP, gsphO, 0, xyBds, ax, fieldNames)
+    result = PlotErrAbs(gsphP, gsphO, 0, xyBds, ax, fieldNames)
    assert isinstance(result, np.ndarray)

-def test_PlotEqErrRel(gamera_pipe):
+def test_PlotErrRel(gamera_pipe):
    fig, ax = plt.subplots()
    gsphP = gamera_pipe
    gsphO = gamera_pipe
    fieldNames = ['Bz']
    xyBds = [-100, 100, -100, 100]
-    result = PlotEqErrRel(gsphP, gsphO, 0, xyBds, ax, fieldNames)
+    result = PlotErrRel(gsphP, gsphO, 0, xyBds, ax, fieldNames)
    assert isinstance(result, np.ndarray)

 def test_PlotLogicalErrAbs(gamera_pipe):