Unify noise-model and fix more axes in plots

dataval/mag2flux.py

@@ -27,14 +27,20 @@ def mag2flux(dval):
 	logger = logging.getLogger('dataval')
 	logger.info('Plotting Magnitude to Flux conversion...')
 
-	zp_grid = np.linspace(19.0, 22.0, 100)
-	mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 200)
+	# Get the maximum median flux to be used to set the plot limits:
+	max_mean_flux = dval.search_database(
+		select=['MAX(mean_flux) AS maxflux'],
+		search=["method_used='aperture'"])
+	max_mean_flux = max_mean_flux[0]['maxflux']
+
+	zp_grid = np.linspace(19.0, 22.0, 300)
+	mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 300)
 
 	# The interpolation is linear *in log*
 	xmin = np.array([0, 1.5, 9, 12.6, 13, 14, 15, 16, 17, 18, 19])
 	ymin = np.array([8e7, 1.8e7, 12500, 250, 59, 5, 1, 1, 1, 1, 1])
 	min_bound_log = InterpolatedUnivariateSpline(xmin, np.log10(ymin), k=1, ext=3)
-	min_bound = lambda x: 10**(min_bound_log(x)) # noqa: E731
+	min_bound = lambda x: np.clip(10**(min_bound_log(x)), 1, None) # noqa: E731
 
 	norm = colors.Normalize(vmin=0, vmax=1)
 	fig2, ax2 = plt.subplots()
@@ -125,13 +131,14 @@ def chi2(c):
 		#ax31.plot(bin_centers2, 1.4826*3*bin_means2, color='g')
 
 		# Add line with best fit, and the minimim bound:
-		ax1.plot(mags, 10**(-0.4*(mags - cc.x)), color='k', ls='--')
+		ax1.plot(mags, np.clip(10**(-0.4*(mags - cc.x)), 0, None), color='k', ls='--')
 		ax1.plot(mags, min_bound(mags), 'r-')
 
 		ax1.set_yscale('log')
 		ax1.set_xlim(dval.tmag_limits[1], dval.tmag_limits[0])
+		ax1.set_ylim(0.5, 2*max_mean_flux)
 		ax1.set_xlabel('TESS magnitude')
-		ax1.set_ylabel('Median flux')
+		ax1.set_ylabel('Median flux (e$^-$/s)')
 		ax1.xaxis.set_major_locator(MultipleLocator(2))
 		ax1.xaxis.set_minor_locator(MultipleLocator(1))
 		colorbar(im, ax=ax1, label='Contamination')

dataval/noise_metrics.py

@@ -28,7 +28,15 @@ def noise_metrics(dval):
 	logger = logging.getLogger('dataval')
 	logger.info('Plotting Noise vs. Magnitude...')
 
-	mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 200)
+	table_noise = 'diagnostics_corr' if dval.corr else 'diagnostics'
+	factor = 1 if dval.corr else 1e6
+	max_vals = dval.search_database(select=[
+		'MAX(' + table_noise + '.rms_hour) AS max_rms',
+		'MAX(' + table_noise + '.ptp) AS max_ptp'])
+	max_rms = factor*max_vals[0]['max_rms']
+	max_ptp = factor*max_vals[0]['max_ptp']
+
+	mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 300)
 
 	# Colors for theoretical lines:
 	# Define the colors directly here to avoid having to import
@@ -45,30 +53,16 @@ def noise_metrics(dval):
 
 	for cadence in dval.cadences:
 
-		if dval.corr:
-			factor = 1
-			star_vals = dval.search_database(
-				select=[
-					'todolist.priority',
-					'todolist.sector',
-					'todolist.tmag',
-					'diagnostics_corr.rms_hour',
-					'diagnostics_corr.ptp',
-					'diagnostics.contamination'
-				],
-				search=f'cadence={cadence:d}')
-		else:
-			factor = 1e6
-			star_vals = dval.search_database(
-				select=[
-					'todolist.priority',
-					'todolist.sector',
-					'todolist.tmag',
-					'diagnostics.rms_hour',
-					'diagnostics.ptp',
-					'diagnostics.contamination'
-				],
-				search=f'cadence={cadence:d}')
+		star_vals = dval.search_database(
+			select=[
+				'todolist.priority',
+				'todolist.sector',
+				'todolist.tmag',
+				table_noise + '.rms_hour',
+				table_noise + '.ptp',
+				'diagnostics.contamination'
+			],
+			search=f'cadence={cadence:d}')
 
 		tmags = np.array([star['tmag'] for star in star_vals], dtype='float64')
 		pri = np.array([star['priority'] for star in star_vals], dtype='int64')
@@ -130,6 +124,9 @@ def noise_metrics(dval):
 			axx.set_yscale('log')
 			#axx.legend(loc='upper left')
 
+		ax1.set_ylim(0.3*np.min(tot_noise_rms), 2*max_rms)
+		ax2.set_ylim(0.3*np.min(tot_noise_ptp), 2*max_ptp)
+
 		colorbar(im1, ax=ax1, label='Contamination')
 		colorbar(im2, ax=ax2, label='Contamination')
 

dataval/noise_model.py

@@ -10,6 +10,7 @@
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 import scipy.interpolate as INT
+from .utilities import mag2flux
 
 #--------------------------------------------------------------------------------------------------
 def ZLnoise(gal_lat):
@@ -55,28 +56,6 @@ def Pixinaperture(Tmag, cad=1800):
 
 	return np.asarray(np.maximum(pixels, 3), dtype='int32')
 
-#--------------------------------------------------------------------------------------------------
-def mean_flux_level(Tmag):
-	"""
-	Mean flux from TESS magnitude
-
-	.. codeauthor:: Mikkel N. Lund <mikkelnl@phys.au.dk>
-	"""
-	# Magnitude system based on Sullivan et al.
-	#collecting_area = np.pi*(10.5/2)**2 # square cm
-	#Teff_list = np.array([2450, 3000, 3200, 3400, 3700, 4100, 4500, 5000, 5777, 6500, 7200, 9700]) # Based on Sullivan
-	#Flux_list = np.array([2.38, 1.43, 1.40, 1.38, 1.39, 1.41, 1.43, 1.45, 1.45, 1.48, 1.48, 1.56])*1e6 # photons per sec; Based on Sullivan
-	#Magn_list = np.array([306, -191, -202, -201, -174, -132, -101, -80, -69.5, -40, -34.1, 35])*1e-3 #Ic-Tmag (mmag)
-
-	#Flux_int = INT.UnivariateSpline(Teff_list, Flux_list, k=1, s=0)
-	#Magn_int = INT.UnivariateSpline(Teff_list, Magn_list, k=1, s=0)
-
-	#Imag = Magn_int(Teff)+Tmag
-	#Flux = 10**(-0.4*Imag) * Flux_int(Teff) * collecting_area
-
-	Flux = 10**(-0.4*(Tmag - 20.54))
-	return Flux
-
 #--------------------------------------------------------------------------------------------------
 def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix=1800):
 	"""
@@ -136,7 +115,7 @@ def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix=
 	Flux_factor = np.sqrt(integrations * pixels)
 
 	# Mean flux level in electrons per cadence
-	mean_level_ppm = mean_flux_level(Tmag) * timescale # electrons (based on measurement) #, Teff
+	mean_level_ppm = mag2flux(Tmag) * timescale # electrons (based on measurement) #, Teff
 
 	# Shot noise
 	shot_noise = 1e6/np.sqrt(mean_level_ppm)
@@ -152,8 +131,10 @@ def phot_noise(Tmag, timescale=3600, coord=None, sysnoise=60, Teff=5775, cadpix=
 
 	# Put individual components together in single table:
 	noise_vals = np.column_stack((shot_noise, zodiacal_noise, read_noise, systematic_noise))
+	noise_vals = np.clip(noise_vals, 0, None)
 
 	# Calculate the total noise model by adding up the individual contributions:
 	total_noise = np.sqrt(np.sum(noise_vals**2, axis=1))
+	total_noise = np.clip(total_noise, 0, None)
 
 	return total_noise, noise_vals # ppm per cadence

dataval/pixinaperture.py

@@ -26,7 +26,13 @@ def pixinaperture(dval):
 	logger = logging.getLogger('dataval')
 	logger.info('Plotting Pixels in aperture vs. Magnitude...')
 
-	mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 200)
+	mags = np.linspace(dval.tmag_limits[0], dval.tmag_limits[1], 500)
+
+	# Find the maximim mask size, used as a constant upper limit on the plots:
+	max_masksize = dval.search_database(
+		select=['MAX(diagnostics.mask_size) AS maxsize'],
+		search=["method_used='aperture'", 'diagnostics.mask_size IS NOT NULL'])
+	max_masksize = max_masksize[0]['maxsize']
 
 	# TODO: Does this need to be separated by cadence as well?
 	#       Are there significant differences beween 120s and 20s mask sizes?
@@ -130,26 +136,33 @@ def pixinaperture(dval):
 
 		# The interpolations are linear *in log space* which is the reason
 		# for the stange lambda wrappers in the following code:
+		# Rounding to 13 decimal places is to avoid numerical noise causing single pixel jumps
 		if datasource == 'ffi':
 			# Minimum bound on FFI data:
 			xmin = np.array([0, 2, 2.7, 3.55, 4.2, 4.8, 5.5, 6.8, 7.6, 8.4, 9.1, 10, 10.5, 11, 11.5, 11.6, 16])
 			ymin = np.array([30, 30, 30, 30, 29, 28, 27, 26, 25, 22, 20, 15, 11.15, 8, 5, 4, 4])
 			min_bound_log = InterpolatedUnivariateSpline(xmin, np.log10(ymin), k=1, ext=3)
-			min_bound = lambda x: np.floor(10**(min_bound_log(x))) # noqa: E731
+			min_bound = lambda x: np.clip(np.floor(np.round(10**(min_bound_log(x)), decimals=13)), 4, None) # noqa: E731
 
 			# Maximum bound on FFI data:
 			xmax = np.array([0, 2, 2.7, 3.55, 4.2, 4.8, 5.5, 6.8, 7.6, 8.4, 9.1, 10, 10.5, 11, 11.5, 12, 12.7, 13.3, 14, 14.5, 15, 16])
 			ymax = np.array([10000, 3200, 2400, 1400, 1200, 900, 800, 470, 260, 200, 170, 130, 120, 100, 94, 86, 76, 67, 59, 54, 50, 50])
 			max_bound_log = InterpolatedUnivariateSpline(xmax, np.log10(ymax), k=1, ext=3)
-			max_bound = lambda x: np.ceil(10**(max_bound_log(x))) # noqa: E731
+			max_bound = lambda x: np.ceil(np.round(10**(max_bound_log(x)), decimals=13)) # noqa: E731
 		else:
 			# Minimum bound on TPF data
 			xmin = np.array([0, 2, 2.7, 3.55, 4.2, 4.8, 5.5, 6.8, 7.6, 8.4, 9.1, 10, 10.5, 11, 11.5, 11.6, 16, 19])
 			ymin = np.array([220, 200, 130, 70, 55, 43, 36, 30, 27, 22, 16, 10, 8, 6, 5, 4, 4, 4])
 			min_bound_log = InterpolatedUnivariateSpline(xmin, np.log10(ymin), k=1, ext=3)
-			min_bound = lambda x: np.floor(10**(min_bound_log(x))) # noqa: E731
+			min_bound = lambda x: np.clip(np.floor(np.round(10**(min_bound_log(x)), decimals=13)), 4, None) # noqa: E731
 			max_bound = None
 
+		# Simple sanity check if the bounds are monotonical decresing:
+		if not np.all(np.diff(min_bound(mags)) <= 0):
+			raise RuntimeError("Minimum bound (" + datasource + ") is not monotonically decreasing") # pragma: no cover
+		if max_bound is not None and not np.all(np.diff(max_bound(mags)) <= 0):
+			raise RuntimeError("Maximum bound (" + datasource + ") is not monotonically decreasing") # pragma: no cover
+
 		# Plot limits:
 		ax1.plot(mags, min_bound(mags), 'r-')
 		if max_bound is not None:
@@ -208,8 +221,7 @@ def pixinaperture(dval):
 		#ax2.plot(mags, pix, color='k', ls='-')
 
 		ax1.set_xlim(dval.tmag_limits)
-		ax1.set_ylim([0.99, np.nanmax(masksizes)+500])
-
+		ax1.set_ylim([0.99, max_masksize+500])
 		ax1.set_xlabel('TESS magnitude')
 		ax1.set_ylabel('Pixels in aperture')
 		xtick_major = np.median(np.diff(ax1.get_xticks()))

dataval/utilities.py

@@ -221,18 +221,20 @@ def mad(x):
 	return mad_to_sigma * nanmedian(np.abs(x - nanmedian(x)))
 
 #--------------------------------------------------------------------------------------------------
-def mag2flux(mag, zp=20.60654144):
+def mag2flux(mag, zp=20.451):
 	"""
 	Convert from magnitude to flux using scaling relation from
 	aperture photometry. This is an estimate.
 
-	The scaling is based on fast-track TESS data from sectors 1 and 2.
+	The default scaling is based on TASOC Data Release 5 from sectors 1-5.
 
 	Parameters:
-		mag (float): Magnitude in TESS band.
+		mag (ndarray): Magnitude in TESS band.
+		zp (float): Zero-point to use in scaling. Default is estimated from
+			TASOC Data Release 5 from TESS sectors 1-5.
 
 	Returns:
-		float: Corresponding flux value
+		ndarray: Corresponding flux value
 	"""
 	return np.clip(10**(-0.4*(mag - zp)), 0, None)