AtHeartEngineer/kaiju
1
0
Fork 0
mirror of https://github.com/JHUAPL/kaiju.git synced 2026-01-09 15:17:56 -05:00
Code Issues Packages Projects Releases Wiki Activity

adding helio ibc python scripts to kaiju

Browse Source
This commit is contained in:
Elena Provornikova
2021-01-04 11:30:21 -07:00
parent 5bd335f661
commit fdf9fd9b34
13 changed files with 1035 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
kaipy/gamhelio/ConfigScripts/startup.config Normal file
View File

@@ -0,0 +1,20 @@
[Gamera]
gameraGridFile = heliogrid_1024x512x1024.h5
Dir = /Users/merkivg1/Git/kaiju/gamera/examples/helio
[WSA]
wsafile = /Users/merkivg1/amun/vgm/work/LFM-helio/Fast_Latitude_Scan/ulys/WSA_SUB2.1/Baseline/1892.fits
density_temperature_infile = yes
gauss_smooth_width = 0 ; 8
plots = yes
normalized = yes
[Constants]
gamma = 1.5 ;1.05
NO2 = 4
[Normalization]
B0 = 1.e-3 ; 100 nT
n0 = 200. ; 200/cc

83
kaipy/gamhelio/lib/ace.py Normal file
View File

@@ -0,0 +1,83 @@
import numpy as np
import datetime
def cdaweb_mfi(ACE_mf_file):
# Mag. field
ACE_mf_time=[]
ACE_bx=[]
ACE_by=[]
ACE_bz=[]
i=1
for line in open(ACE_mf_file,'r'):
i+=1
if 'dd-mm-yyyy' in line: break
start = i
i = 0
for line in open(ACE_mf_file,'r'):
i+=1
if i<start: continue
if line[0]=='#': break
fields = line.strip().split()
ACE_mf_day = int(fields[0].strip().split('-')[0])
ACE_mf_mon = int(fields[0].strip().split('-')[1])
ACE_mf_year = int(fields[0].strip().split('-')[2])
ACE_mf_hr = int(fields[1].strip().split(':')[0])
ACE_mf_min = int(fields[1].strip().split(':')[1])
ACE_mf_sec = int(float(fields[1].strip().split(':')[2]))
ACE_mf_time.append(datetime.datetime(ACE_mf_year, ACE_mf_mon, ACE_mf_day, ACE_mf_hr, ACE_mf_min, ACE_mf_sec))
ACE_bx.append(float(fields[2]))
ACE_by.append(float(fields[3]))
ACE_bz.append(float(fields[4]))
# invalidate bad data
ACE_bx=np.array(ACE_bx)
ACE_by=np.array(ACE_by)
ACE_bz=np.array(ACE_bz)
ACE_bx[np.abs(np.array(ACE_bx)) > 1.e20] = np.nan
ACE_by[np.abs(np.array(ACE_by)) > 1.e20] = np.nan
ACE_bz[np.abs(np.array(ACE_bz)) > 1.e20] = np.nan
return (ACE_mf_time,ACE_bx,ACE_by,ACE_bz)
def cdaweb_swe(ACE_plasma_file):
# Plasma
ACE_plasma_time=[]
ACE_v=[]
ACE_n=[]
i=1
for line in open(ACE_plasma_file,'r'):
i+=1
if 'dd-mm-yyyy' in line: break
start = i
i = 0
for line in open(ACE_plasma_file,'r'):
i+=1
if i<start: continue
if line[0]=='#': break
fields = line.strip().split()
ACE_plasma_day = int(fields[0].strip().split('-')[0])
ACE_plasma_mon = int(fields[0].strip().split('-')[1])
ACE_plasma_year = int(fields[0].strip().split('-')[2])
ACE_plasma_hr = int(fields[1].strip().split(':')[0])
ACE_plasma_min = int(fields[1].strip().split(':')[1])
ACE_plasma_sec = int(float(fields[1].strip().split(':')[2]))
ACE_plasma_time.append(datetime.datetime(ACE_plasma_year, ACE_plasma_mon, ACE_plasma_day, ACE_plasma_hr, ACE_plasma_min, ACE_plasma_sec))
ACE_n.append(float(fields[2]))
ACE_v.append(float(fields[3]))
# invalidate bad data
ACE_v=np.array(ACE_v)
ACE_n=np.array(ACE_n)
ACE_v[np.abs(np.array(ACE_v)) > 1.e20] = np.nan
ACE_n[np.abs(np.array(ACE_n)) > 1.e20] = np.nan
return (ACE_plasma_time,ACE_v,ACE_n)

29
kaipy/gamhelio/lib/cspice.py Normal file
View File

@@ -0,0 +1,29 @@
import numpy as np
import datetime
def cspice(fname):
dati=[]
x=[]
y=[]
z=[]
r=[]
phi=[]
theta=[]
for line in open(fname,'r'):
fields = line.strip().split()
year = int(fields[0].strip().split('-')[0])
doy = int(fields[0].strip().split('-')[1])
hr = int(fields[1].strip().split(':')[0])
mn = int(fields[1].strip().split(':')[1])
dati.append(datetime.datetime(year,1,1,hr,mn)+datetime.timedelta(days=doy-1))
x.append(float(fields[2]))
y.append(float(fields[3]))
z.append(float(fields[4]))
r.append(float(fields[5]))
theta.append(float(fields[6]))
phi.append(float(fields[7]))
return(dati,x,y,z,r,theta,phi)

293
kaipy/gamhelio/lib/lfmhlib.py Normal file
View File

@@ -0,0 +1,293 @@
from pyhdf.SD import SD, SDC
from numpy import sqrt,arccos,arctan2,cos,sin,pi,roll
import time
def read(filename):
"""
Returns R,theta,phi,Rc,thetac,phic,br,btheta,bphi,vr,vtheta,vphi,rho,cs
"""
hdffile = SD(filename,SDC.READ)
x=hdffile.select('X_grid').get()
y=hdffile.select('Y_grid').get()
z=hdffile.select('Z_grid').get()
bx=hdffile.select('bx_').get()[:-1,:-1,:-1]
by=hdffile.select('by_').get()[:-1,:-1,:-1]
bz=hdffile.select('bz_').get()[:-1,:-1,:-1]
vx=hdffile.select('vx_').get()[:-1,:-1,:-1]
vy=hdffile.select('vy_').get()[:-1,:-1,:-1]
vz=hdffile.select('vz_').get()[:-1,:-1,:-1]
rho=hdffile.select('rho_').get()[:-1,:-1,:-1]
cs=hdffile.select('c_').get()[:-1,:-1,:-1]
t=hdffile.time
hdffile.end()
# =========== Cell centers ==============
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:] )
# =======================================
R=sqrt(x**2+y**2+z**2)
theta=arccos(z/R)
phi=arctan2(y,x)
phi[phi<0]+=2*pi
Rc=sqrt(xc**2+yc**2+zc**2)
thetac=arccos(zc/Rc)
phic=arctan2(yc,xc)
phic[phic<0]+=2*pi
br = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
btheta = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
bphi =-bx*sin(phic) + by*cos(phic)
vr = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
vtheta = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
vphi =-vx*sin(phic) + vy*cos(phic)
return(t,R,theta,phi,Rc,thetac,phic,br,btheta,bphi,vr,vtheta,vphi,rho,cs)
def r_theta_phi_uniform(filename):
"""
Return R, theta, phi 1-d arrays, assuming uniform grid
"""
hdffile = SD(filename,SDC.READ)
# note, minimal data are read from file
ni=hdffile.select('X_grid').ni-1
nj=hdffile.select('X_grid').nj-1
nk=hdffile.select('X_grid').nk-1
# first get R; in principle, could just take z along the axis
# but are allowing for the possibility that the axis is cut out
x=hdffile.select('X_grid').get(start=(0,0,0),count=(1,1,ni+1)).squeeze()
y=hdffile.select('Y_grid').get(start=(0,0,0),count=(1,1,ni+1)).squeeze()
z=hdffile.select('Z_grid').get(start=(0,0,0),count=(1,1,ni+1)).squeeze()
R = sqrt(x**2+y**2+z**2)
z=hdffile.select('Z_grid').get(start=(0,0,0),count=(1,nj+1,1)).squeeze()
theta = arccos(z/R[0])
x=hdffile.select('X_grid').get(start=(0,nj/2,0),count=(nk+1,1,1)).squeeze()
y=hdffile.select('Y_grid').get(start=(0,nj/2,0),count=(nk+1,1,1)).squeeze()
phi=arctan2(y,x)
phi[phi<0]+=2*pi
phi[-1]=phi[0]+2*pi
hdffile.end()
R/=6.96e10 # FIX ME HARD CODED UNITS
Rc = 0.5*(R[1:]+R[:-1])
thetac = 0.5*(theta[1:]+theta[:-1])
phic = 0.5*(phi[1:]+phi[:-1])
return R,theta,phi,Rc,thetac,phic
def read_var(filename,var_name):
hdffile = SD(filename,SDC.READ)
if var_name not in ['br','btheta','bphi','vr','vtheta','vphi']:
var=hdffile.select(var_name+'_').get()[:-1,:-1,:-1]
else:
R,theta,phi=r_theta_phi_uniform(filename)
thetac = 0.5*(theta[1:]+theta[:-1])
phic = 0.5*(phi[1:]+phi[:-1])
phic = phic[:,None,None]
thetac = thetac[None,:,None]
if var_name in ['br','btheta','bphi']:
bx=hdffile.select('bx_').get()[:-1,:-1,:-1]
by=hdffile.select('by_').get()[:-1,:-1,:-1]
bz=hdffile.select('bz_').get()[:-1,:-1,:-1]
if var_name=='br':
var = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
elif var_name=='btheta':
var = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
else:
var =-bx*sin(phic) + by*cos(phic)
else:
vx=hdffile.select('vx_').get()[:-1,:-1,:-1]
vy=hdffile.select('vy_').get()[:-1,:-1,:-1]
vz=hdffile.select('vz_').get()[:-1,:-1,:-1]
if var_name=='vr':
var = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
elif var_name=='vtheta':
var = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
else:
var =-vx*sin(phic) + vy*cos(phic)
hdffile.end()
return(var)
def read_var_point(filename,var_name,i,j,k,thetac,phic):
thetac = thetac[j]
phic = phic[k]
hdffile = SD(filename,SDC.READ)
if var_name not in ['br','btheta','bphi','vr','vtheta','vphi']:
var=hdffile.select(var_name+'_').get(start=(k,j,i),count=(1,1,1)).squeeze()
else:
# R,theta,phi=r_theta_phi_uniform(filename)
if var_name in ['br','btheta','bphi']:
bx=hdffile.select('bx_').get(start=(k,j,i),count=(1,1,1)).squeeze()
by=hdffile.select('by_').get(start=(k,j,i),count=(1,1,1)).squeeze()
bz=hdffile.select('bz_').get(start=(k,j,i),count=(1,1,1)).squeeze()
if var_name=='br':
var = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
elif var_name=='btheta':
var = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
else:
var =-bx*sin(phic) + by*cos(phic)
else:
vx=hdffile.select('vx_').get(start=(k,j,i),count=(1,1,1)).squeeze()
vy=hdffile.select('vy_').get(start=(k,j,i),count=(1,1,1)).squeeze()
vz=hdffile.select('vz_').get(start=(k,j,i),count=(1,1,1)).squeeze()
if var_name=='vr':
var = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
elif var_name=='vtheta':
var = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
else:
var =-vx*sin(phic) + vy*cos(phic)
hdffile.end()
return(var)
def read_var_islice(filename,var_name,i,thetac,phic):
nk = phic.size
nj = thetac.size
phic = phic[:,None]
thetac = thetac[None,:]
hdffile = SD(filename,SDC.READ)
if var_name not in ['br','btheta','bphi','vr','vtheta','vphi']:
var=hdffile.select(var_name+'_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
else:
if var_name in ['br','btheta','bphi']:
bx=hdffile.select('bx_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
by=hdffile.select('by_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
bz=hdffile.select('bz_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
if var_name=='br':
var = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
elif var_name=='btheta':
var = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
else:
var =-bx*sin(phic) + by*cos(phic)
else:
vx=hdffile.select('vx_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
vy=hdffile.select('vy_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
vz=hdffile.select('vz_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
if var_name=='vr':
var = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
elif var_name=='vtheta':
var = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
else:
var =-vx*sin(phic) + vy*cos(phic)
hdffile.end()
return(var)
def read_var_jslice(filename,var_name,j,thetac,phic):
hdffile = SD(filename,SDC.READ)
ni = hdffile.select('X_grid').ni-1
nk = hdffile.select('X_grid').nk-1
phic = phic[:,None]
thetac = thetac[j]
if var_name not in ['br','btheta','bphi','vr','vtheta','vphi']:
if var_name in ['X_grid','Y_grid','Z_grid']:
var=hdffile.select(var_name).get(start=(0,j,0),count=(nk,1,ni)).squeeze()
else:
var=hdffile.select(var_name+'_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
else:
if var_name in ['br','btheta','bphi']:
bx=hdffile.select('bx_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
by=hdffile.select('by_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
bz=hdffile.select('bz_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
if var_name=='br':
var = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
elif var_name=='btheta':
var = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
else:
var =-bx*sin(phic) + by*cos(phic)
else:
vx=hdffile.select('vx_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
vy=hdffile.select('vy_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
vz=hdffile.select('vz_').get(start=(0,j,0),count=(nk,1,ni)).squeeze()
if var_name=='vr':
var = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
elif var_name=='vtheta':
var = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
else:
var =-vx*sin(phic) + vy*cos(phic)
hdffile.end()
return(var)
def read_var_ikslice(filename,var_name,i,k,tc,pc):
hdffile = SD(filename,SDC.READ)
ni = hdffile.select('X_grid').ni-1
nj = hdffile.select('X_grid').nj-1
if var_name not in ['br','btheta','bphi','vr','vtheta','vphi','bt','bp','vt','vp']:
var=hdffile.select(var_name+'_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
else:
if var_name in ['br','btheta','bphi','bt','bp']:
bx=hdffile.select('bx_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
by=hdffile.select('by_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
bz=hdffile.select('bz_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
if var_name=='br':
var = bx*cos(pc[k])*sin(tc) + by*sin(pc[k])*sin(tc) + bz*cos(tc)
elif (var_name=='btheta' or var_name=='bt'):
var = bx*cos(pc[k])*cos(tc) + by*sin(pc[k])*cos(tc) - bz*sin(tc)
else:
var =-bx*sin(pc[k]) + by*cos(pc[k])
else:
vx=hdffile.select('vx_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
vy=hdffile.select('vy_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
vz=hdffile.select('vz_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
if var_name=='vr':
var = vx*cos(pc[k])*sin(tc) + vy*sin(pc[k])*sin(tc) + vz*cos(tc)
elif (var_name=='vtheta' or var_name=='vt'):
var = vx*cos(pc[k])*cos(tc) + vy*sin(pc[k])*cos(tc) - vz*sin(tc)
else:
var =-vx*sin(pc[k]) + vy*cos(pc[k])
hdffile.end()
return(var)
def get_time(filename):
hdffile = SD(filename, SDC.READ)
t=hdffile.time
hdffile.end()
return(t)
def time_shift(t,var,Tsolar=25.38*24*3600.):
"""
Shift variables to the beginning of the CR.
Parameters:
t: current time
var: variable -- assumed to be a 2D i-slice (nk,nj) of the LFM data
Tsolar: solar rotation period in units of time used in the code. Default to 25.38 days and seconds as time units.
"""
time_shift= (t % Tsolar)/Tsolar
nshift = int(time_shift*(var.shape[0]+1)) # note, using nk+1 here, which gives us full rotation
return roll(var,-nshift,axis=0)

138
kaipy/gamhelio/lib/mas.py Normal file
View File

@@ -0,0 +1,138 @@
from pyhdf.SD import SD,SDC
import os,sys
import ConfigParser
mas_units = {'length':6.96e10,
'time':1445.87,
'vr':481.3711*1.e5,
'vt':481.3711*1.e5,
'vp':481.3711*1.e5,
'n':1.e8,
'rho':1.6726e-16,
'p':0.3875717,
't':2.807067e7,
'br':2.2068908,
'bt':2.2068908,
'bp':2.2068908}
mas_var_names = ['t','rho','vt','vp','vr','bt','bp','br']
def read_var(fname,varname,normalized=False):
f = SD(fname,SDC.READ)
phi = f.select('fakeDim0')[:]
theta = f.select('fakeDim1')[:]
r = f.select('fakeDim2')[:]
var = f.select('Data-Set-2')[:]
f.end()
if normalized:
return(phi,theta,r,var)
else:
return(phi,theta,r*mas_units['length'],var*mas_units[varname])
def units():
return(mas_units)
def read_all_vars(dir,timeLabel,normalized=False):
vars = {}
for var_name in mas_var_names:
vars[var_name]={}
(vars[var_name]['phi'],
vars[var_name]['theta'],
vars[var_name]['r'],
vars[var_name]['data']) = read_var(os.path.join(dir,var_name+'%03d.hdf'%timeLabel),var_name,normalized)
return(vars)
def set_plot_limits(v,plotConfigFile = 'plot.config'):
plotConfig = ConfigParser.ConfigParser()
plotConfig.read(plotConfigFile)
# assume here v has the necessary keys
# bascially, this needs to be called after mas.read_all_varaibles
try:
v['bt']['lims'] = plotConfig.getfloat('MAS','bt_lim_lo'),plotConfig.getfloat('MAS','bt_lim_hi')
v['bt']['fmt'] = plotConfig.get('MAS','bt_fmt')
v['bp']['lims'] = plotConfig.getfloat('MAS','bp_lim_lo'),plotConfig.getfloat('MAS','bp_lim_hi')
v['bp']['fmt'] = plotConfig.get('MAS','bp_fmt')
v['br']['lims'] = plotConfig.getfloat('MAS','br_lim_lo'),plotConfig.getfloat('MAS','br_lim_hi')
v['br']['fmt'] = plotConfig.get('MAS','br_fmt')
v['vr']['lims'] = plotConfig.getfloat('MAS','vr_lim_lo'),plotConfig.getfloat('MAS','vr_lim_hi')
v['vr']['fmt'] = plotConfig.get('MAS','vr_fmt')
v['vt']['lims'] = plotConfig.getfloat('MAS','vt_lim_lo'),plotConfig.getfloat('MAS','vt_lim_hi')
v['vt']['fmt'] = plotConfig.get('MAS','vt_fmt')
v['vp']['lims'] = plotConfig.getfloat('MAS','vp_lim_lo'),plotConfig.getfloat('MAS','vp_lim_hi')
v['vp']['fmt'] = plotConfig.get('MAS','vp_fmt')
v['rho']['lims'] = plotConfig.getfloat('MAS','rho_lim_lo'),plotConfig.getfloat('MAS','rho_lim_hi')
v['rho']['fmt'] = plotConfig.get('MAS','rho_fmt')
v['t']['lims'] = plotConfig.getfloat('MAS','t_lim_lo'),plotConfig.getfloat('MAS','t_lim_hi')
v['t']['fmt'] = plotConfig.get('MAS','t_fmt')
except KeyError:
print('Some keys are missing from the variable dictionary.')
def mas2h5(fname,vname):
from numpy import meshgrid,sin,cos
if vname not in mas_var_names:
sys.exit('Wrong variable name')
fout = os.path.splitext(fname)[0]+'.h5'
import mas
phi,theta,r,var=mas.read_var(fname,vname)
import h5py
f=h5py.File(fout,'w')
P,T,R = meshgrid(phi,theta,r,indexing='ij')
x = R*sin(T)*cos(P)/mas_units['length']
y = R*sin(T)*sin(P)/mas_units['length']
z = R*cos(T)/mas_units['length']
f.create_dataset('x',data= x.astype(dtype='float32'))
f.create_dataset('y',data= y.astype(dtype='float32'))
f.create_dataset('z',data= z.astype(dtype='float32'))
f.create_dataset(vname,data= var.astype(dtype='float32'))
f.close()
# write xmf file
mas2xdmf(fout,vname,phi.size,theta.size,r.size)
def mas2xdmf(h5name,varname,n1,n2,n3):
fname=os.path.splitext(h5name)[0]+'.xmf'
dsname = os.path.basename(h5name)
f = open(fname,'w')
f.write('<?xml version="1.0" ?>' + '\n'+
'<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>' + '\n'+
'<Xdmf Version="2.0">' + '\n'+
' <Domain>' + '\n'+
' <Grid Name="mesh1" GridType="Uniform">' + '\n'+
' <Topology TopologyType="3DSMesh" NumberOfElements="%d %d %d"/>'%(n1,n2,n3) + '\n'+
' <Geometry GeometryType="X_Y_Z">' + '\n'+
' <DataItem Name="X" Dimensions="%d %d %d" NumberType="Float" Precision="4" Format="HDF">' %(n1,n2,n3)+ '\n'+
' %s:/x'%dsname + '\n'+
' </DataItem>' + '\n'+
' <DataItem Name="Y" Dimensions="%d %d %d" NumberType="Float" Precision="4" Format="HDF">' %(n1,n2,n3)+ '\n'+
' %s:/y'%dsname + '\n'+
' </DataItem>' + '\n'+
' <DataItem Name="Z" Dimensions="%d %d %d" NumberType="Float" Precision="4" Format="HDF">'%(n1,n2,n3) + '\n'+
' %s:/z'%dsname + '\n'+
' </DataItem>' + '\n'+
' </Geometry>' + '\n'+
' <Attribute Name="%s" AttributeType="Scalar" Center="Node">'%varname + '\n'+
' <DataItem Dimensions="%d %d %d" NumberType="Float" Precision="4" Format="HDF">'%(n1,n2,n3) + '\n'+
' %s:/%s'%(dsname,varname) + '\n'+
' </DataItem>' + '\n'+
' </Attribute>' + '\n'+
' </Grid>' + '\n'+
' </Domain>' + '\n'+
'</Xdmf>' )
f.close()

51
kaipy/gamhelio/lib/msgr.py Normal file
View File

@@ -0,0 +1,51 @@
import numpy as np
import datetime
def msgr(msgr_file,ave_1h=True):
r=[]
lat=[]
lon=[]
br=[]
bt=[]
bn=[]
dati=[]
for line in open(msgr_file,'r'):
fields = line.strip().split()
year = int(fields[0])
doy = int(fields[1])
hr = int(fields[2])
mn = int(fields[3])
r.append(float(fields[7]))
lat.append(float(fields[8]))
lon.append(float(fields[9]))
br.append(float(fields[10]))
bt.append(float(fields[11]))
bn.append(float(fields[12]))
dati.append(datetime.datetime(year,1,1,hr,mn)+datetime.timedelta(days=doy-1))
if ave_1h:
# assuming 1-min data. Average to 1-h
n=len(lat)-len(lat)%60
br_1h=np.mean(np.array(br)[:n].reshape(-1,60),1)
bt_1h=np.mean(np.array(bt)[:n].reshape(-1,60),1)
bn_1h=np.mean(np.array(bn)[:n].reshape(-1,60),1)
return (dati[30:n+30:60],br_1h,bt_1h,bn_1h)
else:
return (dati,np.array(br),np.array(bt),np.array(bn))
def msgr_v1(msgr_file,year):
data = np.loadtxt(msgr_file)
start_dati = datetime.datetime(year,1,1,0)
dati = [start_dati+datetime.timedelta(days=d-1) for d in data[:,0]]
br=data[:,4]
bt=data[:,5]
bn=data[:,6]
return (dati,br,bt,bn)

128
kaipy/gamhelio/lib/poisson.py Normal file
View File

@@ -0,0 +1,128 @@
from scipy.optimize import newton_krylov,anderson
from numpy import ediff1d,zeros_like,cos,sin,zeros_like,arange
class poisson():
def __init__(self,theta,phi):
self.theta = theta
self.phi = phi
self.nj = theta.shape[0]
self.nk = phi.shape[0]
self.dtheta = theta[1]-theta[0] # generalize to the case of arbitrary grid spacing later: ediff1d(theta)
self.dphi = phi[1]-phi[0] # ediff1d(phi)
def setRHS(self,RHS):
self.RHS = RHS
def residual_xsin2(self,P,lhs=False):
dphi = self.dphi
dtheta = self.dtheta
theta = self.theta
nk = self.nk
d2x = zeros_like(P)
d2y = zeros_like(P)
dy = zeros_like(P)
d2x[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2]) /dphi**2
d2y[1:-1,:] = (P[2:,:] - 2*P[1:-1,:] + P[:-2,:]) /dtheta**2*sin(theta[1:-1,None])**2
dy[1:-1,:] = 0.5*(P[2:,:] - P[:-2,:]) /dtheta*sin(theta[1:-1,None])*cos(theta[1:-1,None])
# periodic boundary
d2x[:,0] = (P[:,1]-2*P[:,0]+P[:,-1]) /dphi**2
d2x[:,-1] = (P[:,0]-2*P[:,-1]+P[:,-2]) /dphi**2
result = self.RHS*sin(theta[:,None])**2
# Pole boundary condition. Note, the below does not mean that d/dtheta=0 literally.
# We set sin(theta)*dPsi/dtheta=0 at the pole (see notes for details),
# and d2y here represents 1/sin(theta)*d/dtheta (sin(theta)*dPsi/dtheta) in that approxiamtion.
# NOTE, BELOW ASSUMES CONSTANT SPACING IN THETA!!! (0.5 FACTOR IS BASED ON THAT).
dy[0,:]=0.
dy[-1,:]=0.
d2y[0,:] = 0.5*(P[1,:]-P[0,:])
d2y[-1,:] = -0.5*(P[-1,:]-P[-2,:])
###################################################################################################
if not lhs:
return d2x + d2y + dy - result
else:
return d2x + d2y + dy
def residual(self,P,lhs=False):
dphi = self.dphi
dtheta = self.dtheta
theta = self.theta
nk = self.nk
d2x = zeros_like(P)
d2y = zeros_like(P)
dy = zeros_like(P)
d2x[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2]) /dphi**2/sin(theta[:,None])**2
d2y[1:-1,:] = (P[2:,:] - 2*P[1:-1,:] + P[:-2,:]) /dtheta**2
dy[1:-1,:] = 0.5*(P[2:,:] - P[:-2,:]) /dtheta/sin(theta[1:-1,None])*cos(theta[1:-1,None])
# periodic boundary
d2x[:,0] = (P[:,1]-2*P[:,0]+P[:,-1]) /dphi**2/sin(theta)**2
d2x[:,-1] = (P[:,0]-2*P[:,-1]+P[:,-2]) /dphi**2/sin(theta)**2
result = self.RHS
# Pole boundary condition. Note, the below does not mean that d/dtheta=0 literally.
# We set sin(theta)*dPsi/dtheta=0 at the pole (see notes for details),
# and d2y here represents 1/sin(theta)*d/dtheta (sin(theta)*dPsi/dtheta) in that approxiamtion.
# NOTE, BELOW ASSUMES CONSTANT SPACING IN THETA!!! (0.5 FACTOR IS BASED ON THAT).
dy[0,:]=0.
dy[-1,:]=0.
d2y[0,:] = 0.5*(P[1,:]-P[0,:])/sin(theta[0])**2 # Psin=0; 0.5*( 0.5*(P[1,:]+P[0,:]) - Psin)
d2y[-1,:] = -0.5*(P[-1,:]-P[-2,:])/sin(theta[-1])**2 #-0.5*( 0.5*(P[-2,:]+P[-1,:]) - P[-1,:].mean())
###################################################################################################
if not lhs:
return d2x + d2y + dy - result
else:
return d2x + d2y + dy
if __name__=='__main__':
from scipy.optimize import newton_krylov,anderson
from pylab import *
import wsa
(phi_wsa_v,theta_wsa_v,phi,theta,br0,v_wsa,n_wsa,T_wsa)=wsa.read('/Users/merkivg1/work/LFM-helio_2.0/WSA/fits/2010_2008_WSA_ADAPT/vel_201001022300R001_ans.fits',False)
(phi_wsa_v,theta_wsa_v,phi,theta,br1,v_wsa,n_wsa,T_wsa)=wsa.read('/Users/merkivg1/work/LFM-helio_2.0/WSA/fits/2010_2008_WSA_ADAPT/vel_201001032300R001_ans.fits',False)
omega=2*pi/27.27
phi_prime=phi-omega*1.
ind0=where(phi_prime>0)[0][0]
phi_prime=roll(phi_prime,-ind0)
phi_prime[phi_prime<0]+=2*pi
import scipy
from scipy.interpolate import interpolate
br1_rolled=roll(br1,-ind0,axis=1)
fbi = scipy.interpolate.RectBivariateSpline(theta,phi,br1_rolled,kx=1,ky=1)
br1_new=fbi(theta,phi_prime)
import poisson
pois = poisson.poisson(theta,phi)
pois.setRHS(br1_new - br0)
guess=zeros_like(br1_new)
sol = newton_krylov(pois.residual,guess, method='lgmres', verbose=1,f_rtol=1.e-6) #iter=100
print('Residual: %g' % abs(pois.residual(sol)).max())
figure();pcolormesh(sol);colorbar();title('Solution')
# check curl
Etheta = -1/sin(theta[:,None])*diff(sol,axis=1)
Ephi = diff(sol,axis=0)
figure();pcolormesh(Etheta);colorbar();title('Etheta')
figure();pcolormesh(Ephi);colorbar();title('Ephi')

24
kaipy/gamhelio/lib/util.py Normal file
View File

@@ -0,0 +1,24 @@
import numpy as np
def sph2cart(a,theta,phi):
ar,at,ap = a
ax = ar*cos(phi)*sin(theta)+at*cos(phi)*cos(theta)-ap*sin(phi)
ay = ar*sin(phi)*sin(theta)+at*sin(phi)*cos(theta)+ap*cos(phi)
az = ar*cos(theta)-at*sin(theta)
return(ax,ay,az)
def radial_interp(vars_from,var_name,radii_to,coef_label='r_interp_coefs',ind_label='r_interp_above_ind'):
q=[]
ind=[]
for r_to in radii_to: # interpolating to LFM cell centers
r_from = vars_from[var_name]['r']
ind_r_a = np.flatnonzero(r_from>=r_to)[0] # MAS radius above
r_a = r_from[ind_r_a]
r_b = r_from[ind_r_a-1]
q.append( (r_a-r_to)/(r_a-r_b) )
ind.append(ind_r_a)
vars_from[var_name][coef_label]=np.asarray(q,dtype='float')
vars_from[var_name][ind_label]=np.asarray(ind,dtype='int')

88
kaipy/gamhelio/lib/wsa.py Executable file
View File

@@ -0,0 +1,88 @@
#!/usr/bin/env python
from astropy.io import fits
from numpy import linspace,pi,meshgrid,sin,cos,zeros,ones,dstack,diff,sqrt,array,savetxt,flatnonzero,insert,asarray,zeros_like,argmin,unravel_index
def read(wsa_file,densTempInfile,normalized=False,verbose=True):
if verbose: info(wsa_file)
hdul = fits.open(wsa_file)
n_phi_wsa_v = hdul[0].header['NAXIS1']+1 # number of cell vertices
n_phi_wsa_c = hdul[0].header['NAXIS1'] # number of cell centers
phi_wsa_v = linspace(0,360,n_phi_wsa_v)/180.*pi
phi_wsa_c = 0.5*(phi_wsa_v[:-1]+phi_wsa_v[1:])
n_theta_wsa_v = hdul[0].header['NAXIS2']+1 # number of cell vertices
n_theta_wsa_c = hdul[0].header['NAXIS2'] # number of cell centers
theta_wsa_v = linspace(0,180,n_theta_wsa_v)/180.*pi
theta_wsa_c = 0.5*(theta_wsa_v[:-1]+theta_wsa_v[1:])
bi_wsa = hdul[0].data[0,::-1,:] # note the theta reversal to convert from wsa theta to lfm theta definition
v_wsa = hdul[0].data[1,::-1,:] # note the theta reversal to convert from wsa theta to lfm theta definition
if densTempInfile:
n_wsa = hdul[0].data[2,::-1,:]
T_wsa = hdul[0].data[3,::-1,:]
else:
n_wsa = 112.64+9.49e7/v_wsa**2
T0 = 8e5
n0 = 300.
B0 = bi_wsa[unravel_index(argmin(abs(n_wsa-n0)),n_wsa.shape)] # this is in nT
n0 = n_wsa[unravel_index(argmin(abs(n_wsa-n0)),n_wsa.shape)] # this is in cm^-3
T_wsa = n0*T0/n_wsa
# The code below allows total pressure conservation in the non-radial
# direction: nT+B**2/8pi = n0*T0+B0**2/8pi. The unit conversion assumes
# B and B0 in nT (which they should be by now) and botzmann constant in the
# denominator.
# T_wsa = n0*T0/n_wsa + (B0**2-bi_wsa**2)/1.38/8./pi*1.e6/n_wsa
hdul.close()
if normalized:
return(phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa,v_wsa,n_wsa*1.67e-24,T_wsa)
else:
return(phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa*1.e-5,v_wsa*1.e5,n_wsa*1.67e-24,T_wsa)
def info(wsa_file):
hdul = fits.open(wsa_file)
print(repr(hdul[0].header))
hdul.close()
def plot(wsa_file,savefig):
hdul = fits.open(wsa_file)
bi_wsa = hdul[0].data[0,:,:]
v_wsa = hdul[0].data[1,:,:]
hdul.close()
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(16,12))
ax1 = plt.subplot(211)
p1=ax1.pcolormesh(bi_wsa,cmap='RdBu_r',vmin=bi_wsa.min(),vmax=-bi_wsa.min())
plt.colorbar(p1,ax=ax1).set_label('Br')
ax1.set_xlim((0,bi_wsa.shape[1]))
ax1.set_ylim((0,bi_wsa.shape[0]))
ax2 = plt.subplot(212,sharex=ax1)
p2=ax2.pcolormesh(v_wsa)
plt.colorbar(p2,ax=ax2).set_label('V')
ax2.set_xlim((0,bi_wsa.shape[1]))
ax2.set_ylim((0,bi_wsa.shape[0]))
fig.suptitle(wsa_file)
if not savefig:
plt.show()
else:
plt.savefig(wsa_file[:-4]+'png')
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('wsaFile',help='WSA file to use')
parser.add_argument('--savefig',help='Drop figure?',default=False,action='store_true')
args = parser.parse_args()
info(args.wsaFile)
plot(args.wsaFile,args.savefig)

31
kaipy/gamhelio/lib/wsa2h5.py Executable file
View File

@@ -0,0 +1,31 @@
#!/usr/bin/env python
# vgm, Sep 2019: Script to convert WSA fits files into simple h5 files
from astropy.io import fits
import h5py
import glob
wsaFiles = glob.glob('./*.fits')
for wsaFile in wsaFiles:
h5file = wsaFile[:-4]+'h5'
with fits.open(wsaFile) as hdul:
with h5py.File(h5file,'w') as hf:
# take care of the comments
# removing the last ones since we're incorporating them into dataset attributes directly
for key in hdul[0].header[:-4]:
# stores header data as attributes of the file
hf.attrs[key] = hdul[0].header[key]
hf.attrs[key+" COMMENT"] = str(hdul[0].header.comments[key])
# take care of the datasets
br = hf.create_dataset("br",data=hdul[0].data[0])
br.attrs['Units'] = 'nT'
br.attrs['Comment'] = 'Coronal field at RADOUT'
vr = hf.create_dataset("vr",data=hdul[0].data[1])
vr.attrs['Units'] = 'km/s'
vr.attrs['Comment'] = 'Radial velocity at RADOUT'

1
kaipy/gamhelio/wsa2gamera/__init__.py Normal file
View File

@@ -0,0 +1 @@
from . import params

21
kaipy/gamhelio/wsa2gamera/params.py Normal file
View File

@@ -0,0 +1,21 @@
import configparser
class params():
def __init__(self,ConfigFileName):
config = configparser.ConfigParser(inline_comment_prefixes=(';','#'))
config.read(ConfigFileName)
self.gameraGridFile = config['Gamera']['gameraGridFile']
self.dirGameraGridFile = config['Gamera']['Dir']
self.wsaFile = config['WSA']['wsafile']
self.gaussSmoothWidth = config.getint('WSA','gauss_smooth_width')
self.plots = config.getboolean('WSA','plots')
self.densTempInfile = config.getboolean('WSA','density_temperature_infile')
self.normalized = config.getboolean('WSA','normalized')
self.gamma = config.getfloat('Constants','gamma')
self.NO2 = config.getint('Constants','NO2')
self.B0 = config.getfloat('Normalization','B0')
self.n0 = config.getfloat('Normalization','n0')

128
scripts/wsa2gamera.py Executable file
View File

@@ -0,0 +1,128 @@
#! /usr/bin/env python
helioLibLocation='/glade/u/home/elenap/kaiju/kaipy/gamhelio/lib'
import numpy as np
import os,sys,glob
from scipy import interpolate
import time
import h5py
import matplotlib.pyplot as plt
#from wsa2gamera.params import params
#if helioLibLocation not in sys.path: sys.path.append(helioLibLocation)
#import wsa
import kaipy.gamhelio.wsa2gamera.params as params
import kaipy.gamhelio.lib.wsa as wsa
#----------- 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(args.ConfigFileName)
Ng=prm.NO2
gamma = prm.gamma
B0 = prm.B0
n0 = prm.n0
# constants
mp = 1.67e-24
############### WSA STUFF #####################
phi_wsa_v,theta_wsa_v,phi_wsa_c,theta_wsa_c,bi_wsa,v_wsa,n_wsa,T_wsa = wsa.read(prm.wsaFile,prm.densTempInfile,prm.normalized)
# convert the units; remember to use the same units in gamera
# TODO: probably store units in the h5 file?
# B0 = 1.e-3 Gs
# n0 = 200./cc
V0 = B0/np.sqrt(4*np.pi*mp*n0)
bi_wsa /= B0
n_wsa /= (mp*n0)
v_wsa /= V0
# keep temperature in K
############### WSA STUFF #####################
############### GAMERA STUFF #####################
# GAMERA GRID
with h5py.File(os.path.join(prm.dirGameraGridFile,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:])
# 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
P = np.arctan2(y[Ng:-Ng-1,Ng:-Ng-1,:],x[Ng:-Ng-1,Ng:-Ng-1,:])
P[P<0]=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.
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]=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.T,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
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.T,kx=1,ky=1)
vr = fv(Pc[:,0,0],Tc[0,:,0])
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,n_wsa.T,kx=1,ky=1)
rho = f(Pc[:,0,0],Tc[0,:,0])
f = interpolate.RectBivariateSpline(phi_wsa_c,theta_wsa_c,T_wsa.T,kx=1,ky=1)
temp = f(Pc[:,0,0],Tc[0,:,0])
# note, redefining interpolation functions we could also
# interpolate from bi_wsa as above, but then we would have to
# smooth bk, if necessary. The way we're doing it here, bk will be
# smoothed or not, dependent on whether br has been smoothed.
# note also, this has to extrapolate
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])
# 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
with h5py.File('innerbc.h5','w') as hf:
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)