From 9129541bb89950a45b89ff83862a1ad0f8ddd709 Mon Sep 17 00:00:00 2001 From: "Gijs Mulders (Laptop)" Date: Sun, 6 Oct 2019 15:38:21 -0500 Subject: [PATCH] transit duration, impact parameter, correlated spacings --- EPOS/__init__.py | 2 +- EPOS/classes.py | 9 ++++-- EPOS/kepler.py | 9 ++++-- EPOS/plot/multi.py | 2 +- EPOS/run.py | 69 ++++++++++++++++++++++++++++++++++++++-------- 5 files changed, 73 insertions(+), 18 deletions(-) diff --git a/EPOS/__init__.py b/EPOS/__init__.py index 69d3b06..b8e0342 100644 --- a/EPOS/__init__.py +++ b/EPOS/__init__.py @@ -1,5 +1,5 @@ __shortversion__= u'3.0' -__version__= u'3.0.1.dev0' # [N!]N(.N)*[{a|b|rc}N][.postN][.devN] +__version__= u'3.0.1.dev1' # [N!]N(.N)*[{a|b|rc}N][.postN][.devN] __all__ = ['epos','fitparameters','kepler','rv','run','population','plot','occurrence', diff --git a/EPOS/classes.py b/EPOS/classes.py index bd1bde9..4770474 100644 --- a/EPOS/classes.py +++ b/EPOS/classes.py @@ -223,7 +223,7 @@ def __init__(self, name, Debug=False, seed=True, title=None, print ('Survey: None selected') def set_observation(self, xvar, yvar, starID, nstars=1.6862e5, - radiusError=0.1, score=None, Verbose=True): + radiusError=0.1, score=None, tdur=None, Verbose=True): ''' Observed planet population Args: @@ -244,6 +244,9 @@ def set_observation(self, xvar, yvar, starID, nstars=1.6862e5, if score is not None: self.obs_score= np.asarray(score)[order] + + if tdur is not None: + self.obs_tdur= np.asarray(tdur)[order] assert self.obs_xvar.ndim == self.obs_yvar.ndim == self.obs_starID.ndim == 1, 'only 1D arrays' assert self.obs_xvar.size == self.obs_yvar.size == self.obs_starID.size, 'arrays not same length' @@ -634,7 +637,7 @@ def set_parametric(self, func): self.fitpars=self.pdfpars self.summarystatistic= ['N','xvar','yvar'] - def set_multi(self, spacing=None): + def set_multi(self, spacing=None, Correlated=False): if not self.Parametric: raise ValueError('Define a parametric planet population first') self.Multi=True @@ -642,6 +645,8 @@ def set_multi(self, spacing=None): self.RandomPairing= (spacing==None) self.spacing= spacing # None, brokenpowerlaw, dimensionless + self.Correlated=Correlated + # skipping yvar here self.summarystatistic= ['N','xvar','Nk','dP','Pin'] # dR? diff --git a/EPOS/kepler.py b/EPOS/kepler.py index d2ff01e..0be3309 100644 --- a/EPOS/kepler.py +++ b/EPOS/kepler.py @@ -58,7 +58,7 @@ def dr25(subsample='all', score=0.9, Gaia=False, Huber=True, Vetting=False, Verb if Verbose: print (' removed {} planets with no stellar radii'.format(nremove)) koi= {key: np.array(ipac[key][nonzero]) for key in ['kepid','koi_prad','koi_period', - 'koi_steff', 'koi_slogg', 'koi_srad', 'koi_depth', + 'koi_steff', 'koi_slogg', 'koi_srad', 'koi_depth', 'koi_duration', 'koi_pdisposition','koi_score'] } if Gaia: @@ -160,9 +160,10 @@ def dr25(subsample='all', score=0.9, Gaia=False, Huber=True, Vetting=False, Verb #if Huber: slice['subgiants']=issubgiant obs= {'xvar':koi['koi_period'][slice], - 'yvar':koi['koi_prad'][slice], + 'yvar':koi['koi_prad'][slice], 'starID':koi['kepid'][slice], - 'score':koi['koi_score'][slice]} + 'score':koi['koi_score'][slice], + 'tdur':koi['koi_duration'][slice]} ''' Load pre-calculated detection efficiencies ''' # from dr25_epos.py @@ -277,6 +278,8 @@ def vetting_parameters(score=0.9, subsample='all', Gaia=True): else: raise ValueError(errormessage) + return vetpars + ''' NOTE: Kepler files from Q16-epos.py diff --git a/EPOS/plot/multi.py b/EPOS/plot/multi.py index 4ba0dcf..26b97b6 100644 --- a/EPOS/plot/multi.py +++ b/EPOS/plot/multi.py @@ -559,7 +559,7 @@ def radiusratio(epos, MC=False, Input=False, MCMC=False, color='C1', NB=False): # plot extra epos runs for ss in epos.ss_extra: - ax.hist(ss['multi']['Pratio'], bins=bins, + ax.hist(ss['multi']['Rratio'], bins=bins, histtype='step', label=ss['name']) # Observed zoom diff --git a/EPOS/run.py b/EPOS/run.py index ab64270..c693cbc 100644 --- a/EPOS/run.py +++ b/EPOS/run.py @@ -264,8 +264,9 @@ def prep_obs(epos): x= epos.obs_xvar[ix&iy] y= epos.obs_yvar[ix&iy] + t= epos.obs_tdur[ix&iy] - for key, var in zip(['x','y'], [x,y]): + for key, var in zip(['x','y','t'], [x,y,t]): z[key]= np.sort(var) # add (x0,0) and (x1,1)? z[key+' cum']= np.arange(z[key].size,dtype=float) z[key+' cdf']= z[key+' cum']/z[key+' cum'][-1] @@ -370,7 +371,7 @@ def MC(epos, fpara, Store=False, Sample=False, StorePopulation=False, Extra=None ''' Draw multiplanet distributions ''' allX, allY, allI, allN, allID= \ - draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store) + draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store, epos.Correlated) ''' convert to observable parameters ''' allP= allX @@ -469,7 +470,8 @@ def MC(epos, fpara, Store=False, Sample=False, StorePopulation=False, Extra=None itrans= p_trans >= np.random.uniform(0,1,allP.size) else: #multi-transit probability - itrans= istransit(epos, allID, allI, allP, f_iso, f_inc, Verbose=Verbose) + itrans, allB, allTd= \ + istransit(epos, allID, allI, allP, allY, f_iso, f_inc, Verbose=Verbose) # print (multi statistics if Verbose and epos.Multi and not epos.RV: @@ -508,6 +510,9 @@ def MC(epos, fpara, Store=False, Sample=False, StorePopulation=False, Extra=None ''' uncertainty in stellar radius? ''' MC_Y=MC_R * (1.+epos.radiusError*np.random.normal(size=MC_R.size) ) + if epos.Multi and (dInc!=None): + MC_B=allB[itrans] + MC_Td=allTd[itrans] ''' Store (transiting) planet sample for verification plot @@ -676,12 +681,16 @@ def MC(epos, fpara, Store=False, Sample=False, StorePopulation=False, Extra=None pop['ID']= allID pop['k']= allN pop['inc']= allI + pop['b']= allB + pop['tdur']= allTd for key, subset in zip(['system', 'single', 'multi'],[isysdet, isingle, imulti]): pop[key]={} pop[key]['Y']= allY[subset] # R? M? pop[key]['P']= allP[subset] pop[key]['ID']= allID[subset] pop[key]['inc']= allI[subset] + pop[key]['b']= allB[subset] + pop[key]['tdur']= allTd[subset] pop[key]['order']= order[subset] pop[key]['detectable']= itransdet[subset] @@ -706,7 +715,12 @@ def MC(epos, fpara, Store=False, Sample=False, StorePopulation=False, Extra=None ss['multi']['cdf']= multi.cdf(det_ID[ix&iy]) if not epos.RV and epos.Parametric: ss['multi']['PN'],ss['multi']['dPN']= multi.periodratio( - det_ID[ix&iy], det_P[ix&iy], N=det_N[ix&iy]) + det_ID[ix&iy], det_P[ix&iy], N=det_N[ix&iy]) + + # impact, transit duration + if (dInc!=None): + ss['b']= MC_B[idet] + ss['tdur']= MC_Td[idet] epos.prob=prob epos.lnprob=lnprob @@ -902,7 +916,7 @@ def draw_from_2D_distribution(epos, pps, fpara, npl=1): return allX, allY -def draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store): +def draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store, Correlated=False): ''' assign ID to each system ''' sysID= np.arange(sysX.size) #allID= np.repeat(sysID, npl) # array with star ID for each planet @@ -920,7 +934,15 @@ def draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store): _, toplanet, sysnpl= np.unique(allID, return_inverse=True,return_counts=True) allX= sysX[toplanet] allY= sysY[toplanet] - allI= np.random.rayleigh(dInc, allID.size) + + # correlated mutual inc, spacing, size + if Correlated: + xcor= np.random.uniform(0,1,allID.size) + grid_inc= np.geomspace(1e-3,90.) + cdf_inc= 1.-np.exp(-.5*(grid_inc/dInc)**2.) + allI= np.interp(xcor, cdf_inc, grid_inc) + else: + allI= np.random.rayleigh(dInc, allID.size) #allN= np.ones_like(allID) # index to planet in system allN= np.where(allX>=epos.xzoom[0],1,0) # also yzoom? #print (allX[:3] @@ -964,12 +986,26 @@ def draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store): im= i1[sysnpl>=i] # index to ith planet in each system #print (' planet {} in system, {} systems'.format(i,im.size) #dP= draw_from_function(brokenpowerlaw1D,xgrid,im.size, dPbreak, dP1, dP2) - dP= np.interp(np.random.uniform(cdf[0],cdf[-1],im.size), cdf, Pgrid) + + if Correlated: + dP_all= np.interp(cdf[0]+xcor*(cdf[-1]-cdf[0]), cdf, Pgrid) + dP= dP_all[im+(i-1)] + else: + dP= np.interp(np.random.uniform(cdf[0],cdf[-1],im.size), cdf, Pgrid) allX[im+(i-1)]= allX[im+(i-2)]* dP + #np.random.uniform(dP-0.7, dP+0.7, im.size) #print (allX[:3] #np.random.norm() - allY[im+(i-1)]= allY[im+(i-2)]* 10.**np.random.normal(0, dR, im.size) + + if Correlated: + grid_dR= np.linspace(-5.*dR,5*dR) # 5 stdev range + cdf_dR= norm(0,dR).cdf(grid_dR) + dR_all= np.interp(xcor, cdf_dR, grid_dR) + #allY[im+(i-1)]= allY[im]* 10.**dR_all[im+(i-1)] + allY[im+(i-1)]= allY[im+(i-2)]* 10.**dR_all[im+(i-1)] + else: + allY[im+(i-1)]= allY[im+(i-2)]* 10.**np.random.normal(0, dR, im.size) # nth planet in system (zoom range) #allN[im+(i-1)]= allN[im+(i-2)]+1 @@ -992,7 +1028,7 @@ def draw_multi(epos, sysX, sysY, npl, dInc, dR, fpara, Store): return allX, allY, allI, allN, allID -def istransit(epos, allID, allI, allP, f_iso, f_inc, Verbose=False): +def istransit(epos, allID, allI, allP, allR, f_iso, f_inc, Verbose=False): # draw same numbers for multi-planet systems IDsys, toplanet= np.unique(allID, return_inverse=True) # return_counts=True if Verbose: print (' {}/{} systems'.format(IDsys.size, allID.size)) @@ -1003,6 +1039,7 @@ def istransit(epos, allID, allI, allP, f_iso, f_inc, Verbose=False): assert inc_pl.size == allP.size R_a= epos.fgeo_prefac *allP**epos.Pindex # == p_trans + mutual_inc= allI * f_inc #mutual_inc= 0.0 # planar distribution #mutual_inc= 1.0 # fit @@ -1013,16 +1050,26 @@ def istransit(epos, allID, allI, allP, f_iso, f_inc, Verbose=False): delta_inc= mutual_inc *np.cos(np.random.uniform(0,np.pi,allP.size)) *np.pi/180. #rad itrans= np.abs(inc_pl+delta_inc) < np.arcsin(R_a) + # transit duration, impact + impact= np.abs(np.sin((inc_pl+delta_inc))/R_a) + # allow for a fraction of isotropic systems if f_iso > 0: p_trans= epos.fgeo_prefac *allP**epos.Pindex - itrans_iso= p_trans >= np.random.uniform(0,1,allP.size) + rand_fgeo= np.random.uniform(0,1,allP.size) + itrans_iso= p_trans >= rand_fgeo iso_sys= (np.random.uniform(0,1,IDsys.size) < f_iso) iso_pl= iso_sys[toplanet] itrans = np.where(iso_pl, itrans_iso, itrans) + + impact_iso= np.abs(rand_fgeo/R_a) + impact = np.where(iso_pl, impact_iso,impact) + + t_d= allP*(cgs.day/cgs.hour)/np.pi* \ + np.arcsin(R_a*np.sqrt(np.maximum((1.+allR*cgs.Rearth/cgs.Rsun)**2.-impact**2.,0))) - return itrans + return itrans, impact, t_d def storepopulation(allID, allP, det_ID, idetected): # add f_iso?