sun_abs_mags = {'U' : 5.61,
                'B' : 5.48,
                'V' : 4.83,
                'R' : 4.42,
                'I' : 4.08,
                'J' : 3.64,
                'H' : 3.32,
                'K' : 3.28,
                '3.6' : 3.24, # Oh et al. 2008
                'u' : 6.77, #SDSS bands from http://mips.as.arizona.edu/~cnaw/sun.html
                'g' : 5.36,
                'r' : 4.67,
                'i' : 4.48,
                'z' : 4.42
               }

s4g_rawdata = [] 

with open('IpacTableFromSource.tbl') as s4g:
    for line in s4g.readlines():
        if line[0] == '\\' or line[0] == '|':
            pass
        else:
            s4g_rawdata.append(" ".join(line.split()).split(' '))

print 'name    [3.6]   [4.5]   log10(stellar mass)   distance(Mpc)  std(distance)'

s4g_rawdata[0]

"""
Explanation: Сравнение S4G и 2MASS массовых моделей
Идея в том, чтобы сравнить полные массы галактик из S4G и из 2MASS фотометрии, чтобы понять, почему диски так недооцениваются.
Данные из https://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-dd?catalog=s4gcat и https://ned.ipac.caltech.edu/forms/gmd.html
End of explanation
"""


m_to_ls = []
for l in s4g_rawdata:
    if l[1] != 'null' and l[2] != 'null':
        m_to_ls.append(np.power(10, -0.339*(float(l[1]) - float(l[2])) -0.336))

plt.hist(m_to_ls, bins=20, alpha=0.2)
plt.axvline(x=0.6);

"""
Explanation: Разброс M/L по калибровке $log(M/L) = −0.339(±0.057)\times([3.6]−[4.5])−0.336(±0.002)$ https://arxiv.org/abs/1410.0009:
End of explanation
"""


twomass_rawdata = [] 

with open('NED_S4G_2MASS_crossdata.txt') as crossdata:
    for line in crossdata.readlines():
            twomass_rawdata.append(" ".join(line.split()).split('|'))

twomass_rawdata[0]

print 'Name         J         H      K_s'
print twomass_rawdata[0][2], twomass_rawdata[0][10], twomass_rawdata[0][14], twomass_rawdata[0][18]

print s4g_rawdata[522]
print '{:2.2E}'.format(np.power(10., float(s4g_rawdata[522][3])))

for ind, l in enumerate(twomass_rawdata):
    if l[2] == 'NGC4544 ':
        print ind
        print l

M = float(twomass_rawdata[699][14]) - 5*np.log10(float(s4g_rawdata[522][4])) + 5 - 30.
m_to_l = 1.0
print '{:2.2E}'.format(m_to_l*np.power(10., 0.4*(sun_abs_mags['H'] - M)))

s4g_dict = {}
for l in s4g_rawdata:
    if l[3] != 'null' and l[4] != 'null' and l[5] != 'null':
        s4g_dict[l[0]+' '] = (float(l[3]), float(l[4]), float(l[5]))

final = []
for gal in twomass_rawdata:
    if gal[2] in s4g_dict.keys() and gal[10] != ' ' and gal[14] != ' ' and gal[18] != ' ':
        log_stellar, dist, err_dist = s4g_dict[gal[2]]
        stellar_mass_s4g = np.power(10., log_stellar)
        magJ, magH, magK = float(gal[10]), float(gal[14]), float(gal[18])
        photom = []
        for m_to_l, mag, band in zip([1.2, 1., 0.8], [magJ, magH, magK], ['J','H','K']):
            M = mag - 5*np.log10(dist) - 30. + 5.
            twomass_stellar = m_to_l*np.power(10., 0.4*(sun_abs_mags[band] - M))
            photom.append(twomass_stellar)
        final.append((gal[2], stellar_mass_s4g, photom))

len(final)

fig = plt.figure(figsize = [16, 5])
ax = plt.subplot('131')
ax.plot([l[1] for l in final], [l[2][0] for l in final], '.')
ax.set_title('J')
ax.set_xlabel('S4G')
ax.set_ylabel('2MASS')
ax.plot(ax.get_xlim(), ax.get_ylim(), ls="--", c=".3")

ax = plt.subplot('132')
ax.plot([l[1] for l in final], [l[2][1] for l in final], '.')
ax.set_title('H')
ax.set_xlabel('S4G')
ax.set_ylabel('2MASS')
ax.plot(ax.get_xlim(), ax.get_ylim(), ls="--", c=".3")

ax = plt.subplot('133')
ax.plot([l[1] for l in final], [l[2][2] for l in final], '.')
ax.set_title('K_s')
ax.set_xlabel('S4G')
ax.set_ylabel('2MASS')
ax.plot(ax.get_xlim(), ax.get_ylim(), ls="--", c=".3")

plt.show()

"""
Explanation: Разброс достаточно маленький.
End of explanation
"""


fig = plt.figure(figsize = [8, 5])
ax = plt.subplot()
ax.plot([l[1] for l in final], [l[2][0] for l in final], '.', label='J')

ax = plt.subplot()
ax.plot([l[1] for l in final], [l[2][1] for l in final], '.', label='H')

ax = plt.subplot()
ax.plot([l[1] for l in final], [l[2][2] for l in final], '.', label='K')
ax.set_title('J, H, K_s')
ax.set_xlabel('S4G')
ax.set_ylabel('2MASS')
ax.plot(ax.get_xlim(), ax.get_ylim(), ls="--", c=".3")
plt.legend(loc='best')
plt.show()

for l in final:
    if l[0] in ['NGC2985 ', 'NGC4725 ', 'NGC4258 ']:
        print l

for ind, l in enumerate(twomass_rawdata):
    if l[2] in ['NGC2985 ', 'NGC4725 ', 'NGC4258 ']:
        print ind
        print l

"""
Explanation: Как видим, все замечательно согласуется.
End of explanation
"""


for ind, l in enumerate(s4g_rawdata):
    if l[0] in ['NGC2985', 'NGC4725', 'NGC4258']:
        print ind
        print l

print s4g_rawdata[1829]
print '{:2.2E}'.format(np.power(10., float(s4g_rawdata[1829][3])))

"""
Explanation: Что насчет полной массы из калибровок, как она согласуется с табличной?
https://irsa.ipac.caltech.edu/cgi-bin/2MASS/PubGalPS/nph-galps?locstr=147.592417+72.278967+eq&radius=5
End of explanation
"""


M = 7.493 - 5*np.log10(float(s4g_rawdata[1829][4])) + 5 - 30.
m_to_l = 1.0
print '{:2.2E}'.format(m_to_l*np.power(10., 0.4*(sun_abs_mags['H'] - M)))

M = 8.307 - 5*np.log10(float(s4g_rawdata[1829][4])) + 5 - 30.
m_to_l = 1.2
print '{:2.2E}'.format(m_to_l*np.power(10., 0.4*(sun_abs_mags['J'] - M)))

M = 7.360 - 5*np.log10(float(s4g_rawdata[1829][4])) + 5 - 30.
m_to_l = 0.8
print '{:2.2E}'.format(m_to_l*np.power(10., 0.4*(sun_abs_mags['K'] - M)))

final_s4g = []
for l in s4g_rawdata:
    if l[1] != 'null' and l[2] != 'null' and l[3] != 'null':
        stellar_mass_s4g = np.power(10., float(l[3]))
        
        m_to_l = np.power(10, -0.339*(float(l[1]) - float(l[2])) -0.336)
        photom_stellar = m_to_l*np.power(10., 0.4*(sun_abs_mags['3.6'] - float(l[1])))
        final_s4g.append((l[0], stellar_mass_s4g, photom_stellar))

len(final_s4g)

fig = plt.figure(figsize = [5, 5])
ax = plt.subplot()
ax.plot(zip(*final_s4g)[1], zip(*final_s4g)[2], '.')

ax.set_title('S4G calibration')
ax.set_xlabel('stellar mass from table')
ax.set_ylabel('from calibration')
ax.plot(ax.get_xlim(), ax.get_ylim(), ls="--", c=".3")
plt.show()

"""
Explanation: j_m_ext   8.307
h_m_ext   7.493
k_m_ext   7.360
End of explanation
"""


for l in final_s4g:
    if l[0] in ['NGC2985', 'NGC4725', 'NGC4258']:
        print l

"""
Explanation: Хорошо согласуется.
End of explanation
"""


