hst_to_wfi.py

#! /usr/bin/env python
'''
Reads in HST F606/F814/F125/F160 AB mags and produces STIPS
simulation input files for WFI filters
'''
from __future__ import (absolute_import, division, print_function, unicode_literals)

import numpy as np
from astropy.io import ascii
from astropy.table import Table
from astropy import constants as const

c = const.c.cgs.value
h = const.h.cgs.value
fo = 3631 # Jy

hst_filt = ['F606W','F814W','F125W','F160W']                                   # HST, ACS/ACS/WFC3/WFC3 filters
wfi_filt = ['X606','X625','X670','Z087','Y106','J129','H158','F184']           # WFI, 3 blue filter options

wv1 = np.array([0.599548,0.811524,1.236464,1.537033])*1e-4                     # lambda_pivot (cm), HST
wv2 = np.array([0.5965,0.6254,0.6725,0.8709,1.0611,1.2945,1.5809,1.8452])*1e-4 # same ... WFI
awp = np.array([0.7458,0.9712,0.8779,0.5419,0.5353,0.5054,0.5222,0.3658])*1e4  # Aeff*Weff/lambda_pivot (cm^2)

def DoAll(wv1,wv2):
    glxs, mags = [],[]
    for i,filt in enumerate(hst_filt):
        glxs.append(ascii.read(filt+'.dat'))
        mags.append(glxs[i]['mag'])
    m1  = np.array(mags).T

    in1 = np.arange(m1.shape[0])\
          [~(np.isnan(m1[:,0])|np.isnan(m1[:,1])|np.isnan(m1[:,2])|np.isnan(m1[:,3]))]

    # temporary cut at H=20
    in1 = in1[m1[in1,3]>20]

    m1 = m1[in1,:]

    f_nu1 = 10**(m1/(-2.5)) * fo * 1e-23                                       # ergs/s/cm^2/Hz
    N_lm1 = np.log10(f_nu1/(wv1*h))                                            # Photons/s/cm^2/cm
    N_lm2 = np.zeros((N_lm1.shape[0],wv2.shape[0]))
    wv1, wv2 = np.log10(wv1),np.log10(wv2)

    tmp = np.array([np.polyfit(wv1[0:2], [N_lm1[i,0], N_lm1[i,1]], 1) \
                    for i in range(N_lm1.shape[0])])
    A, B = tmp[:,0], tmp[:,1]
    N_lm2[:,0], N_lm2[:,1], N_lm2[:,2] = A*wv2[0]+B, A*wv2[1]+B, A*wv2[2]+B    # F606/814 >>> X606/625/670
    del A, B, tmp

    tmp = np.array([np.polyfit(wv1[1:3], [N_lm1[i,1], N_lm1[i,2]], 1) \
                    for i in range(N_lm1.shape[0])])
    A, B = tmp[:,0], tmp[:,1]
    N_lm2[:,3], N_lm2[:,4] = A*wv2[3]+B, A*wv2[4]+B                            # F814/110 >>> Z087/Y106
    del A, B, tmp

    tmp = np.array([np.polyfit(wv1[2:4], [N_lm1[i,2], N_lm1[i,3]], 1) \
                    for i in range(N_lm1.shape[0])])
    A, B = tmp[:,0], tmp[:,1]
    N_lm2[:,5], N_lm2[:,6], N_lm2[:,7] = A*wv2[5]+B, A*wv2[6]+B, A*wv2[7]+B    # F110/160 >>> J129/H158/F184
    del A, B, tmp

    N = 10**N_lm2 * awp * 10**wv2                                              # Counts/s; WFI
    t = [0,0,0,1,2,2,3,3]

    return [write_stips_table(glxs[t[i]],N[:,i],wfi_filt[i],in1) for i in range(len(t))]

def write_stips_table(gals,flux,filt,in1):
    tab = [gals['id'][in1], gals['ra'][in1], gals['dec'][in1], flux,\
           gals['type'][in1], gals['n'][in1], gals['re'][in1],\
           gals['phi'][in1], gals['ratio'][in1], gals['notes'][in1]]
    nms  = ('id', 'ra', 'dec', 'flux', 'type', 'n', 're', 'phi', 'ratio', 'notes')
    fmt  = {'id':'%10d', 'ra':'%15.7f', 'dec':'%15.7f',  'flux':'%15.7f',  'type':'%8s',\
            'n':'%10.3f', 're':'%15.7f', 'phi':'%15.7f', 'ratio':'%15.7f', 'notes':'%8s'}
    t    = Table(tab, names=nms)
    outfile = '.'.join([filt,'txt'])
    radec  = np.mean(np.array([gals['ra'][in1],gals['dec'][in1]]).T, axis=0)
    ascii.write(t, outfile, format='fixed_width', delimiter='', formats=fmt)
    return print('\nWrote out:', outfile) #, ';\tCenter = ', radec)

if __name__ == '__main__':
    import sys, time
    tic = time.time()
    print('\n\nPython check: This should be 1.500 = %.3f \n' % (3/2))
    DoAll(wv1,wv2)
    print('\n\nCompleted in %.3f seconds \n' % (time.time()-tic))