In [22]:
%autosave 20
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
import astropy.units as u
from scipy import optimize
import pandas as pd
Autosaving every 20 seconds
In [20]:
coord = SkyCoord(ra=5*u.hourangle, dec=-5*u.deg-23*u.arcmin,
                 radial_velocity=5*u.km/u.s,
                 distance=100*u.pc,
                 pm_ra_cosdec=1*u.mas/u.yr, pm_dec=-1*u.mas/u.yr)
display(coord.transform_to('galactic'))
<SkyCoord (Galactic): (l, b, distance) in (deg, deg, pc)
    (204.61668338, -27.18344969, 100.)
 (pm_l_cosb, pm_b, radial_velocity) in (mas / yr, mas / yr, km / s)
    (1.35101031, 0.41805639, 5.)>
<SkyCoord (CIRS: obstime=J2000.000): (ra, dec, distance) in (deg, deg, pc)
    (75.00537301, -5.38592396, 99.99999628)
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    (-56064.66875155, -64131.82653705, 17.96993101)>
In [37]:
data = np.genfromtxt('diagram.csv', delimiter=',', names=('d', 'v'))
plt.plot(data['d'], data['v'], 'x')

def v_d(d, a, b):
    return a + b * d

def v_d_jac(d, a, b):
    jac = np.empty(d.shape + (2,))
    jac[:,0] = 1
    jac[:,1] = d
    return jac

(v0, H0), _ = optimize.curve_fit(v_d, data['d'], data['v'],
                                 (0, 500), jac=v_d_jac)

d_ = np.r_[data['d'].min():data['d'].max():100j]
plt.plot(d_, v_d(d_, v0, H0))
Out[37]:
[<matplotlib.lines.Line2D at 0x7f538a359208>]
In [52]:
from astroquery import vizier
from astroquery import simbad
# conda install -c astropy astroquery
# pip3 install astroquery
m87 = simbad.Simbad.query_object('M87')
display(m87)
m87_coord = SkyCoord(ra=m87['RA'], dec=m87['DEC'],
                     unit=(u.hourangle, u.deg))
display(m87_coord)

viz = vizier.Vizier(
    column_filters={'Jmag': '<10'},
    row_limit=100,
)
tables = viz.query_region(m87_coord[0], radius=1*u.deg,
                          catalog='2MASS')
tables[0]
Table masked=True length=1
MAIN_IDRADECRA_PRECDEC_PRECCOO_ERR_MAJACOO_ERR_MINACOO_ERR_ANGLECOO_QUALCOO_WAVELENGTHCOO_BIBCODE
"h:m:s""d:m:s"masmasdeg
objectstr13str13int16int16float32float32int16str1str1object
M 8712 30 49.4233+12 23 28.043990.2500.01790AR2009A&A...493..317L
<SkyCoord (ICRS): (ra, dec) in deg
    [(187.70593042, 12.39112306)]>
Out[52]:
Table masked=True length=48
RAJ2000DEJ2000_2MASSJmage_JmagHmage_HmagKmage_KmagQflgRflgBflgCflgXflgAflg
degdegmagmagmagmagmagmag
float64float64bytes17float32float32float32float32float32float32bytes3bytes3bytes3bytes3uint8uint8
187.38466811.49256112293232+11293327.7340.0247.5500.0277.5250.024AAA11111100000
187.03457211.69636412280829+11414699.3240.0278.9520.0228.9060.023AAA22211100000
187.45246811.72599712294859+11433358.3990.0268.0270.0407.9510.017AAA11111100000
187.47985111.72515712295516+11433059.8020.0239.3010.0239.2050.021AAA22211100000
187.63028111.69170712303126+11413019.4510.0369.0640.0448.8680.027EEE222111c0000
187.12560011.76738212283014+11460258.8840.0308.6110.0348.5160.021AAA11211100000
187.85192711.75506012312446+11451829.2120.0229.0350.0238.9380.023AAA22211100000
188.07578811.79057212321818+11472609.5840.0239.2960.0269.2020.021AAA22211100000
187.06193411.82652712281486+11493545.5600.0234.9310.0154.6560.026AAA11111100000
.............................................
186.87817612.75253712273076+12450918.5880.0218.0220.0427.8640.020AAA11111100000
187.46644412.80791712295194+12482859.9290.0209.7820.0199.7810.020AAA22211100000
187.51749512.81025612300419+12483698.7630.0348.4010.0558.2030.017AAA11111100000
187.27803613.04375612290672+13023759.9230.0209.4320.0179.3380.016AAA22211100000
187.19948313.08793312284787+13051658.4340.0218.2420.0178.1510.021AAA11111100000
187.97930813.10705812315503+13062549.4300.0228.6810.0238.4720.019AAA21211100000
187.27665813.24975312290639+13145919.6850.0219.3750.0179.3210.016AAA22211100000
187.47753113.24913212295460+13145689.3220.0208.7610.0178.6170.016AAA22211100000
187.53685313.23347212300884+13140048.6650.0397.9310.0617.7430.027AAA11111100000
188.05693713.24672212321366+13144819.8320.0219.3050.0229.2010.019AAA22211100000
In [81]:
from scipy import integrate

# JLA
viz = vizier.Vizier(row_limit=1000)
tables = viz.get_catalogs('J/A+A/568/A22/tablef3')
t = tables[0]

def distance(z, Omega):
    """Luminosity distance, pc"""
    z = np.asarray(z, dtype=float)
    H0 = 70 / 1e6
    c = 3e5
    d = np.empty_like(z)
    for i, z_i in np.ndenumerate(z):
        d[i] = c / H0 * (1+z_i) * integrate.quad(
            lambda z: 1 / np.sqrt(Omega + (1-Omega) * (1+z)**3),
            0, z_i
        )[0]
    return d

def res(params, m, z, x1, c):
    Omega, M0, alpha, beta = params
    d = distance(z, Omega)
    m_th = M0 - 5 + 5*np.log10(d) - alpha * x1 + beta * c
    return m_th - m

result = optimize.least_squares(
    res,
    x0=(0.7, -19, 0.15, 2),
    args=(t['mb'], t['zcmb'], t['x1'], t['c']),
    bounds=([0, -22, 0, 0], [1, -18, 1, 3]),
)
print(result)
result.x

# import lmfit
 active_mask: array([0, 0, 0])
        cost: 10.570052101742341
         fun: <MaskedColumn name='x1' dtype='float64' length=740>
0.021282058144056037
-0.17635845722222498
 -0.0182790885464037
-0.14872565447865682
-0.27182356955313125
 0.27798067612144806
 0.06314672613605055
 0.13666835411755685
-0.12107218644429807
 -0.1098988633869844
-0.08689218861948333
 0.12286884771465267
                 ...
-0.16798386349589478
0.008762473695025363
-0.16625423469364264
 0.12894200899473063
0.019600891127151243
 0.08367782935281198
0.007105040712996669
 0.11298463625138666
   0.194534279932693
-0.12394845300535806
 0.01920676416535727
-0.06461656038566232
        grad: <MaskedColumn name='x1' dtype='float64' length=3>
-0.006843629535302376
  0.03138103254929803
  -0.2826112195655166
         jac: array([[ 0.77539682,  1.        , -1.        ],
       [ 0.87925053,  1.        , -1.        ],
       [ 0.76443148,  1.        ,  1.        ],
       ...,
       [ 0.03911877,  1.        , -0.5       ],
       [ 0.03911877,  1.        , -1.        ],
       [ 0.0358572 ,  1.        , -0.25      ]])
     message: '`xtol` termination condition is satisfied.'
        nfev: 11
        njev: 3
  optimality: 0.2496852159836657
      status: 3
     success: True
           x: array([  0.72298385, -19.08708444,   0.11650636])
Out[81]:
array([  0.72298385, -19.08708444,   0.11650636])
In [ ]:
?optimize.least_squares