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

# # pyphot - A tool for computing photometry from spectra
# 
# 
# Some examples are provided in this notebook
# 
# Full documentation available at http://mfouesneau.github.io/docs/pyphot/

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')

import pylab as plt
import numpy as np

import sys
sys.path.append('../')
from pyphot import sandbox as pyphot


# ## Quick Start

# Quick start example to access the library and it's content

# In[2]:


# get the internal default library of passbands filters
lib = pyphot.get_library()
print("Library contains: ", len(lib), " filters")
# find all filter names that relates to IRAC
# and print some info
f = lib.find('irac')
for name in f:
    lib[name].info(show_zeropoints=True)


# Suppose one has a calibrated spectrum and wants to compute the vega magnitude throug the HST WFC3 F110W passband,

# In[3]:


# convert to magnitudes
import numpy as np

# We'll use Vega spectrum as example
from pyphot.vega import Vega
vega = Vega()
f = lib['HST_WFC3_F110W']
# compute the integrated flux through the filter f
# note that it work on many spectra at once
fluxes = f.get_flux(vega.wavelength, vega.flux, axis=-1)
# Note that fluxes is now with units of erg/s/cm2/AA
# pyphot gives Vega in flam and can convert between flux density units. 
fluxes, vega.wavelength, vega.flux


# In[4]:


# convert to vega magnitudes
mags = -2.5 * np.log10(fluxes.magnitude) - f.Vega_zero_mag
print("Vega magnitude of Vega in {0:s} is : {1:f} mag".format(f.name, mags))
mags = -2.5 * np.log10(fluxes.magnitude) - f.AB_zero_mag
print("AB magnitude of Vega in {0:s} is : {1:f} mag".format(f.name, mags))
mags = -2.5 * np.log10(fluxes.magnitude) - f.ST_zero_mag
print("ST magnitude of Vega in {0:s} is : {1:f} mag".format(f.name, mags))


# ## Provided Filter library
# This section shows the content of the provided library with respective properties of the passband filters. The code to generate the table is also provided in the documentation.

# In[5]:


import pyphot

# define header and table format (as csv)
hdr = ("name", "detector type", "wavelength units",
       "central wavelength", "pivot wavelength", "effective wavelength",
       "Vega mag", "Vega flux", "Vega Jy",
       "AB mag", "AB flux", "AB Jy",
       "ST mag", "ST flux", "ST Jy")
fmt = "{0:s},{1:s},{2:s},{3:.3f},{4:.3f},{5:.3f},{6:.5f},{7:.5g},{8:.5g},{9:.5f},{10:.5g},{11:.5g},{12:.5f},{13:.5g},{14:.5g}\n"

l = pyphot.get_library()

with open('table.csv', 'w') as output:
    output.write(','.join(hdr) + '\n')

    for k in sorted(l.content):
        fk = l[k]
        rec = (fk.name, fk.dtype, fk.wavelength_unit,
               fk.cl.magnitude, fk.lpivot.magnitude, fk.leff.magnitude,
               fk.Vega_zero_mag, fk.Vega_zero_flux.magnitude, fk.Vega_zero_Jy.magnitude,
               fk.AB_zero_mag, fk.AB_zero_flux.magnitude, fk.AB_zero_Jy.magnitude,
               fk.ST_zero_mag, fk.ST_zero_flux.magnitude, fk.ST_zero_Jy.magnitude)
        output.write(fmt.format(*rec)) 


# **Table description**
# 
# * name: the identification name of the filter in the library.
# * detector type: energy or photon counter.
# * wavelength units: filter defined with these units and all wavelength properties: central wavelength, pivot wavelength, and effective wavelength.
# * <X> mag: magnitude in Vega, AB or ST system (w.r.t. the detector type)
# * <X> flux: flux in $erg/s/cm^2/AA $ in the X system
# * <X> Jy: flux in $Jy$ (Jansky) in the X system

# In[6]:


import pandas as pd
df = pd.read_csv('./table.csv')
df.head()


# ## Extention to Lick indices
# 
# We also include functions to compute lick indices and provide a series of commonly use ones.
# 
# The Lick system of spectral line indices is one of the most commonly used methods of determining ages and metallicities of unresolved (integrated light) stellar populations.

# In[7]:


# convert to magnitudes
import numpy as np
from pyphot.sandbox import UnitLickLibrary as LickLibrary
from pyphot.vega import Vega

vega = Vega()
# using the internal collection of indices
lib = LickLibrary()
f = lib['CN_1']
# work on many spectra at once
index = f.get(vega.wavelength, vega.flux, axis=-1)
print("The index of Vega in {0:s} is {1:f} {2:s}".format(f.name, index, f.index_unit))


# Similarly, we show the content of the provided library with respective properties of the passband filters. 
# The table below is also part of the documentation.

# In[9]:


# define header and table format (as csv)
hdr = ("name", "wavelength units", "index units", "min", "max" "min blue", "max blue", "min red", "max red")
fmt = "{0:s},{1:s},{2:s},{3:.3f},{4:.3f},{5:.3f},{6:.5f},{7:.3f},{8:.3f}\n"

l = pyphot.LickLibrary()

with open('licks_table.csv', 'w') as output:
    output.write(','.join(hdr) + '\n')

    for k in sorted(l.content):
        fk = l[k]
        # wavelength have units
        band = fk.band.magnitude
        blue = fk.blue.magnitude
        red = fk.red.magnitude
        rec = (fk.name, fk.wavelength_unit, fk.index_unit, band[0], band[1],
               blue[0], blue[1], red[0], red[1])
        output.write(fmt.format(*rec))


# In[19]:


import pandas as pd
df = pd.read_csv('./licks_table.csv', index_col=False)
df.head()


# In[ ]: