<jupyter_start><jupyter_text>Note:This notebook is aimed at getting a handle on how the FFT spectral resolution convolution works, and how we can modify it to have wavelength-dependent spectral resolution.  See Issue 8 in the `jammer` repository.<jupyter_code>import Starfish
import os
from Starfish.grid_tools import HDF5Creator
instrument = eval("Starfish.grid_tools." + Starfish.data["instruments"][0])()
instrument.FWHM
instrument.name
try:
    air = instrument.air
    print("New in v0.3: Using explicit air/vacuum state from Instrument class.")
except AttributeError:
    air = True
    if (Starfish.data["grid_name"] == "PHOENIX") & (len(Starfish.grid['parname']) == 3):
        mygrid = eval("Starfish.grid_tools." + Starfish.data["grid_name"]+ "GridInterfaceNoAlpha")(air=air)
    else:
        mygrid = eval("Starfish.grid_tools." + Starfish.data["grid_name"]+ "GridInterface")(air=air)
hdf5_path = os.path.expandvars(Starfish.grid["hdf5_path"])
Starfish.grid["parrange"]
creator = HDF5Creator(mygrid, hdf5_path, instrument,
                      ranges=Starfish.grid["parrange"])
plt.plot(creator.wl_native[:-1], np.diff(creator.wl_native), label='Native Sampling')
plt.plot(creator.wl_FFT[:-1], np.diff(creator.wl_FFT), label='FFT Sampling')
plt.plot(creator.wl_final[:-1], np.diff(creator.wl_final), label='Final Sampling')
plt.ylabel('$\delta \lambda (\AA)$')
plt.xlabel('$\lambda \; (\AA)$')
plt.ylim(0,30)
plt.legend(loc='best');<jupyter_output><empty_output><jupyter_text>Starfish's `grid_tools` slightly oversamples the spectrum with linear interpolation between the pixels.<jupyter_code>creator.dv_native, creator.dv_FFT, creator.dv_final
creator.wl_final
creator.points
import itertools
param_list = []

# use itertools.product to create permutations of all possible values
for i in itertools.product(*creator.points):
    param_list.append(np.array(i))

all_params = np.array(param_list)
test_par = [800.0, 5.0]
creator.Instrument.FWHM
#creator.Instrument.FWHM = 9200.0/2.0
sigma = creator.Instrument.FWHM / 2.35 # in km/s
print(creator.Instrument)
#creator.taper = np.exp(-2 * (np.pi ** 2) * (sigma ** 2) * (creator.ss ** 2))
fl, header = creator.process_flux(test_par)
raw, hdr = mygrid.load_flux(test_par)
sns.set_context('poster', font_scale=1.7)
wl = creator.wl_final<jupyter_output><empty_output><jupyter_text>Read the data and previous models<jupyter_code>import pandas as pd
dat1 = pd.read_csv('../../sf/Gl570D/output/marley_grid/run01/spec_config.csv')
#dat2 = pd.read_csv('../sf/Gl570D/output/marley_grid/run01/models_draw.csv')<jupyter_output><empty_output><jupyter_text>Read .flx files from Mark on April 20, 2017<jupyter_code>dat_MM_flx = pd.read_csv('./t800g1000nc.flx', delim_whitespace=True, 
                         names = ['id', 'wl_um', 'junk', 'flux_density_per_hz'])
dat_MM_flx2 = pd.read_csv('./t800g1000nc_16.flx', delim_whitespace=True, 
                          names = ['id', 'wl_um', 'junk', 'flux_density_per_hz'])
import astropy.units as u
x = dat_MM_flx.wl_um.values*u.micron
f = dat_MM_flx.flux_density_per_hz
f = (f.values*u.erg/u.cm**2/u.s/u.Hz).to(
             u.erg/u.cm**2/u.s/u.Angstrom, 
             equivalencies=u.spectral_density(x))
x2 = dat_MM_flx2.wl_um.values*u.micron
f2 = dat_MM_flx2.flux_density_per_hz
f2 = (f2.values*u.erg/u.cm**2/u.s/u.Hz).to(
             u.erg/u.cm**2/u.s/u.Angstrom, 
             equivalencies=u.spectral_density(x))
plt.plot(creator.wl_native, raw/3, 'k-', ms=0.5, alpha=1, label='Native resolution')
plt.step(wl, fl/3, 'r', label='New model ~SpeX res')
plt.step(dat1.wl, dat1.data*10**18, 'b-', ms=0.5, alpha=1, label='Gl570D Line et al.')
plt.step(x*10000.0, f/9, 'k', label = '2010 models')
plt.step(x2*10000.0, f2/9, 'g', label = '2016 models')
plt.xlim(10000, 25000)
plt.ylim(-100, 1000)
#plt.yscale('log')
plt.xlabel('$\lambda\; (\AA)$')
plt.ylabel('$f_\lambda$')
plt.legend(loc='best')
#plt.savefig('../../results/ex_model_spectrum.pdf', bbox_inches='tight');<jupyter_output><empty_output><jupyter_text>Key FFT programming idea-- The `ss[0]` term is the DC component-- it controls the average flux level.  The hack in `grid_tools.py` setting `ss[0]` to a non-zero value corrupted the average flux level when the taper was re-computed after updating `ss[0]`.  **I have simply turned off the ss[0] = 0.01 hack.**<jupyter_code>len(wl), len(fl)
plt.plot(creator.taper[0:2000])
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>The FFT step.<jupyter_code>from numpy.fft import fft, ifft, fftfreq, rfftfreq
from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
# Interpolate the native spectrum to a log-lam FFT grid
interp = InterpolatedUnivariateSpline(creator.wl_native, raw, k=5)
fl = interp(creator.wl_FFT)
del interp
# Do the FFT
FF = np.fft.rfft(fl)
FF_tap = FF * creator.taper
FF.shape
n_pix, = creator.wl_final.shape
fl_final = creator.wl_final*0.0
n_pix//8<jupyter_output><empty_output><jupyter_text>FWHM = 299792.46/R<jupyter_code>for j in range(n_pix//8):
    #R_est = np.mean(creator.Instrument.res_gradient(creator.wl_final[j:j+8]))
    #sigma = 299792.46/R_est / 2.35 # in km/s
    # Instrumentally broaden the spectrum by multiplying with a Gaussian in Fourier space
    #creator.taper = np.exp(-2 * (np.pi ** 2) * (sigma ** 2) * (self.ss ** 2))
    FF_tap = FF * creator.taper
    fl_tapered = np.fft.irfft(FF_tap)
    interp = InterpolatedUnivariateSpline(creator.wl_FFT, fl_tapered, k=5)
    fl_final[8*j:8*j+8] = interp(creator.wl_final[8*j:8*j+8])
fl_final_orig = fl_final.copy()
import gc
for j in range(n_pix//8):
    R_est = np.mean(creator.Instrument.res_gradient()(creator.wl_final[8*j:8*j+8]))
    sigma = 299792.46/R_est / 2.35 # in km/s
    # Instrumentally broaden the spectrum by multiplying with a Gaussian in Fourier space
    creator.ss[0] = 0.00
    taper = np.exp(-2 * (np.pi ** 2) * (sigma ** 2) * (creator.ss ** 2))
    FF_tap = FF * taper
    fl_tapered = np.fft.irfft(FF_tap)
    interp = InterpolatedUnivariateSpline(creator.wl_FFT, fl_tapered, k=5)
    fl_final[8*j:8*j+8] = interp(creator.wl_final[8*j:8*j+8])
    del interp
    gc.collect()
fl_final.shape
fl_fin_multi = np.empty((2048, 256))
R_trial = np.arange(50, 50+256)
for j in range(256):
    R_est = R_trial[j]
    sigma = 299792.46/R_est / 2.35 # in km/s
    # Instrumentally broaden the spectrum by multiplying with a Gaussian in Fourier space
    taper = np.exp(-2 * (np.pi ** 2) * (sigma ** 2) * (creator.ss ** 2))
    creator.ss[0] = 0.00
    FF_tap = FF * taper
    fl_tapered = np.fft.irfft(FF_tap)
    interp = InterpolatedUnivariateSpline(creator.wl_FFT, fl_tapered, k=5)
    fl_fin_multi[:, j] = interp(creator.wl_final)
    del interp
    gc.collect()
sns.set_context('notebook')
plt.plot(creator.wl_final, fl_final_orig)
plt.plot(creator.wl_final, fl_final)
#for j in range(256):
#    plt.plot(creator.wl_final, fl_fin_multi[:, j], alpha=0.1)

plt.plot(creator.wl_final, (fl_final_orig-fl_final))<jupyter_output><empty_output><jupyter_text>WTF is with that??<jupyter_code>plt.plot(creator.wl_final, fl_final_orig)
plt.plot(creator.wl_final, fl_final-np.min(fl_final-4.5))
plt.yscale('log')<jupyter_output><empty_output><jupyter_text>OK, so this looks good-- matching resolution at long wavelengths (expected behavior), with lower resolution at shorter wavelengths (expected).<jupyter_code>plt.plot(creator.ss, np.abs(FF**2), '.')
plt.plot(creator.ss, np.abs(FF_tap**2), '.')
plt.plot(creator.ss, creator.taper, '.')
plt.xscale('log')
plt.yscale('log')
plt.ylim(1.0e-20, 1.0e20);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15606783253842077, 'fraction_numerical': 0.036301006889242184, 'mean_word_length': 4.740684410646388, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '493299', 'license': 'MIT', 'n_tokens_mistral': 3219, 'n_tokens_neox': 2967, 'n_words': 561, 'path': 'sf/grid_experiments/fft_visualization.ipynb', 'repo_name': 'gully/jammer-Gl570D', 'size': '7548'}