frms_comparison.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from astropy.io import fits
from astropy.stats import sigma_clip
import matplotlib.pyplot as plt

# Nicer plotting if possible, requires matplotlib>=1.5.0
try:
    plt.style.use(['ggplot', 'seaborn-talk', 'seaborn-ticks'])
except:
    pass


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)
    axes[0].legend(loc='best')

    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')
    fig.tight_layout()
    plt.show()