New KeplerQ9v3-long training set

docs/starclass.rst

@@ -13,6 +13,7 @@ Training Sets
 
 	starclass.training_sets.TrainingSet
 	starclass.training_sets.keplerq9v3
+	starclass.training_sets.keplerq9v3_long
 	starclass.training_sets.keplerq9v2
 	starclass.training_sets.keplerq9
 	starclass.training_sets.tdasim

docs/starclass.training_sets.keplerq9v3_long.rst

@@ -0,0 +1,14 @@
+Kepler Q9 (version 3), long, Training Set (``starclass.training_sets.keplerq9v3_long``)
+=======================================================================================
+
+.. autoclass:: starclass.training_sets.keplerq9v3_long
+	:show-inheritance:
+	:members:
+	:undoc-members:
+	:inherited-members:
+
+.. autoclass:: starclass.training_sets.keplerq9v3_long_instr
+	:show-inheritance:
+	:members:
+	:undoc-members:
+	:inherited-members:

notes/constant_stars/extract_constant.py

@@ -0,0 +1,37 @@
+
+import numpy as np
+import os.path
+import glob
+from tqdm import tqdm
+import matplotlib.pyplot as plt
+
+files = glob.glob(r'G:\keplerq9\keplerq9v3-long\CONSTANT\constant_*.txt')
+files = sorted(files, key=lambda x: int(os.path.basename(x).replace('constant_', '').replace('.txt', '')))
+
+with open('keplerq9v3_targets_constant.txt', 'w') as fid:
+
+	fid.write("# Kepler Q9 Training Set Targets (version 3)\n")
+	fid.write("# Column 1: Constant star identifier\n")
+	fid.write("# Column 2: Standard deviation (sigma)\n")
+	fid.write("#-------------------------------------------\n")
+
+	for fname in tqdm(files):
+
+		data = np.loadtxt(fname, dtype='float64')
+
+		#print(data)
+		#sigma = np.std(data[:,1])
+		flux_err = np.median(data[:,2])
+
+		#print(sigma, flux_err)
+
+		#plt.figure()
+		#plt.scatter(data[:,0], data[:,1])
+		#plt.show()
+
+		fid.write("{0:s},{1:.16e}\n".format(
+			os.path.basename(fname).replace('.txt', ''),
+			flux_err
+		))
+
+	fid.write("#-------------------------------------------\n")

notes/constant_stars/noisedistribution.py

@@ -0,0 +1,101 @@
+import astropy.constants as const
+import astropy.units as u
+from astroquery.mast import Catalogs
+import noiseestimation as ne
+import numpy as np
+import glob, os
+
+#TESS white noise lightcurve simulator
+def TIC_byID(ID):
+    """
+    """
+    catTable = Catalogs.query_criteria(ID=ID, catalog="Tic")
+
+    return catTable['ID','ra','dec','Tmag','Teff']
+
+
+#generate random integer
+num = 0
+cadence = 30*60. #30 min cadence
+whitenoise = []
+teffdist = []
+tmagdist = []
+radist = []
+decdist = []
+#
+while num < 1000:
+    id = np.random.randint(0,472000000.)
+    catdata = TIC_byID(id)
+    if len(catdata)>1:
+        print(catdata)
+    if len(catdata==1):
+        if catdata['Tmag'] < 15 and catdata['Teff'] < 9700. and catdata['Teff'] > 2450:
+            PARAM = {}
+            PARAM['RA'] = catdata['ra']
+            PARAM['DEC'] = catdata['dec']
+            noise, PARAM = ne.phot_noise(catdata['Tmag'], catdata['Teff'], cadence, PARAM)
+            totnoise = np.sqrt(noise[0]**2+noise[1]**2+noise[2]**2+noise[3]**2)
+            whitenoise.append(totnoise[0])
+            tmagdist.append(catdata['Tmag'][0])
+            teffdist.append(catdata['Teff'][0])
+            radist.append(catdata['ra'][0])
+            decdist.append(catdata['dec'][0])
+            num += 1      
+            print(num)  
+    else:
+        print('ID not in TIC')
+
+whitenoise = np.array(whitenoise)
+tmagdist = np.array(tmagdist)
+teffdist = np.array(teffdist)
+radist = np.array(radist)
+decdist = np.array(decdist)
+
+
+np.savetxt('whitenoisedist.txt',whitenoise)
+np.savetxt('tmagdist.txt',tmagdist)
+np.savetxt('teffdist.txt',teffdist)
+np.savetxt('radist.txt',radist)
+np.savetxt('decdist.txt',decdist)
+
+kepdata = '/Users/davidarmstrong/Software/Python/TESS/batch06_noisy/'
+keplist = glob.glob(os.path.join(kepdata,'*.noisy'))
+import pylab as p
+p.ion()
+
+for i,lcjit in enumerate(whitenoise):
+    print(i)
+    #load random kepler lc.
+    kepfile = np.genfromtxt(np.random.choice(keplist))
+    timestamps = kepfile[:,0]
+    lcjit_relflux = lcjit*1e-6
+    flux = np.random.normal(1,lcjit_relflux,len(timestamps))
+    flux_err = np.ones(len(timestamps)) * lcjit_relflux
+    output = np.array([timestamps,flux,flux_err]).T
+    #p.figure(1)
+    #p.clf()
+    #p.errorbar(timestamps,flux,yerr=flux_err,fmt='b.')
+    #p.pause(1)
+    #input('a')
+    np.savetxt('constantlcs/constant_'+str(i)+'.txt',output)
+
+#get Tmag, Teff, RA, DEC
+
+
+
+#add noise sources
+
+#define time baseline
+
+#make curve
+
+
+
+
+
+
+#nsamples = 1000
+#Tmagdist = np.random.uniform(5.,15.,nsamples)
+#Teffdist = np.random.uniform(2450.,9700.,nsamples)
+#RAdist = np.random.uniform(0.,360.,nsamples)
+#DECdist = np.random.uniform(-89.,-20.,nsamples)

notes/constant_stars/noiseestimation.py

@@ -0,0 +1,86 @@
+import scipy.interpolate as INT
+from astropy.coordinates import SkyCoord
+import astropy.units as u
+import numpy as np
+
+def ZLnoise(gal_lat):
+    # RMS noise from Zodiacal background
+    rms = (16-10)*(gal_lat/90 -1)**2 + 10 # e-1 / pix in 2sec integration
+    return rms
+
+def Pixinaperture(Tmag):
+    # Approximate relation for pixels in aperture (based on plot in Sullivan et al.)
+    pixels = (30 + (((3-30)/(14-7)) * (Tmag-7)))*(Tmag<14) + 3*(Tmag>=14) 
+    return int(np.max([pixels, 3]))
+
+def mean_flux_level(Tmag, Teff):
+    # 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
+
+    return Flux
+
+
+def phot_noise(Tmag, Teff, cad, PARAM, verbose=False, sysnoise=60):
+
+	# Calculate galactic latitude for Zodiacal noise
+	gc= SkyCoord(PARAM['RA']*u.degree, PARAM['DEC']*u.degree, frame='icrs')
+#    gc = SkyCoord(lon=PARAM['ELON']*u.degree, lat=PARAM['ELAT']*u.degree, frame='barycentrictrueecliptic')
+	gc_gal = gc.transform_to('galactic')
+	gal_lat0 = gc_gal.b.deg
+
+	gal_lat = np.arcsin(np.abs(np.sin(gal_lat0*np.pi/180)))*180/np.pi
+
+	# Number of 2 sec integrations in cadence
+	integrations = cad/2
+
+	# Number of pixels in aperture given Tmag
+	pixels = int(Pixinaperture(Tmag))
+
+	# noise values are in rms, so square-root should be used when factoring up
+	Flux_factor = np.sqrt(integrations * pixels)
+
+	# Mean flux level in electrons per cadence
+	mean_level_ppm = mean_flux_level(Tmag, Teff) * cad # electrons
+
+	# Shot noise
+	shot_noise = 1e6/np.sqrt(mean_level_ppm)
+
+	# Read noise
+	read_noise = 10 * Flux_factor *1e6/mean_level_ppm # ppm
+
+	# Zodiacal noise
+	zodiacal_noise = ZLnoise(gal_lat) * Flux_factor *1e6/mean_level_ppm # ppm
+
+	# Systematic noise in ppm
+	systematic_noise_ppm = sysnoise / np.sqrt(cad/(60*60)) # ppm / sqrt(hr)
+
+
+	if verbose:
+		print('Galactic latitude', gal_lat)
+		print('Systematic noise in ppm', systematic_noise_ppm)
+		print('Integrations', integrations)
+		print('Pixels', pixels)
+		print('Flux factor', Flux_factor)
+		print('Mean level ppm', mean_level_ppm)
+		print('Shot noise', shot_noise)
+		print('Read noise', read_noise)
+		print('Zodiacal noise', zodiacal_noise)
+
+
+	PARAM['Galactic_lat'] = gal_lat
+	PARAM['Pixels_in_aper'] = pixels
+
+	noise_vals = np.array([shot_noise, zodiacal_noise, read_noise, systematic_noise_ppm])
+	return noise_vals, PARAM # ppm per cadence
+
+