In [ ]:

    

import numpy as np
import os
import matplotlib
matplotlib.rcParams['figure.dpi'] = 150
import matplotlib.pyplot as plt
from tools import setup_cl


    

In [ ]:

    

def get_dat(dirname):
    datfnames = ['MHD_Jet.dat', 'MHD_Jet_10Myr.dat', 'MHD_Jet_nojiggle.dat']
    for datfname in datfnames:
        if os.path.isfile(os.path.join(dirname, datfname)):
            data = np.genfromtxt(os.path.join(dirname, datfname))
            break
    Emags = data[:,8]
    Ekins = data[:,6]
    mask = np.logical_and(Emags < 7E58, Ekins < 2E59)

    tt = data[mask,0]
    Masses = data[mask,1]
    Etots = data[mask,5]
    Ekins = data[mask,6]
    Eints = data[mask,7] 
    Emags = Emags[mask]

    
    Ptots = (Etots[1:]-Etots[:-1])/(tt[1:]-tt[:-1])
    Pkins = (Ekins[1:]-Ekins[:-1])/(tt[1:]-tt[:-1])
    Pints = (Eints[1:]-Eints[:-1])/(tt[1:]-tt[:-1])
    
    #print max(Ptots), max(Pkins), max(Pints)
    return tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints


    

In [ ]:

    

dirnames = ['/home/ychen/data/0only_0605_hinf',
            '/home/ychen/data/0only_0529_h1',
            '/home/ychen/data/0only_0605_h0',
            
            '/home/ychen/data/0only_0204_hinf_10Myr',
            '/home/ychen/data/0only_1022_h1_10Myr',
            '/home/ychen/data/0only_0204_h0_10Myr',
            ]

table = []
for dirname in dirnames:
    table.append(get_dat(dirname))


    

In [ ]:

    

colors, labels = setup_cl(dirnames)
plt.figure(figsize=(6,3))
for dir, data in zip(dirnames, table):
    tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = get_dat(dir)
    if '10Myr' not in dir:
        plt.plot(tt[:]/3.154E13, Ekins, c=colors[dir], label=labels[dir], lw=1, alpha=1)
    else:
        plt.plot(tt[:]/3.154E13, Ekins, c=colors[dir], lw=1, alpha=1)
plt.ylim(0, 1.5E59)
#plt.semilogx()
plt.xlim(0, 100)
plt.legend()
plt.xlabel('time (Myr)')
plt.ylabel('Total Kinetic Energy')


    

In [ ]:

    

colors, labels = setup_cl(dirnames)
plt.figure(figsize=(6,3))
for dir, data in zip(dirnames, table):
    tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = get_dat(dir)
    if '10Myr' not in dir:
        Eint0 = Eints[0]
        plt.plot(tt[:]/3.154E13, Eints-Eint0, c=colors[dir], label=labels[dir], lw=1, alpha=0.7)
    else:
        plt.plot(tt[:]/3.154E13, Eints-Eint0, c=colors[dir], lw=1, alpha=0.7)
plt.ylim(0, 6E59)
#plt.semilogx()
plt.xlim(0, 100)
plt.legend()
plt.xlabel('time (Myr)')
plt.ylabel('Total Internal Energy')


    

In [ ]:

    

colors, labels = setup_cl(dirnames)
plt.figure(figsize=(6,3))
for dir, data in zip(dirnames, table):
    tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = data
    if '10Myr' not in dir:
        plt.plot(tt[:]/3.154E13, Emags, c=colors[dir], label=labels[dir], lw=1, markersize=0, alpha=1)
    else:
        plt.plot(tt[:]/3.154E13, Emags, c=colors[dir], lw=1, markersize=0, alpha=1)
plt.ylim(0, 4E57)
#plt.semilogx()
plt.xlim(0, 20)
plt.legend()
plt.xlabel('time (Myr)')
plt.ylabel('Total Magnetic Energy')


    

In [ ]:

    

dirnames = [
    '/home/ychen/data/0only_0525_hinf_nojiggle/',
    '/home/ychen/data/0only_0314_h1_nojiggle/',
    '/home/ychen/data/0only_0330_h0_nojiggle/',
]
table = []
for dirname in dirnames:
    table.append(get_dat(dirname))


    

In [ ]:

    

colors, labels = setup_cl(dirnames)
plt.figure(figsize=(6,3))
for dir, data in zip(dirnames, table):
    tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = data
    if '10Myr' not in dir:
        plt.plot(tt[:]/3.154E13, Emags, c=colors[dir], label=labels[dir], lw=0, marker='.', markersize=1, alpha=1)
    else:
        plt.plot(tt[:]/3.154E13, Emags, c=colors[dir], lw=0, marker='.', markersize=1, alpha=1)
plt.ylim(0, 2E57)
#plt.semilogx()
plt.xlim(0, 5)
plt.legend()
plt.xlabel('time (Myr)')
plt.ylabel('Total Magnetic Energy')


    

In [ ]:

    

colors, labels = setup_cl(dirnames)
plt.figure(figsize=(6,3))
for dir, data in zip(dirnames, table):
    tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = data
    if '10Myr' not in dir:
        plt.plot(tt[:]/3.154E13, Ekins, c=colors[dir], label=labels[dir], lw=0, marker='.', markersize=1, alpha=0.7)
    else:
        plt.plot(tt[:]/3.154E13, Ekins, c=colors[dir], lw=0, marker='.', markersize=1, alpha=0.7)
plt.ylim(0, 1E59)
#plt.semilogx()
plt.xlim(0, 8)
plt.legend()
plt.xlabel('time (Myr)')
plt.ylabel('Total Kinetic Energy')


    

In [ ]:

    

colors, labels = setup_cl(dirnames)
plt.figure(figsize=(6,3))
for dir, data in zip(dirnames, table):
    tt, Masses, Etots, Ekins, Eints, Emags, Ptots, Pkins, Pints = data
    if '10Myr' not in dir:
        plt.plot(tt[:]/3.154E13, Eints-Eints[0], c=colors[dir], label=labels[dir], lw=0, marker='.', markersize=1, alpha=0.7)
    else:
        plt.plot(tt[:]/3.154E13, Eints-Eints[0], c=colors[dir], lw=0, marker='.', markersize=1, alpha=0.7)
#plt.ylim(0, 1E59)
#plt.semilogx()
plt.xlim(0, 8)
plt.legend()
plt.xlabel('time (Myr)')
plt.ylabel('Total Internal Energy')


    

In [ ]: