<commit_before>#! /usr/bin/python
import numpy as np
from astropy import units as u
from astropy import constants as c

'''
Radial Velocity uncertainty based off of signal in a somewhat idealized spectrograph, and stellar type. Equations and data from Beatty and Gaudi (2015). DOI: 10.1086/684264

Also uses http://www.aanda.org/articles/aa/full/2001/29/aa1316/aa1316.right.html for reference.

Current limitations: Detector properties not closely based on actual hardware, need to double check vsini, microturbulence fixed at ~1 km/s, granulation/starspots/other jitter not considered.
'''

# Beatty spectra (simulated): 100 A chunks; Solar metallicity;
# R = 100k (4000-6500 A), 165000 (6500-10000 A), 88000 (10000A - 25000 A)
# In general need to specify wavelength range, as corrections vary with the one chosen (optical, red, near IR).
# Table from http://www.personal.psu.edu/tgb15/beattygaudi/table1.dat
beatty = np.genfromtxt("table1.dat", dtype=None)
beatty.dtype.names = ('Angstrom', 'Teff', 'Uncertaintykms')

# Currently assuming a spectrometer entirely in the 'optical' range.
# Optical is basically Johnson V-band, so m = 0 means 3460 Jy
λ_min = 4000
λ_max = 6500
Δλ = 100
BeattyWaves = np.arange(4050, 6550, Δλ) #Optical range wavelength bits

Teff = 5800 #may need to interpolate between 2 later, as Beatty's table goes in 200 K chunks.
FeH = 0.0 #solar metallicity, can load from elsewhere
logg = 4.5 #solar surface gravity, can load from elsewhere
R = 100000 #actually varies with wavelength band, but...
vsini = 2.0 #solarish placeholder, will use actual values if possible
theta_rot = 1.13*vsini # Rough approximation, but rotational effects are near-linear, no matter limb-darkening.
# need stellar diameter/distance (or apparent magnitude in wavelength range)
# for correct brightness
rstar = 1.0 * u.solRad
dstar = 10.0 * u.pc

exptime = 300 * u.s
efficiency = 0.1 # total optical system efficiency, currently includes CCD quantum efficiency
telescope = 3.5 * u.m #diameter
area = np.pi * (telescope/2)**2

# Alpha Cen B test
# How much do R and wavelength range matter?
exptime = 10 * u.minute
telescope = 3.6 * u.m
area = np.pi * (telescope/2)**2
efficiency = 0.1
λ_min = 3800
λ_max = 6800
BeattyWaves = np.arange(3850, 6850, Δλ)
R = 110000
Teff = 5200 # 5214 K
logg = 4.37
FeH = 0.23
rstar = 0.865 * u.solRad
dstar = 1.3 * u.pc
vsini = np.pi*rstar/(38.7*u.day) * u.s/u.km

BTSettl = np.genfromtxt(str(Teff), dtype=float)
# BT Settl spectra are labled by Teff
# These spectra have a wavelength (Angstroms), and a flux (1e8 erg/s/cm^2/Angstrom) column
# sum of flux(λ)*Δλ(λ)/1e8 == total flux (in power/area) emitted. 
# Multiply bt stellar_radius^2/distance^2 for recieved.
# Spectra downloaded from other sources will use different units!

#De-duping, this is not necessary if spectra are downloaded from another source.
model = [np.array([0,0])]
for x in BTSettl:
	if np.array_equal(x,model[-1]) == False:
		model.append(x)

#I_0 = np.zeros((len(BeattyWaves), 4)) #Wavelength bin, intensity at that velocity/wavelength element in terms of power and photons, overall SNR.
I_0 = np.zeros(len(BeattyWaves), dtype=[('wavelength','f4'),('intensity','f8'),('photons','f8'),('SNR','f8'),('real photons','f8')])
for λ in np.arange(0, len(BeattyWaves)): #need some C-style array traversal.
	for i in np.arange(0, len(model)-1):
		if ((model[i][0] >= BeattyWaves[λ]) and (model[i][0] < BeattyWaves[λ]+Δλ) and (model[i][0] >= λ_min) and (model[i][0] <= λ_max)):
			# alternative rescaling, sum up power over eg: V-band -> PV
			# PV/(k*10^(-0.4mv)) = flux rescaling needed
			# k is relevant conversion factor because Vega mags are terrible
			power = (model[i][1]*1e-8*u.erg/u.cm**2/u.s/u.angstrom) * ((model[i+1][0]-model[i][0]) * u.angstrom)
			power *= rstar**2/dstar**2 #rescaling emitted spectrum based on stellar surface area and distance from us
			I_0[λ][1] += power.si * u.m**2/u.watt
			photons = power * model[i][0] * u.angstrom / (c.h * c.c)
			I_0[λ][2] += photons * area * efficiency * exptime
			I_0[λ][0] = BeattyWaves[λ]

gain = 1 # electrons generated per photon
read_noise = 2.5 #RMS of ± spurious electrons
dark_current = 4 / u.hour
#convert photon intensities to electrons and therefore signal
#Resolution elements: approximately 5 pix wide by 5 pix long for first guess -- 25 in total
n_pix = 25
# But a resolution element is not 100 A long
# R = l/dl, so 1 resolution element at wavelength lambda is lambda/R wide
# Or, a bin contains 100 A/(lambda/R) = 100A*R/lambda resolution elements
for λ in np.arange(0, len(BeattyWaves)):
	pix = n_pix*Δλ*R/I_0[λ][0]
	# For now, assume equal signal per pixel. A gaussian would be better.
	I_0[λ][3] = I_0[λ][2]/pix*gain / np.sqrt(I_0[λ][2]/pix*gain + (gain*read_noise*2.2*2)**2 + (dark_current*exptime)**2)
	I_0[λ][4] = I_0[λ][3]**2 * pix / (c.c/(u.km/u.s) * Δλ/I_0[λ][0])

print(I_0)
#print(sum(I_0[:,1]), "W/m^2")
print(sum(I_0['intensity']), "W/m^2")
I_SNR = sum(I_0['photons'])
print("Total SNR", I_SNR/pix*gain / np.sqrt(I_SNR/pix*gain + (gain*read_noise*2.2*2)**2 + (dark_current*exptime)**2), "Average SNR", np.mean(I_0['SNR']))

<commit_msg>Fix units (some versions)<commit_after>#! /usr/bin/python
import numpy as np
from astropy import units as u
from astropy import constants as c

'''
Radial Velocity uncertainty based off of signal in a somewhat idealized spectrograph, and stellar type. Equations and data from Beatty and Gaudi (2015). DOI: 10.1086/684264

Also uses http://www.aanda.org/articles/aa/full/2001/29/aa1316/aa1316.right.html for reference.

Current limitations: Detector properties not closely based on actual hardware, need to double check vsini, microturbulence fixed at ~1 km/s, granulation/starspots/other jitter not considered.
'''

# Beatty spectra (simulated): 100 A chunks; Solar metallicity;
# R = 100k (4000-6500 A), 165000 (6500-10000 A), 88000 (10000A - 25000 A)
# In general need to specify wavelength range, as corrections vary with the one chosen (optical, red, near IR).
# Table from http://www.personal.psu.edu/tgb15/beattygaudi/table1.dat
beatty = np.genfromtxt("table1.dat", dtype=None)
beatty.dtype.names = ('Angstrom', 'Teff', 'Uncertaintykms')

# Currently assuming a spectrometer entirely in the 'optical' range.
# Optical is basically Johnson V-band, so m = 0 means 3460 Jy
λ_min = 4000
λ_max = 6500
Δλ = 100
BeattyWaves = np.arange(4050, 6550, Δλ) #Optical range wavelength bits

Teff = 5800 #may need to interpolate between 2 later, as Beatty's table goes in 200 K chunks.
FeH = 0.0 #solar metallicity, can load from elsewhere
logg = 4.5 #solar surface gravity, can load from elsewhere
R = 100000 #actually varies with wavelength band, but...
vsini = 2.0 #solarish placeholder, will use actual values if possible
theta_rot = 1.13*vsini # Rough approximation, but rotational effects are near-linear, no matter limb-darkening.
# need stellar diameter/distance (or apparent magnitude in wavelength range)
# for correct brightness
rstar = 1.0 * u.solRad
dstar = 10.0 * u.pc

exptime = 300 * u.s
efficiency = 0.1 # total optical system efficiency, currently includes CCD quantum efficiency
telescope = 3.5 * u.m #diameter
area = np.pi * (telescope/2)**2

# Alpha Cen B test
# How much do R and wavelength range matter?
exptime = 10 * u.minute
telescope = 3.6 * u.m
area = np.pi * (telescope/2)**2
efficiency = 0.1
λ_min = 3800
λ_max = 6800
BeattyWaves = np.arange(3850, 6850, Δλ)
R = 110000
Teff = 5200 # 5214 K
logg = 4.37
FeH = 0.23
rstar = 0.865 * u.solRad
dstar = 1.3 * u.pc
vsini = np.pi*rstar/(38.7*u.day) * u.s/u.km

BTSettl = np.genfromtxt(str(Teff), dtype=float)
# BT Settl spectra are labled by Teff
# These spectra have a wavelength (Angstroms), and a flux (1e8 erg/s/cm^2/Angstrom) column
# sum of flux(λ)*Δλ(λ)/1e8 == total flux (in power/area) emitted. 
# Multiply bt stellar_radius^2/distance^2 for recieved.
# Spectra downloaded from other sources will use different units!

#De-duping, this is not necessary if spectra are downloaded from another source.
model = [np.array([0,0])]
for x in BTSettl:
	if np.array_equal(x,model[-1]) == False:
		model.append(x)

#I_0 = np.zeros((len(BeattyWaves), 4)) #Wavelength bin, intensity at that velocity/wavelength element in terms of power and photons, overall SNR.
I_0 = np.zeros(len(BeattyWaves), dtype=[('wavelength','f4'),('intensity','f8'),('photons','f8'),('SNR','f8'),('real photons','f8')])
for λ in np.arange(0, len(BeattyWaves)): #need some C-style array traversal.
	for i in np.arange(0, len(model)-1):
		if ((model[i][0] >= BeattyWaves[λ]) and (model[i][0] < BeattyWaves[λ]+Δλ) and (model[i][0] >= λ_min) and (model[i][0] <= λ_max)):
			# alternative rescaling, sum up power over eg: V-band -> PV
			# PV/(k*10^(-0.4mv)) = flux rescaling needed
			# k is relevant conversion factor because Vega mags are terrible
			power = (model[i][1]*1e-8*u.erg/u.cm**2/u.s/u.angstrom) * ((model[i+1][0]-model[i][0]) * u.angstrom)
			power *= rstar**2/dstar**2 #rescaling emitted spectrum based on stellar surface area and distance from us
			I_0[λ][1] += power.si * u.m**2/u.watt
			photons = power * model[i][0] * u.angstrom / (c.h * c.c)
			I_0[λ][2] += photons * area * efficiency * exptime
			I_0[λ][0] = BeattyWaves[λ]

gain = 1 # electrons generated per photon
read_noise = 2.5 #RMS of ± spurious electrons
dark_current = 4 / u.hour
#convert photon intensities to electrons and therefore signal
#Resolution elements: approximately 5 pix wide by 5 pix long for first guess -- 25 in total
n_pix = 25
# But a resolution element is not 100 A long
# R = l/dl, so 1 resolution element at wavelength lambda is lambda/R wide
# Or, a bin contains 100 A/(lambda/R) = 100A*R/lambda resolution elements
for λ in np.arange(0, len(BeattyWaves)):
	pix = n_pix*Δλ*R/I_0[λ][0]
	# For now, assume equal signal per pixel. A gaussian would be better.
	I_0[λ][3] = I_0[λ][2]/pix*gain / np.sqrt(I_0[λ][2]/pix*gain + (gain*read_noise*2.2*2)**2 + (dark_current*exptime)**2)
	I_0[λ][4] = I_0[λ][3]**2 * pix / ((c.c/(u.km/u.s)).si * Δλ/I_0[λ][0])

print(I_0)
#print(sum(I_0[:,1]), "W/m^2")
print(sum(I_0['intensity']), "W/m^2")
I_SNR = sum(I_0['photons'])
print("Total SNR", I_SNR/pix*gain / np.sqrt(I_SNR/pix*gain + (gain*read_noise*2.2*2)**2 + (dark_current*exptime)**2), "Average SNR", np.mean(I_0['SNR']))


{'fraction_non_alphanumeric': 0.1414111922141119, 'fraction_numerical': 0.0570316301703163, 'mean_word_length': 4.4572490706319705, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '992966', 'n_tokens_mistral': 4229, 'n_tokens_neox': 3712, 'n_words': 1193}