In [1]:

    

%%file output/custom.css
body {
    font-size: 160%;
    font-family: Lato, Ariel, sans-serif !important;
}

div.slides {
    margin-top: -7%;
}

.left {
    width: 49%;
    float: left;
}

.right {
    width: 49%;
    float: right;
}

.centre {
    text-align: center;
}

h1 {
    text-align: center;
}

h2 {
    text-align: center;
}

div.prompt {
    display: none;
}

div.highlight {
    font-size: 85%; /* This Seems to give an approximate 80char line at 1024 */
}

div.output_html {
    font-size: 85%;
}

div.output_subarea {
    max-width: 100% !important;
}

div.output_png {
    text-align: center;
}

li {
    padding-top: 1em !important;
}

ul.logos {
    margin: 0px;
    padding: 0px;
    width: 100%;
}

ul.logos li {
    list-style: none;
    height:150px;

}

.output_stderr {
    display: none;
}

div.output_area .rendered_html img, div.output_area .rendered_html table {
    margin-left: auto;
    margin-right: auto;
}

.data-table {
    border-collapse: collapse;
}

.border-bottom {
    border-bottom: 1px solid #000;
}

.align-center {
    text-align: center;
}

.reveal i {
    font-size: inherit !important;
    font-style: italic !important;
}

.cite2c-dialog .twitter-typeahead {
    margin: 8px;
    width: 90%;
}

.cite2c-dialog .tt-input {
    width: 100%;
}

.cite2c-dialog .tt-dropdown-menu {
    background-color: #fff;
    border: 1px solid #ccc;
    border-radius: 4px;
}

.cite2c-dialog .tt-suggestion {
    padding: 4px;
}
.cite2c-dialog .tt-suggestion.tt-cursor {
    background-color: #eee;
}

.csl-entry {
    padding: 0.4em 0 !important;
}

.rendered_html blockquote {
    font-family: serif;
    font-size: 150%;
    line-height: 130%;
    margin: 5px auto;
    width: 90%;
    text-align: center;
}


    

    




Overwriting output/custom.css

In [2]:

    

import sys

sys.path.append('./python/')

from mayavi import mlab
from mayavi.tools.sources import vector_field, scalar_field
mlab.options.offscreen = True

import numpy as np

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import astropy.units as u
import sunpy.map

import yt
import pysac.yt

%matplotlib inline

font_size = 20
pgf_with_latex = {
    "font.size": font_size,
    "axes.labelsize": font_size,  # LaTeX default is 10pt font.
    "legend.fontsize":font_size + 5,
    "xtick.labelsize": font_size,
    "ytick.labelsize": font_size,
    "ytick.labelsize": font_size,
    "savefig.transparent": True,
    "animation.html": "html5"
    }

matplotlib.rcParams.update(pgf_with_latex)


    

    




WARNING:traits.has_traits:DEPRECATED: traits.has_traits.wrapped_class, 'the 'implements' class advisor has been deprecated. Use the 'provides' class decorator.

In [3]:

    

x_range, y_range = [-300, 300]*u.arcsec, [-250, 250]*u.arcsec

plt.ioff()
fig = plt.figure(dpi=50, figsize=(11,8))
mm = sunpy.map.Map('/home/stuart/VivaData/gband_data/Gband_cospatial_cotemporal_00000.fits').submap(x_range, y_range)
mm = mm.submap([-440,440]*u.arcsec, [-440,440]*u.arcsec)
mm.plot_settings['cmap'] = 'gray'
mm.plot_settings['title'] = ''
im = mm.plot()
fig.savefig('./images/gband-plot.png', transparent=True)


    

In [4]:

    

import pysac.mhs_atmosphere as atm

#Read in the VAL3C model
empirical_data = atm.hs_atmosphere.read_VAL3c_MTW(MTW_file=False)


# Create a Z array at the interpolated resolution and interpolate.
ZZ = u.Quantity(np.linspace(empirical_data['Z'][0], empirical_data['Z'][-1], 128), unit=empirical_data['Z'].unit)
table = atm.hs_atmosphere.interpolate_atmosphere(empirical_data, ZZ, s=0)


# Create a figure and make space for the axes.
fig, ax = plt.subplots(gridspec_kw={'right':0.77, 'left':0.16, 'bottom':0.13}, figsize=(13,8))

# Shortcut all the Mm conversion.
Z = empirical_data['Z'].to(u.Mm)

lrho, = ax.plot(Z, empirical_data['rho'].quantity.si, 'x', color='blue')
lrho_i, = ax.plot(ZZ.to(u.Mm), table['rho'].quantity.si, color='blue')

ax2 = ax.twinx()
lp, = ax2.plot(Z, empirical_data['p'].to(u.Pa), 'x', color='green')
lp_i, = ax2.plot(ZZ.to(u.Mm), table['p'].to(u.Pa), color='green')


ax3 = ax.twinx()
ax3.spines["right"].set_position(("axes", 1.2))
lt, = ax3.plot(Z, empirical_data['T'].to(u.K), 'x', color='red')
lt_i, = ax3.plot(ZZ.to(u.Mm), table['T'].to(u.K), color='red')


# Set primary axes properties and labels
ax.semilogy()
ax.set_ylabel(r"Density [{}]".format(lrho._yorig.unit._repr_latex_()))
ax.set_xlabel(r"Height [{}]".format(lrho._xorig.unit._repr_latex_()))
ax.set_xlim(Z[0].value-0.1, Z[-1].value+0.1)


# Pressure Axis
ax2.semilogy()
ax2.set_ylabel(r"Pressure [{}]".format(lp._yorig.unit._repr_latex_()))


# Temp axis
ax3.set_ylabel(r"Temperature [{}]".format(lt._yorig.unit._repr_latex_()))

ax.set_xlim([-0.02,1.62])
ax3.set_ylim([3500,7500])
# Set the colours for the ticks and the labels.
ax.tick_params(axis='y', colors=lrho.get_color())
ax2.tick_params(axis='y', colors=lp.get_color())
ax3.tick_params(axis='y', colors=lt.get_color())

ax.yaxis.label.set_color(lrho.get_color())
ax2.yaxis.label.set_color(lp.get_color())
ax3.yaxis.label.set_color(lt.get_color())
fig


    

    
Out[4]:

In [5]:

    

from pysac.mhs_atmosphere.parameters.model_pars import mfe_setup as model_pars
import pysac.mhs_atmosphere as atm

# Cheeky Reset to Photosphere
model_pars['xyz'][4] = 0*u.Mm
#==============================================================================
# Build the MFE flux tube model using pysac
#==============================================================================
# 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'])
table['rho'] += table['rho'].min()*3.6

# 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'])

# Make the vertical profile 3D
pressure_z, rho_z, Rgas_z = atm.hs_atmosphere.vertical_profile(coords['Z'], table, magp_meanz,
                                                               physical_constants, coords['dz'])

# Generate 3D coordinate arrays
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)

# Get default MFE flux tube parameters out of pysac
xi, yi, Si = atm.flux_tubes.get_flux_tubes(model_pars, coords, option_pars)

# Generate the 3D magnetic field
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

# 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)

#### YT STUFF ####

magnetic_field_x = lambda field, data: data['mag_field_x']
yt.add_field(("gas","magnetic_field_x"), function=magnetic_field_x, units=yt.units.T.units)
magnetic_field_y = lambda field, data: data['mag_field_y']
yt.add_field(("gas","magnetic_field_y"), function=magnetic_field_y, units=yt.units.T.units)
magnetic_field_z = lambda field, data: data['mag_field_z']
yt.add_field(("gas","magnetic_field_z"), function=magnetic_field_z, units=yt.units.T.units)

# Add derived Fields
def magnetic_field_strength(field, data):
    return np.sqrt(data["mag_field_x"]**2 + data["mag_field_y"]**2 + data["mag_field_z"]**2)
yt.add_field(("gas","magnetic_field_strength"), function=magnetic_field_strength, units=yt.units.T.units)

#def alfven_speed(field, data):
#    return np.sqrt(2.*data['magnetic_pressure']/data['density'])
#yt.add_field(("gas","alfven_speed"), function=alfven_speed, units=(yt.units.m/yt.units.s).units)

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

# Now build a yt DataSet with the generated data:
data = {'mag_field_x':yt.YTQuantity.from_astropy(Bx.decompose()),
        'mag_field_y':yt.YTQuantity.from_astropy(By.decompose()),
        'mag_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)


    

    




/home/stuart/GitHub/SWAT/pysac/pysac/mhs_atmosphere/mhs_model/flux_tubes.py:312: Warning: pbbal.max() = -0.111891563752 kg / (m s2)
  warnings.warn("pbbal.max() = {}".format(pbbal.max().decompose()), Warning)
yt : [INFO     ] 2015-12-04 14:48:45,830 Parameters: current_time              = 0.0
yt : [INFO     ] 2015-12-04 14:48:45,831 Parameters: domain_dimensions         = [129 129 128]
yt : [INFO     ] 2015-12-04 14:48:45,832 Parameters: domain_left_edge          = [-1000000. -1000000.        0.]
yt : [INFO     ] 2015-12-04 14:48:45,834 Parameters: domain_right_edge         = [ 1000000.   1000000.   1587786.3]
yt : [INFO     ] 2015-12-04 14:48:45,835 Parameters: cosmological_simulation   = 0.0

In [6]:

    

slc = yt.SlicePlot(ds, 'x', 'density', origin='lower-center-domain', axes_unit='Mm')
slc.set_figure_size(9)
slc.set_cmap('all', 'viridis')
slc.set_font_size(20)
slc.set_zlim('all', 0, 2.5e-4)


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':'Blues', 'linewidth':2,
'norm': matplotlib.colors.LogNorm(*ds.all_data().quantities.extrema("magnetic_field_strength"))})
#force render
slc.save('/tmp/test.png')
#use the raw Figure for transparent bg
slc.plots['density'].figure


    

    




yt : [INFO     ] 2015-12-04 14:48:45,929 Loading field plugins.
yt : [INFO     ] 2015-12-04 14:48:45,931 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,932 Loaded astro (15 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,932 Loaded cosmology (22 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,934 Loaded fluid (63 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,935 Loaded fluid_vector (95 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,936 Loaded geometric (111 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,937 Loaded local (115 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,937 Loaded magnetic_field (121 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,938 Loaded my_plugins (121 new fields)
yt : [INFO     ] 2015-12-04 14:48:45,939 Loaded species (123 new fields)
yt : [INFO     ] 2015-12-04 14:48:46,276 xlim = -1000000.000000 1000000.000000
yt : [INFO     ] 2015-12-04 14:48:46,276 ylim = 0.000000 1587786.300000
yt : [INFO     ] 2015-12-04 14:48:46,279 xlim = -1000000.000000 1000000.000000
yt : [INFO     ] 2015-12-04 14:48:46,280 ylim = 0.000000 1587786.300000
yt : [INFO     ] 2015-12-04 14:48:46,290 Making a fixed resolution buffer of (('gas', 'density')) 800 by 800
yt : [WARNING  ] 2015-12-04 14:48:46,343 Plot image for field ('gas', 'density') has both positive and negative values. Min = -0.000032, Max = 0.000267.
yt : [WARNING  ] 2015-12-04 14:48:46,344 Switching to symlog colorbar scaling unless linear scaling is specified later
yt : [INFO     ] 2015-12-04 14:48:47,927 Saving plot /tmp/test.png






    
Out[6]:

In [7]:

    

#Define tvtk notebook viewer
from IPython.core.display import Image
import subprocess
def mlab_view(scene, azimuth=153, elevation=62, distance=400, focalpoint=np.array([  25.,   63.,  60.]), aa=16):
    scene.anti_aliasing_frames = aa
    mlab.view(azimuth=azimuth, elevation=elevation, distance=distance, focalpoint=focalpoint)
    scene.save('offscreen.png', size=(750, 750))
    subprocess.call(["convert", "offscreen.png", "-transparent", "white", "offscreen.png"])
    return Image(filename='offscreen.png')


    

In [8]:

    

ds = pysac.yt.SACGDFDataset('/home/stuart/VivaData/Slog_p240-0_A10_B005_00001.gdf')


    

    




yt : [WARNING  ] 2015-12-04 14:48:49,348 'field_units' was overridden by 'dataset_units/density_bg'
yt : [WARNING  ] 2015-12-04 14:48:49,350 'field_units' was overridden by 'dataset_units/density_pert'
yt : [WARNING  ] 2015-12-04 14:48:49,354 'field_units' was overridden by 'dataset_units/internal_energy_bg'
yt : [WARNING  ] 2015-12-04 14:48:49,356 'field_units' was overridden by 'dataset_units/internal_energy_pert'
yt : [WARNING  ] 2015-12-04 14:48:49,360 'field_units' was overridden by 'dataset_units/mag_field_x_bg'
yt : [WARNING  ] 2015-12-04 14:48:49,362 'field_units' was overridden by 'dataset_units/mag_field_x_pert'
yt : [WARNING  ] 2015-12-04 14:48:49,364 'field_units' was overridden by 'dataset_units/mag_field_y_bg'
yt : [WARNING  ] 2015-12-04 14:48:49,366 'field_units' was overridden by 'dataset_units/mag_field_y_pert'
yt : [WARNING  ] 2015-12-04 14:48:49,368 'field_units' was overridden by 'dataset_units/mag_field_z_bg'
yt : [WARNING  ] 2015-12-04 14:48:49,370 'field_units' was overridden by 'dataset_units/mag_field_z_pert'
yt : [WARNING  ] 2015-12-04 14:48:49,377 'field_units' was overridden by 'dataset_units/velocity_x'
yt : [WARNING  ] 2015-12-04 14:48:49,379 'field_units' was overridden by 'dataset_units/velocity_y'
yt : [WARNING  ] 2015-12-04 14:48:49,381 'field_units' was overridden by 'dataset_units/velocity_z'
yt : [INFO     ] 2015-12-04 14:48:49,420 Parameters: current_time              = 1.00339395144
yt : [INFO     ] 2015-12-04 14:48:49,421 Parameters: domain_dimensions         = [128 128 128]
yt : [INFO     ] 2015-12-04 14:48:49,423 Parameters: domain_left_edge          = [  781250.    781250.   3664122.1]
yt : [INFO     ] 2015-12-04 14:48:49,424 Parameters: domain_right_edge         = [  1.99218750e+08   1.99218750e+08   1.58778630e+08]
yt : [INFO     ] 2015-12-04 14:48:49,426 Parameters: cosmological_simulation   = 0.0

In [9]:

    

from tvtk.api import tvtk

#pysac imports
import pysac.yt
import pysac.analysis.tube3D.tvtk_tube_functions as ttf
import pysac.plot.mayavi_plotting_functions as mpf
from pysac.plot.mayavi_seed_streamlines import SeedStreamline
from pysac.plot.divergingcolourmaps import get_mayavi_colourmap
from pysac.analysis.tube3D.process_utils import get_yt_mlab

### Load in and Config ###

# loaded above
ds = pysac.yt.SACGDFDataset('/home/stuart/VivaData/Slog_p240-0_A10_B005_00001.gdf')
tube_r = 60

#if running this creates a persistant window just get it out of the way!
mlab.options.offscreen = True
fig = mlab.figure(bgcolor=(1, 1, 1))

cg = ds.index.grids[0]

#Slices
cube_slice = np.s_[:,:,:-5]
x_slice = np.s_[:,:,:,:-5]

#Define the size of the domain
xmax, ymax, zmax = np.array(cg['density'].to_ndarray()[cube_slice].shape) - 1
domain = {'xmax':xmax, 'ymax':ymax, 'zmax':zmax}

bfield, vfield = get_yt_mlab(ds, cube_slice, flux=False)

#Create a scalar field of the magntiude of the vector field
bmag = mlab.pipeline.extract_vector_norm(bfield, name="Field line Normals")


    

    




yt : [INFO     ] 2015-12-04 14:48:49,672 Loading field plugins.
yt : [INFO     ] 2015-12-04 14:48:49,673 Loaded angular_momentum (8 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,674 Loaded astro (15 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,675 Loaded cosmology (22 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,678 Loaded fluid (63 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,680 Loaded fluid_vector (95 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,682 Loaded geometric (111 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,683 Loaded local (115 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,685 Loaded magnetic_field (120 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,687 Loaded my_plugins (120 new fields)
yt : [INFO     ] 2015-12-04 14:48:49,688 Loaded species (122 new fields)

In [10]:

    

xc = domain['xmax']/2
yc = domain['ymax']/2
ti = 0
n = 100

surf_seeds = []
for theta in np.linspace(0, 2 * np.pi, n, endpoint=False):
    surf_seeds.append([tube_r * np.cos(theta + 0.5 * ti) + xc,
    tube_r * np.sin(theta + 0.5 * ti) + yc, domain['zmax']])

seeds = np.array(surf_seeds)
#Add axes:
axes, outline = mpf.add_axes(np.array(zip(ds.domain_left_edge,ds.domain_right_edge)).flatten()/1e8)

#Add seed points to plot:
seed_points = mlab.points3d(seeds[:,0], seeds[:,1], seeds[:,2],
color=(0.231, 0.298, 0.752), scale_mode='none',
scale_factor=1.5)

mlab_view(fig.scene)


    

    




/opt/miniconda/envs/thesis/lib/python2.7/site-packages/mayavi/tools/camera.py:288: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if focalpoint is not None and not focalpoint == 'auto':






    
Out[10]:

In [11]:

    

field_lines = SeedStreamline(seed_points = np.array(seeds))
bmag.add_child(field_lines)
field_lines.actor.mapper.scalar_visibility = False
field_lines.actor.property.color = (0,0,0)
field_lines.actor.property.line_width = 1.5

mlab_view(fig.scene)


    

    
Out[11]:

In [12]:

    

pd_seeds = ttf.make_circle_seeds(100, 60, **domain)
fieldlines, surface = ttf.create_flux_surface(bfield.outputs[0], pd_seeds)
surface.output.lines = None
flux_surface = mlab.pipeline.surface(surface.output)
flux_surface.actor.mapper.scalar_visibility = False
flux_surface.actor.property.color = (0.8,0.8,0.8)
#flux_surface.actor.property.line_width = 0

mlab_view(fig.scene)


    

    
Out[12]:

In [13]:

    

axes.visible = False
outline.visible = False
flux_surface.actor.property.edge_visibility = True
mlab_view(fig.scene, azimuth = 90, elevation = 75, distance=80, focalpoint=[63, 120, 110], aa=20)


    

    
Out[13]:

In [14]:

    

poly_norms = ttf.make_poly_norms(surface.output)
normvec = mlab.pipeline.glyph(poly_norms.output)
normvec.glyph.glyph_source.glyph_source = normvec.glyph.glyph_source.glyph_dict['arrow_source']
normvec.glyph.glyph.scale_mode = 'data_scaling_off'
normvec.glyph.glyph.color_mode = 'color_by_scale'
normvec.glyph.glyph.scale_factor = 5
normvec.glyph.glyph_source.glyph_position = 'tail'

mlab_view(fig.scene, azimuth=85, elevation=80, distance=50, focalpoint=[63, 120, 110], aa=20)


    

    
Out[14]:

In [33]:

    

import td_plotting_helpers as ph
import time_distance_plots as tdp
from matplotlib.image import NonUniformImage
import texfigure
ch4 = texfigure.Manager(None, number=4, base_path='/home/stuart/GitHub/Thesis/thesis/Chapter4/')

figsize = (17,9)

pvel_labels = {'par_label':r'$V_\parallel[$ ms$^{-1}$]', 
                'perp_label':r'$V_\perp$ [ms$^{-1}$]',
                'phi_label':r'$V_\phi$ [ms$^{-1}$]'}


pflux_labels = {'par_label':r'$F^2_\parallel / F^2$ ', 
                'perp_label':r'$F^2_\perp / F^2$',
                'phi_label':r'$F^2_\phi / F^2$'}

post_amp = "A10"
period = "p240"
tube_r = 'r30'
drivers = ['horiz', 'vert', 'Suni', 'Sarch', 'Slog']
exp_facs = [None, None, 'B0', 'B0005', 'B005']
captions = ['Horizontal', 'Vertical', 'Circular', 'Archemedian Spiral', 'Logarithmic Spiral']

figures = {}

for j, (driver, exp_fac, caption) in enumerate(zip(drivers, exp_facs, captions)):
    
    all_times, y, all_spoints = tdp.get_xy(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)
    data, beta_line = tdp.get_data(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)
    va_line, cs_line = tdp.get_speeds(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)

    Fdata, beta_line, avgs = tdp.get_flux(ch4.data_dir, driver, period, post_amp, tube_r, exp_fac)

    
    fd = lambda args: [a.T for a in args]
       
    ph.xxlim = -150
    
    fig, axes = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=figsize)
    ph.triple_plot(axes[:,0], all_times, y, *fd(data), beta_line=1./beta_line,
                   levels=[1.,3.,5.,7.,10.,100.], manual_locations=False, **pvel_labels)
    
    
    
    ph.triple_plot(axes[:,1], all_times, y, *fd(Fdata), beta_line=1./beta_line,
                   levels=[1.,3.,5.,7.,10.,100.], manual_locations=False, cmap='PRGn', **pflux_labels)
    
    for ax in axes[:,0]:
        ph.add_phase_speeds(ax, all_times, y, va_line, cs_line, dx_scale=1e6, color='g')
    
    for ax in axes[:,1]:
        ph.add_phase_speeds(ax, all_times, y, va_line, cs_line, dx_scale=1e6, color='c')
    
    axes[0,1].legend(bbox_to_anchor=(1.60, 1.05))
        
    #Remove the top two x labels
    axes[0,0].set_xlabel('')
    axes[1,0].set_xlabel('')
    axes[0,1].set_xlabel('')
    axes[1,1].set_xlabel('')
    fig.tight_layout(h_pad=0.05)
    figures[driver] = fig


    

In [34]:

    

figures['vert']


    

    
Out[34]:

In [35]:

    

figures['horiz']


    

    
Out[35]:

In [36]:

    

figures['Slog']


    

    
Out[36]:

In [19]:

    

from flux_comparison import make_flux_bar_chart, get_averages

Favgs = get_averages(ch4.data_dir)

fig = make_flux_bar_chart((12,10), ch4.data_dir)
fig


    

    
Out[19]:

In [38]:

    

from sacconfig import SACConfig
cfg = SACConfig(cfg_file="./python/ch5_config.cfg")

BL = np.array([0.015, 0.05, 0.15, 0.45, 1.5])


    

In [39]:

    

fig, ax = plt.subplots(figsize=(9,2), dpi=600)
ax.plot(BL, np.ones(BL.size), 'x', markersize=10, mew=2)
ax.errorbar([0.15], [1], xerr=np.array([[-1*(0.15-1/(6.4-1.6)), 0.15+1/(6.4+1.6)]]).T, mew=2, elinewidth=2)
ax.semilogx()
ax.get_yaxis().set_visible(False)
ax.set_frame_on(False)
ax.get_xaxis().tick_bottom()
ax.xaxis.set_tick_params(width=2)
ax.xaxis.set_tick_params(width=2, which='minor')
ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax.xaxis.set_ticks(BL)
xmin, xmax, ymin, ymax = ax.axis()
ax.add_artist(plt.Line2D((xmin, xmax), (ymin, ymin), color='black', linewidth=1.4))
l = ax.set_xlim([0.01, 2.0])
l = ax.set_xlabel(r'$B_L$', fontsize=18)

fig.tight_layout(h_pad=0.01)
fig.savefig('./images/ch5_space.png', transparent=True)


    

In [40]:

    

beta = False
def add_triple_phase(ax, tube_r, color='g'):
    ps = ph.get_phase_speeds(cfg, tube_r)
    for ax0 in ax:
        ph.add_phase_speeds(ax0, color=color, **ps)

bl_figures = {}
for j, bl in enumerate(BL):
    cfg.exp_fac = bl
    
    fig, ax = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=figsize)
    
    kwargs = ph.get_single_velocity(cfg, 'r30', beta=beta)
    kwargs.update(pvel_labels)
    
    ph.triple_plot(ax[:,0], **kwargs)
    
    kwargs = ph.get_single_percentage_flux(cfg, 'r30', beta=beta)
    kwargs.update(pflux_labels)
    kwargs.update({'cmap': 'PRGn'})
    
    ph.triple_plot(ax[:,1], **kwargs)
    
    #Remove the top two x labels
    ax[0,0].set_xlabel('')
    ax[1,0].set_xlabel('')
    add_triple_phase(ax[:,0], 'r30')
    ax[0,1].set_xlabel('')
    ax[1,1].set_xlabel('')
    add_triple_phase(ax[:,1], 'r30', color='c')
    #add_all_bpert(ax, 'r30')
    fig.tight_layout(h_pad=0.05)
    ax[0,1].legend(bbox_to_anchor=(1.60, 1.05))
    bl_figures['B{}'.format(str(bl).replace('.',''))] = fig


    

In [41]:

    

bl_figures['B0015']


    

    
Out[41]:

In [42]:

    

bl_figures['B015']


    

    
Out[42]:

In [43]:

    

bl_figures['B15']


    

    
Out[43]:

In [44]:

    

import os
int_periods = np.floor(600./cfg.period)*180

def calc_avgs(tube_r):
    AvgsP = np.zeros([3,len(BL)])
    for i, bl in enumerate(BL):
        cfg.exp_fac = bl
        
        times = np.load(os.path.join(cfg.data_dir, 'Times_{}.npy'.format(cfg.get_identifier())))
        max_index = np.argmin(np.abs(int_periods - times))
        
        Fpar, Fperp, Fphi = map(np.load, ph.glob_files(cfg, tube_r, 'LineFlux*Fp*npy'))
        #Fpar, Fperp, Fphi = map(np.load, ph.glob_files(cfg, tube_r, '*vp*npy'))
        Fpar[np.abs(Fpar)<1e-5], Fperp[np.abs(Fperp)<1e-5], Fphi[np.abs(Fphi)<1e-5] = 0., 0., 0.
        Fpar, Fperp, Fphi = Fpar[:max_index,:], Fperp[:max_index,:], Fphi[:max_index,:]
        
        Ftot2 = (Fpar**2 + Fperp**2 + Fphi**2)
        Fpar2, Fperp2, Fphi2 = np.array([Fpar, Fperp, Fphi])**2
        FparP, FperpP, FphiP = (Fpar2/Ftot2)*100, (Fperp2/Ftot2)*100, (Fphi2/Ftot2)*100
        
        FparP = np.mean(np.ma.masked_array(FparP, np.isnan(FparP)))
        FperpP = np.mean(np.ma.masked_array(FperpP, np.isnan(FperpP)))
        FphiP = np.mean(np.ma.masked_array(FphiP, np.isnan(FphiP)))
        
        AvgsP[:, i] = FparP, FperpP, FphiP
    return AvgsP


fig, axs = plt.subplots(nrows=3, figsize=(11,11), sharex=True)
titles = ["Flux Surface Radius $=156$ km", "Flux Surface Radius $=468$ km", "Flux Surface Radius $=936$ km"]
tubes = []
for i, ax in enumerate(axs):
    AvgsP = calc_avgs(cfg.tube_radii[i])
    tubes.append(AvgsP)
    ax.semilogx()
    ax.plot(BL, AvgsP[0], 'o', label=r"$F_\parallel^2$", mew=0, ms=10)
    ax.plot(BL, AvgsP[1], '_', label=r"$F_\perp^2$", mew=2, ms=10)
    ax.plot(BL, AvgsP[2], 'x', label=r"$F_\phi^2$", mew=2, ms=10)
    ax.set_ylabel("% Square Wave Flux")
    ax.set_title(titles[i])
    ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    ax.xaxis.set_ticks(BL)
    ax.set_xlim([0.01, 2.01])
    ax.set_ylim([0, 85])
    err = np.array([-1*(0.15-1/(6.4-1.6)), 0.15+1/(6.4+1.6)])
    ax.fill_betweenx(np.linspace(-5,105), err[0], err[1], alpha=0.3, color='green', linewidth=0)

axs[0].legend(bbox_to_anchor=(1.20, 1.05))
axs[-1].set_xlabel("Logarithmic Spiral Expansion Factor ($B_L$)")

plt.tight_layout()
fig


    

    
Out[44]:

In [46]:

    

from period_amps import periods, str_amps, sim_params
all_periods = sim_params[:10]
periods = periods[:10]

cfg = SACConfig(cfg_file='./python/ch6_config.cfg')
cfg.data_dir = '/home/stuart/GitHub/Thesis/thesis/Chapter6/Data/'


    

In [47]:

    

ph.xxlim = 600
p_figures= {}
for i, paf in enumerate(all_periods):
    [setattr(cfg, f, getattr(paf, f)) for f in paf._fields]
    
    fig, ax = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=figsize)
    
    kwargs = ph.get_single_velocity(cfg, 'r30', beta=beta)
    kwargs.update(pvel_labels)
    
    ph.triple_plot(ax[:,0], **kwargs)
    
    kwargs = ph.get_single_percentage_flux(cfg, 'r30', beta=beta)
    kwargs.update(pflux_labels)
    kwargs.update({'cmap': 'PRGn'})
    
    ph.triple_plot(ax[:,1], **kwargs)
    
    #Remove the top two x labels
    ax[0,0].set_xlabel('')
    ax[1,0].set_xlabel('')
    add_triple_phase(ax[:,0], 'r30')
    ax[0,1].set_xlabel('')
    ax[1,1].set_xlabel('')
    add_triple_phase(ax[:,1], 'r30', color='c')
    #add_all_bpert(ax, 'r30')
    fig.tight_layout(h_pad=0.05)
    ax[0,1].legend(bbox_to_anchor=(1.6, 1.05))
    p_figures['{}'.format(str(paf.period))] = fig


    

In [48]:

    

p_figures['30.0']


    

    
Out[48]:

In [49]:

    

p_figures['150.0']


    

    
Out[49]:

In [50]:

    

p_figures['300.0']


    

    
Out[50]:

In [51]:

    

from period_amps import periods, sim_params
sim_params = sim_params[:10]
periods = np.array(periods[:10])

fig, axs = plt.subplots(nrows=3, figsize=(11,11), sharex=True)
titles = ["Flux Surface Radius $=156$ km", "Flux Surface Radius $=468$ km", "Flux Surface Radius $=936$ km"]
tubes = []
for i, ax in enumerate(axs):
    AvgsP = ph.get_all_avgs(cfg, cfg.tube_radii[i], sim_params)
    tubes.append(AvgsP)
    ax.plot(periods, AvgsP[0], 'o', label=r"$F_\parallel^2$", mew=0, ms=12)
    ax.plot(periods, AvgsP[1], '_', label=r"$F_\perp^2$", mew=2, ms=12)
    ax.plot(periods, AvgsP[2], 'x', label=r"$F_\phi^2$", mew=2, ms=12)
    ax.set_ylabel("% Square Wave Flux")
    ax.set_title(titles[i])
    ax.xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    ax.xaxis.set_ticks(periods)
    ax.set_xticklabels(["{}\n[{:n}]".format(p, 600//p) for p in periods])
    ax.set_ylim([10, 75])
    ax.set_xlim([25, 305])

axs[-1].set_xlabel("Period [s] \n [Number of periods averaged]")
#axs[0].legend(loc=2)

axs[0].legend(bbox_to_anchor=(1.21, 1.05))
plt.tight_layout(h_pad=0.1)
fig


    

    
Out[51]: