In [95]:
import matplotlib.pyplot as plt
In [101]:
import scipy
import numpy as np
def subtractContinuum(s, n=3, plot=False):
'''Take a 1D array, use spline to subtract off continuum.
subtractContinuum(s, n=3)
required:
s = the array
optional:
n = 3, the number of spline points to use
'''
x = np.arange(len(s))
points = (np.arange(n)+1)*len(s)/(n+1)
spline = scipy.interpolate.LSQUnivariateSpline(x,s,points)
if plot:
plt.ion()
plt.figure()
plt.plot(x, s)
plt.plot(x, spline(x), linewidth=5, alpha=0.5)
return s - spline(x)
def ccf(f, g, scale=1.0):
'''Calculate the normalized cross-correlation function of two identically-size arrays.
[required]:
f = an N-element array (for example, spectrum of target star)
g = an N-element array (for example, spectrum of template star)
scale = a scalar indicating what the indices of f and g (one unit of "lag") correspond to
'''
# how long are our arrays
N = len(f)
# define the x-axis, if not supplied
assert(len(f) == len(g))
x = np.arange(-N+1, N, 1.0)*scale
# calculation the normalized cross-correlation function
sigma_f = np.sqrt(np.sum(f**2)/N)
sigma_g = np.sqrt(np.sum(g**2)/N)
C_fg = np.correlate(f, g, 'full')/N/sigma_f/sigma_g
# WILL THIS WORK?
return scipy.interpolate.interp1d(x,C_fg, fill_value=0.0, bounds_error=False)
def todcor(f, g1, g2, scale=1.0, luminosity_ratio=None):
'''Calculate the 2D correlation of a 1D array with two template arrays.
'''
assert(len(f) == len(g1))
assert(len(f) == len(g2))
C_1 = ccf(f, g1, scale=scale)
C_2 = ccf(f, g2, scale=scale)
C_12 = ccf(g1, g2, scale=scale)
N = len(f)
sigma_g1 = np.sqrt(np.sum(g1**2)/N)
sigma_g2 = np.sqrt(np.sum(g2**2)/N)
def bestalphaprime(s1, s2):
return sigma_g1/sigma_g2*(C_1(s1)*C_12(s2 - s1) - C_2(s2))/(C_2(s2)*C_12(s2-s1) - C_1(s1))
def R(s1, s2):
#a =
if luminosity_ratio is None:
a = np.maximum(np.minimum(bestalphaprime(s1,s2), sigma_g2/sigma_g1), 0.0)
flexiblecorrelation = (C_1(s1) + a*C_2(s2))/np.sqrt(1.0 + 2*a*C_12(s2 - s1) + a**2)
ok = np.isfinite(a)
peak = np.argmax(flexiblecorrelation[ok].flatten())
a = a[ok].flatten()[peak]
else:
a = luminosity_ratio*sigma_g2/sigma_g1
# print("alpha spans", np.min(a), np.max(a))
return (C_1(s1) + a*C_2(s2))/np.sqrt(1.0 + 2*a*C_12(s2 - s1) + a**2), a*sigma_g1/sigma_g2
return R
In [108]:
from astropy.io import fits
spec1 = fits.getdata("../tests/testdata/lte05200-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits")
spec2 = fits.getdata("../tests/testdata/lte03400-4.50-0.0.PHOENIX-ACES-AGSS-COND-2011-HiRes.fits")
wav = fits.getdata("../tests/testdata/WAVE_PHOENIX-ACES-AGSS-COND-2011.fits")
wav = wav/10
from PyAstronomy import pyasl
spec1_shifted, wlprime = pyasl.dopplerShift(wav, spec1, -5)
spec2_shifted, wlprime = pyasl.dopplerShift(wav, spec2, 15)
spec_3 = spec1_shifted + spec2_shifted
spec_1 = spec1[(wav > 2100) * (wav< 2200)]
spec_2 = spec2[(wav > 2100) * (wav< 2200)]
spec_3 = spec_3[(wav > 2100) * (wav< 2200)]
wav_3 = wav[(wav > 2100) * (wav< 2200)]
# need to log resample so that evenly spaced in rv
def log_resample(x, drv=0.01):
logx = np.log(x)
logx = np.arange(logx[0], logx[-1], drv)
return np.exp(logx)
new_full_wav = log_resample(wav, 0.001)
new_wav = log_resample(wav_3)
new_s1 = np.interp(new_full_wav, wav, spec1)
new_s2 = np.interp(new_full_wav, wav, spec2)
new_s3 = np.interp(new_wav, wav_3, spec_3)
In [107]:
plt.plot(wav, spec1, label="1")
plt.plot(new_full_wav, new_s1, label="1")
plt.show()
plt.plot(wav_3, spec2, label="2")
plt.show()
plt.plot(wav_3, spec_3, label="3")
plt.legend()
plt.show()
In [ ]:
s1 = subtractContinuum(spec1, n=5)
s2 = subtractContinuum(spec2,n=5)
s3 = subtractContinuum(spec_3,n=5)
plt.plot(s1)
plt.plot(s2)
plt.plot(s3)
plt.show()
In [93]:
R = todcor(s3, s1, s2)
In [97]:
rv = np.arange(-50, 50, 0.5)
gamma = np.arange(-50, 50, 1)
RV, GAM = np.meshgrid(rv, gamma, indexing="ij")
c_n1_n2 = np.empty((len(rv), len(gamma)))
a_n1_n2 = np.empty((len(rv), len(gamma)))
for jj, y in enumerate(gamma):
for ii, x in enumerate(rv):
res = R(x, y)
c_n1_n2[ii, jj] = res[0]
a_n1_n2[ii, jj] = res[1]
In [ ]:
In [109]:
plt.contourf(RV, GAM, c_n1_n2)
plt.title("correlation value")
plt.show()
a_n1_n2
plt.contourf(RV, GAM, a_n1_n2)
plt.title("flux ratio")
plt.show()
In [ ]:
In [ ]:
In [ ]: