In [1]:

    

%pylab inline
from astropy.io import fits


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

print('## Local normalisation, compare to Czekaj+14 table 7 ##')
outputDir = '../output/more_bulge_100pc'
#outputDir = '../output/mockcat_new_SFR_new_IMF_100pc_more_thick'
x = fits.getdata(outputDir + '/' + 'GDR2mock_20.7Gmag.fits')
smaller_sphere = 50
sphere = (4/3)*np.pi*smaller_sphere**3
for popid in np.arange(9):
    selection = (x.popid==popid) & (x.logg < 7) & (np.divide(1,x.parallax)<smaller_sphere/1000)
    if (popid!=8):
        print('popid=',popid,' %.1f Msun/pc^3*10^-3' %(sum(x.smass[selection])/sphere*1000))
    else:
        print('popid=',popid,' %.1f Msun/pc^3*10^-6' %(sum(x.smass[selection])/sphere*1e6))

selection = (x.popid<=6) & (x.logg > 7)& (np.divide(1,x.parallax)<smaller_sphere/1000)
print('thindisc WD',' %.1f Msun/pc^3*10^-3 %.1f mact' %(sum(x.smass[selection])/sphere*1000,sum(x.mact[selection])/sphere*1000))
selection = (x.popid==7) & (x.logg > 7)& (np.divide(1,x.parallax)<smaller_sphere/1000)
print('thick disc WD',' %.1f Msun/pc^3*10^-3 %.1f mact' %(sum(x.smass[selection])/sphere*1000,sum(x.mact[selection])/sphere*1000))


    

    




## Local normalisation, compare to Czekaj+14 table 7 ##






    




---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-2-de2d2df0e57d> in <module>()
      2 outputDir = '../output/more_bulge_100pc'
      3 #outputDir = '../output/mockcat_new_SFR_new_IMF_100pc_more_thick'
----> 4 x = fits.getdata(outputDir + '/' + 'GDR2mock_20.7Gmag.fits')
      5 smaller_sphere = 50
      6 sphere = (4/3)*np.pi*smaller_sphere**3

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/convenience.py in getdata(filename, header, lower, upper, view, *args, **kwargs)
    190     mode, closed = _get_file_mode(filename)
    191 
--> 192     hdulist, extidx = _getext(filename, mode, *args, **kwargs)
    193     try:
    194         hdu = hdulist[extidx]

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/convenience.py in _getext(filename, mode, ext, extname, extver, *args, **kwargs)
   1009         raise TypeError('extver alone cannot specify an extension.')
   1010 
-> 1011     hdulist = fitsopen(filename, mode=mode, **kwargs)
   1012 
   1013     return hdulist, ext

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in fitsopen(name, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    149 
    150     return HDUList.fromfile(name, mode, memmap, save_backup, cache,
--> 151                             lazy_load_hdus, **kwargs)
    152 
    153 

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in fromfile(cls, fileobj, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    385         return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,
    386                              save_backup=save_backup, cache=cache,
--> 387                              lazy_load_hdus=lazy_load_hdus, **kwargs)
    388 
    389     @classmethod

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in _readfrom(cls, fileobj, data, mode, memmap, save_backup, cache, lazy_load_hdus, **kwargs)
    972             if not isinstance(fileobj, _File):
    973                 # instantiate a FITS file object (ffo)
--> 974                 fileobj = _File(fileobj, mode=mode, memmap=memmap, cache=cache)
    975             # The Astropy mode is determined by the _File initializer if the
    976             # supplied mode was None

~/anaconda3/lib/python3.6/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    486                         # one with the name of the new argument to the function
    487                         kwargs[new_name[i]] = value
--> 488             return function(*args, **kwargs)
    489 
    490         return wrapper

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in __init__(self, fileobj, mode, memmap, overwrite, cache)
    173             self._open_fileobj(fileobj, mode, overwrite)
    174         elif isinstance(fileobj, str):
--> 175             self._open_filename(fileobj, mode, overwrite)
    176         else:
    177             self._open_filelike(fileobj, mode, overwrite)

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in _open_filename(self, filename, mode, overwrite)
    529 
    530         if not self._try_read_compressed(self.name, magic, mode, ext=ext):
--> 531             self._file = fileobj_open(self.name, IO_FITS_MODES[mode])
    532             self.close_on_error = True
    533 

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/util.py in fileobj_open(filename, mode)
    386     """
    387 
--> 388     return open(filename, mode, buffering=0)
    389 
    390 

FileNotFoundError: [Errno 2] No such file or directory: '../output/more_bulge_100pc/GDR2mock_20.7Gmag.fits'

In [3]:

    

# Galaxia parametrization
sfrv = np.array([1.60552043, 1.79400035, 2.02272975, 2.42164718,3.07851692, 4.15556317])
#thindisk0
hr1=5.0;
hr2=3.0;
rho0=1.5120222/(5.5683282*(hr1*hr1*hr1-hr2*hr2*hr2))
print(rho0)
rho0=2.37/(5.5683282*(hr1*hr1*hr1-hr2*hr2*hr2))
print(rho0)

#thindisk1        
hr1=2.400#the value now is used in Czekaj+2014, 2.200; used in Reyle+2009//2.530; used in Robin+ 2003, Galaxia original 
hr2=1.320#;
for i,sfr in enumerate(sfrv):
    rho0=sfr/(23.719602*(hr1*hr1*hr1-hr2*hr2*hr2))
    print(i+1, rho0)
#thickdisc        
hr1=2.5;
hr2=0.8;
xl=0.4;
c1=1/((hr2*xl)*(2+xl/hr2));
c2=np.exp(xl/hr2)/(1+0.5*xl/hr2);
rho0=1.34e6*4.4634;
print(rho0)
#Halo:
ac=0.5;
rho0=9.32e3*6.4159*1.1;

        
        
#Bulge:
rho0=0.255*13.7e9;
a=78.9 ; b=3.5; g=91.3;
x0=1.59; y0=0.424; z0=0.424;
Rc=2.54;


    

    




0.0027708133222434112
0.004343076162318837
1 0.005873594806356947
2 0.006563124916674242
3 0.007399902693398728
4 0.008859291998717789
5 0.011262367425990327
6 0.015202605832178794
5980956.0

In [4]:

    

agev = np.array([0.075e9, 0.575e9 ,1.5e9, 2.5e9, 4e9,   6e9,   8.5e9])
agev = np.divide(agev,1e9)
dagev = np.array([0.075e9, 0.425e9, 0.5e9, 0.5e9, 1e9,   1e9,   1.5e9])
dagev = np.divide(dagev,1e9)


    

In [5]:

    

print('## Local normalisation, compare to Czekaj+14 table 7 ##')
outputDir = '../output/GDR2mock_207'
#outputDir = '../output/mockcat_new_SFR_new_IMF_100pc_more_thick'
x = fits.getdata(outputDir + '/' + 'GDR2mock_50pc.fits')
smaller_sphere = 50
sphere = (4/3)*np.pi*smaller_sphere**3
for i in np.arange(7):
    selection = (x.age<agev[i]+dagev[i]) & (x.age > agev[i]-dagev[i]) 
    print('popid=',i,' %.1f Msun/pc^3*10^-3' %(sum(x.mass[selection])/sphere*1e3))

selection = (x.age>=10) & (x.age <= 11)
print('thick disc',' %.1f Msun/pc^3*10^-3' %(sum(x.mass[selection])/sphere*1000))
selection = (x.age>11)
print('halo',' %.1f Msun/pc^3*10^-6' %(sum(x.mass[selection])/sphere*1e6))


    

    




## Local normalisation, compare to Czekaj+14 table 7 ##
popid= 0  2.9 Msun/pc^3*10^-3
popid= 1  7.6 Msun/pc^3*10^-3
popid= 2  5.3 Msun/pc^3*10^-3
popid= 3  3.8 Msun/pc^3*10^-3
popid= 4  5.7 Msun/pc^3*10^-3
popid= 5  4.6 Msun/pc^3*10^-3
popid= 6  5.0 Msun/pc^3*10^-3
thick disc  1.2 Msun/pc^3*10^-3
halo  10.2 Msun/pc^3*10^-6

In [ ]: