

In [2]:

    
import numpy as np
#import pyPLUTO as pp
from astropy.io import ascii
import os
import sys
from ipywidgets import interactive, widgets,fixed
from IPython.display import Audio, display
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.ticker as ticker
from matplotlib.animation import FuncAnimation,FFMpegWriter
from matplotlib import rc,rcParams
from scipy.integrate import quad
rc('text', usetex=True)
rcParams['figure.figsize'] = (15., 6.0)
rcParams['ytick.labelsize'],rcParams['xtick.labelsize'] = 17.,17.
rcParams['axes.labelsize']=19.
rcParams['legend.fontsize']=17.
rcParams['text.latex.preamble'] = ['\\usepackage{siunitx}']
import seaborn
seaborn.despine()
seaborn.set_style('white', {'axes.linewidth': 0.5, 'axes.edgecolor':'black'})
seaborn.despine(left=True)
%load_ext autoreload









    



The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload






    





<matplotlib.figure.Figure at 0x7fe16c358a50>



In [9]:

    
%autoreload 1



In [10]:

    
%aimport f



In [5]:

    
def quadruple(d,VAR,tdk='Myrs',Save_Figure='',cl='',nn=0,mspeed='km',rows=2,cols=2,xlim=[None,None],
              ylim=[None,None],tlim=None,VARlim=[None,None],datafolder='../Document/DataImages/'):
    """
    Plot a rows(=2) x cols(=2) Variable 
    """
    X,Y=d['X'],d['Y']
    Vx=d['Vx'] if nn>0 else 0
    Vy=d['Vy'] if nn>0 else 0
    T=np.linspace(0,d['T'].shape[0]-1,rows*cols,dtype=int) if tlim==None else np.linspace(0,tlim,rows*cols,dtype=int)
    fig, axes = plt.subplots(nrows=rows, ncols=cols, sharex=True, sharey=True,
                            figsize=(cols*5,rows*5))
    i=0
    td=1e3 if tdk=='kyrs' else 1e6
    vmin=VARlim[0] if VARlim[0]==None else VAR.min()
    vmax=VARlim[1] if VARlim[1]==None else VAR.max()
    
    for ax in axes.flat:
        ext=[X.min(),X.max(),Y.min(),Y.max()]
        ax.get_yaxis().get_major_formatter().set_useOffset(False)
        #ax.add_artist(plt.Circle((0, 0), 1.0, color='r',fill=False,linestyle='--'))
        label = '{:.1f} {}'.format(d['T'][T[i]]/td,tdk)
        ax.set_title(label,fontsize=20)
        ax.grid(False)
        pc = ax.imshow(VAR[:,:,T[i]].T,cmap='viridis',origin='lower',aspect='equal',
                       extent=ext,vmin=vmin,vmax=vmax)
        if nn>0:
            k=nn #distance from boundaries for first/last arrows
            sc=2. if mspeed =='max' else 5. if mspeed == 'c' else 1e-4
            q=pc.axes.quiver(X[k:-k:nn],Y[k:-k:nn],
                            Vx[:,:,T[i]][k:-k:nn,k:-k:nn].T,
                            Vy[:,:,T[i]][k:-k:nn,k:-k:nn].T,
                             scale=sc,alpha=0.5,width=0.002)
            if mspeed == 'c':
                pc.axes.quiverkey(q,0.05,1.02,1.,r'$1\si{c}$',labelpos='E',fontproperties={'weight': 'bold'})
            elif mspeed == 'max':
                mV=np.max(np.sqrt(Vx[np.argmin((d['Y']-ylim[0])**2):np.argmin((d['Y']-ylim[1])**2),
                                     np.argmin((d['X']-xlim[0])**2):np.argmin((d['X']-xlim[1])**2),T[i]]**2+
                                  Vy[np.argmin((d['Y']-ylim[0])**2):np.argmin((d['Y']-ylim[1])**2),
                                     np.argmin((d['X']-xlim[0])**2):np.argmin((d['X']-xlim[1])**2),T[i]]**2))
                pc.axes.quiverkey(q,0.05,1.02,mV,'{:.2f} c'.format(mV),labelpos='E',
                                  fontproperties={'weight': 'bold'})
            else:
                pc.axes.quiverkey(q,0.02,1.02,3.36e-6,r'$1\si{km.s^{-1}}$',labelpos='E',fontproperties={'weight': 'bold'})
            
        i=i+1
    ax.set_xlim(xlim[0],xlim[1])
    ax.set_ylim(ylim[0],ylim[1])
    plt.tight_layout()
    cbar_ax = fig.add_axes([0., 1.015, 1., 0.025*(np.float(cols)/rows)])#*(np.float(cols)/rows)
    cb=fig.colorbar(pc, cax=cbar_ax,orientation="horizontal",label=cl)
    cb.ax.tick_params(labelsize=17)
    cb.ax.xaxis.offsetText.set(size=20)
    cb.ax.xaxis.set_ticks_position('top')
    cb.ax.xaxis.set_label_position('top')
    form='.png'
    if Save_Figure <> '': plt.savefig(datafolder+Save_Figure+form,bbox_inches='tight',format='png', dpi=100)



In [3]:

    
#JC=np.load('../Data/jet-3c.npz')
#aJC=np.load('../Data/afterJet.npz')
hJC=np.load('../Data/jet-3c-n.npz')
sJC=np.load('../Data/jet-sneq.npz')
nJC=np.load('../Data/jet-3c-nn.npz')
tJC=np.load('../Data/jet-3c-tab.npz')
atJC=np.load('../Data/afterjet-tab.npz')

No Cooling Jet



In [21]:

    
plt.figure(figsize=(10,10))
plt.imshow(np.log10(nJC['RHO'][:,:,0].T),origin='low',extent=[-150,150,-150,150])
Save_Figure='Jet0'
datafolder='../Document/DataImages/'
form='.png'
#plt.savefig(datafolder+Save_Figure+form,bbox_inches='tight',format='png', dpi=100)



In [50]:

    
f.quadruple(nJC,np.sqrt(nJC['Vx']**2+nJC['Vy']**2),tdk='kyrs',tlim=16,xlim=[-4,4],ylim=[-14,-8],Save_Figure='VnoCool')



In [49]:

    
f.quadruple(nJC,np.log10(nJC['RHO']),tdk='kyrs',tlim=16,xlim=[-4,4],ylim=[-14,-8],Save_Figure='RHOnoCool')



In [79]:

    
f.quadruple(nJC,np.log10(nJC['RHO']),tlim=200,tdk='kyrs')#,tlim=8,xlim=[-2,2],ylim=[-14,-12],VARlim=[0.,10.],Save_Figure='RHOnoCool0-10')



In [33]:

    
f.quadruple(nJC,np.sqrt(nJC['Vx']**2+nJC['Vy']**2),tdk='kyrs',tlim=8,xlim=[-2,2],ylim=[-14,-12])

Tabulated Cooling



In [65]:

    
f.quadruple(tJC,np.sqrt(tJC['Vx']**2+tJC['Vy']**2),tdk='kyrs',tlim=100)



In [80]:

    
f.quadruple(nJC,np.log10(nJC['RHO']),tdk='kyrs',tlim=100)



In [63]:

    
f.quadruple(tJC,np.log10(tJC['RHO']),tdk='kyrs',tlim=100)



In [17]:

    
def radial_profile(data, center):
    y, x = np.indices((data.shape))
    r = np.sqrt((x - center[0])**2 + (y - center[1])**2)
    r = r.astype(np.int)

    tbin = np.bincount(r.ravel(), data.ravel())
    nr = np.bincount(r.ravel())
    radialprofile = tbin / nr
    return radialprofile



In [20]:

    
V=np.sqrt(tJC['Vx'][:,:,-1]**2+tJC['Vy'][:,:,-1]**2)



In [86]:

    
fV=np.fft.fft2(V)
fsV=np.abs(np.fft.fftshift(fV))
kx=np.fft.fftshift(np.fft.fftfreq(fsV.shape[0],10))
ky=np.fft.fftshift(np.fft.fftfreq(fsV.shape[1],10))



In [87]:

    
plt.imshow(np.log10(fsV**2),extent=[kx[0],kx[-1],ky[0],ky[-1]])









    Out[87]:





<matplotlib.image.AxesImage at 0x7fd3cd9d1d50>



In [95]:

    
k=kx[256:]**2+ky[256:]**2
radprof=radial_profile(fsV[:256,:256]**2,[256,256])[:256]
plt.loglog(k,radprof)









    



/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in divide
  






    Out[95]:





[<matplotlib.lines.Line2D at 0x7fd3cd691050>]



In [93]:

    
from scipy.optimize import curve_fit



In [94]:

    
def fit(x,A,B): return A-B*x



In [100]:

    
y=np.log10(radprof[1:])
x=np.log10(k[1:])
s,sd2=curve_fit(fit,x,y,[1e8,5./3])



In [106]:

    
plt.loglog(k,radprof)
plt.loglog(k,10**s[0]*k**-s[1],label='{:.2f} Cascade'.format(-s[1]))
plt.legend()









    



/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in power
  






    Out[106]:





<matplotlib.legend.Legend at 0x7fd3cd4e3150>



In [ ]:

    
mm=max(JC['T'].shape[0],aJC['T'].shape[0])
TT=np.arange(mm)
times=JC['T'] if JC['T'].shape[0]>aJC['T'].shape[0] else aJC['T']



In [69]:

    
def get_gas(RHO,aV,Temp,TT):
    #return np.array([RHO[:,:,t][np.logical_and(np.logical_and(aV[:,:,t]>6.7e-4,aV[:,:,t]<2.3e-3),Temp[:,:,t]<1e7)].sum()/RHO[:,:,t].sum() for t in TT])
    return np.array([RHO[:,:,t][aV[:,:,t]>6.7e-4].sum()/RHO[:,:,t].sum() for t in TT])
    #return np.array([RHO[:,:,t][aV[:,:,t]>2e-3].sum()/RHO[:,:,t].sum() for t in TT])



In [70]:

    
%%time
Mftj7=get_gas(tJC['RHO'],np.sqrt(tJC['Vx']**2+tJC['Vy']**2),tJC['PRS']*f.Temp0/tJC['RHO'],np.arange(tJC['T'].shape[0]))
Mfnj7=get_gas(nJC['RHO'],np.sqrt(nJC['Vx']**2+nJC['Vy']**2),nJC['PRS']*f.Temp0/nJC['RHO'],np.arange(nJC['T'].shape[0]))









    



CPU times: user 1min 2s, sys: 5.99 s, total: 1min 8s
Wall time: 1min 28s



In [97]:

    
ttlim=150
plt.plot(tJC['T'][:tlim]/1e3,Mftj[:ttlim],label='Escaped Gas ($V>\SI{600}{km.s^{-1}}$)')
#plt.plot(tJC['T'][:tlim]/1e3,Mfnj[:tlim],label='No Cooling')
plt.plot(tJC['T'][:tlim]/1e3,Mftj7[:ttlim],label='Gas in Wind ($V>\SI{200}{km.s^{-1}}$)')
#plt.plot(tJC['T'][:tlim]/1e3,Mfnj[:tlim],label='No Cooling')
plt.xlabel('Time (kyrs)')
plt.legend()
plt.savefig('/home/astromix/astro/MasterThesis/Document/DataImages/RatioEscapedGas.png',bbox_inches='tight')
f.quadruple(tJC,np.log10(tJC['RHO']),rows=3,tdk='kyrs',tlim=ttlim,Save_Figure='JetCloudRHO')
f.quadruple(tJC,np.sqrt(tJC['Vx']**2+tJC['Vy']**2),rows=3,tdk='kyrs',tlim=ttlim,Save_Figure='JetCloudV')



In [93]:

    
f.quadruple(tJC,np.sqrt(tJC['Vx']**2+tJC['Vy']**2),tdk='kyrs',tlim=8,xlim=[-2,2],ylim=[-14,-12],Save_Figure='JetISMV')
f.quadruple(tJC,np.log10(tJC['RHO']),tdk='kyrs',tlim=8,xlim=[-2,2],ylim=[-14,-12],Save_Figure='JetISMRHO')



In [95]:

    
f.quadruple(tJC,np.sqrt(tJC['Vx']**2+tJC['Vy']**2),tdk='kyrs',tlim=30,xlim=[-2,2],ylim=[-11,0])#,Save_Figure='JetISMV')
f.quadruple(tJC,np.log10(tJC['RHO']),tdk='kyrs',tlim=30,xlim=[-2,2],ylim=[-11,0])#,Save_Figure='JetISMRHO')



In [ ]:

    
p=85
plt.plot(times[TT][0:min(mm,Mfj.shape[0])]/1e3,Mfj)
plt.plot(times[TT][p+1:min(mm,Mfaj.shape[0]+p)]/1e3,Mfaj[1:mm-p])



In [ ]:

    
Temp=JC['PRS']*f.Temp0/JC['RHO']
Tempm=np.ma.masked_where(JC['RHO']<10,Temp)
plt.imshow(Tempm[:,:,-1].T,cmap='viridis',origin='lower')
#plt.contour(np.log10(JC['RHO'][:,:,-1]).T,origin='lower',levels=[1.,2.,3.])
plt.colorbar()



In [ ]:

    
aTemp=aJC['PRS']*f.Temp0/aJC['RHO']
aTempm=np.ma.masked_where(aJC['RHO']<10,aTemp)
plt.imshow(aTempm[:,:,-1].T,cmap='viridis',origin='lower')
#plt.contour(np.log10(JC['RHO'][:,:,-1]).T,origin='lower',levels=[1.,2.,3.])
plt.colorbar()



In [ ]:

    
plt.plot(times[TT][0:min(mm,Mfj.shape[0])]/1e3,Tempm.mean(axis=(0,1)))
plt.plot(times[TT][p+1:min(mm,Mfaj.shape[0]+p)]/1e3,aTempm.mean(axis=(0,1))[1:mm-p])



In [ ]:

    
plt.plot(times[TT][0:min(mm,Mfj.shape[0])]/1e3,Temp.mean(axis=(0,1)))
plt.plot(times[TT][p+1:min(mm,Mfaj.shape[0]+p)]/1e3,aTemp.mean(axis=(0,1))[1:mm-p])
plt.yscale('log')



In [ ]:

    
plt.plot(times[TT][0:min(mm,Mfj.shape[0])]/1e3,(np.sqrt(JC['Vx']**2+JC['Vy']**2)).mean(axis=(0,1)))
plt.plot(times[TT][p+1:min(mm,Mfaj.shape[0]+p)]/1e3,(np.sqrt(aJC['Vx']**2+aJC['Vy']**2)).mean(axis=(0,1))[1:mm-p])



In [ ]:

    
plt.plot(times[TT][0:min(mm,Mfj.shape[0])]/1e3,(np.sqrt(JC['Vx']**2+JC['Vy']**2)).std(axis=(0,1)))
plt.plot(times[TT][p+1:min(mm,Mfaj.shape[0]+p)]/1e3,(np.sqrt(aJC['Vx']**2+aJC['Vy']**2)).std(axis=(0,1))[1:mm-p])
#plt.yscale('log')

TT=np.linspace(0,JC['RHO'].shape[2]-1,6,dtype=int) fig, axes = plt.subplots(nrows=1,ncols=len(TT)) for i,t in enumerate(TT): axes[i].imshow(np.log10(JC['RHO'][:,:,t].T),origin='lower',extent=[-10,10,-10,10]) axes[i].set_xlim(-2,2) axes[i].set_ylim(-10,-3) axes[i].vlines(-1,-10,10,linewidth=0.5)

for t in TT: plt.plot(NC['Y']*10.+100.,NC['RHO'][256-26,:,t],label='{:.2f} kyrs'.format(NC['T'][t]/1e3)) plt.yscale('log') plt.xlabel('Y $(\si{pc})$') plt.ylabel('n $(\si{cm^{-3}})$') plt.xlim(10,50.) plt.ylim(1e-1,None) plt.legend()



In [46]:

    
%%time
FF=tJC
vfile='TabulatedV'#'test'
step=1

VAR=np.sqrt(FF['Vx']**2+FF['Vy']**2)####
#VAR=np.log10(FF['RHO'])
T = FF['T'][::step]
fig=plt.figure(figsize=(10,10))
fig.set_tight_layout(True)
ax1=plt.subplot()
ax1.get_yaxis().get_major_formatter().set_useOffset(False)
ext=[FF['X'].min(),FF['X'].max(),FF['Y'].min(),FF['Y'].max()]

pc = ax1.imshow(VAR[:,:,0].T,cmap='viridis',origin='lower',aspect='equal',extent=ext,
                            vmin=VAR.min(axis=2).min(),vmax=VAR.max(axis=2).max())
def update(i):
    ax1.cla()
    pc = ax1.imshow(VAR[:,:,i].T,cmap='viridis',origin='lower',aspect='equal',extent=ext,
                    vmin=VAR.min(axis=2).min(),vmax=VAR.max(axis=2).max())
    label = 'Time = {0:.1f} kyrs'.format(T[i]/1000.)
    ax1.set_title(label)
    return ax1
anim = FuncAnimation(fig, update, frames=range(T.shape[0]), interval=150)
anim.save(vfile+'.gif',writer='imagemagic')









    



CPU times: user 1min 44s, sys: 1min 42s, total: 3min 27s
Wall time: 2min 30s



In [ ]: