In [ ]:
%matplotlib inline

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

import astropy.units as u

import yt
import pysac.yt

import h5py

from astropy.modeling import models, fitting

In [ ]:
ds = yt.load('/nocrypt/Slog_p240-0_A10_B005_00001.gdf')
type(ds)

In [ ]:
cg = ds.index.grids[0]

xmin = ds.domain_left_edge[0].convert_to_units('Mm')
xmax = ds.domain_right_edge[0].convert_to_units('Mm')
nx = ds.domain_dimensions[0]


f5 = h5py.File(ds.filename)
x_lin = f5['x'][0][:,63,0]
x_lin = x_lin*u.m
x_lin = x_lin.to(u.Mm)
#x_lin = np.linspace(xmin, xmax, nx)
bmag = cg['magnetic_field_strength'][:,63,0].to_astropy().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 [ ]:
fig = plt.figure(figsize=(15,8))
plt.plot(x_lin, bmag.to(u.mT), 'x')
X = np.linspace(xmin, xmax, 1000)
plt.plot(X, bmag_fit(X.value))

plt.xlabel("X [{}]".format(x_lin.unit))
plt.ylabel("Magnetic Field Strength [{}]".format(bmag.unit))
plt.axhline(y=bmag_fit.amplitude.value/10)

In [ ]:
f5['x'][2][0,0,:][-1]

In [ ]:
_.to(u.Mm)

In [ ]:
x_lin

In [ ]:
x_lin[1:] - x_lin[:-1]

In [ ]:
np.save('mfe_x_coords_Mm.npy', x_lin.value)

In [ ]:
np.save('mfe_bmag_x_mT.npy', bmag.value)

In [ ]:
import numpy as np

from astropy.modeling import models, fitting
import astropy.units as u


x_lin = np.load('mfe_x_coords_Mm.npy')*u.Mm
bmag = np.load('mfe_bmag_x_mT.npy')*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 [ ]:
cg['plasma_beta'].max(), cg['plasma_beta'].min()

In [ ]:
cg['density'].max(), cg['density'].min()

In [ ]:
cg['thermal_pressure'].max(), cg['thermal_pressure'].min()

In [ ]:
slc = yt.SlicePlot(ds, 'x', 'thermal_pressure', 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('mag_field_y', 'mag_field_z', field_color='magnetic_field_strength',
                         plot_args={'start_points':seed_points, 'density':15, 'cmap':'winter'})

slc.annotate_contour('plasma_beta', clim=(0.1,50),label=True, take_log=False)

In [ ]:
cg['mag_field_z'][63,63,:]

In [ ]:
np.save('mfe_Bz_2015.npy', _)

In [ ]: