import numpy as np import multiprocessing as mp import argparse import emcee """ Modeling rotation curve of the Milky Way based on proper motions and radial velocities of Classical Cepheids. P. Mroz @ OAUW, 12 Sep 2018 """ def calculate_tangential_velocity (ra,dec,l,b,dist,pm_ra,pm_ra_err,pm_dec,pm_dec_err,pm_corr,ra_G,dec_G): """ Calculating components of tangential velocity and their uncertainties Parameters: ----------- ra : right ascension (radians) dec : declination (radians) l : Galactic longitude (radians) b : Galactic latitude (radians) dist : Heliocentric distance (kpc) pm_ra : proper motion in right ascension (mas/yr) pm_ra_err : uncertainty of proper motion in right ascension (mas/yr) pm_dec : proper motion in declination (mas/yr) pm_dec_err : uncertainty of proper motion in declination (mas/yr) pm_corr : correlation coefficient between pm_ra and pm_dec v_rad : Heliocentric radial velocity (km/s) v_rad_err : uncertainty of Heliocentric radial velocity (km/s) ra_G : right ascension of the North Galactic pole (radians) dec_G : declination of the North Galactic pole (radians) Returns: -------- v_l : tangential velocity along l (km/s) v_l_err : uncertainty of tangential velocity along l (km/s) v_b : tangential velocity along b (km/s) v_b_err : uncertainty of tangential velocity along b (km/s) pm_lb_corr : correlation coefficient between v_l and v_b """ # Transforming equatorial proper motions to the Galactic system # see Poleski 2013 (arXiv:1306.2945) C1 = np.sin(dec_G)*np.cos(dec)-np.cos(dec_G)*np.sin(dec)*np.cos(ra-ra_G) C2 = np.cos(dec_G)*np.sin(ra-ra_G) norm = np.sqrt(C1*C1+C2*C2) C1 /= norm C2 /= norm pm_l = C1*pm_ra+C2*pm_dec pm_b = C1*pm_dec-C2*pm_ra # Transforming error bars of proper motions J1 = np.array([[C1,C2],[-C2,C1]]) J2 = np.transpose(J1) CORR = np.array([[pm_ra_err**2,pm_corr*pm_ra_err*pm_dec_err], [pm_corr*pm_ra_err*pm_dec_err,pm_dec_err**2]]) COV = np.dot(np.dot(J1,CORR),J2) pm_l_err = np.sqrt(COV[0,0]) pm_b_err = np.sqrt(COV[1,1]) pm_lb_corr = COV[0,1] / pm_l_err / pm_b_err # Calculating tangential components of velocity vector and their error bars v_l = 4.74 * dist * pm_l v_b = 4.74 * dist * pm_b v_l_err = np.absolute(v_l*pm_l_err/pm_l) v_b_err = np.absolute(v_b*pm_b_err/pm_b) return v_l, v_l_err, v_b, v_b_err, pm_lb_corr def lnprob_v1 (pars,l,b,dist,v_l,v_l_err,v_b,v_b_err,v_r,v_r_err): """ Likelihood function for model 1 Parameters: ----------- pars : Model parameters l : Galactic longitude (radians) b : Galactic latitude (radians) dist : heliocentric distance (kpc) v_l : Tangential velocity (l component) (km/s) v_l_err : Uncertainty of tangential velocity (l component) (km/s) v_b : Tangential velocity (b component) (km/s) v_b_err : Uncertainty of tangential velocity (b component) (km/s) v_rad : Radial velocity (km/s) v_rad_err : Uncertainty of radial velocity (km/s) Returns: -------- lnprob : log-likelihood """ R0,theta_0,Us,Vs,Ws,U0,V0,W0,eps_l,eps_b,eps_r = pars # PRIORS: if Us < -100.0: return -np.inf if Us > 100.0: return -np.inf if Vs < -100.0: return -np.inf if Vs > 100.0: return -np.inf if Ws < -100.0: return -np.inf if Ws > 100.0: return -np.inf if R0 < 7.0: return -np.inf if R0 > 9.0: return -np.inf if eps_l < 0.0: return -np.inf if eps_b < 0.0: return -np.inf if eps_r < 0.0: return -np.inf # priors on U0,V0,W0 from Schonrich et al. (2010) lprior = 0.0 lprior += -0.5 * (U0-11.1)**2 / (1.3*1.3) lprior += -0.5 * (V0-12.2)**2 / (2.1*2.1) lprior += -0.5 * (W0-7.3)**2 / (0.7*0.7) lprior += -0.5 * (R0-8.122)**2 / (0.031*0.031) # calculating projected distance from the Galactic center and angle beta # see Appendix of Reid et al. (2009) Dp = dist * np.cos(b) Rp = R0*R0 + Dp*Dp - 2.0*R0*Dp*np.cos(l) Rp = np.sqrt(Rp) sinbeta = Dp * np.sin(l) / Rp cosbeta = (R0 - Dp * np.cos(l)) / Rp beta = np.arctan2(sinbeta,cosbeta) # rotating velocity vector through the angle -beta and subtracting # velocity of the Sun (equations 1--3) theta_S = theta_0 U2 = Us * np.cos(beta) + (Vs + theta_S) * np.sin(beta) V2 = -Us * np.sin(beta) + (Vs + theta_S) * np.cos(beta) W2 = Ws U1 = U2 - U0 V1 = V2 - V0 - theta_0 W1 = W2 - W0 # calculating predicted tangential and radial velocities # (equations 4--6): v_l_model = V1 * np.cos(l) - U1 * np.sin(l) tmp = U1 * np.cos(l) + V1 * np.sin(l) v_b_model = W1 * np.cos(b) - tmp * np.sin(b) v_r_model = tmp * np.cos(b) + W1 * np.sin(b) # correcting the error bars v_l_err = np.sqrt(v_l_err**2+(eps_l)**2) v_b_err = np.sqrt(v_b_err**2+(eps_b)**2) v_r_err = np.sqrt(v_r_err**2+(eps_r)**2) # calculating the residuals from the model and log-likelihood: R_r = (v_r-v_r_model)/v_r_err R_l = (v_l-v_l_model)/v_l_err R_b = (v_b-v_b_model)/v_b_err lnprob = np.sum(-0.5*R_l*R_l-np.log(v_l_err)) lnprob += np.sum(-0.5*R_b*R_b-np.log(v_b_err)) lnprob += np.sum(-0.5*R_r*R_r-np.log(v_r_err)) return (lnprob+lprior) def lnprob_v2 (pars,l,b,dist,v_l,v_l_err,v_b,v_b_err,v_r,v_r_err): """ Likelihood function for model 2 Parameters: ----------- pars : Model parameters l : Galactic longitude (radians) b : Galactic latitude (radians) dist : heliocentric distance (kpc) v_l : Tangential velocity (l component) (km/s) v_l_err : Uncertainty of tangential velocity (l component) (km/s) v_b : Tangential velocity (b component) (km/s) v_b_err : Uncertainty of tangential velocity (b component) (km/s) v_rad : Radial velocity (km/s) v_rad_err : Uncertainty of radial velocity (km/s) Returns: -------- lnprob : log-likelihood """ R0,theta_0,dtheta,Us,Vs,Ws,U0,V0,W0,eps_l,eps_b,eps_r = pars if Us < -100.0: return -np.inf if Us > 100.0: return -np.inf if Vs < -100.0: return -np.inf if Vs > 100.0: return -np.inf if Ws < -100.0: return -np.inf if Ws > 100.0: return -np.inf if R0 < 7.0: return -np.inf if R0 > 9.0: return -np.inf if eps_l < 0.0: return -np.inf if eps_b < 0.0: return -np.inf if eps_r < 0.0: return -np.inf # priors on U0,V0,W0 from Schonrich et al. (2010) lprior = 0.0 lprior += -0.5 * (U0-11.1)**2 / (1.3*1.3) lprior += -0.5 * (V0-12.2)**2 / (2.1*2.1) lprior += -0.5 * (W0-7.3)**2 / (0.7*0.7) lprior += -0.5 * (R0-8.122)**2 / (0.031*0.031) # calculating projected distance from the Galactic center and angle beta # see Appendix of Reid et al. (2009) Dp = dist * np.cos(b) Rp = R0*R0 + Dp*Dp - 2.0*R0*Dp*np.cos(l) Rp = np.sqrt(Rp) sinbeta = Dp * np.sin(l) / Rp cosbeta = (R0 - Dp * np.cos(l)) / Rp beta = np.arctan2(sinbeta,cosbeta) # rotating velocity vector through the angle -beta and subtracting # velocity of the Sun (equations 1--3) theta_S = theta_0 + dtheta * (Rp - R0) U2 = Us * np.cos(beta) + (Vs + theta_S) * np.sin(beta) V2 = -Us * np.sin(beta) + (Vs + theta_S) * np.cos(beta) W2 = Ws U1 = U2 - U0 V1 = V2 - V0 - theta_0 W1 = W2 - W0 # calculating predicted tangential and radial velocities # (equations 4--6): v_l_model = V1 * np.cos(l) - U1 * np.sin(l) tmp = U1 * np.cos(l) + V1 * np.sin(l) v_b_model = W1 * np.cos(b) - tmp * np.sin(b) v_r_model = tmp * np.cos(b) + W1 * np.sin(b) # correcting the error bars v_l_err = np.sqrt(v_l_err**2+(eps_l)**2) v_b_err = np.sqrt(v_b_err**2+(eps_b)**2) v_r_err = np.sqrt(v_r_err**2+(eps_r)**2) # calculating the residuals from the model and log-likelihood: R_r = (v_r-v_r_model)/v_r_err R_l = (v_l-v_l_model)/v_l_err R_b = (v_b-v_b_model)/v_b_err lnprob = np.sum(-0.5*R_l*R_l-np.log(v_l_err)) lnprob += np.sum(-0.5*R_b*R_b-np.log(v_b_err)) lnprob += np.sum(-0.5*R_r*R_r-np.log(v_r_err)) return (lnprob+lprior) def lnprob_v3 (pars,l,b,dist,v_l,v_l_err,v_b,v_b_err,v_r,v_r_err): """ Likelihood function for model 2 Parameters: ----------- pars : Model parameters l : Galactic longitude (radians) b : Galactic latitude (radians) dist : heliocentric distance (kpc) v_l : Tangential velocity (l component) (km/s) v_l_err : Uncertainty of tangential velocity (l component) (km/s) v_b : Tangential velocity (b component) (km/s) v_b_err : Uncertainty of tangential velocity (b component) (km/s) v_rad : Radial velocity (km/s) v_rad_err : Uncertainty of radial velocity (km/s) Returns: -------- lnprob : log-likelihood """ R0,a1,a2,a3,Us,Vs,Ws,U0,V0,W0,eps_l,eps_b,eps_r = pars if Us < -100.0: return -np.inf if Us > 100.0: return -np.inf if Vs < -100.0: return -np.inf if Vs > 100.0: return -np.inf if Ws < -100.0: return -np.inf if Ws > 100.0: return -np.inf if R0 < 7.0: return -np.inf if R0 > 9.0: return -np.inf if eps_l < 0.0: return -np.inf if eps_b < 0.0: return -np.inf if eps_r < 0.0: return -np.inf # priors on U0,V0,W0 from Schonrich et al. (2010) lprior = 0.0 lprior += -0.5 * (U0-11.1)**2 / (1.3*1.3) lprior += -0.5 * (V0-12.2)**2 / (2.1*2.1) lprior += -0.5 * (W0-7.3)**2 / (0.7*0.7) lprior += -0.5 * (R0-8.122)**2 / (0.031*0.031) # calculating projected distance from the Galactic center and angle beta # see Appendix of Reid et al. (2009) Dp = dist * np.cos(b) Rp = R0*R0 + Dp*Dp - 2.0*R0*Dp*np.cos(l) Rp = np.sqrt(Rp) sinbeta = Dp * np.sin(l) / Rp cosbeta = (R0 - Dp * np.cos(l)) / Rp beta = np.arctan2(sinbeta,cosbeta) # rotating velocity vector through the angle -beta and subtracting # velocity of the Sun (equations 1--3) x = (Rp/R0) / a2 bet = 0.72 Vd = bet * 1.97 * pow(x,1.22) / pow(x**2+0.78**2,1.43) Vh = (1.0-bet) * (1.0+a3*a3) * x * x / (x*x + a3*a3) theta_S = a1 * np.sqrt(Vd+Vh) x0 = 1.0 / a2 Vd0 = bet * 1.97 * pow(x0,1.22) / pow(x0**2+0.78**2,1.43) Vh0 = (1.0-bet) * (1.0+a3*a3) * x0 * x0 / (x0*x0 + a3*a3) theta_0 = a1 * np.sqrt(Vd0+Vh0) U2 = Us * np.cos(beta) + (Vs + theta_S) * np.sin(beta) V2 = -Us * np.sin(beta) + (Vs + theta_S) * np.cos(beta) W2 = Ws U1 = U2 - U0 V1 = V2 - V0 - theta_0 W1 = W2 - W0 # calculating predicted tangential and radial velocities # (equations 4--6): v_l_model = V1 * np.cos(l) - U1 * np.sin(l) tmp = U1 * np.cos(l) + V1 * np.sin(l) v_b_model = W1 * np.cos(b) - tmp * np.sin(b) v_r_model = tmp * np.cos(b) + W1 * np.sin(b) # correcting the error bars v_l_err = np.sqrt(v_l_err**2+(eps_l)**2) v_b_err = np.sqrt(v_b_err**2+(eps_b)**2) v_r_err = np.sqrt(v_r_err**2+(eps_r)**2) # calculating the residuals from the model and log-likelihood: R_r = (v_r-v_r_model)/v_r_err R_l = (v_l-v_l_model)/v_l_err R_b = (v_b-v_b_model)/v_b_err lnprob = np.sum(-0.5*R_l*R_l-np.log(v_l_err)) lnprob += np.sum(-0.5*R_b*R_b-np.log(v_b_err)) lnprob += np.sum(-0.5*R_r*R_r-np.log(v_r_err)) return (lnprob+lprior) def professional_mcmc (best, lnprob, N, l, b, dist, v_l, v_l_err, v_b, v_b_err, v_r, v_r_err, ndim=3, nwalkers=100, burnin=100, threads=1, wersor=[0.01,0.01,0.01], save_to_file=True, root='./', name='results'): """ MCMC fitting using the Affine Invariant Markov chain Monte Carlo Ensemble sampler: http://adsabs.harvard.edu/abs/2013PASP..125..306F Parameters: ----------- best : starting parameters for the MCMC lnprob : log-likelihood function N : number of MCMC steps dist, v_l, v_l_err, v_b, v_b_err, v_r, v_r_err : input data ndim : number of dimensions nwalkers : number of MCMC walkers burnin : number of burn-in steps threads : number of threads (for parallelization) wersor : defines the shape of initial distributions save_to_file : if True, saves the MCMC chain to file $root$name.mcmc """ if len(best) != ndim: raise ValueError('Error: Vector best has incorrect lenght.') if len(wersor) != ndim: raise ValueError('Error: Vector wersor has incorrect lenght.') # Starting points for each walker are drawn from a multivariate # normal distribution with variances defined by the wersor vector. wersor = np.array(wersor) scale = np.diag(wersor**2) start = np.array([np.random.multivariate_normal(best, scale) for i in xrange(nwalkers)]) maxcpu = mp.cpu_count() if threads <=0 or threads > maxcpu: threads = maxcpu - 2 sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(l, b, dist, v_l, v_l_err, v_b, v_b_err, v_r, v_r_err), threads=threads) # Burn-in steps pos, prob, state = sampler.run_mcmc(start, burnin) sampler.reset() # Run MCMC chain sampler.run_mcmc(pos, N) samples = sampler.flatchain chi2_vals = (-2.0)*sampler.flatlnprobability # Save MCMC chain if save_to_file: output = root+name+'.mcmc' f = open(output, 'w') for ss, cc in zip(samples, chi2_vals): for x in ss: f.write('%f ' % x) f.write('%f\n' % cc) f.close() # Print output to STD chi2_min = min(chi2_vals) filtr = (chi2_vals <= (chi2_min + 25.0)) # print chi2_min for i in xrange(ndim): lo, med, up = np.percentile(samples[:,i][filtr], (16,50,84)) print 'Parameter %d: %f +%f / -%f' % (i+1,med,up-med,med-lo) return 0 def main(): desc = """ Modeling rotation curve of the Milky Way based on proper motions and radial velocities of Classical Cepheids. Usage: fit_model_public.py model model = 1 : flat rotation curve \Theta(R) = \Theta_0 model = 2 : polynomial rotation curve \Theta(R) = \Theta_0 + a1*(R-R_0) model = 3 : universal rotation curve (Persic et al. 1996) """ parser = argparse.ArgumentParser(description=desc) parser.add_argument('model', default=0, type=int, help='model number') args = parser.parse_args() if args.model < 1 or args.model > 3: print desc #-----Reading data from file---------------------------------------# filename = 'Cepheids.dat' dtype = 'S20,int,float,float,float,float,float,float,float,float, float,float,float,float,float' try: ident, flag, ra_ogle, dec_ogle, glon_ogle, glat_ogle, dist, dist_err, pm_ra, pm_ra_err, pm_dec, pm_dec_err, pm_corr, v_r, v_r_err = np.loadtxt(filename, unpack=True, dtype=dtype) except IOError: exit('Error while reading %s. Exiting.'%filename) N = len(ident) #-----Preparing data-----------------------------------------------# deg = np.pi / 180.0 # degrees to radians dec_G = 27.12825 * deg # declination of North Galactic pole ra_G = 192.85948 * deg # right ascension of North Galactic pole l = glon_ogle * deg b = glat_ogle * deg ra = ra_ogle * deg dec = dec_ogle * deg #-----Calculating tangential velocities----------------------------# v_l, v_l_err, v_b, v_b_err, v_lb_corr = [], [], [], [], [] for i in xrange(N): _v_l, _v_l_err, _v_b, _v_b_err, _v_lb_corr = calculate_tangential_velocity (ra[i],dec[i],l[i],b[i],dist[i],pm_ra[i],pm_ra_err[i],pm_dec[i],pm_dec_err[i],pm_corr[i],ra_G,dec_G) v_l.append(_v_l) v_l_err.append(_v_l_err) v_b.append(_v_b) v_b_err.append(_v_b_err) v_lb_corr.append(_v_lb_corr) v_l_err = np.array(v_l_err) v_b_err = np.array(v_b_err) v_b = np.array(v_b) v_l = np.array(v_l) v_lb_corr = np.array(v_lb_corr) #-----Running MCMC-------------------------------------------------# np.random.seed() filtr = (flag == 1) model = args.model if model == 1: # model 1: Theta(R) = const pars0 = (8.1,230.0,1.97,-4.99,-0.05,11.1,12.2,7.3,10.0,6.0,10.0) sigma = (0.1,10.0,2.0,2.0,2.0,2.0,2.0,2.0,0.5,0.5,0.5) num = 5000 burnin = 2000 ndim = 11 professional_mcmc(pars0,lnprob_v1,num,l[filtr],b[filtr],dist[filtr],v_l[filtr],v_l_err[filtr],v_b[filtr],v_b_err[filtr],v_r[filtr],v_r_err[filtr],burnin=burnin,wersor=sigma,ndim=ndim,save_to_file=True) elif model == 2: # model 2: Theta(R) = Theta_0 + dTheta/dR (R-R_0) pars0 = (8.1,230.0,0.0,1.97,-4.99,-0.05,11.1,12.2,7.3,10.0,6.0,10.0) sigma = (0.1,10.0,5.0,2.0,2.0,2.0,2.0,2.0,2.0,0.5,0.5,0.5) num = 1000 burnin = 1000 ndim = 12 professional_mcmc(pars0,lnprob_v2,num,l[filtr],b[filtr],dist[filtr],v_l[filtr],v_l_err[filtr],v_b[filtr],v_b_err[filtr],v_r[filtr],v_r_err[filtr],burnin=burnin,wersor=sigma,ndim=ndim,save_to_file=True) elif model == 3: # model 3: universal rotation curve (Persic et al. 1996) pars0 = (8.1,230.0,1.0,1.5,1.97,-4.99,-0.05,11.1,12.2,7.3,10.0,6.0,10.0) sigma = (0.1,10.0,0.1,0.1,2.0,2.0,2.0,2.0,2.0,2.0,0.5,0.5,0.5) num = 5000 burnin = 2000 ndim = 13 professional_mcmc(pars0,lnprob_v3,num,l[filtr],b[filtr],dist[filtr],v_l[filtr],v_l_err[filtr],v_b[filtr],v_b_err[filtr],v_r[filtr],v_r_err[filtr],burnin=burnin,wersor=sigma,ndim=ndim,save_to_file=True) if __name__ == '__main__': main()