notebook.community

Edit and run 





In [1]:



    


import sys
import pickle
import numpy
from galpy.util import bovy_plot
from matplotlib import cm
import define_rcsample

Tri-broken exp flare





In [2]:



    


with open('../mapfits/tribrokenexpflare.sav','rb') as savefile:
    bf= numpy.array(pickle.load(savefile))
    samples= numpy.array(pickle.load(savefile))
    bf_g15= numpy.array(pickle.load(savefile))
    samples_g15= numpy.array(pickle.load(savefile))
    bf_zero= numpy.array(pickle.load(savefile))
    samples_zero= numpy.array(pickle.load(savefile))
maps= define_rcsample.MAPs()



    











In [5]:



    


plotthis= numpy.zeros(len(bf))+numpy.nan
plotthise= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(numpy.exp(samples[ii,3]))
    terr= 1.4826*numpy.nanmedian(numpy.fabs(numpy.exp(samples[ii,3])-tmed))
    if tmed < 5.:
        tmed= 0.
        terr= numpy.nan
    plotthis[ii]= tmed
    plotthise[ii]= terr
maps.plot(plotthis,vmin=5.,vmax=13.,minnstar=15,zlabel=r'$R_{\mathrm{peak}}\,(\mathrm{kpc})$')
bovy_plot.bovy_text(-0.55,0.225,r'$\mathrm{single}$',
                         size=16.,color='w')
bovy_plot.bovy_text(-0.495,0.19,r'$\mathrm{exponential}$',
                     size=16.,color='w')
maps.plot(plotthise,vmin=0.,vmax=1.,minnstar=15,zlabel=r'$R_{\mathrm{peak}}\,(\mathrm{kpc})$')
print numpy.nanmedian(plotthise)



    


















    






0.310307233657










    

















    
























In [4]:



    


plotthisf= numpy.zeros(len(bf))+numpy.nan
plotthisfe= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(samples[ii,4])
    terr= 1.4826*numpy.nanmedian(numpy.fabs(samples[ii,4])-tmed)
    if terr > 0.6: tmed= numpy.nan
    plotthisf[ii]= tmed
    plotthisfe[ii]= terr
#plotthisf[plotthis == 0.]= numpy.nan
maps.plot(plotthisf,vmin=-0.3,vmax=0.1,minnstar=15)
maps.plot(plotthisfe,vmin=0.,vmax=.6,minnstar=15)
print numpy.nanmedian(plotthisf[plotthis > 0.]), 1.4826*numpy.nanmedian(numpy.fabs(plotthisf[plotthis > 0.]-numpy.nanmedian(plotthisf[plotthis > 0.])))/numpy.sqrt(numpy.sum(plotthis > 0.))
print numpy.nanmedian(plotthisf[plotthis == 0.]), 1.4826*numpy.nanmedian(numpy.fabs(plotthisf[plotthis == 0.]-numpy.nanmedian(plotthisf[plotthis == 0.])))/numpy.sqrt(numpy.sum(plotthis == 0.))



    


















    






-0.0986811478888 0.00950197448604
-0.0065494242185 0.0179334389286










    

















    
























In [5]:



    


overplot=False
for ii, map in enumerate(maps.map()):
    if plotthis[ii] == 0.: continue
    tmed= numpy.nanmedian(samples[ii,4])
    terr= 1.4826*numpy.nanmedian(numpy.fabs(samples[ii,4])-tmed)
    if terr > 0.6: continue
    bovy_plot.bovy_hist(samples[ii,4],bins=51,range=[-0.3,0.1],histtype='step',
                        color=cm.coolwarm(((numpy.nanmedian(map['ALPHAFE'])+0.0)/0.1)),
                        overplot=overplot)
    overplot=True



    


















    
























In [6]:



    


overplot=False
for ii, map in enumerate(maps.map()):
    if plotthis[ii] > 0.: continue
    tmed= numpy.nanmedian(samples[ii,4])
    terr= 1.4826*numpy.nanmedian(numpy.fabs(samples[ii,4])-tmed)
    if terr > 0.3: continue
    bovy_plot.bovy_hist(samples[ii,4],bins=51,range=[-0.3,0.1],histtype='step',
                        color=cm.coolwarm(((numpy.nanmedian(map['ALPHAFE'])-0.05)/0.2)),
                        overplot=overplot)
    overplot=True



    


















    
























In [7]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    trp= numpy.nanmedian(numpy.exp(samples[ii,3]))
    terr= numpy.std(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]) < 0.03)])
    if trp < 5. and terr > 0.1:
        tmed= numpy.median(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]) < 0.03)])
        if tmed > 1.: 
            tmed= sorted(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]) < 0.03)])
            tmed= tmed[int(round(0.05*len(tmed)))]
    elif terr > 0.1:
        tmed= numpy.median(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]+0.1) < 0.03)])
    else:
        tmed= numpy.median(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))])
    #if numpy.sum((True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]+0.1) < 0.03)) < 1000: tmed= numpy.nan
    plotthisz[ii]= tmed
maps.plot(plotthisz*1000.,
          vmin=200.,vmax=1000.,
          minnstar=15,
          zlabel=r'$h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    






/Library/Python/2.7/site-packages/numpy/core/_methods.py:83: RuntimeWarning: Degrees of freedom <= 0 for slice
  warnings.warn("Degrees of freedom <= 0 for slice", RuntimeWarning)











    
Out[7]:








<matplotlib.image.AxesImage at 0x1089d5850>










    
























In [8]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    trp= numpy.nanmedian(numpy.exp(samples[ii,3]))
    terr= numpy.std(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]) < 0.03)])
    if trp < 5. and terr > 0.1:
        tmed= numpy.median(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]) < 0.03)])
        if tmed > 1.: 
            tmed= sorted(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]) < 0.03)])
            tmed= tmed[int(round(0.05*len(tmed)))]
    elif terr > 0.1:
        tmed= numpy.median(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]+0.1) < 0.03)])
    else:
        tmed= numpy.median(1./samples[ii,1,(True-numpy.isnan(1./samples[ii,1]))])
    #if numpy.sum((True-numpy.isnan(1./samples[ii,1]))*(numpy.fabs(samples[ii,4]+0.1) < 0.03)) < 1000: tmed= numpy.nan
    plotthisz[ii]= tmed
maps.plot(1./plotthisz,
          vmin=1.,vmax=4.,
          minnstar=15,
          zlabel=r'$1/h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    
Out[8]:








<matplotlib.image.AxesImage at 0x108964f90>










    
























In [9]:



    


plotthisa= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(map['AVG_ALPHAFE'])-(maps.ymin+0.05*(ii%7)+0.05/2.)
    plotthisa[ii]= tmed
maps.plot(plotthisa,vmin=-0.01,vmax=0.01,minnstar=15)



    


















    






/Library/Python/2.7/site-packages/numpy/lib/nanfunctions.py:598: RuntimeWarning: Mean of empty slice
  warnings.warn("Mean of empty slice", RuntimeWarning)











    
Out[9]:








<matplotlib.image.AxesImage at 0x108a97310>

Tri broken exp





In [10]:



    


with open('../mapfits/tribrokenexp.sav','rb') as savefile:
    bf= numpy.array(pickle.load(savefile))
    samplesnf= numpy.array(pickle.load(savefile))



    











In [11]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samplesnf[ii,1])
    plotthisz[ii]= tmed
maps.plot(plotthisz*1000.,
          vmin=200.,vmax=1000.,
          minnstar=15,
          zlabel=r'$h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    
Out[11]:








<matplotlib.image.AxesImage at 0x108927d90>

Tri broken exp, fixed flare





In [12]:



    


with open('../mapfits/tribrokenexpfixedflare.sav','rb') as savefile:
    bf= numpy.array(pickle.load(savefile))
    samplesff= numpy.array(pickle.load(savefile))



    











In [13]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samplesff[ii,1])
    plotthisz[ii]= tmed
maps.plot(plotthisz*1000.,
          vmin=200.,vmax=1000.,
          minnstar=15,
          zlabel=r'$h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    
Out[13]:








<matplotlib.image.AxesImage at 0x108afba50>










    
























In [14]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samplesff[ii,1])/numpy.nanmedian(1./samples[ii,1])
    plotthisz[ii]= tmed
maps.plot(plotthisz,
          vmin=0.7,vmax=1.3,
          minnstar=15,
          zlabel=r'$h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    
Out[14]:








<matplotlib.image.AxesImage at 0x102793910>

Combination





In [15]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
plotthisze= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 5.:
        tmed= numpy.nanmedian(1./samplesnf[ii,1])
        terr= numpy.nanstd(1./samplesnf[ii,1])
    else:
        tmed= numpy.nanmedian(1./samplesff[ii,1])
        terr= numpy.nanstd(1./samplesff[ii,1])
    plotthisz[ii]= tmed
    plotthisze[ii]= terr
plotthisz[plotthisze/plotthisz > 0.2]= numpy.nan
maps.plot(plotthisz*1000.,
          vmin=200.,vmax=1000.,
          minnstar=15,
          zlabel=r'$h_Z\,(\mathrm{pc})$',
          shrink=1.)
maps.plot(plotthisz/plotthisze,
          vmin=5.,vmax=20.,
          minnstar=15,
          zlabel=r'$\sigma h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    
Out[15]:








<matplotlib.image.AxesImage at 0x108b329d0>










    

















    
























In [15]:

Plots of the radial profile's best-fit





In [35]:



    


plotthis= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samples[ii,2])
    terr= numpy.nanstd(1./samples[ii,2])
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 5.: tmed= numpy.nan
    #    terr= 1.4826*numpy.nanmedian(numpy.fabs(numpy.exp(samples[ii,3])-tmed))
    #if tmed < 5.:
    #    tmed= 0.
    plotthis[ii]= tmed
maps.plot(plotthis,vmin=1.,vmax=5.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}\,(\mathrm{kpc}^{-1})$')



    


















    
Out[35]:








<matplotlib.image.AxesImage at 0x10c1e5dd0>










    
























In [132]:



    


plotthis= numpy.zeros(len(bf))+numpy.nan
plotthise= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samples[ii,0])
    terr= numpy.nanstd(samples[ii,0])
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 6.: tmed= numpy.nan
    #    terr= 1.4826*numpy.nanmedian(numpy.fabs(numpy.exp(samples[ii,3])-tmed))
    #if tmed < 5.:
    #    tmed= 0.
    plotthis[ii]= tmed
    plotthise[ii]= terr  
#plotthis[plotthise > 0.1]= numpy.nan
maps.plot(plotthis,vmax=5.,vmin=1./.6,minnstar=15,zlabel=r'$h_{R,{\mathrm{in}}}\,(\mathrm{kpc})$')



    


















    
Out[132]:








<matplotlib.image.AxesImage at 0x118361bd0>










    
























In [51]:



    


plotthis= numpy.zeros(len(bf))+numpy.nan
plotthise= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(samples[ii,0]/samples[ii,2])
    terr= numpy.nanstd(samples[ii,0]/samples[ii,2])
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 6.: 
        tmed= numpy.nan
        terr= numpy.nan
    #    terr= 1.4826*numpy.nanmedian(numpy.fabs(numpy.exp(samples[ii,3])-tmed))
    #if tmed < 5.:
    #    tmed= 0.
    plotthis[ii]= tmed
    plotthise[ii]= terr  
#plotthis[plotthise > 0.1]= numpy.nan
maps.plot(plotthis,vmin=0.,vmax=2.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}/h_{R,{\mathrm{in}}}$')
maps.plot(plotthise,vmin=0.,vmax=2.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}/h_{R,{\mathrm{in}}}$')
maps.plot((1.-plotthis)/plotthise,vmin=-4.,vmax=4.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}/h_{R,{\mathrm{in}}}$')



    


















    
Out[51]:








<matplotlib.image.AxesImage at 0x1122c3d10>










    

















    

















    
























In [128]:



    


figsize(8,8)
plotthis= numpy.zeros(len(bf))+numpy.nan
plotthise= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(samples[ii,2]-samples[ii,0])
    terr= numpy.nanstd(samples[ii,2]-samples[ii,0])
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 6.: 
        tmed= numpy.nan
        terr= numpy.nan
    #    terr= 1.4826*numpy.nanmedian(numpy.fabs(numpy.exp(samples[ii,3])-tmed))
    #if tmed < 5.:
    #    tmed= 0.
    plotthis[ii]= tmed
    plotthise[ii]= terr  
#plotthis[plotthise > 0.1]= numpy.nan
maps.plot(plotthis,vmin=0.,vmax=.6,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}-h_{R,{\mathrm{in}}}$')
maps.plot(plotthis/plotthise,vmin=0.,vmax=5.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}-h_{R,{\mathrm{in}}}$')



    


















    
Out[128]:








<matplotlib.image.AxesImage at 0x1176fefd0>










    

















    
























In [151]:



    


figsize(8,8)
plotthis= numpy.zeros(len(bf))+numpy.nan
plotthise= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    diff= numpy.amin(numpy.array([numpy.log(2.)/samples[ii,0],numpy.exp(samples[ii,3])]),axis=1)
    diff+= numpy.amin(numpy.array([numpy.log(2.)/samples[ii,2],14.-numpy.exp(samples[ii,3])]),axis=1)
    tmed= numpy.nanmedian(diff)
    terr= numpy.nanstd(diff)
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 6.: 
        tmed= numpy.nan
        terr= numpy.nan
    #    terr= 1.4826*numpy.nanmedian(numpy.fabs(numpy.exp(samples[ii,3])-tmed))
    #if tmed < 5.:
    #    tmed= 0.
    plotthis[ii]= tmed
    plotthise[ii]= terr  
#plotthis[plotthise > 0.1]= numpy.nan
maps.plot(plotthis,vmin=4.,vmax=6.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}+h_{R,{\mathrm{in}}}$')
maps.plot(plotthise,vmin=0.,vmax=5.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}+h_{R,{\mathrm{in}}}$')
print numpy.nanmedian(plotthis)



    


















    






4.64182521842










    

















    
























In [62]:



    


import densprofiles
import copy



    











In [90]:



    


Rs= numpy.linspace(4.,12.,1001)
ldp= []
ldpe= []
tsamples= copy.deepcopy(samples)
tsamples[:,3,:]= numpy.log(8.)
tsamples[:,2,:]-= tsamples[:,0,:]
tsamples[:,0,:]= 0.
for ii, map in enumerate(maps.map()):
    print ii
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 6.:  continue
    if len(map) < 16: continue
    tdp= numpy.ones((len(Rs),len(tsamples[ii,0])//10))
    for jj in range(len(tsamples[ii,0])//10):
        tdp[Rs >= 8.,jj]= numpy.exp(-tsamples[ii,2,jj*10]*(Rs[Rs >= 8.]-8.))
        #densprofiles.tribrokenexpflaredisk(Rs,0.*Rs,0.*Rs,params=tsamples[ii,:,jj*10])/(tsamples[ii,1,jj*10]*numpy.exp((Rs-8.0)*tsamples[ii,4,jj]))
    ldp.append(numpy.nanmedian(tdp,axis=-1))
    ldpe.append(numpy.nanstd(tdp,axis=-1))



    


















    






0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

















In [113]:



    


figsize(16,30)
subplot(3,1,1)
ii= 0
for jj, map in enumerate(maps.map()):
    if numpy.nanmedian(numpy.exp(samples[jj,3])) < 6.:  continue
    if len(map) < 16: continue
    ii+= 1
    if jj%7 != 0: continue
    if ii >= len(ldp): break
    semilogy(Rs,ldp[ii]/ldp[ii][501],lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
    semilogy(Rs,(ldp[ii]+ldpe[ii])/(ldp[ii][501]+ldpe[ii][501]),ls=':',lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
    semilogy(Rs,(ldp[ii]-ldpe[ii])/(ldp[ii][501]-ldpe[ii][501]),ls='--',lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
semilogy(Rs,numpy.median(numpy.array([l/l[501] for l in ldp]),axis=0),lw=3.,color='k')
ylim(0.05,20.)
subplot(3,1,2)
ii= 0
for jj, map in enumerate(maps.map()):
    if numpy.nanmedian(numpy.exp(samples[jj,3])) < 6.:  continue
    if len(map) < 16: continue
    ii+= 1
    if jj%7 != 1: continue
    if ii >= len(ldp): break
    semilogy(Rs,ldp[ii]/ldp[ii][501],lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
    semilogy(Rs,(ldp[ii]+ldpe[ii])/(ldp[ii][501]+ldpe[ii][501]),ls=':',lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
    semilogy(Rs,(ldp[ii]-ldpe[ii])/(ldp[ii][501]-ldpe[ii][501]),ls='--',lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
semilogy(Rs,numpy.median(numpy.array([l/l[501] for l in ldp]),axis=0),lw=3.,color='k')
ylim(0.05,20.)
subplot(3,1,3)
ii= 0
for jj, map in enumerate(maps.map()):
    if numpy.nanmedian(numpy.exp(samples[jj,3])) < 6.:  continue
    if len(map) < 16: continue
    ii+= 1
    if jj%7 != 2: continue
    if ii >= len(ldp): break
    semilogy(Rs,ldp[ii]/ldp[ii][501],lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
    semilogy(Rs,(ldp[ii]+ldpe[ii])/(ldp[ii][501]+ldpe[ii][501]),ls=':',lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
    semilogy(Rs,(ldp[ii]-ldpe[ii])/(ldp[ii][501]-ldpe[ii][501]),ls='--',lw=3.,color=cm.coolwarm((numpy.median(map['FE_H'])+0.5)/.7))
semilogy(Rs,numpy.median(numpy.array([l/l[501] for l in ldp]),axis=0),lw=3.,color='k')
ylim(0.05,20.)



    


















    
Out[113]:








(0.05, 20.0)










    
























In [83]:



    


plot(densprofiles.tribrokenexpflaredisk(Rs,0.*Rs,0.*Rs,params=[0.,1.,0.5,numpy.log(8.),-0.1]))



    


















    
Out[83]:








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










    
























In [96]:



    


len(ldp)



    


















    
Out[96]:








26

















In [116]:



    


-numpy.log(0.3)/4.



    


















    
Out[116]:








0.30099320108148403

















In [12]:



    


plotthis= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samples[ii,2])
    terr= numpy.nanstd(1./samples[ii,2])
    if numpy.nanmedian(numpy.exp(samples[ii,3])) >= 5.: 
        tmed= numpy.nan
    if terr/tmed > 0.5:
        tmed= numpy.nan
    plotthis[ii]= tmed
maps.plot(plotthis,vmin=0.,vmax=4.,minnstar=15,zlabel=r'$h_{R,{\mathrm{out}}}$')
print numpy.nanmedian(plotthis), 1.4826*numpy.nanmedian(numpy.fabs(plotthis-numpy.nanmedian(plotthis)))/numpy.sqrt(numpy.sum(True-numpy.isnan(plotthis)))



    


















    






2.16984397169 0.149473395962

Two vertical exponentials?





In [152]:



    


with open('../mapfits/tribrokentwoexp.sav','rb') as savefile:
    bf_twoexp= numpy.array(pickle.load(savefile))
    samples_twoexp= numpy.array(pickle.load(savefile))



    











In [180]:



    


figsize(5,5)
plotthis= numpy.zeros(len(bf))+numpy.nan
plotthishz= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    hzindx= (numpy.fabs((1./samples_twoexp[ii,1]-1./samples_twoexp[ii,5])*samples_twoexp[ii,1]) > 0.5)
    #hzindx= numpy.ones(len(samples[ii,4]),dtype='bool')
    #print numpy.sum(hzindx)
    tup= sorted(densprofiles.ilogit(samples_twoexp[ii,4,hzindx]))[int(round(0.95*numpy.sum(hzindx)))]
    tmed= numpy.nanmedian(1./samples_twoexp[ii,5,hzindx])
    plotthis[ii]= tup
    plotthishz[ii]= tmed  
#plotthis[plotthise > 0.1]= numpy.nan
maps.plot(plotthis,vmin=0.,vmax=0.5,minnstar=15,zlabel=r'$\beta_2$')
maps.plot(plotthishz*1000.,vmin=0.,vmax=1000.,minnstar=15,zlabel=r'$h_{Z,2}\,(\mathrm{pc})$')
print numpy.nanmedian(plotthis)



    


















    






0.133979409307










    

















    
























In [174]:



    


plotthisz= numpy.zeros(len(bf))+numpy.nan
plotthisze= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    tmed= numpy.nanmedian(1./samples_twoexp[ii,1])
    terr= numpy.nanstd(1./samples_twoexp[ii,1])
    plotthisz[ii]= tmed
    plotthisze[ii]= terr
plotthisz[plotthisze/plotthisz > 0.2]= numpy.nan
maps.plot(plotthisz*1000.,
          vmin=200.,vmax=1000.,
          minnstar=15,
          zlabel=r'$h_Z\,(\mathrm{pc})$',
          shrink=1.)



    


















    
Out[174]:








<matplotlib.image.AxesImage at 0x11b72f7d0>










    
























In [207]:



    


figsize(5,5)
plotthisfrac= numpy.zeros(len(bf))+numpy.nan
plotthisfrace= numpy.zeros(len(bf))+numpy.nan
plotthisnhz= numpy.zeros(len(bf))+numpy.nan
plotthisnhze= numpy.zeros(len(bf))+numpy.nan
plotthishz= numpy.zeros(len(bf))+numpy.nan
plotthishze= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    hzindx= densprofiles.ilogit(samples_twoexp[ii,4]) > 0.1
    tmed= numpy.nanmedian(samples_twoexp[ii,1,hzindx]/samples_twoexp[ii,5,hzindx])
    terr= numpy.nanstd(samples_twoexp[ii,1,hzindx]/samples_twoexp[ii,5,hzindx])
    plotthisfrac[ii]= tmed
    plotthisfrace[ii]= terr
    if numpy.nanmedian(numpy.exp(samples[ii,3])) < 5.:
        tmed= numpy.nanmedian(1./samplesnf[ii,1])
        terr= numpy.nanstd(1./samplesnf[ii,1])
    else:
        tmed= numpy.nanmedian(1./samplesff[ii,1])
        terr= numpy.nanstd(1./samplesff[ii,1])
    plotthishz[ii]= tmed
    plotthishze[ii]= terr
plotthisfrac[plotthishze/plotthishz > 0.2]= numpy.nan
#plotthisfrac[plotthisfrace/plotthisfrac > 2.]= numpy.nan
maps.plot(plotthisfrac,vmin=0.,vmax=2.,minnstar=15,zlabel=r'$h_{Z,2}\,(\mathrm{pc})$')
maps.plot(plotthisfrace,vmin=0.,vmax=2.,minnstar=15,zlabel=r'$h_{Z,2}\,(\mathrm{pc})$')
print numpy.nanmedian(plotthisfrac)



    


















    






0.945832673875










    

















    
























In [252]:



    


figsize(6,6)
plotthisz1= numpy.zeros(len(bf))+numpy.nan
plotthisz1e= numpy.zeros(len(bf))+numpy.nan
plotthisz2= numpy.zeros(len(bf))+numpy.nan
plotthisz2e= numpy.zeros(len(bf))+numpy.nan
for ii, map in enumerate(maps.map()):
    hzindx= densprofiles.ilogit(samples_twoexp[ii,4]) > 0.15
    tmed= numpy.nanmedian(1./samples_twoexp[ii,1,hzindx])
    terr= numpy.nanstd(1./samples_twoexp[ii,1,hzindx])
    plotthisz1[ii]= tmed
    plotthisz1e[ii]= terr
    tmed= numpy.nanmedian(1./samples_twoexp[ii,5,hzindx])
    terr= numpy.nanstd(1./samples_twoexp[ii,5,hzindx])
    plotthisz2[ii]= tmed
    plotthisz2e[ii]= terr
plotthisz1[plotthisz1e/plotthisz1 > 0.5]= numpy.nan
plot(plotthisz2*1000.,plotthisz1*1000.,'ko')
errorbar(plotthisz2*1000.,plotthisz1*1000.,
         yerr=plotthisz1e*1000.,
         marker='o',color='k',ls='none')
plot((0,2000),(0,2000),'k-')
xlim(0.,1200.)
ylim(0.,1200.)
maps.plot(plotthisz1*1000.,vmin=200.,vmax=1000.,minnstar=15,zlabel=r'$h_{Z,2}\,(\mathrm{pc})$')



    


















    
Out[252]:








<matplotlib.image.AxesImage at 0x1186ddd50>

Flaring PDFs

Low alpha





In [1]:



    


import sys
import pickle
import numpy
from galpy.util import bovy_plot
from matplotlib import cm
import define_rcsample
import os
from matplotlib import cm
import extreme_deconvolution as XD
os.environ['OMP_NUM_THREADS']= '1'



    











In [2]:



    


with open('../mapfits/tribrokenexpflare.sav','rb') as savefile:
    bf= numpy.array(pickle.load(savefile))
    samples= numpy.array(pickle.load(savefile))
maps= define_rcsample.MAPs()



    











In [4]:



    


plotmaps= [9,16,23,29,36,43,50,57,64,71]
cmap= cm.coolwarm
xrange= [-0.4,0.2]
overplot= False
plotXDFit= True
ngauss= 2
allxamp= numpy.empty((len(plotmaps),ngauss))
allxmean= numpy.empty((len(plotmaps),ngauss,1))
allxcovar= numpy.empty((len(plotmaps),ngauss,1,1))
cnt= 0
for ii, map in enumerate(maps.map()):
    if not ii in plotmaps: continue
    tfeh= round(numpy.median(map['FE_H'])*20.)/20.
    if tfeh == 0.35: tfeh= 0.4
    if tfeh == -0.0: tfeh= 0.0
    print tfeh
    bovy_plot.bovy_hist(samples[ii,4],range=xrange,bins=51,overplot=overplot,
                        histtype='step',normed=True,zorder=2,
                        color=cmap((tfeh+0.5)*0.95/0.9+0.05))
    # Fit PDFs with XD
    xamp= numpy.array([0.45,0.5])
    xmean= numpy.array([numpy.mean(samples[ii,4])+numpy.random.normal()*numpy.std(samples[ii,4]),
                         numpy.mean(samples[ii,4])+numpy.random.normal()*numpy.std(samples[ii,4])])[:,numpy.newaxis]
    xcovar= numpy.reshape(numpy.tile(numpy.var(samples[ii,4]),(2,1)),(2,1,1))
    XD.extreme_deconvolution(samples[ii,4][:,numpy.newaxis],numpy.zeros((len(samples[ii,4]),1)),
                             xamp,xmean,xcovar)
    allxamp[cnt]= xamp
    allxmean[cnt]= xmean
    allxcovar[cnt]= xcovar
    if plotXDFit:
        txs= numpy.linspace(xrange[0],xrange[1],1001)
        plot(txs,1./numpy.sqrt(2.*numpy.pi)*(xamp[0]/numpy.sqrt(xcovar[0,0,0])*numpy.exp(-0.5*(txs-xmean[0,0])**2./xcovar[0,0,0])
                                             +xamp[1]/numpy.sqrt(xcovar[1,0,0])*numpy.exp(-0.5*(txs-xmean[1,0])**2./xcovar[1,0,0])),
             color=cmap((tfeh+0.5)*0.95/0.9+0.05),zorder=1)
    overplot=True
    cnt+= 1
txs= numpy.linspace(xrange[0],xrange[1],1001)
comb= numpy.ones_like(txs)
for ii in range(len(plotmaps)):
    comb*= 1./numpy.sqrt(2.*numpy.pi)*(allxamp[ii,0]/numpy.sqrt(allxcovar[ii,0,0,0])*numpy.exp(-0.5*(txs-allxmean[ii,0,0])**2./allxcovar[ii,0,0,0])
                                       +allxamp[ii,1]/numpy.sqrt(allxcovar[ii,1,0,0])*numpy.exp(-0.5*(txs-allxmean[ii,1,0])**2./allxcovar[ii,1,0,0]))
comb/= numpy.sum(comb)*(txs[1]-txs[0])
plot(txs,comb/2.,'k-',lw=2.,zorder=20)



    


















    






-0.5
-0.4
-0.3
-0.2
-0.1
0.0
0.1
0.2
0.3
0.4










    
Out[4]:








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










    
























In [6]:



    


numpy.sum(comb*txs)/numpy.sum(comb), numpy.sqrt(numpy.sum(comb*txs**2.)/numpy.sum(comb)-(numpy.sum(comb*txs)/numpy.sum(comb))**2.)



    


















    
Out[6]:








(-0.12016519159285134, 0.0082092724756257177)

High alpha





In [8]:



    


plotmaps= [12,19,26,32,39,45]#,51]
cmap= cm.coolwarm
xrange= [-0.3,0.3]
overplot= False
plotXDFit= True
ngauss= 2
allxamp= numpy.empty((len(plotmaps),ngauss))
allxmean= numpy.empty((len(plotmaps),ngauss,1))
allxcovar= numpy.empty((len(plotmaps),ngauss,1,1))
cnt= 0
for ii, map in enumerate(maps.map()):
    if not ii in plotmaps: continue
    tfeh= round(numpy.median(map['FE_H'])*20.)/20.
    if tfeh == 0.35: tfeh= 0.4
    if tfeh == -0.0: tfeh= 0.0
    print tfeh
    bovy_plot.bovy_hist(samples[ii,4],range=xrange,bins=51,overplot=overplot,
                        histtype='step',normed=True,zorder=2,
                        color=cmap((tfeh+0.2)*0.65/0.6+0.05))
    # Fit PDFs with XD
    xamp= numpy.array([0.45,0.5])
    xmean= numpy.array([numpy.mean(samples[ii,4])+numpy.random.normal()*numpy.std(samples[ii,4]),
                         numpy.mean(samples[ii,4])+numpy.random.normal()*numpy.std(samples[ii,4])])[:,numpy.newaxis]
    xcovar= numpy.reshape(numpy.tile(numpy.var(samples[ii,4]),(2,1)),(2,1,1))
    XD.extreme_deconvolution(samples[ii,4][:,numpy.newaxis],numpy.zeros((len(samples[ii,4]),1)),
                             xamp,xmean,xcovar)
    allxamp[cnt]= xamp
    allxmean[cnt]= xmean
    allxcovar[cnt]= xcovar
    if plotXDFit:
        txs= numpy.linspace(xrange[0],xrange[1],1001)
        plot(txs,1./numpy.sqrt(2.*numpy.pi)*(xamp[0]/numpy.sqrt(xcovar[0,0,0])*numpy.exp(-0.5*(txs-xmean[0,0])**2./xcovar[0,0,0])
                                             +xamp[1]/numpy.sqrt(xcovar[1,0,0])*numpy.exp(-0.5*(txs-xmean[1,0])**2./xcovar[1,0,0])),
             color=cmap((tfeh+0.5)*0.95/0.9+0.05),zorder=1)
    overplot=True
    cnt+= 1
txs= numpy.linspace(xrange[0],xrange[1],1001)
comb= numpy.ones_like(txs)
for ii in range(len(plotmaps)):
    comb*= 1./numpy.sqrt(2.*numpy.pi)*(allxamp[ii,0]/numpy.sqrt(allxcovar[ii,0,0,0])*numpy.exp(-0.5*(txs-allxmean[ii,0,0])**2./allxcovar[ii,0,0,0])
                                       +allxamp[ii,1]/numpy.sqrt(allxcovar[ii,1,0,0])*numpy.exp(-0.5*(txs-allxmean[ii,1,0])**2./allxcovar[ii,1,0,0]))
comb/= numpy.sum(comb)*(txs[1]-txs[0])
plot(txs,comb/2.,'k-',lw=2.,zorder=20)



    


















    






-0.5
-0.4
-0.3
-0.2
-0.1
0.0
0.1










    
Out[8]:








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










    
























In [9]:



    


numpy.sum(comb*txs)/numpy.sum(comb), numpy.sqrt(numpy.sum(comb*txs**2.)/numpy.sum(comb)-(numpy.sum(comb*txs)/numpy.sum(comb))**2.)



    


















    
Out[9]:








(-0.033220558878686614, 0.021213741080536291)

















In [13]:



    


txs= numpy.linspace(xrange[0],xrange[1],1001)
comb= numpy.ones_like(txs)
for ii in range(len(plotmaps)-1):
    comb*= 1./numpy.sqrt(2.*numpy.pi)*(allxamp[ii,0]/numpy.sqrt(allxcovar[ii,0,0,0])*numpy.exp(-0.5*(txs-allxmean[ii,0,0])**2./allxcovar[ii,0,0,0])
                                       +allxamp[ii,1]/numpy.sqrt(allxcovar[ii,1,0,0])*numpy.exp(-0.5*(txs-allxmean[ii,1,0])**2./allxcovar[ii,1,0,0]))
comb/= numpy.sum(comb)*(txs[1]-txs[0])
plot(txs,comb/2.,'k-',lw=2.,zorder=20)



    


















    
Out[13]:








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










    
























In [14]:



    


numpy.sum(comb*txs)/numpy.sum(comb), numpy.sqrt(numpy.sum(comb*txs**2.)/numpy.sum(comb)-(numpy.sum(comb*txs)/numpy.sum(comb))**2.)



    


















    
Out[14]:








(-0.0047754633718367539, 0.024751623039118684)

















In [ ]:

Content source: jobovy/apogee-maps

Similar notebooks:

notebook.community | gallery | about