Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
from import fits
from astropy.stats import sigma_clip
import matplotlib.pyplot as plt
# Nicer plotting if possible, requires matplotlib>=1.5.0
try:['ggplot', 'seaborn-talk', 'seaborn-ticks'])
def compute_frms_srw(flux):
'''Fractional rms computation vs Simon W's method'''
sc = sigma_clip(flux, axis=1, iters=1, sig=3.)
med_flux = np.median(sc, axis=1)
mad_flux = np.median(np.abs(sc - med_flux[:, np.newaxis]), axis=1)
std_flux = 1.48 * mad_flux
frms = std_flux / med_flux
return med_flux, frms
def compute_frms_eg(flux):
'''Fractional rms computation vs Ed G's method'''
med_flux = np.median(flux, axis=1)
relative_flux = flux / med_flux[:, np.newaxis]
relative_mag = -2.5 * np.log10(relative_flux)
N = relative_mag.shape[1]
frms = np.sqrt(np.sum(relative_mag**2, axis=1) / N)
return med_flux, frms
if __name__ == '__main__':
# Example file to use
fname = '/ngts/pipedev/ParanalOutput/nightly_data/20150911-ng2000-802-custom-flat-high-quality.fits'
flux = fits.getdata(fname, 'tamflux')
med_srw, frms_srw = compute_frms_srw(flux)
_, frms_eg = compute_frms_eg(flux)
# Plot
fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].loglog(med_srw, frms_srw, '.', label='SRW', ms=5)
axes[0].loglog(med_srw, frms_eg, '.', label='EG', ms=5)
ratio = frms_eg / frms_srw
axes[1].semilogx(med_srw, ratio, '.', ms=5)
axes[0].set(xlim=(1E1, 1E7),
xlabel='Median flux / ADU',
ylim=(1E-3, 1E0),
ylabel='Fractional rms')
axes[1].set(xlim=(1E1, 1E7),
xlabel='Median flux / ADU',
ylim=(0, 2),
ylabel='EG / SRW')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment