planet formation model updates and notebook

EPOS/classes.py

@@ -63,10 +63,10 @@ def add(self, key, value, fixed=False, min=-np.inf, max=np.inf,
 			dx=0.1*value if dx is None else dx
 			fp['dx']=abs(dx) if (dx!=0) else 0.1
 
-	def default(self, key, value, Verbose=True):
+	def default(self, key, value, isnorm=False, Verbose=True):
 		if not key in self.keysall: 
 			if Verbose: print ('  Set {} to default {}'.format(key, value))
-			self.add(key, value, fixed=True)
+			self.add(key, value, fixed=True, isnorm=isnorm)
 
 	def set(self, key, value):
 		self.fitpars[key]['value_init']=value
@@ -209,7 +209,7 @@ def __init__(self, name, Debug=False, seed=True, title=None,
 
 		''' Load Default Settings for each Survey '''
 		if survey is 'Kepler':
-			print ('\nSurvey: Kepler-Gaia all dwarfs')
+			print ('\nSurvey: Kepler-Gaia all dwarfs, reliability > 0.9')
 			obs, survey= EPOS.kepler.dr25(Huber=False, Gaia=True, Vetting=True, score=0.9, 
 				Verbose=False)
 			self.set_observation(Verbose=False, radiusError=0.02, **obs)
@@ -314,7 +314,8 @@ def set_survey(self, xvar, yvar, eff_2D, Rstar=1.0, Mstar=1.0, vet_2D=None):
 			self.fgeo_prefac= self.Rstar*cgs.Rsun * fourpi2_GM**(1./3.) / cgs.day**(2./3.)
 			P, R= np.meshgrid(self.eff_xvar, self.eff_yvar, indexing='ij')
 
-			self.completeness= self.eff_2D * self.fgeo_prefac*P**self.Pindex
+			self.fgeo= self.fgeo_prefac*P**self.Pindex
+			self.completeness= self.eff_2D * self.fgeo
 
 		if vet_2D is not None:
 			self.vetting= np.asarray(vet_2D)
@@ -718,6 +719,7 @@ def set_population(self, name, sma, mass,
 		if pfm['np'] > pfm['ns']:
 
 			order= np.lexsort((pfm['sma'],pfm['ID'])) # sort by ID, then sma
+			pfm['order']=order
 			for key in ['ID','sma','M','P','inc','ecc','tag','R']:
 				if key in pfm: pfm[key]=pfm[key][order]
 
@@ -761,6 +763,8 @@ def set_population(self, name, sma, mass,
 			if 'tag' in pfm:
 				pfm['system tag']= pfm['tag']
 
+			pfm['order']= np.arange(len(sma))
+
 		pfm['M limits']=[np.min(pfm['M']),np.max(pfm['M'])]
 		pfm['P limits']=[np.min(pfm['P']),np.max(pfm['P'])]
 

EPOS/occurrence.py

@@ -11,7 +11,7 @@
 from EPOS.population import periodradius
 import EPOS.analytics
 
-def all(epos, BinnedMetric=False, MLE=False):
+def all(epos, BinnedMetric=False, MLE=False, Verbose=True):
 	if hasattr(epos,'occurrence'):
 		planets(epos)
 
@@ -86,12 +86,17 @@ def models(epos, Log=False):
 		#print completeness
 	else:
 		# if zeros in detection efficiency?
-		pl_comp= interpolate.RectBivariateSpline(epos.eff_xvar, 
-			epos.eff_yvar, epos.completeness)
-		completeness= pl_comp(pfm['P'], pfm['R'], grid=False)
+		pl_ftot= interpolate.RectBivariateSpline(epos.eff_xvar, epos.eff_yvar, epos.completeness)
+		pl_fgeo= interpolate.RectBivariateSpline(epos.eff_xvar, epos.eff_yvar, epos.fgeo)
+		pl_fdet= interpolate.RectBivariateSpline(epos.eff_xvar, epos.eff_yvar, epos.eff_2D)
+		pl_fvet= interpolate.RectBivariateSpline(epos.eff_xvar, epos.eff_yvar, epos.vetting)
 
 	focc['model']={}
-	focc['model']['completeness']= completeness
+	focc['model']['completeness']= pl_ftot(pfm['P'], pfm['R'], grid=False)
+	focc['model']['fgeo']= pl_fgeo(pfm['P'], pfm['R'], grid=False)
+	focc['model']['eff_2D']= pl_fdet(pfm['P'], pfm['R'], grid=False)
+	focc['model']['vetting']= pl_fvet(pfm['P'], pfm['R'], grid=False)
+
 	focc['model']['eta']= epos.fitpars.getpps() if hasattr(epos, 'fitpars') else 1.
 	#focc['model']['occ']= focc['model']['eta']/completeness/epos.nstars
 	#print epos.planet_occurrence

EPOS/pfmodel.py

@@ -8,7 +8,7 @@
 ''' Helper functions to read in planet formation models'''
 
 #def symba(name='HMSim1', dir='hdf5/Sim1', plts_mass=0, istep=None, Verbose=False):
-def symba(name, fname, plts_mass=0, cut=-np.inf, smacut=np.inf, istep=None, 
+def symba(name, fname, plts_mass=0, cut=None, smacut=None, istep=None, 
 	Verbose=False, Saved=True):
 	''' 
 	returns a list of planetary systems
@@ -17,11 +17,20 @@ def symba(name, fname, plts_mass=0, cut=-np.inf, smacut=np.inf, istep=None,
 	remove planetesimals (m < plt_mass)
 	cut planets that are still forming (m < cut*au)
 	'''
-	
-	dir= 'npz/{}'.format(name)
+
+	dir= 'npz'
 	if not os.path.exists(dir): os.makedirs(dir)
-	fnpz= '{}/{}.npz'.format(dir, name)	
-
+
+	print (cut, smacut)
+	if cut is None and smacut is None:
+		cut=-np.inf
+		smacut=np.inf
+		fnpz= '{}/{}.npz'.format(dir, name)	
+	else:
+		if cut is None: cut=-np.inf
+		if smacut is None: smacut=np.inf
+		fnpz= '{}/{}-cut.npz'.format(dir, name)	
+
 	''' Load hdf5 file or npz dictionary for quicker access'''
 	if os.path.isfile(fnpz) and Saved:
 		print ('\nLoading saved status from {}'.format(fnpz))

EPOS/run.py

@@ -63,9 +63,9 @@ def once(epos, fac=1.0, Extra=None, goftype='KS', Verbose=False):
 
 		else:
 			# set defaults for planet formation models here
-			epos.fitpars.default('eta',0.5)
+			epos.fitpars.default('eta',0.3, isnorm=True)
 			epos.fitpars.default('f_cor',0.5)
-			epos.fitpars.default('f_iso',0.5)
+			epos.fitpars.default('f_iso',0.5) # 1.0
 			epos.fitpars.default('f_inc',1.0)
 			epos.fitpars.default('f_dP',1.0)
 

docs/tutorials/generate_population.ipynb

@@ -419,7 +419,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.7.4"
   }
  },
  "nbformat": 4,