import h5py as h5 import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec import numpy as np import datetime import kaipy.kdefs as kd import kaipy.kaiTools as kt import kaipy.kaijson as kj import kaipy.kaiViz as kv import fileIO isotfmt = '%Y-%m-%dT%H:%M:%S' #--- # Controls #--- # File names fname_slice = "slice.h5" fname_moments = "imagMoments.h5" fname_scTrack = "RBSPB.sc.h5" fname_imagTrack = "rcmTrack_RBSP-B.json" fname_dsts = "Dsts.json" # Calculations r_DPS = 6 # Radial distance of volume to calculate DPS-Dst within # Plotting utPlotStr = "2013-03-17T10:09:00" # UT to display plot data for int_cmap = "CMRmap" # Colormap name for intensity plot int_norm = kv.genNorm(1E4, 5E6, doLog=True) # Intensity normalization in energy_bounds = [10, 200] # Energy bounds for intensity plot in keV eq_xbounds = [-15, 10] # Bounds in X direction for equatorial pressure plot #--- # Helper functions #--- def dtIdx(dtList,dt): # Find index of closest datetime value within a datetime array return np.array([np.abs((dtList[n]-dt).total_seconds()) for n in range(len(dtList))]).argmin() def calcSlcDPSDst(fname_slice, slcInfo, r_DPS, utStart, utEnd): # Open file slc5 = h5.File(fname_slice,'r') # Calculate area of each equatorial grid cell dA_arr = np.zeros(slc5['Bz0'].shape) Bz0 = slc5['Bz0'][:] dpsdst = np.zeros(slcInfo.Nt) dDps_dt = np.zeros(slcInfo.Nt) X = slc5['X'][:] Y = slc5['Y'][:] r_cc = np.sqrt(kt.to_center2D(X)**2 + kt.to_center2D(Y)**2) for i in range(X.shape[0]-1): for j in range(X.shape[1]-1): if r_cc[i,j] < r_DPS: # shoelace theorem dA_arr[i,j] = 0.5*np.abs(X[i ,j ]*Y[i+1,j ] - Y[i ,j ]*X[i+1,j ] + X[i+1,j ]*Y[i+1,j+1] - Y[i+1,j ]*X[i+1,j+1] + X[i+1,j+1]*Y[i ,j+1] - Y[i+1,j+1]*X[i ,j+1] + X[i ,j+1]*Y[i ,j ] - Y[i ,j+1]*X[i ,j ]) # Now calc dpsdst for each time step dpsdst = np.zeros(slcInfo.Nt) for s in range(slcInfo.Nt): s5 = slc5[slcInfo.stepStrs[s]] bVol = s5['dvB'][:]*kd.REarth/1e-9 * 4.581 # [Re/CODE_B] -> [m/T] # bVol will be zero on cells that are on non-closed flux tubes KE = s5['bP'][:]*1e-9*bVol*(s5['dBz'][:]+Bz0)*1e-9*dA_arr*kd.REarth**2 / kd.kev2J # p[J/m^3] * bVol[m/T] * B[T] * Re^2[m^2] = [J] * keV/J = [keV] # Clean up NaNs KE = np.nan_to_num(KE) # [keV] dpsdst[s] += -4.2*(1.0e-30)*np.sum(KE) # Calculate time derivative # End points, just do side-by-side derivative dt = slcInfo.times[1:]-slcInfo.times[:-1] dDps_dt[0] = (dpsdst[1]-dpsdst[0])/dt[0] dDps_dt[-1] = (dpsdst[-1]-dpsdst[-2])/dt[-1] #The rest, dDps_dt[i] = (dps[i+1]-dps[i-1])/(2h) dDps_dt[1:-1] = (dpsdst[2:] - dpsdst[:-2])/(dt[:-1]+dt[1:]) # Save to file if we want one dpsDstDict = {'UT' : slcInfo.UTs, 'dpsdst' : dpsdst, 'dDps_dt' : dDps_dt} # Chop down to the provided ut range sStart = dtIdx(slcInfo.UTs, utStart) sEnd = dtIdx(slcInfo.UTs, utEnd) for k in dpsDstDict.keys(): dpsDstDict[k] = dpsDstDict[k][sStart:sEnd] return dpsDstDict def getMomentData(fname_moments, momentInfo, utPlot): # Get index closest to desired plot time utIdx = dtIdx(momentInfo.UTs, utPlot) f5 = h5.File(fname_moments,'r') step5 = f5[momentInfo.stepStrs[utIdx]] IOpen = step5['IOpen'][:] # Information about whether a cell is on an open/closed flux tube P = step5['P'][:] # RCM pressure Pmask = np.ma.masked_where(IOpen>-1, P) return { 'xMin' : step5['xMin'], 'yMin' : step5['yMin'], 'P' : Pmask } def plotDsts(Ax, bsDstDict, dpsDstDict): Ax.plot(bsDstDict['UT'] , bsDstDict['symh'] , label="Observed Sym-H") Ax.plot(bsDstDict['UT'] , bsDstDict['bsdst'] , label="Model Biot-Savart Dst") Ax.plot(dpsDstDict['UT'], dpsDstDict['dpsdst'], label="Model DPS Dst within {} R$_E$".format(r_DPS)) Ax.legend() Ax.grid(True,alpha=0.6) Ax.set_xlim([bsDstDict['UT'][0], bsDstDict['UT'][-1]]) Ax.set_xlabel("UT") Ax.set_ylabel("Sym-H [nT]") def plotPressEq(Ax, momentData, utPlot): Ax.pcolormesh(momentData['xMin'], momentData['yMin'], momentData['P']) Ax.set_xlabel("X [R$_E$]") Ax.set_ylabel("Y [R$_E$]") Ax.yaxis.set_ticks_position('right') Ax.yaxis.set_label_position('right') kv.addEarth2D(ax=AxPressEq) # Add UT label utTextStr = datetime.datetime.strftime(utPlot, "2013-03-17 10:09:00") tLbl = Ax.text(0.98,0.97,utTextStr, horizontalalignment='right', \ verticalalignment='top',transform=Ax.transAxes,fontsize=14,zorder=100) tLbl.set_bbox(dict(facecolor='white',edgecolor='black',alpha=1)) # Make sure plot aspect ratio is 1 boxRatio = kv.get_aspect(AxPressEq) dx = eq_xbounds[1] - eq_xbounds[0] dy = boxRatio*dx ybMin = -dy/2 ybMax = dy/2 AxPressEq.set_xlim(eq_xbounds) AxPressEq.set_ylim([ybMin, ybMax]) def plotSatTrack(Ax, imagTrack, utPlot): # Plot whole satellite track, and current satellite position at time utPlot # Get index closest to desired plot time utIdx = dtIdx(imagTrack['UT'], utPlot) # Get sposition info x_sc = imagTrack['xeq'] y_sc = imagTrack['yeq'] req_sc = imagTrack['Req'] xm = np.ma.masked_where(req_sc < 0.2, x_sc) ym = np.ma.masked_where(req_sc < 0.2, y_sc) # Whole trajectory Ax.plot(xm, ym, '-',color='#f7a116', alpha=0.5,zorder=10, linewidth=3) # Current location if req_sc[utIdx] > 1E-8: satCircleB = plt.Circle((x_sc[utIdx], y_sc[utIdx]), 0.18, color='#f7e4a3',zorder=12) satCircleW = plt.Circle((x_sc[utIdx], y_sc[utIdx]), 0.26, color='#f7a116',zorder=11) Ax.add_patch(satCircleW) Ax.add_patch(satCircleB) def plotIntensity(Ax, imagTrack): energyGrid = imagTrack['ions']['energies'] UT = np.array(imagTrack['UT']) utGrid = np.repeat(UT[:,np.newaxis], energyGrid.shape[1], axis=1) Ax.pcolormesh(utGrid, energyGrid, imagTrack['ions']['diffFlux'], norm=int_norm, cmap=int_cmap, shading='auto') Ax.set_yscale('log') Ax.set_ylim(energy_bounds) Ax.set_ylabel('Energy [keV]') Ax.set_xlabel('UT') def plotTimeLine(Ax, utPlot): ylim = Ax.get_ylim() Ax.plot([utPlot, utPlot], ylim, '-k',linewidth=2) Ax.plot([utPlot, utPlot], ylim, '-w',linewidth=1) Ax.set_ylim(ylim) #--- # Main #--- if __name__=="__main__": utPlot = datetime.datetime.strptime(utPlotStr, isotfmt) #--- # Retrieve Data #--- # Read imagTrack data from json print("Reading imagTrack info") imagTrack = kj.load(fname_imagTrack) # Dictionary containing RCM quantities along RBSP-B satellite orbit # Read Dst data print("Reading Dst info") bsDstDict = kj.load(fname_dsts) # Dictionary containing 'UT' and 'Dst' # Dst is calculated from magnetosphere domain using Biot Savart, evaluated at the center of the Earth # Retrieve info from the slice file for easy reference later print("Gathering slice info") slcInfo = fileIO.SlcInfo.getInfo(fname_slice) # Calculate DPS-Dst within integration volume with radius r_DPS print("Calculating DPS-Dst within {} Re from slice file".format(r_DPS)) dpsDstDict = calcSlcDPSDst(fname_slice, slcInfo, # Slice filename and collected information r_DPS, # Radius to calculate DPS-Dst within bsDstDict['UT'][0], # Starting UT, used for trimming bsDstDict['UT'][-1] # Ending UT, used for trimming ) # Retrieve equatorial pressure from RCM print("Retrieving equatorial pressure for " + utPlotStr) momentInfo = fileIO.FileInfo.getInfo(fname_moments) momentData = getMomentData(fname_moments, momentInfo, utPlot) #--- # Plotting #--- # Prep figures fig = plt.figure(1,figsize=(12,9)) gs = gridspec.GridSpec(2,20, wspace=0.5, hspace=0.3) AxDsts = fig.add_subplot(gs[0,:10]) AxPressEq = fig.add_subplot(gs[0,10:]) AxIntensity = fig.add_subplot(gs[1,:-1]) AxInt_CB = fig.add_subplot(gs[1,-1]) # Plot Dst plotDsts(AxDsts, bsDstDict, dpsDstDict) plotTimeLine(AxDsts, utPlot) # Plot Pressure plotPressEq(AxPressEq, momentData, utPlot) # Plot satellite trajectory on top of pressure plot plotSatTrack(AxPressEq, imagTrack, utPlot) # Plot Intensity kv.genCB(AxInt_CB,int_norm,'Intensity\n[$cm^{-2} sr^{-1} s^{-1} keV^{-1}$]',cM=int_cmap,doVert=True,cbSz=14) plotIntensity(AxIntensity, imagTrack) plotTimeLine(AxIntensity, utPlot) plt.show()