In [1]:
from __future__ import print_function
In [2]:
from functools import partial
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import astropy.units as u
import numpy as np
from astropy.modeling import models, fitting
Now we try and build the mfe single tube model.
In [3]:
import pysac.mhs_atmosphere as atm
from pysac.mhs_atmosphere.parameters.model_pars import mfe_setup as model_pars
In [4]:
# Cheeky Reset to Photosphere
#model_pars['xyz'][4] = 36641.22*u.m
#model_pars['xyz'][5] = 1587786.3*u.m
#model_pars['B_corona'] = 0.00055*u.T
#model_pars['radial_scale'] = 0.03938*u.Mm
##model_pars['radial_scale'] = 0.05*u.Mm
model_pars['pBplus'] = 1.2e-3*u.T
#model_pars['chrom_scale'] = 0.4*u.Mm
In [5]:
keys = ['pBplus', 'B_corona', 'coratio', 'chratio', 'radial_scale', 'chrom_scale', 'corona_scale', 'xyz']
for key in keys:
print(key, model_pars[key])
In [6]:
# model setup
scales, physical_constants = atm.units_const.get_parameters()
option_pars = atm.options.set_options(model_pars, False, l_gdf=True)
coords = atm.model_pars.get_coords(model_pars['Nxyz'], u.Quantity(model_pars['xyz']))
#interpolate the hs 1D profiles from empirical data source[s]
empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(mu=physical_constants['mu'])
table = atm.hs_atmosphere.interpolate_atmosphere(empirical_data, coords['Zext'])
In [ ]:
table['rho'] += table['rho'].min()*3.6
In [ ]:
#==============================================================================
#calculate 1d hydrostatic balance from empirical density profile
#==============================================================================
# the hs pressure balance is enhanced by pressure equivalent to the
# residual mean coronal magnetic pressure contribution once the magnetic
# field has been applied
magp_meanz = np.ones(len(coords['Z'])) * u.one
magp_meanz *= model_pars['pBplus']**2/(2*physical_constants['mu0'])
pressure_z, rho_z, Rgas_z = atm.hs_atmosphere.vertical_profile(coords['Z'], table, magp_meanz,
physical_constants, coords['dz'])
In [ ]:
x, y, z = u.Quantity(np.mgrid[coords['xmin']:coords['xmax']:1j*model_pars['Nxyz'][0],
coords['ymin']:coords['ymax']:1j*model_pars['Nxyz'][1],
coords['zmin']:coords['zmax']:1j*model_pars['Nxyz'][2]], unit=coords['xmin'].unit)
In [ ]:
xi, yi, Si = atm.flux_tubes.get_flux_tubes(model_pars, coords, option_pars)
In [ ]:
Si[0][0] = 0.1436*u.T
In [ ]:
allmag = atm.flux_tubes.construct_magnetic_field(x, y, z,
xi[0], yi[0], Si[0],
model_pars, option_pars,
physical_constants,
scales)
pressure_m, rho_m, Bx, By ,Bz, Btensx, Btensy = allmag
In [ ]:
# local proc 3D mhs arrays
pressure, rho = atm.mhs_3D.mhs_3D_profile(z,
pressure_z,
rho_z,
pressure_m,
rho_m)
magp = (Bx**2 + By**2 + Bz**2) / (2.*physical_constants['mu0'])
energy = atm.mhs_3D.get_internal_energy(pressure, magp, physical_constants)
In [ ]:
x_lin = np.linspace(coords['xmin'], coords['xmax'], model_pars['Nxyz'][0])
bmag = np.sqrt((Bx**2 + By**2 + Bz**2))[:,64,0].to(u.mT)
gaussian = models.Gaussian1D()
bmag_fit = fitting.LevMarLSQFitter()(gaussian, x_lin, bmag)
fwhm = 2.*np.sqrt(2*np.log(2))*bmag_fit.stddev.value
fwhm = np.abs(fwhm) * u.Mm
fwtm = 2.*np.sqrt(2*np.log(10))*bmag_fit.stddev.value
fwtm = np.abs(fwtm) * u.Mm
u.Quantity([fwhm, fwtm]).to(u.km)
In [ ]:
%matplotlib inline
fig = plt.figure(figsize=(15,8))
plt.plot(x_lin, bmag.to(u.mT), 'x')
X = np.linspace(coords['xmin'], coords['xmax'], 1000)
plt.plot(X, bmag_fit(X.value))
plt.xlabel("X [{}]".format(x.unit))
plt.ylabel("Magnetic Field Strength [{}]".format(bmag.unit))
plt.axhline(y=bmag_fit.amplitude.value/10)
In [ ]:
import yt
In [ ]:
def magnetic_field_strength(field, data):
return np.sqrt(data["magnetic_field_x"]**2 + data["magnetic_field_y"]**2 + data["magnetic_field_z"]**2)
yt.add_field(("gas","magnetic_field_strength"),
function=magnetic_field_strength,
units = "T")
In [ ]:
def plasma_beta(field, data):
return data['pressure'] / data['magnetic_pressure']
yt.add_field(("gas","plasma_beta"),
function=plasma_beta,
units = "")
In [ ]:
def alfven_speed(field, data):
return np.sqrt(2.*data['magnetic_pressure']/data['density'])
yt.add_field(("gas","alfven_speed"),
function=alfven_speed,
units = "m/s")
In [ ]:
def sound_speed(field, data):
cspeed = np.sqrt(physical_constants['gamma']*pressure/rho)
return cspeed
yt.add_field(("gas","sound_speed"),
function=sound_speed,
units = "m/s", force_override=True)
In [ ]:
bbox = u.Quantity([u.Quantity([coords['xmin'], coords['xmax']]),
u.Quantity([coords['ymin'], coords['ymax']]),
u.Quantity([coords['zmin'], coords['zmax']])]).to(u.m).value
In [ ]:
data = {'magnetic_field_x':yt.YTQuantity.from_astropy(Bx.decompose()),
'magnetic_field_y':yt.YTQuantity.from_astropy(By.decompose()),
'magnetic_field_z':yt.YTQuantity.from_astropy(Bz.decompose()),
'pressure': yt.YTQuantity.from_astropy(pressure.decompose()),
'magnetic_pressure': yt.YTQuantity.from_astropy(magp.decompose()),
'density': yt.YTQuantity.from_astropy(rho.decompose())}
ds = yt.load_uniform_grid(data, x.shape,
length_unit='m', magnetic_unit='T', mass_unit='kg',
periodicity=[False]*3, bbox=bbox)
In [ ]:
ds.index.grids[0]['density'].min(), ds.index.grids[0]['magnetic_field_z'].max(),
In [ ]:
ds.index.grids[0]['pressure'].min()
In [ ]:
ds.index.grids[0]['magnetic_field_z'][64,64,:]
In [ ]:
slc = yt.SlicePlot(ds, 'x', 'alfven_speed', origin='lower-center-domain', axes_unit='Mm')
slc.set_cmap('all', 'OrRd')
seed_points = np.zeros([11,2]) + 1.52
seed_points[:,0] = np.linspace(-0.99, 0.95, seed_points.shape[0], endpoint=True)
slc.annotate_streamlines('magnetic_field_y', 'magnetic_field_z', field_color='magnetic_field_strength',
plot_args={'start_points':seed_points, 'density':15, 'cmap':'winter',
'norm': mpl.colors.LogNorm(*ds.all_data().quantities.extrema("magnetic_field_strength"))})
slc.annotate_contour('plasma_beta', clim=(1,50), label=True, take_log=False)
In [ ]:
slc = yt.SlicePlot(ds, 'x', ['pressure','magnetic_field_strength'], origin='lower-center-domain',
axes_unit='Mm')
slc.set_cmap('all', 'BuPu')
#slc.annotate_contour('plasma_beta')
seed_points = np.zeros([11,2]) + 1.52
seed_points[:,0] = np.linspace(-0.99, 0.95, seed_points.shape[0], endpoint=True)
slc.annotate_streamlines('magnetic_field_y', 'magnetic_field_z', field_color='magnetic_field_strength',
plot_args={'start_points':seed_points, 'density':15, 'cmap':'winter', 'linewidth': 2,
'norm': slc.plots['magnetic_field_strength'].image.norm})
slc.annotate_contour('plasma_beta', label=True, take_log=False, clim=(1,50))
In [ ]:
slc = yt.SlicePlot(ds, 'z', 'magnetic_field_strength', origin='center-domain',
axes_unit='km', center=("max", "magnetic_field_strength"), width=(0.5, "Mm"))
slc.set_cmap('all', 'Oranges')
slc.annotate_contour('magnetic_field_strength', take_log=False,
plot_args={'levels':[ds.all_data().quantities.extrema("magnetic_field_strength")[1]/2.],
'linewidths':3,
'colors':'black'})
In [ ]:
In [ ]: