AtHeartEngineer/kaipy
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaipy.git synced 2026-01-10 06:47:56 -05:00
Code Issues Packages Projects Releases Wiki Activity

Revert "Enabling pip install to add scripts to PATH (pull request #64)"

Browse Source
This commit is contained in:
Nikhil Rao
2025-02-20 21:40:35 +00:00
parent e9f8540c44
commit ec31f03b36
54 changed files with 557 additions and 790 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
kaipy/satcomp/test_cdas.py
View File

@@ -38,10 +38,7 @@ def test_getCdasData():
pass
def main():
if __name__ == "__main__":
test_getscIds()
test_getCdasData()
if __name__ == "__main__":
main()

0
kaipy/scripts/OHelio/__init__.py
View File

352
kaipy/scripts/OHelio/ih2oh.py
View File

@@ -12,226 +12,222 @@ import kaipy.gamhelio.lib.wsa as wsa
import kaipy.gamera.gamGrids as gg
def main():
#----------- PARSE ARGUMENTS ---------#
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('ConfigFileName',help='The name of the configuration file to use',default='startup.config')
args = parser.parse_args()
#----------- PARSE ARGUMENTS ---------#
#----------- PARSE ARGUMENTS ---------#
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('ConfigFileName',help='The name of the configuration file to use',default='startup.config')
args = parser.parse_args()
#----------- PARSE ARGUMENTS ---------#
# Read params from config file
prm = params.params(args.ConfigFileName)
Ng=prm.NO2
gamma = prm.gamma
# constants
mp = 1.67e-24
kb = 1.38e-16
# Read params from config file
prm = params.params(args.ConfigFileName)
Ng=prm.NO2
gamma = prm.gamma
# constants
mp = 1.67e-24
kb = 1.38e-16
#grid parameters
tMin = prm.tMin
tMax = prm.tMax
Rin = prm.Rin
Rout = prm.Rout
Ni = prm.Ni
Nj = prm.Nj
Nk = prm.Nk
#grid parameters
tMin = prm.tMin
tMax = prm.tMax
Rin = prm.Rin
Rout = prm.Rout
Ni = prm.Ni
Nj = prm.Nj
Nk = prm.Nk
#normalization in IH
B0 = prm.B0
n0 = prm.n0
V0 = B0/np.sqrt(4*np.pi*mp*n0)
T0 = B0*B0/4/np.pi/n0/kb #in K p = nkT
#normalization in IH
B0 = prm.B0
n0 = prm.n0
V0 = B0/np.sqrt(4*np.pi*mp*n0)
T0 = B0*B0/4/np.pi/n0/kb #in K p = nkT
print ("inner helio normalization")
print (B0, n0, V0, T0)
print ("inner helio normalization")
print (B0, n0, V0, T0)
#normalization in OH
B0OH = 5.e-5 # [Gs] 5 nT = 5.e-5 Gs
n0OH = 10 # [cm-3]
V0OH = B0OH/np.sqrt(4*np.pi*mp*n0OH) #Alfven speed at 1 AU 34.5 [km/s]
T0OH = B0OH*B0OH/4/np.pi/n0OH/kb #in K p = nkT
#normalization in OH
B0OH = 5.e-5 # [Gs] 5 nT = 5.e-5 Gs
n0OH = 10 # [cm-3]
V0OH = B0OH/np.sqrt(4*np.pi*mp*n0OH) #Alfven speed at 1 AU 34.5 [km/s]
T0OH = B0OH*B0OH/4/np.pi/n0OH/kb #in K p = nkT
print ("outer helio units")
print (B0OH, n0OH, V0OH, T0OH)
print ("outer helio units")
print (B0OH, n0OH, V0OH, T0OH)
#----------GENERATE HELIO GRID------
#----------GENERATE HELIO GRID------
print("Generating gamera-Ohelio grid ...")
print("Generating gamera-Ohelio grid ...")
X3,Y3,Z3 = gg.GenKSph(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
X3,Y3,Z3 = gg.GenKSph(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
#to generate non-uniform grid for GL cme (more fine in region 0.1-0.3 AU)
#X3,Y3,Z3 = gg.GenKSphNonUGL(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
gg.WriteGrid(X3,Y3,Z3,fOut=os.path.join(prm.GridDir,prm.gameraGridFile))
#to generate non-uniform grid for GL cme (more fine in region 0.1-0.3 AU)
#X3,Y3,Z3 = gg.GenKSphNonUGL(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
gg.WriteGrid(X3,Y3,Z3,fOut=os.path.join(prm.GridDir,prm.gameraGridFile))
print("Gamera-Ohelio grid ready!")
print("Gamera-Ohelio grid ready!")
#----------GENERATE HELIO GRID------
#----------GENERATE HELIO GRID------
############### READ GAMERA solution at 1 AU #####################
f = h5py.File(prm.wsaFile,'r')
#the latest Step saved in inner helio solution wsa.h5
step = 'Step#2'
############### READ GAMERA solution at 1 AU #####################
f = h5py.File(prm.wsaFile,'r')
#the latest Step saved in inner helio solution wsa.h5
step = 'Step#2'
f[step].attrs.keys()
Nphi, Nth, Nr = np.shape(f[step]['Vx'])
f[step].attrs.keys()
Nphi, Nth, Nr = np.shape(f[step]['Vx'])
#coordinates of cell centers
x = 0.125*(f['X'][:-1,:-1,:-1]+f['X'][:-1,:-1,1:]+f['X'][:-1,1:,:-1]+f['X'][:-1,1:,1:]+
f['X'][1:,:-1,:-1]+f['X'][1:,:-1,1:]+f['X'][1:,1:,:-1]+f['X'][1:,1:,1:])
y = 0.125*(f['Y'][:-1,:-1,:-1]+f['Y'][:-1,:-1,1:]+f['Y'][:-1,1:,:-1]+f['Y'][:-1,1:,1:]+
f['Y'][1:,:-1,:-1]+f['Y'][1:,:-1,1:]+f['Y'][1:,1:,:-1]+f['Y'][1:,1:,1:])
z = 0.125*(f['Z'][:-1,:-1,:-1]+f['Z'][:-1,:-1,1:]+f['Z'][:-1,1:,:-1]+f['Z'][:-1,1:,1:]+
f['Z'][1:,:-1,:-1]+f['Z'][1:,:-1,1:]+f['Z'][1:,1:,:-1]+f['Z'][1:,1:,1:])
#coordinates of cell centers
x = 0.125*(f['X'][:-1,:-1,:-1]+f['X'][:-1,:-1,1:]+f['X'][:-1,1:,:-1]+f['X'][:-1,1:,1:]+
f['X'][1:,:-1,:-1]+f['X'][1:,:-1,1:]+f['X'][1:,1:,:-1]+f['X'][1:,1:,1:])
y = 0.125*(f['Y'][:-1,:-1,:-1]+f['Y'][:-1,:-1,1:]+f['Y'][:-1,1:,:-1]+f['Y'][:-1,1:,1:]+
f['Y'][1:,:-1,:-1]+f['Y'][1:,:-1,1:]+f['Y'][1:,1:,:-1]+f['Y'][1:,1:,1:])
z = 0.125*(f['Z'][:-1,:-1,:-1]+f['Z'][:-1,:-1,1:]+f['Z'][:-1,1:,:-1]+f['Z'][:-1,1:,1:]+
f['Z'][1:,:-1,:-1]+f['Z'][1:,:-1,1:]+f['Z'][1:,1:,:-1]+f['Z'][1:,1:,1:])
r = np.sqrt(x[:]**2 + y[:]**2 + z[:]**2)
rxy = np.sqrt(x[:]**2 + y[:]**2)
r = np.sqrt(x[:]**2 + y[:]**2 + z[:]**2)
rxy = np.sqrt(x[:]**2 + y[:]**2)
#phi and theta of centers
theta = np.arccos(z/r)
phi = np.arctan2(y[:], x[:])
phi[phi<0]=phi[phi<0]+2*np.pi
#phi and theta of centers
theta = np.arccos(z/r)
phi = np.arctan2(y[:], x[:])
phi[phi<0]=phi[phi<0]+2*np.pi
theta_wsa_c = theta[0,:,0]
phi_wsa_c = phi[:,0,0]
theta_wsa_c = theta[0,:,0]
phi_wsa_c = phi[:,0,0]
print ("grid dimensions from 1 AU input solution")
print (theta_wsa_c.shape, phi_wsa_c.shape)
print ("grid dimensions from 1 AU input solution")
print (theta_wsa_c.shape, phi_wsa_c.shape)
#these are normilized according to inner helio normalization
Vr = (f[step]['Vx'][:]*x[:] + f[step]['Vy'][:]*y[:] + f[step]['Vz'][:]*z[:])/r[:]
#Br = f[step]['Br'][:]
Br = (f[step]['Bx'][:]*x[:] + f[step]['By'][:]*y[:] + f[step]['Bz'][:]*z[:])/r[:]
Rho = f[step]['D'][:]
T = f[step]['P'][:]/f[step]['D'][:]
#these are normilized according to inner helio normalization
Vr = (f[step]['Vx'][:]*x[:] + f[step]['Vy'][:]*y[:] + f[step]['Vz'][:]*z[:])/r[:]
#Br = f[step]['Br'][:]
Br = (f[step]['Bx'][:]*x[:] + f[step]['By'][:]*y[:] + f[step]['Bz'][:]*z[:])/r[:]
Rho = f[step]['D'][:]
T = f[step]['P'][:]/f[step]['D'][:]
#take solution from the last cell in i, already normilized
#use wsa variable names for now
bi_wsa = Br[:,:,Nr-1]
v_wsa = Vr[:,:,Nr-1]
n_wsa = Rho[:,:,Nr-1]
T_wsa = T[:,:,Nr-1]
#take solution from the last cell in i, already normilized
#use wsa variable names for now
bi_wsa = Br[:,:,Nr-1]
v_wsa = Vr[:,:,Nr-1]
n_wsa = Rho[:,:,Nr-1]
T_wsa = T[:,:,Nr-1]
print ("1AU arrays")
print (bi_wsa.shape, v_wsa.shape, n_wsa.shape, T_wsa.shape)
print ("1AU arrays")
print (bi_wsa.shape, v_wsa.shape, n_wsa.shape, T_wsa.shape)
#renormalize inner helio solution
bi_wsa = bi_wsa * B0/B0OH
n_wsa = n_wsa * n0/n0OH
v_wsa = v_wsa * V0/V0OH
T_wsa = T_wsa * T0
# keep temperature in K
#renormalize inner helio solution
bi_wsa = bi_wsa * B0/B0OH
n_wsa = n_wsa * n0/n0OH
v_wsa = v_wsa * V0/V0OH
T_wsa = T_wsa * T0
# keep temperature in K
#######Interpolate to GAMERA grid###########
# GAMERA GRID
with h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile),'r') as f:
x=f['X'][:]
y=f['Y'][:]
z=f['Z'][:]
#######Interpolate to GAMERA grid###########
# GAMERA GRID
with h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile),'r') as f:
x=f['X'][:]
y=f['Y'][:]
z=f['Z'][:]
xc = 0.125*(x[:-1,:-1,:-1]+x[:-1,1:,:-1]+x[:-1,:-1,1:]+x[:-1,1:,1:]
+x[1:,:-1,:-1]+x[1:,1:,:-1]+x[1:,:-1,1:]+x[1:,1:,1:])
yc = 0.125*(y[:-1,:-1,:-1]+y[:-1,1:,:-1]+y[:-1,:-1,1:]+y[:-1,1:,1:]
+y[1:,:-1,:-1]+y[1:,1:,:-1]+y[1:,:-1,1:]+y[1:,1:,1:])
zc = 0.125*(z[:-1,:-1,:-1]+z[:-1,1:,:-1]+z[:-1,:-1,1:]+z[:-1,1:,1:]
+z[1:,:-1,:-1]+z[1:,1:,:-1]+z[1:,:-1,1:]+z[1:,1:,1:])
xc = 0.125*(x[:-1,:-1,:-1]+x[:-1,1:,:-1]+x[:-1,:-1,1:]+x[:-1,1:,1:]
+x[1:,:-1,:-1]+x[1:,1:,:-1]+x[1:,:-1,1:]+x[1:,1:,1:])
yc = 0.125*(y[:-1,:-1,:-1]+y[:-1,1:,:-1]+y[:-1,:-1,1:]+y[:-1,1:,1:]
+y[1:,:-1,:-1]+y[1:,1:,:-1]+y[1:,:-1,1:]+y[1:,1:,1:])
zc = 0.125*(z[:-1,:-1,:-1]+z[:-1,1:,:-1]+z[:-1,:-1,1:]+z[:-1,1:,1:]
+z[1:,:-1,:-1]+z[1:,1:,:-1]+z[1:,:-1,1:]+z[1:,1:,1:])
# remove the ghosts from angular dimensions
R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2) # radius of the inner boundary
# remove the ghosts from angular dimensions
R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2) # radius of the inner boundary
#cell corners including ghost cells
r = np.sqrt(x[:]**2+y[:]**2+z[:]**2)
#cell corners including ghost cells
r = np.sqrt(x[:]**2+y[:]**2+z[:]**2)
#corners of physical cells
P = np.arctan2(y[Ng:-Ng,Ng:-Ng,:],x[Ng:-Ng,Ng:-Ng,:])
P[P<0] += 2*np.pi
#P = P % (2*np.pi) # sometimes the very first point may be a very
# small negative number, which the above call sets
# to 2*pi. This takes care of it.
T = np.arccos(z[Ng:-Ng,Ng:-Ng,:]/r[Ng:-Ng,Ng:-Ng,:])
#corners of physical cells
P = np.arctan2(y[Ng:-Ng,Ng:-Ng,:],x[Ng:-Ng,Ng:-Ng,:])
P[P<0] += 2*np.pi
#P = P % (2*np.pi) # sometimes the very first point may be a very
# small negative number, which the above call sets
# to 2*pi. This takes care of it.
T = np.arccos(z[Ng:-Ng,Ng:-Ng,:]/r[Ng:-Ng,Ng:-Ng,:])
#grid (corners) for output into innerbc.h5
P_out = P[:,:,0:Ng+1]
T_out = T[:,:,0:Ng+1]
R_out = r[Ng:-Ng,Ng:-Ng,0:Ng+1]
print ("shapes of output phi and theta ", P_out.shape, T_out.shape, R_out.shape)
#grid (corners) for output into innerbc.h5
P_out = P[:,:,0:Ng+1]
T_out = T[:,:,0:Ng+1]
R_out = r[Ng:-Ng,Ng:-Ng,0:Ng+1]
print ("shapes of output phi and theta ", P_out.shape, T_out.shape, R_out.shape)
#centers
Rc = np.sqrt(xc[Ng:-Ng,Ng:-Ng,:]**2+yc[Ng:-Ng,Ng:-Ng,:]**2+zc[Ng:-Ng,Ng:-Ng,:]**2)
Pc = np.arctan2(yc[Ng:-Ng,Ng:-Ng,:],xc[Ng:-Ng,Ng:-Ng,:])
Pc[Pc<0] += 2*np.pi
Tc = np.arccos(zc[Ng:-Ng,Ng:-Ng,:]/Rc)
#centers
Rc = np.sqrt(xc[Ng:-Ng,Ng:-Ng,:]**2+yc[Ng:-Ng,Ng:-Ng,:]**2+zc[Ng:-Ng,Ng:-Ng,:]**2)
Pc = np.arctan2(yc[Ng:-Ng,Ng:-Ng,:],xc[Ng:-Ng,Ng:-Ng,:])
Pc[Pc<0] += 2*np.pi
Tc = np.arccos(zc[Ng:-Ng,Ng:-Ng,:]/Rc)
# this is fast and better than griddata in that it nicely extrapolates boundaries:
fbi = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,bi_wsa,kx=1,ky=1)
br = fbi(Pc[:,0,0],Tc[0,:,0])
# this is fast and better than griddata in that it nicely extrapolates boundaries:
fbi = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,bi_wsa,kx=1,ky=1)
br = fbi(Pc[:,0,0],Tc[0,:,0])
############### SMOOTHING #####################
if not prm.gaussSmoothWidth==0:
import astropy
from astropy.convolution import convolve,Gaussian2DKernel
############### SMOOTHING #####################
if not prm.gaussSmoothWidth==0:
import astropy
from astropy.convolution import convolve,Gaussian2DKernel
gauss=Gaussian2DKernel(width=prm.gaussSmoothWidth)
br =astropy.convolution.convolve(br,gauss,boundary='extend')
gauss=Gaussian2DKernel(width=prm.gaussSmoothWidth)
br =astropy.convolution.convolve(br,gauss,boundary='extend')
############### INTERPOLATE AND DUMP #####################
fv = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,v_wsa,kx=1,ky=1)
vr = fv(Pc[:,0,0],Tc[0,:,0])
############### INTERPOLATE AND DUMP #####################
fv = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,v_wsa,kx=1,ky=1)
vr = fv(Pc[:,0,0],Tc[0,:,0])
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa,kx=1,ky=1)
rho = f(Pc[:,0,0],Tc[0,:,0])
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa,kx=1,ky=1)
rho = f(Pc[:,0,0],Tc[0,:,0])
fT = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,T_wsa,kx=1,ky=1)
temp = fT(Pc[:,0,0],Tc[0,:,0])
fT = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,T_wsa,kx=1,ky=1)
temp = fT(Pc[:,0,0],Tc[0,:,0])
fbi = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],br,kx=1,ky=1)
fv = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],vr,kx=1,ky=1)
fbi = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],br,kx=1,ky=1)
fv = interpolate.RectBivariateSpline(Pc[:,0,0],Tc[0,:,0],vr,kx=1,ky=1)
br_kface = fbi(P[:,0,0],Tc[0,:,0])
vr_kface = fv (P[:,0,0],Tc[0,:,0])
br_kface = fbi(P[:,0,0],Tc[0,:,0])
vr_kface = fv (P[:,0,0],Tc[0,:,0])
# Scale inside ghost region
(vr,vr_kface,rho,temp,br,br_kface) = [np.dstack(prm.NO2*[var]) for var in (vr,vr_kface,rho,temp,br,br_kface)]
rho *= (R0/Rc[0,0,:Ng])**2
br *= (R0/Rc[0,0,:Ng])**2
br_kface *= (R0/Rc[0,0,:Ng])**2
# Scale inside ghost region
(vr,vr_kface,rho,temp,br,br_kface) = [np.dstack(prm.NO2*[var]) for var in (vr,vr_kface,rho,temp,br,br_kface)]
rho *= (R0/Rc[0,0,:Ng])**2
br *= (R0/Rc[0,0,:Ng])**2
br_kface *= (R0/Rc[0,0,:Ng])**2
#FIX:
#For now I use wsa.h5 which do not have mjd inside
#so hardcoded using WSA fits value + 200*4637/60/60/24 [number of days]
mjd_c = 58005.83415
#FIX:
#For now I use wsa.h5 which do not have mjd inside
#so hardcoded using WSA fits value + 200*4637/60/60/24 [number of days]
mjd_c = 58005.83415
print ("writing out innerbc.h5...")
with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
hf.attrs["MJD"] = mjd_c
hf.create_dataset("vr",data=vr)
hf.create_dataset("vr_kface",data=vr_kface)
hf.create_dataset("rho",data=rho)
hf.create_dataset("temp",data=temp)
hf.create_dataset("br",data=br)
hf.create_dataset("br_kface",data=br_kface)
hf.close()
print ("writing out innerbc.h5...")
with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
hf.attrs["MJD"] = mjd_c
hf.create_dataset("vr",data=vr)
hf.create_dataset("vr_kface",data=vr_kface)
hf.create_dataset("rho",data=rho)
hf.create_dataset("temp",data=temp)
hf.create_dataset("br",data=br)
hf.create_dataset("br_kface",data=br_kface)
hf.close()
#innerbc to plot in Paraview
with h5py.File(os.path.join(prm.IbcDir,'innerbc_OHighostgr.h5'),'w') as hfg:
hfg.create_dataset("X", data=P_out)
hfg.create_dataset("Y", data=T_out)
hfg.create_dataset("Z", data=R_out)
grname = "Step#0"
grp = hfg.create_group(grname)
grp.attrs.create("MJD", mjd_c)
#grp.attrs.create("time", time_sec)
grp.create_dataset("vr",data=vr)
grp.create_dataset("vr_kface",data=vr_kface)
grp.create_dataset("rho",data=rho)
grp.create_dataset("temp",data=temp)
grp.create_dataset("br",data=br)
grp.create_dataset("br_kface",data=br_kface)
hfg.close
if __name__ == "__main__":
main()
#innerbc to plot in Paraview
with h5py.File(os.path.join(prm.IbcDir,'innerbc_OHighostgr.h5'),'w') as hfg:
hfg.create_dataset("X", data=P_out)
hfg.create_dataset("Y", data=T_out)
hfg.create_dataset("Z", data=R_out)
grname = "Step#0"
grp = hfg.create_group(grname)
grp.attrs.create("MJD", mjd_c)
#grp.attrs.create("time", time_sec)
grp.create_dataset("vr",data=vr)
grp.create_dataset("vr_kface",data=vr_kface)
grp.create_dataset("rho",data=rho)
grp.create_dataset("temp",data=temp)
grp.create_dataset("br",data=br)
grp.create_dataset("br_kface",data=br_kface)
hfg.close

0
kaipy/scripts/__init__.py
View File

0
kaipy/scripts/datamodel/__init__.py
View File

5
kaipy/scripts/datamodel/helioSatComp.py
View File

@@ -115,7 +115,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
"""Begin main program."""
# Set up the command-line parser.
@@ -341,6 +341,3 @@ def main():
# if verbose:
# print("Plotting %s trajectory in spacecraft frame to %s." % (sc_id, plot_file_path))
# kv.helioTrajPlot(plot_file_path, sc_id, sc_data)
if __name__ == "__main__":
main()

5
kaipy/scripts/datamodel/msphParallelSatComp.py
View File

@@ -34,7 +34,7 @@ import kaipy.kaiTools as kaiTools
import kaipy.kaijson as kj
import kaipy.satcomp.scutils as scutils
def main():
if __name__ == '__main__':
MainS = """Extracts information from satellite trajectory for various
spacecraft. Space craft data is pulled from CDAWeb. Output CDF files
contain data pulled from CDAWeb along with data extracted from GAMERA.
@@ -111,6 +111,3 @@ def main():
subprocess.run(['rm',os.path.join(fdir,'err.*.txt')])
print('All done!')
if __name__ == '__main__':
main()

4
kaipy/scripts/datamodel/msphPbsSatComp.py
View File

@@ -24,7 +24,7 @@ import spacepy.datamodel as dm
def main():
if __name__ == '__main__':
MainS = """Checks the run for satellites with data avilable and then
sets up PBS job scripts for running interpolation in parallel."""
@@ -179,5 +179,3 @@ def main():
if __name__ == '__main__':
main()

5
kaipy/scripts/datamodel/msphSatComp.py
View File

@@ -19,7 +19,7 @@ import kaipy.kaijson as kj
import kaipy.satcomp.scutils as scutils
import spacepy.datamodel as dm
def main():
if __name__ == '__main__':
MainS = """Extracts information from satellite trajectory for various
spacecraft. Space craft data is pulled from CDAWeb. Output CDF files
contain data pulled from CDAWeb along with data extracted from GAMERA.
@@ -118,6 +118,3 @@ def main():
print('Plotting trajectory to',plotname)
kv.trajPlot(plotname,scId,data,toRe)
if __name__ == '__main__':
main()

5
kaipy/scripts/datamodel/rbspSCcomp.py
View File

@@ -64,7 +64,7 @@ def getJScl(Bmag,Beq,en=2.0):
It[n] = Ic.sum()/I0
return It
def main():
if __name__ == "__main__":
#Defaults
fdir = os.getcwd()
ftag = "eRBpsd.ps.h5"
@@ -314,5 +314,4 @@ def main():
kv.savePic('1dJComp.png')
if __name__ == "__main__":
main()

5
kaipy/scripts/datamodel/rcm_rbsp_satcomp.py → kaipy/scripts/datamodel/rcm-rbsp-satcomp.py
View File

@@ -23,7 +23,7 @@ def fmtTKL(AxTKL):
AxTKL.xaxis.labelpad = -1
AxTKL.title.set_text(str(ut_tkl[0]))
def main():
if __name__=="__main__":
fdir = os.getcwd()
ftag = "msphere"
trtag = "RBSP-%s_MAGNETOMETER_1SEC-GSM_EMFISIS-L3.sc.h5" # Spacecraft trajectory and values along track
@@ -269,6 +269,3 @@ def main():
plt.show()
if __name__ == "__main__":
main()

0
kaipy/scripts/postproc/__init__.py
View File

7
kaipy/scripts/postproc/block_genmpiXDMF.py
View File

@@ -4,7 +4,7 @@
import argparse
from argparse import RawTextHelpFormatter
import numpy as np
import kaipy.gamera.gampp as gampp
import kaipy.gamera.block_gampp as gampp
import xml.etree.ElementTree as et
import xml.dom.minidom
import kaipy.kaiH5 as kh5
@@ -52,7 +52,7 @@ def genName(bStr, i, j, k, Ri, Rj, Rk,isOld):
return fID
def main():
if __name__ == "__main__":
# Defaults
fdir = os.getcwd()
ftag = "msphere"
@@ -228,6 +228,3 @@ def main():
# Prep for next step
tOut = tOut+1
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/embiggen.py
View File

@@ -9,7 +9,7 @@ import h5py
import os
import kaipy.kaiH5 as kh5
def main():
if __name__ == "__main__":
dIn = os.getcwd()
#Input tiling
@@ -191,8 +191,7 @@ def main():
os.remove(fTmp)
os.remove(fTmp2X)
if __name__ == "__main__":
main()
# #!/usr/bin/env python

5
kaipy/scripts/postproc/embiggenMIX.py
View File

@@ -9,7 +9,7 @@ import h5py
import os
import kaipy.kaiH5 as kh5
def main():
if __name__ == "__main__":
dIn = os.getcwd()
@@ -68,6 +68,3 @@ def main():
#Now get
iH5.close()
oH5.close()
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/embiggenRCM.py
View File

@@ -9,7 +9,7 @@ import h5py
import os
import kaipy.kaiH5 as kh5
def main():
if __name__ == "__main__":
dIn = os.getcwd()
@@ -109,6 +109,3 @@ def main():
print("\tRCMSIZEI = %d"%(Nri))
print("\tRCMSIZEJ = %d"%(Nrj))
print("\tRCMSIZEK = %d"%(Nk))
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/embiggenVOLT.py
View File

@@ -9,7 +9,7 @@ import h5py
import os
import kaipy.kaiH5 as kh5
def main():
if __name__ == "__main__":
dIn = os.getcwd()
@@ -68,6 +68,3 @@ def main():
#Now get
iH5.close()
oH5.close()
if __name__ == "__main__":
main()

4
kaipy/scripts/postproc/genXDMF.py
View File

@@ -131,7 +131,7 @@ def addRCMVars(Grid, dimInfo, rcmInfo, sID):
kxmf.addHyperslab(Grid,vName_k,mr_vDimStr,dimStr,startStr,strideStr,numStr,r_vDimStr,text)
def main():
if __name__ == "__main__":
outfname = ''
@@ -300,5 +300,3 @@ def main():
with open(fOutXML,"w") as f:
f.write(xmlStr)
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/genXLine.py
View File

@@ -64,7 +64,7 @@ def getNum(fIn,n,m):
lId = "Line#%d"%(m)
Np = hf[gId][lId].attrs["Np"]
return Np
def main():
if __name__ == "__main__":
#Set defaults
parser = argparse.ArgumentParser(description="Generates XDMF file from CHIMP tracer HDF5 output")
parser.add_argument('h5F',nargs=1,type=str,metavar='tracer.h5',help="Filename of CHIMP tracer HDF5 Output")
@@ -186,6 +186,3 @@ def main():
xmlStr = xml.dom.minidom.parseString(et.tostring(Xdmf)).toprettyxml(indent=" ")
with open(fOutXML,"w") as f:
f.write(xmlStr)
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/genmpiXDMF.py
View File

@@ -12,7 +12,7 @@ import kaipy.kaiH5 as kh5
import os
def main():
if __name__ == "__main__":
#Defaults
fdir = os.getcwd()
ftag = "msphere"
@@ -177,6 +177,3 @@ def main():
#Prep for next step
tOut = tOut+1
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/numSteps.py
View File

@@ -11,7 +11,7 @@ def MJD2Str(m0):
dtObj = Time(m0,format='mjd').datetime
tStr = dtObj.strftime("%H:%M:%S") + " " + dtObj.strftime("%m/%d/%Y")
return tStr
def main():
if __name__ == "__main__":
#Defaults
MainS = """Identifies the domain (in steps and time) of a Kaiju HDF-5 file"""
@@ -48,6 +48,3 @@ def main():
print("\tGit: Hash = %s / Branch = %s"%(hStr,bStr))
#---------------------
if __name__ == "__main__":
main()

6
kaipy/scripts/postproc/pitmerge.py
View File

@@ -44,7 +44,7 @@ def createfile(fIn,fOut,doLink=False):
return oH5
def main():
if __name__ == "__main__":
dIn = os.getcwd()
runid = "msphere"
@@ -157,7 +157,3 @@ def main():
#Done
oH5.close()
if __name__ == "__main__":
main()

4
kaipy/scripts/postproc/printResTimes.py
View File

@@ -39,7 +39,7 @@ def getKVsFromFile(fName):
attrs[k] = v
return attrs
def main():
if __name__=='__main__':
idStr_noMPI = ".gam.Res.*.h5"
idStrMPI = "_0*_0*_0*_0000_0000_0000.gam.Res.*.h5"
@@ -102,5 +102,3 @@ def main():
print(formatString)
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/run_supermag_comparison.py
View File

@@ -220,7 +220,7 @@ def compute_ground_delta_B(runid, mpixml=None):
return delta_B_file
def main():
if __name__ == "__main__":
"""Begin main program."""
# Set up the command-line parser.
@@ -319,6 +319,3 @@ def main():
sm.MakeContourPlots(SM, SMinterp, maxx = 1000, fignumber=2)
contour_plot_file = runid + "_contours.png"
plt.savefig(contour_plot_file)
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/slimFL.py
View File

@@ -24,7 +24,7 @@ def createfile(iH5,fOut):
oH5.create_dataset(sQ,data=iH5[sQ])
return oH5
def main():
if __name__ == "__main__":
#Set defaults
ns = 0
ne = -1 #Proxy for last entry
@@ -118,5 +118,4 @@ def main():
iH5.close()
oH5.close()
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/slimh5.py
View File

@@ -59,7 +59,7 @@ def createfile(iH5,fOut):
return oH5
def main():
if __name__ == "__main__":
#Set defaults
ns = -1
ne = -1 #Proxy for last entry
@@ -183,6 +183,3 @@ def main():
#Close up
iH5.close()
oH5.close()
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/slimh5_classic.py
View File

@@ -59,7 +59,7 @@ def createfile(iH5,fOut):
return oH5
def main():
if __name__ == "__main__":
#Set defaults
ns = -1
ne = -1 #Proxy for last entry
@@ -175,6 +175,3 @@ def main():
#Close up
iH5.close()
oH5.close()
if __name__ == "__main__":
main()

5
kaipy/scripts/postproc/ut2mjd.py
View File

@@ -6,7 +6,7 @@ import datetime
fmt='%m/%d/%Y, %H:%M:%S'
def main():
if __name__ == "__main__":
t0="2010-01-01T00:00:00"
fmt='%Y-%m-%dT%H:%M:%S'
@@ -26,6 +26,3 @@ def main():
mjd = Time(ut).mjd
print("%s (UT) => %f (MJD)"%(utStr,mjd))
if __name__ == "__main__":
main()

4
kaipy/scripts/preproc/XMLGenerator.py
View File

@@ -367,7 +367,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
"""Convert a .ini file to a .xml file."""
# Set up the command-line parser.
parser = create_command_line_parser()
@@ -385,5 +385,3 @@ def main():
print("Converting %s to XML output %s." % (args.ini_file, args.xml_file))
create_xml_template(args.ini_file, args.xml_file)
if __name__ == "__main__":
main()

0
kaipy/scripts/preproc/__init__.py
View File

7
kaipy/scripts/preproc/cda2wind.py
View File

@@ -148,7 +148,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
maxf107 = 300.0
minMfast = 1.5
@@ -578,8 +578,3 @@ def main():
else:
raise Exception('Error: Misunderstood output file format.')
if __name__ == '__main__':
main()

8
kaipy/scripts/preproc/genLFM.py
View File

@@ -52,10 +52,10 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults asdf
Nc0 = 8 #Number of outer i cells to cut out from LFM grid (OCT)
fIn = "lfmG"
fIn = "./lfmG"
doEpsY = True
TINY = 1.0e-8
gLabs = ["Double","Quad","Oct","Hex"]
@@ -120,7 +120,3 @@ def main():
if (doViz):
gg.VizGrid(XX,YY,xxG,yyG,fOut=fOut,doGhost=doVizG)
#gg.genRing(XX,YY,Nk=Nk,Tol=1.0,doVerb=True)
if __name__ == "__main__":
main()

5
kaipy/scripts/preproc/genRCM.py
View File

@@ -68,7 +68,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
# Set up the command-line parser.
parser = create_command_line_parser()
@@ -142,6 +142,3 @@ def main():
if __name__ == "__main__":
main()

660
kaipy/scripts/preproc/wsa2TDgamera.py
View File

@@ -84,405 +84,403 @@ def plotBc(wsa_file, phi, theta, var1, var2, var3, var4):
plt.suptitle(year + ':' + month + ':' + day, y=0.85)
plt.savefig(wsaFile[:-5]+'_bc.png', bbox_inches='tight')
def main():
#----------- PARSE ARGUMENTS ---------#
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('ConfigFileName',help='The name of the configuration file to use',default='startup.config')
args = parser.parse_args()
#----------- PARSE ARGUMENTS ---------#
#----------- PARSE ARGUMENTS ---------#
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('ConfigFileName',help='The name of the configuration file to use',default='startup.config')
args = parser.parse_args()
#----------- PARSE ARGUMENTS ---------#
# Read params from config file
prm = params.params(args.ConfigFileName)
(ni,nj,nk) = (prm.Ni,prm.Nj,prm.Nk)
Ng = prm.NO2
# Read params from config file
prm = params.params(args.ConfigFileName)
(ni,nj,nk) = (prm.Ni,prm.Nj,prm.Nk)
Ng = prm.NO2
#grid parameters
tMin = prm.tMin
tMax = prm.tMax
Rin = prm.Rin
Rout = prm.Rout
Ni = prm.Ni
Nj = prm.Nj
Nk = prm.Nk
#grid parameters
tMin = prm.tMin
tMax = prm.tMax
Rin = prm.Rin
Rout = prm.Rout
Ni = prm.Ni
Nj = prm.Nj
Nk = prm.Nk
#----------GENERATE HELIO GRID------
#----------GENERATE HELIO GRID------
print("Generating gamera-helio grid ...")
print("Generating gamera-helio grid ...")
X3,Y3,Z3 = gg.GenKSph(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
X3,Y3,Z3 = gg.GenKSph(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
#to generate non-uniform grid for GL cme (more fine in region 0.1-0.3 AU)
#X3,Y3,Z3 = gg.GenKSphNonUGL(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
gg.WriteGrid(X3,Y3,Z3,fOut=os.path.join(prm.GridDir,prm.gameraGridFile))
#to generate non-uniform grid for GL cme (more fine in region 0.1-0.3 AU)
#X3,Y3,Z3 = gg.GenKSphNonUGL(Ni=Ni,Nj=Nj,Nk=Nk,Rin=Rin,Rout=Rout,tMin=tMin,tMax=tMax)
gg.WriteGrid(X3,Y3,Z3,fOut=os.path.join(prm.GridDir,prm.gameraGridFile))
print("Gamera-helio grid ready!")
print("Gamera-helio grid ready!")
#----------GENERATE HELIO GRID------
#----------GENERATE HELIO GRID------
# [EP] sorted list of WSA files
wsaFiles = sorted(glob.glob(os.path.join(prm.adaptdir,prm.adaptWildCard)))
# [EP] sorted list of WSA files
wsaFiles = sorted(glob.glob(os.path.join(prm.adaptdir,prm.adaptWildCard)))
print(wsaFiles)
print(wsaFiles)
# [EP] electric fields on edges
#+2 in j directions, two ghost cells at start and end
et_save = np.zeros( (nj+2,nk+1) )
ep_save = np.zeros( (nj+1+2,nk) )
# [EP] electric fields on edges
#+2 in j directions, two ghost cells at start and end
et_save = np.zeros( (nj+2,nk+1) )
ep_save = np.zeros( (nj+1+2,nk) )
#Normalization
Vnorm = 1.e5 #cm/s => km/s
Bnorm = 1.e-5 #Gs => nT
mp = 1.67e-24
kblts = 1.38e-16
#Normalization
Vnorm = 1.e5 #cm/s => km/s
Bnorm = 1.e-5 #Gs => nT
mp = 1.67e-24
kblts = 1.38e-16
#open innerbcTD.h5 for ouput
with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
#open innerbcTD.h5 for ouput
with h5py.File(os.path.join(prm.IbcDir,prm.gameraIbcFile),'w') as hf:
#[EP] going through the list of WSA files
for (fcount,wsaFile) in enumerate(wsaFiles):
#print(fcount)
############### WSA STUFF #####################
isFirstFile = (wsaFile == wsaFiles[0])
#[EP] reading WSA file
jd_c,phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa,v_wsa,n_wsa,T_wsa = wsa.read(wsaFile,prm.densTempInfile,prm.normalized, verbose = isFirstFile)
#bi_wsa in Gs CGS units
#v_wsa in cm/s
#n_wsa in g/cm-3
#T_wsa in K
#[EP] going through the list of WSA files
for (fcount,wsaFile) in enumerate(wsaFiles):
#print(fcount)
############### WSA STUFF #####################
isFirstFile = (wsaFile == wsaFiles[0])
#[EP] reading WSA file
jd_c,phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa,v_wsa,n_wsa,T_wsa = wsa.read(wsaFile,prm.densTempInfile,prm.normalized, verbose = isFirstFile)
#bi_wsa in Gs CGS units
#v_wsa in cm/s
#n_wsa in g/cm-3
#T_wsa in K
#convert julian date from wsa fits into modified julian date
mjd_c = jd_c - 2400000.5
#convert julian date from wsa fits into modified julian date
mjd_c = jd_c - 2400000.5
if isFirstFile:
#take JD from the first wsa file
jd0 = jd_c
if isFirstFile:
#take JD from the first wsa file
jd0 = jd_c
# GAMERA GRID
# read GAMERA grid from innerbc.h5
print ('reading heliogrid.h5 ...')
f = h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile), 'r')
#Nphi, Nth, Nr = np.shape(f['X'])
#corners
x = f['X'][:]
y = f['Y'][:]
z = f['Z'][:]
#centers
xc = 0.125*(f['X'][:-1,:-1,:-1]+f['X'][:-1,:-1,1:]+f['X'][:-1,1:,:-1]+f['X'][:-1,1:,1:]+
f['X'][1:,:-1,:-1]+f['X'][1:,:-1,1:]+f['X'][1:,1:,:-1]+f['X'][1:,1:,1:])
yc = 0.125*(f['Y'][:-1,:-1,:-1]+f['Y'][:-1,:-1,1:]+f['Y'][:-1,1:,:-1]+f['Y'][:-1,1:,1:]+
f['Y'][1:,:-1,:-1]+f['Y'][1:,:-1,1:]+f['Y'][1:,1:,:-1]+f['Y'][1:,1:,1:])
zc = 0.125*(f['Z'][:-1,:-1,:-1]+f['Z'][:-1,:-1,1:]+f['Z'][:-1,1:,:-1]+f['Z'][:-1,1:,1:]+
f['Z'][1:,:-1,:-1]+f['Z'][1:,:-1,1:]+f['Z'][1:,1:,:-1]+f['Z'][1:,1:,1:])
#radius of inner boundary. Index order [k,j,i]
R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2)
#[EP for testing]
#cell corners including ghost cells
r = np.sqrt(x[:]**2+y[:]**2+z[:]**2)
rxy = np.sqrt(x[:]**2+y[:]**2)
# GAMERA GRID
# read GAMERA grid from innerbc.h5
# remove the ghosts from angular dimensions (corners)
P = np.arctan2(y[Ng:-Ng,Ng:-Ng,:],x[Ng:-Ng,Ng:-Ng,:])
P [ P < 0] += 2*np.pi
T = np.arccos(z[Ng:-Ng,Ng:-Ng,:]/r[Ng:-Ng,Ng:-Ng,:])
print ('reading heliogrid.h5 ...')
f = h5py.File(os.path.join(prm.GridDir,prm.gameraGridFile), 'r')
#Nphi, Nth, Nr = np.shape(f['X'])
#corners
x = f['X'][:]
y = f['Y'][:]
z = f['Z'][:]
#grid for output into innerbc.h5
P_out = P[:,:,0:Ng+1]
T_out = T[:,:,0:Ng+1]
R_out = r[Ng:-Ng,Ng:-Ng,0:Ng+1]
print ("shapes of output phi and theta ", P_out.shape, T_out.shape, R_out.shape)
#centers
xc = 0.125*(f['X'][:-1,:-1,:-1]+f['X'][:-1,:-1,1:]+f['X'][:-1,1:,:-1]+f['X'][:-1,1:,1:]+
f['X'][1:,:-1,:-1]+f['X'][1:,:-1,1:]+f['X'][1:,1:,:-1]+f['X'][1:,1:,1:])
yc = 0.125*(f['Y'][:-1,:-1,:-1]+f['Y'][:-1,:-1,1:]+f['Y'][:-1,1:,:-1]+f['Y'][:-1,1:,1:]+
f['Y'][1:,:-1,:-1]+f['Y'][1:,:-1,1:]+f['Y'][1:,1:,:-1]+f['Y'][1:,1:,1:])
zc = 0.125*(f['Z'][:-1,:-1,:-1]+f['Z'][:-1,:-1,1:]+f['Z'][:-1,1:,:-1]+f['Z'][:-1,1:,1:]+
f['Z'][1:,:-1,:-1]+f['Z'][1:,:-1,1:]+f['Z'][1:,1:,:-1]+f['Z'][1:,1:,1:])
#centers spherical grid excluding ghosts in angular directions
#Rc = np.sqrt(xc[Ng:-Ng, Ng:-Ng,:]**2 + yc[Ng:-Ng, Ng:-Ng,:]**2 + zc[Ng:-Ng, Ng:-Ng,:]**2)
#Pc = np.arctan2(yc[Ng:-Ng, Ng:-Ng,:], xc[Ng:-Ng, Ng:-Ng,:])
#Tc = np.arccos(zc[Ng:-Ng,Ng:-Ng,:]/Rc)
#radius of inner boundary. Index order [k,j,i]
R0 = np.sqrt(x[0,0,Ng]**2+y[0,0,Ng]**2+z[0,0,Ng]**2)
#include one extra cell in j direction at start and end
Pg = Ng-1
Rc = np.sqrt(xc[Ng:-Ng, Pg:-Pg,:]**2 + yc[Ng:-Ng, Pg:-Pg,:]**2 + zc[Ng:-Ng, Pg:-Pg,:]**2)
Pc = np.arctan2(yc[Ng:-Ng, Pg:-Pg,:], xc[Ng:-Ng, Pg:-Pg,:])
Tc = np.arccos(zc[Ng:-Ng,Pg:-Pg,:]/Rc)
#[EP for testing]
#cell corners including ghost cells
r = np.sqrt(x[:]**2+y[:]**2+z[:]**2)
rxy = np.sqrt(x[:]**2+y[:]**2)
# remove the ghosts from angular dimensions (corners)
P = np.arctan2(y[Ng:-Ng,Ng:-Ng,:],x[Ng:-Ng,Ng:-Ng,:])
P [ P < 0] += 2*np.pi
T = np.arccos(z[Ng:-Ng,Ng:-Ng,:]/r[Ng:-Ng,Ng:-Ng,:])
Pc [Pc < 0] += 2*np.pi
#GAMERA grid centers at the inner boundary, 1D array
phi = Pc[:,0,0]
theta = Tc[0,:,0]
#grid for output into innerbc.h5
P_out = P[:,:,0:Ng+1]
T_out = T[:,:,0:Ng+1]
R_out = r[Ng:-Ng,Ng:-Ng,0:Ng+1]
print ("shapes of output phi and theta ", P_out.shape, T_out.shape, R_out.shape)
#debug
#print (phi)
#print (theta)
# what exactly does this do???
pois = poisson.poisson(theta,phi)
#time from the 1st wsa map in seconds
time_sec = (jd_c - jd0)*24.*60.*60.
omega=2*np.pi/prm.Tsolar*(25.38/27.27)
#shift of wsa maps needed if wsa solutions are provided in inertial frame (folder UPDATED)
#if wsa solutions are provided in rotating carrington frame (folder CARR), no need to shift.
#shift phi coordinate in wsa data according to the shift of the wsa map relative to the first one
#wsa maps move to the right with cadence 1 day
phi_prime=(phi_wsa_c-omega*prm.adaptCadence*fcount)%(2*np.pi)
#looking for index of shift
if np.where(np.ediff1d(phi_prime)<0)[0].size!=0: #for the first map size =0, for other maps size=1
ind0=np.where(np.ediff1d(phi_prime)<0)[0][0]+1
#print 'ind = ', ind0
else:
ind0=0 # this is for the first map
#centers spherical grid excluding ghosts in angular directions
#shifting phi_prime to the left
phi_prime=np.roll(phi_prime,-ind0)
bi_wsa_rolled=np.roll(bi_wsa,-ind0,axis=1)
v_wsa_rolled=np.roll(v_wsa,-ind0,axis=1)
n_wsa_rolled=np.roll(n_wsa,-ind0,axis=1)
T_wsa_rolled=np.roll(T_wsa,-ind0,axis=1)
#Rc = np.sqrt(xc[Ng:-Ng, Ng:-Ng,:]**2 + yc[Ng:-Ng, Ng:-Ng,:]**2 + zc[Ng:-Ng, Ng:-Ng,:]**2)
#Pc = np.arctan2(yc[Ng:-Ng, Ng:-Ng,:], xc[Ng:-Ng, Ng:-Ng,:])
#Tc = np.arccos(zc[Ng:-Ng,Ng:-Ng,:]/Rc)
#plot br from original wsa map (top plot) and shifted to the origin map(bottom plot)
#changes in time in the bottom plot are purely due to time-dependent variations of B_r (rotation is eliminted)
#include one extra cell in j direction at start and end
Pg = Ng-1
Rc = np.sqrt(xc[Ng:-Ng, Pg:-Pg,:]**2 + yc[Ng:-Ng, Pg:-Pg,:]**2 + zc[Ng:-Ng, Pg:-Pg,:]**2)
Pc = np.arctan2(yc[Ng:-Ng, Pg:-Pg,:], xc[Ng:-Ng, Pg:-Pg,:])
Tc = np.arccos(zc[Ng:-Ng,Pg:-Pg,:]/Rc)
plot(wsaFile, bi_wsa, bi_wsa_rolled)
##plot(wsaFile, v_wsa, v_wsa_rolled)
###INTERPOLATION OF ROLLED WSA MAPS TO GAMERA GRID phi-theta####
Pc [Pc < 0] += 2*np.pi
#GAMERA grid centers at the inner boundary, 1D array
phi = Pc[:,0,0]
theta = Tc[0,:,0]
# bivariate spline approximation over a rectangular mesh
fbi = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,bi_wsa_rolled.T,kx=1,ky=1)
# interpolation to Gamera grid
br = fbi(phi,theta)
#debug
#print (phi)
#print (theta)
#Next Slava used SMOOTHING for br for the paper, we do not need it for now
# what exactly does this do???
pois = poisson.poisson(theta,phi)
fv = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,v_wsa_rolled.T,kx=1,ky=1)
vr = fv(phi,theta)
#time from the 1st wsa map in seconds
time_sec = (jd_c - jd0)*24.*60.*60.
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa_rolled.T,kx=1,ky=1)
rho = f(phi,theta)
omega=2*np.pi/prm.Tsolar*(25.38/27.27)
#shift of wsa maps needed if wsa solutions are provided in inertial frame (folder UPDATED)
#if wsa solutions are provided in rotating carrington frame (folder CARR), no need to shift.
#shift phi coordinate in wsa data according to the shift of the wsa map relative to the first one
#wsa maps move to the right with cadence 1 day
#not interpolating temperature, calculating sound speed cs
#assuming uniform total pressure Rho_max*k*T0 = p+Br^2/8pi
phi_prime=(phi_wsa_c-omega*prm.adaptCadence*fcount)%(2*np.pi)
#looking for index of shift
if np.where(np.ediff1d(phi_prime)<0)[0].size!=0: #for the first map size =0, for other maps size=1
ind0=np.where(np.ediff1d(phi_prime)<0)[0][0]+1
#print 'ind = ', ind0
else:
ind0=0 # this is for the first map
#shifting phi_prime to the left
phi_prime=np.roll(phi_prime,-ind0)
bi_wsa_rolled=np.roll(bi_wsa,-ind0,axis=1)
v_wsa_rolled=np.roll(v_wsa,-ind0,axis=1)
n_wsa_rolled=np.roll(n_wsa,-ind0,axis=1)
T_wsa_rolled=np.roll(T_wsa,-ind0,axis=1)
#TODO: Check Temp calculation
T0 = 0.9e6
Rho0 = 1100.*mp #density in the HCS
#cs = np.sqrt(prm.gamma/rho*(rho.max()*1.38e-16*T0/1.67e-24-br**2/8/np.pi))
Temp = mp/rho/kblts*(Rho0*kblts*T0/mp-br**2/8./np.pi/2.)
#plot br from original wsa map (top plot) and shifted to the origin map(bottom plot)
#changes in time in the bottom plot are purely due to time-dependent variations of B_r (rotation is eliminted)
plot(wsaFile, bi_wsa, bi_wsa_rolled)
##plot(wsaFile, v_wsa, v_wsa_rolled)
# Poisson solver after interpolation onto GAMERA grid
if fcount>0:
print ('fcount = ', fcount)
#right-hand side of laplacian equation
pois.setRHS( (br-br_save).T) #after transponding it becomes (nj,nk)
guess=np.zeros_like(br.T)
#electric field potential psi: (Laplacian(Psi) = dB_r/dt)
Psi = newton_krylov(pois.residual,guess, method='lgmres',verbose=True,iter=100)#,f_rtol=1.e-6) #iter=100
print('Residual: %g' % abs(pois.residual(Psi)).max())
###INTERPOLATION OF ROLLED WSA MAPS TO GAMERA GRID phi-theta####
print ('Psi.shape = ', Psi.shape) # (nj, nk) =(128, 256)
#Psi is defined in cell centers
# bivariate spline approximation over a rectangular mesh
fbi = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,bi_wsa_rolled.T,kx=1,ky=1)
# interpolation to Gamera grid
br = fbi(phi,theta)
#calculate electric field componenet
#suffix _a denotes that this is adapt field
#E_theta = dPsi/dphi/sin(theta)
#E_phi = -dPsi/dtheta/r
#Next Slava used SMOOTHING for br for the paper, we do not need it for now
et_a = np.zeros( (Psi.shape[0],Psi.shape[1]+1) ) #(128, 257)
et_a[:,1:-1] = np.diff(Psi,axis=1)/np.diff(phi) #except first and last cells in k
et_a[:,0] = (Psi[:,0] - Psi[:,-1])/(phi[0]-phi[-1]+2*np.pi) #k=0
et_a[:,-1]=et_a[:,0] #k=Nk+1
et_a /= np.sin(theta[:,None])
print ('E_theta.shape = ', et_a.shape)
"""
note, here we assume theta constant along phi and same theta on the
boundary and in the center of the GAMERA cell - Elena: we should probably fix that
"""
ep_a=np.zeros((Psi.shape[0]+1,Psi.shape[1])) #(129, 256)
ep_a[1:-1,:] = -np.diff(Psi,axis=0)/np.diff(theta)[:,None] #except first and last cells in j
#for j=0 and j=N_j we set nearby values
ep_a[0,:]=ep_a[1,:] # used to set these to zero, but more appropriate to repeat from next theta, since Ephi does not depend on theta at the pole.
ep_a[-1,:]=ep_a[-2,:]
print ('E_phi.shape = ', ep_a.shape)
#[EP]: two lines above are from LFM where theta went from pole to pole
fv = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,v_wsa_rolled.T,kx=1,ky=1)
vr = fv(phi,theta)
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa_rolled.T,kx=1,ky=1)
rho = f(phi,theta)
# I do not understand all that business with interpolation of electric fields in time (see adapt2lfm.py)
# Convert to CGS. FIX ME!!! UNITS HARD CODED
et_a*= prm.Rin*prm.scale/prm.adaptCadence/24./3600.
ep_a*= prm.Rin*prm.scale/prm.adaptCadence/24./3600.
#not interpolating temperature, calculating sound speed cs
#assuming uniform total pressure Rho_max*k*T0 = p+Br^2/8pi
et_save = et_a
ep_save = ep_a
#[EP] for debug
dbr = br - br_save
#et_save and ep_save are defined at times of adapt and on cell edges
#TODO: Check Temp calculation
T0 = 0.9e6
Rho0 = 1100.*mp #density in the HCS
#cs = np.sqrt(prm.gamma/rho*(rho.max()*1.38e-16*T0/1.67e-24-br**2/8/np.pi))
Temp = mp/rho/kblts*(Rho0*kblts*T0/mp-br**2/8./np.pi/2.)
br_save = br
# Poisson solver after interpolation onto GAMERA grid
if fcount>0:
print ('fcount = ', fcount)
#right-hand side of laplacian equation
pois.setRHS( (br-br_save).T) #after transponding it becomes (nj,nk)
guess=np.zeros_like(br.T)
#electric field potential psi: (Laplacian(Psi) = dB_r/dt)
Psi = newton_krylov(pois.residual,guess, method='lgmres',verbose=True,iter=100)#,f_rtol=1.e-6) #iter=100
print('Residual: %g' % abs(pois.residual(Psi)).max())
print ('Psi.shape = ', Psi.shape) # (nj, nk) =(128, 256)
#Psi is defined in cell centers
#calculate electric field componenet
#suffix _a denotes that this is adapt field
#E_theta = dPsi/dphi/sin(theta)
#E_phi = -dPsi/dtheta/r
et_a = np.zeros( (Psi.shape[0],Psi.shape[1]+1) ) #(128, 257)
et_a[:,1:-1] = np.diff(Psi,axis=1)/np.diff(phi) #except first and last cells in k
et_a[:,0] = (Psi[:,0] - Psi[:,-1])/(phi[0]-phi[-1]+2*np.pi) #k=0
et_a[:,-1]=et_a[:,0] #k=Nk+1
et_a /= np.sin(theta[:,None])
print ('E_theta.shape = ', et_a.shape)
"""
After we obtained et_save ep_save at cell edges we calculate B_theta and B_phi
at cell centers and faces
note, here we assume theta constant along phi and same theta on the
boundary and in the center of the GAMERA cell - Elena: we should probably fix that
"""
vrt = vr.T #(nj,nk) in cell centers
bp_a = np.zeros_like(vrt) #B_phi in cell centers
bt_a = np.zeros_like(vrt) #B_theta in cell centers
bp_kface_a = np.zeros( (vrt.shape[0],vrt.shape[1]+1) ) #(nj,nk+1)
bt_jface_a = np.zeros( (vrt.shape[0]+1,vrt.shape[1]) ) #(nj+1,nk)
vrt_kface = np.zeros( (vrt.shape[0],vrt.shape[1]+1) )
vrt_jface = np.zeros( (vrt.shape[0]+1,vrt.shape[1]) )
if fcount >0:
# B_phi and B_theta defined at cell centers
bp_a = 0.5*(et_save[:,:-1]+et_save[:,1:])/vrt
bt_a = -0.5*(ep_save[:-1,:]+ep_save[1:,:])/vrt
# the above are at cell centers, also need at the
# corresponding faces, see below
# First interpolate velocity to faces
vrt_kface[:,1:-1] = 0.5*(vrt[:,:-1]+vrt[:,1:]); vrt_kface[:,0] = 0.5*(vrt[:,-1]+vrt[:,0]); vrt_kface[:,-1] = vrt_kface[:,0]
vrt_jface[1:-1,:] = 0.5*(vrt[1:,:]+vrt[:-1,:]) ; vrt_jface[0,:]=vrt[1,:].mean(); vrt_jface[-1,:]=vrt[-2,:].mean();
#B_phi and B_theta at faces
bp_kface_a = et_save/vrt_kface
bt_jface_a = -ep_save/vrt_jface
#transponse again to agree with GAMERA indexing nk,nj,ni
# Note, these are defined at cell centers on the boundary (at rmin)
bp_a = bp_a.T
bt_a = bt_a.T
#in kaiju we do not to save B-components at cell centers, so we do not need bp_a and bt_a
# at faces; change shapes to match order in gamera nk, nj
bt_jface_a = bt_jface_a.T
bp_kface_a = bp_kface_a.T
et_save = et_save.T
ep_save = ep_save.T
ep_a=np.zeros((Psi.shape[0]+1,Psi.shape[1])) #(129, 256)
ep_a[1:-1,:] = -np.diff(Psi,axis=0)/np.diff(theta)[:,None] #except first and last cells in j
#for j=0 and j=N_j we set nearby values
ep_a[0,:]=ep_a[1,:] # used to set these to zero, but more appropriate to repeat from next theta, since Ephi does not depend on theta at the pole.
ep_a[-1,:]=ep_a[-2,:]
print ('E_phi.shape = ', ep_a.shape)
#[EP]: two lines above are from LFM where theta went from pole to pole
# Scale inside ghost region
#print(rho.shape)
(vr,rho,Temp,br,bp_kface_a,bt_jface_a,et_save,ep_save) = [np.dstack(prm.NO2*[var]) for var in (vr,rho,Temp,br,bp_kface_a,bt_jface_a,et_save,ep_save)]
rho*=(R0/Rc[0,0,:Ng])**2
Temp*=(R0/Rc[0,0,:Ng])
br*=(R0/Rc[0,0,:Ng])**2
bp_kface_a*=(R0/Rc[0,0,:Ng])
et_save*=(R0/Rc[0,0,:Ng])
# I do not understand all that business with interpolation of electric fields in time (see adapt2lfm.py)
# Convert to CGS. FIX ME!!! UNITS HARD CODED
et_a*= prm.Rin*prm.scale/prm.adaptCadence/24./3600.
ep_a*= prm.Rin*prm.scale/prm.adaptCadence/24./3600.
#tangential velocities are set to zero
#vp = zeros_like(vr)
#vt = zeros_like(vr)
et_save = et_a
ep_save = ep_a
#[EP] for debug
dbr = br - br_save
#et_save and ep_save are defined at times of adapt and on cell edges
#print vr.shape, rho.shape, cs.shape, br.shape, bt_jface_a.shape, bp_kface_a.shape
#print et_save.shape, ep_save.shape
br_save = br
#Agreement. For innerbcTD.h5 the output units are V[km/s], Rho[cm-3], T[K], B[nT]
#v_wsa /= Vnorm
#n_wsa /= mp
#bi_wsa /= Bnorm
"""
After we obtained et_save ep_save at cell edges we calculate B_theta and B_phi
at cell centers and faces
"""
vrt = vr.T #(nj,nk) in cell centers
bp_a = np.zeros_like(vrt) #B_phi in cell centers
bt_a = np.zeros_like(vrt) #B_theta in cell centers
bp_kface_a = np.zeros( (vrt.shape[0],vrt.shape[1]+1) ) #(nj,nk+1)
bt_jface_a = np.zeros( (vrt.shape[0]+1,vrt.shape[1]) ) #(nj+1,nk)
vrt_kface = np.zeros( (vrt.shape[0],vrt.shape[1]+1) )
vrt_jface = np.zeros( (vrt.shape[0]+1,vrt.shape[1]) )
print (wsaFile)
#removing two bounding cells in theta and normalizing
vrp = vr[:,1:-1,:]/Vnorm
vp = np.zeros_like(vrp)/Vnorm
vt = np.zeros_like(vrp)/Vnorm
rhop = rho[:,1:-1,:]/mp
Tempp = Temp[:,1:-1,:]
brp = br[:,1:-1,:]/Bnorm
bt_jface_a_p = bt_jface_a[:,1:-1,:]/Bnorm
if fcount >0:
# B_phi and B_theta defined at cell centers
bp_a = 0.5*(et_save[:,:-1]+et_save[:,1:])/vrt
bt_a = -0.5*(ep_save[:-1,:]+ep_save[1:,:])/vrt
bp_kface_a_p = bp_kface_a[:,1:-1,:]/Bnorm
et_save_p = et_save[:,1:-1,:]
ep_save_p = ep_save[:,1:-1,:]
#print vrp.shape, rhop.shape, csp.shape, brp.shape, bt_jface_a_p.shape, bp_kface_a_p.shape
#print et_save_p.shape, ep_save_p.shape
# the above are at cell centers, also need at the
# corresponding faces, see below
#V in cm/s B in Gs n in gcm-3
print (fcount, time_sec, mjd_c)
if prm.dumpBC:
if fcount == 0:
#write out phi and th coords of corners at inner boundary grid
hf.create_dataset("X", data=P_out)
hf.create_dataset("Y", data=T_out)
hf.create_dataset("Z", data=R_out)
grname = "Step#"+str(fcount)
grp = hf.create_group(grname)
grp.attrs.create("time", time_sec)
grp.attrs.create("MJD", mjd_c)
grp.create_dataset("vr",data=vrp) #cc
grp.create_dataset("vp",data=vp) #cc !zeros
grp.create_dataset("vt",data=vt) #cc !zeros
#hf.create_dataset("vr_kface",data=vr_kface) #kface
grp.create_dataset("rho",data=rhop) #cc
grp.create_dataset("T",data=Tempp) #cc
grp.create_dataset("br",data=brp) #cc
#hf.create_dataset("br_kface",data=br_kface) #kface
#hf.create_dataset("bp",data=bp_a) #cc
#hf.create_dataset("bt",data=bt_a) #cc
grp.create_dataset("bt_jface",data=bt_jface_a_p) #jface
grp.create_dataset("bp_kface",data=bp_kface_a_p) #kface
grp.create_dataset("et",data=et_save_p) #k-edges
grp.create_dataset("ep",data=ep_save_p) #j-edges
plotBc(wsaFile,phi, theta[1:-1], vrp[:,:,Ng-1], brp[:,:,Ng-1], rhop[:,:,Ng-1], Tempp[:,:,Ng-1])
# First interpolate velocity to faces
vrt_kface[:,1:-1] = 0.5*(vrt[:,:-1]+vrt[:,1:]); vrt_kface[:,0] = 0.5*(vrt[:,-1]+vrt[:,0]); vrt_kface[:,-1] = vrt_kface[:,0]
vrt_jface[1:-1,:] = 0.5*(vrt[1:,:]+vrt[:-1,:]) ; vrt_jface[0,:]=vrt[1,:].mean(); vrt_jface[-1,:]=vrt[-2,:].mean();
#B_phi and B_theta at faces
bp_kface_a = et_save/vrt_kface
bt_jface_a = -ep_save/vrt_jface
# [EP] test if calculated tengential electric fields give Br from wsa
if fcount > 30:
print ('Elena debug')
print (fcount)
#transponse again to agree with GAMERA indexing nk,nj,ni
# Note, these are defined at cell centers on the boundary (at rmin)
bp_a = bp_a.T
bt_a = bt_a.T
#in kaiju we do not to save B-components at cell centers, so we do not need bp_a and bt_a
dphi = phi[2]-phi[1]
dtheta = theta[2]-theta[1]
# at faces; change shapes to match order in gamera nk, nj
bt_jface_a = bt_jface_a.T
bp_kface_a = bp_kface_a.T
dbrp = dbr[:,1:-1]
et_save = et_save.T
ep_save = ep_save.T
#cell edges dphi and dtheta at inner boundary face
dlp = dphi*rxy[Ng:-Ng-1,Ng:-Ng,Ng] #(256,129) dlp change with nj
#dlp = dphi*r[Ng:-Ng-1,Ng:-Ng,Ng]*sin(T[:,:])
dlt = dtheta*R0 #dlt is same for all cells
et_use = et_save_p.T #(257,128)
ep_use = ep_save_p.T #(256,129)
# Scale inside ghost region
#print(rho.shape)
(vr,rho,Temp,br,bp_kface_a,bt_jface_a,et_save,ep_save) = [np.dstack(prm.NO2*[var]) for var in (vr,rho,Temp,br,bp_kface_a,bt_jface_a,et_save,ep_save)]
rho*=(R0/Rc[0,0,:Ng])**2
Temp*=(R0/Rc[0,0,:Ng])
br*=(R0/Rc[0,0,:Ng])**2
bp_kface_a*=(R0/Rc[0,0,:Ng])
et_save*=(R0/Rc[0,0,:Ng])
#rotE of the cell face
circE = zeros((nk,nj))
#tangential velocities are set to zero
#vp = zeros_like(vr)
#vt = zeros_like(vr)
#circE1 = - (ep_use[:,:-1]*dlp + et_use[1:,:]*dlt - ep_use[:,1:]*dlp - et_use[:-1,:]*dlt)
#print vr.shape, rho.shape, cs.shape, br.shape, bt_jface_a.shape, bp_kface_a.shape
#print et_save.shape, ep_save.shape
for k in range(256):
for j in range(128):
circE[k,j] = - (ep_use[k,j+1]*dlp[k,j+1] - ep_use[k,j]*dlp[k,j] + et_use[k,j]*dlt - et_use[k+1,j]*dlt)
#Agreement. For innerbcTD.h5 the output units are V[km/s], Rho[cm-3], T[K], B[nT]
#v_wsa /= Vnorm
#n_wsa /= mp
#bi_wsa /= Bnorm
print (wsaFile)
#removing two bounding cells in theta and normalizing
vrp = vr[:,1:-1,:]/Vnorm
vp = np.zeros_like(vrp)/Vnorm
vt = np.zeros_like(vrp)/Vnorm
rhop = rho[:,1:-1,:]/mp
Tempp = Temp[:,1:-1,:]
brp = br[:,1:-1,:]/Bnorm
bt_jface_a_p = bt_jface_a[:,1:-1,:]/Bnorm
bp_kface_a_p = bp_kface_a[:,1:-1,:]/Bnorm
et_save_p = et_save[:,1:-1,:]
ep_save_p = ep_save[:,1:-1,:]
#print vrp.shape, rhop.shape, csp.shape, brp.shape, bt_jface_a_p.shape, bp_kface_a_p.shape
#print et_save_p.shape, ep_save_p.shape
#V in cm/s B in Gs n in gcm-3
print (fcount, time_sec, mjd_c)
if prm.dumpBC:
if fcount == 0:
#write out phi and th coords of corners at inner boundary grid
hf.create_dataset("X", data=P_out)
hf.create_dataset("Y", data=T_out)
hf.create_dataset("Z", data=R_out)
grname = "Step#"+str(fcount)
grp = hf.create_group(grname)
grp.attrs.create("time", time_sec)
grp.attrs.create("MJD", mjd_c)
grp.create_dataset("vr",data=vrp) #cc
grp.create_dataset("vp",data=vp) #cc !zeros
grp.create_dataset("vt",data=vt) #cc !zeros
#hf.create_dataset("vr_kface",data=vr_kface) #kface
grp.create_dataset("rho",data=rhop) #cc
grp.create_dataset("T",data=Tempp) #cc
grp.create_dataset("br",data=brp) #cc
#hf.create_dataset("br_kface",data=br_kface) #kface
#hf.create_dataset("bp",data=bp_a) #cc
#hf.create_dataset("bt",data=bt_a) #cc
grp.create_dataset("bt_jface",data=bt_jface_a_p) #jface
grp.create_dataset("bp_kface",data=bp_kface_a_p) #kface
grp.create_dataset("et",data=et_save_p) #k-edges
grp.create_dataset("ep",data=ep_save_p) #j-edges
plotBc(wsaFile,phi, theta[1:-1], vrp[:,:,Ng-1], brp[:,:,Ng-1], rhop[:,:,Ng-1], Tempp[:,:,Ng-1])
# [EP] test if calculated tengential electric fields give Br from wsa
if fcount > 30:
print ('Elena debug')
print (fcount)
dphi = phi[2]-phi[1]
dtheta = theta[2]-theta[1]
dbrp = dbr[:,1:-1]
dt = prm.adaptCadence*24.*3600.
dbrdt = dlp[:,:-1]*dlt*prm.scale*dbrp/dt
resid = dbrdt - circE
fig1 = plt.figure(); plt.pcolormesh(np.log10(np.abs(resid.T))); plt.colorbar()
fig1.suptitle(wsaFile)
plt.savefig(wsaFile[:-5]+'_testFL.png')
#cell edges dphi and dtheta at inner boundary face
dlp = dphi*rxy[Ng:-Ng-1,Ng:-Ng,Ng] #(256,129) dlp change with nj
#dlp = dphi*r[Ng:-Ng-1,Ng:-Ng,Ng]*sin(T[:,:])
dlt = dtheta*R0 #dlt is same for all cells
et_use = et_save_p.T #(257,128)
ep_use = ep_save_p.T #(256,129)
#rotE of the cell face
circE = zeros((nk,nj))
#circE1 = - (ep_use[:,:-1]*dlp + et_use[1:,:]*dlt - ep_use[:,1:]*dlp - et_use[:-1,:]*dlt)
for k in range(256):
for j in range(128):
circE[k,j] = - (ep_use[k,j+1]*dlp[k,j+1] - ep_use[k,j]*dlp[k,j] + et_use[k,j]*dlt - et_use[k+1,j]*dlt)
dt = prm.adaptCadence*24.*3600.
dbrdt = dlp[:,:-1]*dlt*prm.scale*dbrp/dt
resid = dbrdt - circE
fig1 = plt.figure(); plt.pcolormesh(np.log10(np.abs(resid.T))); plt.colorbar()
fig1.suptitle(wsaFile)
plt.savefig(wsaFile[:-5]+'_testFL.png')
#if fcount==2:
# sys.exit("passed two wsa files")
#if fcount==2:
# sys.exit("passed two wsa files")
if __name__ == '__main__':
main()

0
kaipy/scripts/quicklook/__init__.py
View File

8
kaipy/scripts/quicklook/dbVid.py
View File

@@ -55,7 +55,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
rad2deg = 180.0/np.pi
bMag = dbViz.dbMag
@@ -163,8 +163,4 @@ def main():
#Save
npl = vO[i]-nS
fOut = oDir+"/vid.%04d.png"%(npl)
kv.savePic(fOut,bLenX=45)
if __name__ == "__main__":
main()
kv.savePic(fOut,bLenX=45)

5
kaipy/scripts/quicklook/dbpic.py
View File

@@ -103,7 +103,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
"""Plot the ground magnetic field perturbations."""
# Set up the command-line parser.
@@ -277,6 +277,3 @@ def main():
if debug:
print("fOut = %s" % fOut)
kv.savePic(fOut)
if __name__ == "__main__":
main()

5
kaipy/scripts/quicklook/dstpic.py
View File

@@ -45,7 +45,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
iMax = -1
@@ -128,6 +128,3 @@ def main():
ax.set_xlim(xMin,xMax)
kv.savePic("qkdstpic.png")
if __name__ == "__main__":
main()

5
kaipy/scripts/quicklook/gamerrVid.py
View File

@@ -129,7 +129,7 @@ def makeImage(i,gsph1,gsph2,tOut,doVerb,xyBds,fnList,oDir,errTimes,errListRel,er
fOut = oDir+"/vid.%04d.png"%(npl)
kv.savePic(fOut,bLenX=45,saveFigure=fig,doClose=True)
def main():
if __name__ == "__main__":
#Defaults
fdir1 = os.getcwd()
ftag1 = "msphere"
@@ -254,6 +254,3 @@ def main():
makeMovie(oDir,oSub)
if __name__ == "__main__":
main()

5
kaipy/scripts/quicklook/gamerrpic.py
View File

@@ -52,7 +52,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
doMPI = False #[Add MPI tiling]
# Set up the command-line parser.
@@ -109,6 +109,3 @@ def main():
kv.savePic(oName,bLenX=45)
if __name__ == "__main__":
main()

6
kaipy/scripts/quicklook/gamsphVid.py
View File

@@ -67,7 +67,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
doDen = False
doMPI = False #[Add MPI tiling]
@@ -205,7 +205,3 @@ def main():
fOut = oDir+"/vid.%04d.png"%(npl)
kv.savePic(fOut,bLenX=45)
if __name__ == "__main__":
main()

5
kaipy/scripts/quicklook/mixpic.py
View File

@@ -277,7 +277,7 @@ def makePlot(i, remixFile, nStp):
#plt.close(fig)
def main():
if __name__ == "__main__":
"""Plot remix data, either a single time step or as a movie"""
# Set up the command-line parser.
@@ -380,6 +380,3 @@ def main():
with Pool(processes=ncpus) as pl:
pl.starmap(makePlot,ag)
print("Done making all the images. Go to mixVid folder")
if __name__ == "__main__":
main()

5
kaipy/scripts/quicklook/msphpic.py
View File

@@ -321,7 +321,7 @@ def makePlot(i,spacecraft,nStp):
kv.savePic(outPath, bLenX=45)
def main():
if __name__ == "__main__":
"""Make a quick figure of a Gamera magnetosphere run."""
# Set up the command-line parser.
@@ -440,6 +440,3 @@ def main():
with Pool(processes=ncpus) as pl:
pl.starmap(makePlot,ag)
print("Done making all the images. Go to mixVid folder")
if __name__ == "__main__":
main()

6
kaipy/scripts/quicklook/rcmDataProbe.py
View File

@@ -62,7 +62,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
cmap = plt.cm.plasma
# Set up the command-line parser.
@@ -210,7 +210,3 @@ def main():
#Ax.set_ylim(1.0e+10,1.0e+19)
sTag = "e" if doElectrons else "i"
kv.savePic("qkrcm%s_%s.png"%(filetag, sTag))
if __name__ == "__main__":
main()

4
kaipy/scripts/quicklook/rcmpic.py
View File

@@ -474,7 +474,7 @@ def makePlot(i,rcmdata,nStp):
# Save the figure to a file.
kv.savePic(outPath, dpiQ=300)
def main():
if __name__ == "__main__":
"""Make a quick figure of a Gamera magnetosphere run."""
# Set up the command-line parser.
@@ -597,5 +597,3 @@ def main():
pl.starmap(makePlot,ag)
print("Done making all the images. Go to mixVid folder")
if __name__ == "__main__":
main()

5
kaipy/scripts/quicklook/remixTimeSeries.py
View File

@@ -45,7 +45,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
parser = create_command_line_parser()
args = parser.parse_args()
#Open the file and read the time information
@@ -138,7 +138,6 @@ def main():
dset.attrs['name'] = ipfac['name']
if __name__ == "__main__":
main()

4
kaipy/scripts/quicklook/swpic.py
View File

@@ -32,7 +32,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
# Set up the command-line parser.
@@ -103,5 +103,3 @@ def main():
swBCplots.swQuickPlot(UTall,D,Temp,Vx,Vy,Vz,Bx,By,Bz,SYMH,pltInterp,fOut,doEps=doEps,doTrim=doTrim,t0fmt=t0Fmt)
if __name__ == "__main__":
main()

7
kaipy/scripts/quicklook/vizTrj.py
View File

@@ -67,7 +67,7 @@ def create_command_line_parser():
return parser
def main():
if __name__ == "__main__":
#Defaults
LW = 0.5
fs = 16
@@ -107,7 +107,4 @@ def main():
titS = "Particle trajectory, pID = %d"%(id0)
ax.set_title(titS,fontsize=fs)
axEqual3d(ax)
plt.show()
if __name__ == "__main__":
main()
plt.show()

5
kaipy/solarWind/OMNI.py
View File

@@ -762,9 +762,6 @@ class OMNIW(OMNI):
return (dates,dataArray, hasBeenInterpolated, dataOrigin)
def main():
if __name__ == '__main__':
import doctest
doctest.testmod()
if __name__ == '__main__':
main()

5
kaipy/solarWind/SWPC.py
View File

@@ -528,9 +528,6 @@ class ACESWPC(DSCOVRNC):
units='n/a',
data=metadata)
def main():
if __name__ == '__main__':
import doctest
doctest.testmod()
if __name__ == "__main__":
main()

5
kaipy/solarWind/TimeSeries.py
View File

@@ -131,9 +131,6 @@ class TimeSeries(dict):
except TypeError:
self[key]['data'] = data
def main():
if __name__ == '__main__':
import doctest
doctest.testmod()
if __name__ == '__main__':
main()

5
kaipy/solarWind/WIND.py
View File

@@ -407,9 +407,6 @@ class WIND(OMNI):
return (diff.days*24.0*60.0 + diff.seconds/60.0)
def main():
if __name__ == '__main__':
import doctest
doctest.testmod()
if __name__ == '__main__':
main()

52
pyproject.toml
View File

@@ -7,7 +7,7 @@ include-package-data = true
[project]
name = "kaipy"
version = "1.0.6"
version = "1.0.5"
authors = [
{name = "Eric Winter", email = "eric.winter@jhuapl.edu"},
]
@@ -34,53 +34,3 @@ dependencies = [
"spacepy",
"sunpy",
]
[project.scripts]
test_cdas = "kaipy.satcomp.test_cdas:main"
ih2oh = "kaipy.scripts.OHelio.ih2oh:main"
helioSatComp = "kaipy.scripts.datamodel.helioSatComp:main"
msphParallelSatComp = "kaipy.scripts.datamodel.msphParallelSatComp:main"
msphPbsSatComp = "kaipy.scripts.datamodel.msphPbsSatComp:main"
msphSatComp = "kaipy.scripts.datamodel.msphSatComp:main"
rbspSCcomp = "kaipy.scripts.datamodel.rbspSCcomp:main"
block_genmpiXDMF = "kaipy.scripts.postproc.block_genmpiXDMF:main"
embiggen = "kaipy.scripts.postproc.embiggen:main"
embiggenMIX = "kaipy.scripts.postproc.embiggenMIX:main"
embiggenRCM = "kaipy.scripts.postproc.embiggenRCM:main"
embiggenVOLT = "kaipy.scripts.postproc.embiggenVOLT:main"
genXDMF = "kaipy.scripts.postproc.genXDMF:main"
genXLine = "kaipy.scripts.postproc.genXLine:main"
genmpiXDMF = "kaipy.scripts.postproc.genmpiXDMF:main"
numSteps = "kaipy.scripts.postproc.numSteps:main"
pitmerge = "kaipy.scripts.postproc.pitmerge:main"
printResTimes = "kaipy.scripts.postproc.printResTimes:main"
run_supermag_comparison = "kaipy.scripts.postproc.run_supermag_comparison:main"
slimFL = "kaipy.scripts.postproc.slimFL:main"
slimh5 = "kaipy.scripts.postproc.slimh5:main"
slimh5_classic = "kaipy.scripts.postproc.slimh5_classic:main"
ut2mjd = "kaipy.scripts.postproc.ut2mjd:main"
cda2wind = "kaipy.scripts.preproc.cda2wind:main"
genLFM = "kaipy.scripts.preproc.genLFM:main"
genRCM = "kaipy.scripts.preproc.genRCM:main"
wsa2TDgamera = "kaipy.scripts.preproc.wsa2TDgamera:main"
XMLGenerator = "kaipy.scripts.preproc.XMLGenerator:main"
dbVid = "kaipy.scripts.quicklook.dbVid:main"
dbpic = "kaipy.scripts.quicklook.dbpic:main"
dstpic = "kaipy.scripts.quicklook.dstpic:main"
gamerrVid = "kaipy.scripts.quicklook.gamerrVid:main"
gamerrpic = "kaipy.scripts.quicklook.gamerrpic:main"
gamsphVid = "kaipy.scripts.quicklook.gamsphVid:main"
mixpic = "kaipy.scripts.quicklook.mixpic:main"
msphpic = "kaipy.scripts.quicklook.msphpic:main"
rcmDataProbe = "kaipy.scripts.quicklook.rcmDataProbe:main"
rcmPrecipSpecFlux = "kaipy.scripts.quicklook.rcmPrecipSpecFlux:main"
rcmSpecFlux = "kaipy.scripts.quicklook.rcmSpecFlux:main"
rcmpic = "kaipy.scripts.quicklook.rcmpic:main"
remixTimeSeries = "kaipy.scripts.quicklook.remixTimeSeries:main"
swpic = "kaipy.scripts.quicklook.swpic:main"
vizTrj = "kaipy.scripts.quicklook.vizTrj:main"
OMNI = "kaipy.solarWind.OMNI:main"
SWPC = "kaipy.solarWind.SWPC:main"
TimeSeries = "kaipy.solarWind.TimeSeries:main"
WIND = "kaipy.solarWind.WIND:main"
ols = "kaipy.solarWind.ols:main"

54
setup.py
View File

@@ -2,7 +2,7 @@ from setuptools import setup, find_packages
setup(
name='kaipy',
version='1.0.6',
version='1.0.5',
description='Python software for CGS MAGE and other Kaiju models',
author='Kaiju team',
author_email='wiltbemj@ucar.edu',
@@ -37,56 +37,4 @@ setup(
'Programming Language :: Python :: 3.9',
'Programming Language :: Python :: 3.10',
],
entry_points={
'console_scripts': [
'test_cdas=kaipy.satcomp.test_cdas:main',
'ih2oh=kaipy.scripts.OHelio.ih2oh:main',
'helioSatComp=kaipy.scripts.datamodel.helioSatComp:main',
'msphParallelSatComp=kaipy.scripts.datamodel.msphParallelSatComp:main',
'msphPbsSatComp=kaipy.scripts.datamodel.msphPbsSatComp:main',
'msphSatComp=kaipy.scripts.datamodel.msphSatComp:main',
'rbspSCcomp=kaipy.scripts.datamodel.rbspSCcomp:main',
'block_genmpiXDMF=kaipy.scripts.postproc.block_genmpiXDMF:main',
'embiggen=kaipy.scripts.postproc.embiggen:main',
'embiggenMIX=kaipy.scripts.postproc.embiggenMIX:main',
'embiggenRCM=kaipy.scripts.postproc.embiggenRCM:main',
'embiggenVOLT=kaipy.scripts.postproc.embiggenVOLT:main',
'genXDMF=kaipy.scripts.postproc.genXDMF:main',
'genXLine=kaipy.scripts.postproc.genXLine:main',
'genmpiXDMF=kaipy.scripts.postproc.genmpiXDMF:main',
'numSteps=kaipy.scripts.postproc.numSteps:main',
'pitmerge=kaipy.scripts.postproc.pitmerge:main',
'printResTimes=kaipy.scripts.postproc.printResTimes:main',
'run_supermag_comparison=kaipy.scripts.postproc.run_supermag_comparison:main',
'slimFL=kaipy.scripts.postproc.slimFL:main',
'slimh5=kaipy.scripts.postproc.slimh5:main',
'slimh5_classic=kaipy.scripts.postproc.slimh5_classic:main',
'ut2mjd=kaipy.scripts.postproc.ut2mjd:main',
'cda2wind=kaipy.scripts.preproc.cda2wind:main',
'genLFM=kaipy.scripts.preproc.genLFM:main',
'genRCM=kaipy.scripts.preproc.genRCM:main',
'wsa2TDgamera=kaipy.scripts.preproc.wsa2TDgamera:main',
'XMLGenerator=kaipy.scripts.preproc.XMLGenerator:main',
'dbVid=kaipy.scripts.quicklook.dbVid:main',
'dbpic=kaipy.scripts.quicklook.dbpic:main',
'dstpic=kaipy.scripts.quicklook.dstpic:main',
'gamerrVid=kaipy.scripts.quicklook.gamerrVid:main',
'gamerrpic=kaipy.scripts.quicklook.gamerrpic:main',
'gamsphVid=kaipy.scripts.quicklook.gamsphVid:main',
'mixpic=kaipy.scripts.quicklook.mixpic:main',
'msphpic=kaipy.scripts.quicklook.msphpic:main',
'rcmDataProbe=kaipy.scripts.quicklook.rcmDataProbe:main',
'rcmPrecipSpecFlux=kaipy.scripts.quicklook.rcmPrecipSpecFlux:main',
'rcmSpecFlux=kaipy.scripts.quicklook.rcmSpecFlux:main',
'rcmpic=kaipy.scripts.quicklook.rcmpic:main',
'remixTimeSeries=kaipy.scripts.quicklook.remixTimeSeries:main',
'swpic=kaipy.scripts.quicklook.swpic:main',
'vizTrj=kaipy.scripts.quicklook.vizTrj:main',
'OMNI=kaipy.solarWind.OMNI:main',
'SWPC=kaipy.solarWind.SWPC:main',
'TimeSeries=kaipy.solarWind.TimeSeries:main',
'WIND=kaipy.solarWind.WIND:main',
'ols=kaipy.solarWind.ols:main'
]
}
)