file_listcal = "alma_sourcecat_searchresults_20180419.csv"

q = databaseQuery()

"""
Explanation: New function to make a list and to select calibrator
I add a function to retrieve all the flux from the ALMA Calibrator list with its frequency and observing date, and to retrieve redshift (z) from NED.
End of explanation
"""


listcal = q.read_calibratorlist(file_listcal, fluxrange=[0.1, 999999])

len(listcal)

print("Name: ", listcal[0][0])
print("J2000 RA, dec: ", listcal[0][1], listcal[0][2])
print("Alias: ", listcal[0][3])
print("Flux density: ", listcal[0][4])
print("Band: ", listcal[0][5])
print("Freq: ", listcal[0][6])
print("Obs date: ", listcal[0][4])

"""
Explanation: Example, retrieve all the calibrator with a flux > 0.1 Jy:
End of explanation
"""


report, resume = q.select_object_from_sqldb("calibrators_brighterthan_0.1Jy_20180419.db", \
                                    maxFreqRes=999999999, array='12m', \
                                    excludeCycle0=True, \
                                    selectPol=False, \
                                    minTimeBand={3:0., 6:0., 7:0.}, \
                                    silent=True)

"""
Explanation: Select all calibrators that heve been observed at least in 3 Bands [ >60s in B3, B6, B7]

already queried and convert it to SQL
exclude Cycle 0, array 12m
End of explanation
"""


print("Name: ", resume[0][0])
print("From NED: ")
print("Name: ", resume[0][3])
print("J2000 RA, dec: ", resume[0][4], resume[0][5])
print("z: ", resume[0][6])
print("Total # of projects: ", resume[0][7])
print("Total # of UIDs: ", resume[0][8])
print("Gal lon: ", resume[0][9])
print("Gal lat: ", resume[0][10])

"""
Explanation: We can write a "report file" or only use the "resume data", some will have redshift data retrieved from NED.
End of explanation
"""


for i, obj in enumerate(resume):
    for j, cal in enumerate(listcal):
        if obj[0] == cal[0]: # same name
            obj.append(cal[4:]) # add [flux, band, flux obsdate] in the "resume"

"""
Explanation: Sometimes there is no redshift information found in NED
Combining listcal and resume information.
End of explanation
"""


def collect_z_and_flux(Band):
    z = []
    flux = []
    for idata in resume:
        if idata[6] is not None: # select object which has redshift information
            fluxnya = idata[11][0]
            bandnya = idata[11][1]
            freqnya = idata[11][2]
            datenya  = idata[11][3]
            for i, band in enumerate(bandnya):
                if band == str(Band): # take only first data
                    flux.append(fluxnya[i])
                    z.append(idata[6])
                    break
                    
    return z, flux

z3, f3 = collect_z_and_flux(Band=3)
print("Number of seleted source in B3: ", len(z3))
z6, f6 = collect_z_and_flux(6)
print("Number of seleted source in B6: ", len(z6))
z7, f7 = collect_z_and_flux(7)
print("Number of seleted source in B7: ", len(z7))

"""
Explanation: Select objects which has redshift

collect the flux, band, freq, and obsdate
plot based on the Band
End of explanation
"""


plt.figure(figsize=(15,10))

plt.subplot(221)
plt.plot(z3, f3, 'ro')
plt.xlabel("z")
plt.ylabel("Flux density (Jy)")
plt.title("B3")

plt.subplot(222)
plt.plot(z6, f6, 'go')
plt.xlabel("z")
plt.ylabel("Flux density (Jy)")
plt.title("B6")

plt.subplot(223)
plt.plot(z7, f7, 'bo')
plt.xlabel("z")
plt.ylabel("Flux density (Jy)")
plt.title("B7")

plt.subplot(224)
plt.plot(z3, f3, 'ro', z6, f6, 'go', z7, f7, 'bo', alpha=0.3)
plt.xlabel("z")
plt.ylabel("Flux density (Jy)")
plt.title("B3, B6, B7")

"""
Explanation: Plot
-- same object will located in the same z, some of them will not have flux in all 3 bands.
End of explanation
"""


from astropy.cosmology import FlatLambdaCDM
cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)

"""
Explanation: Plot log(Luminosity) vs redshift
End of explanation
"""


def calc_power(z, flux):
    """
    z = redshift
    flux in Jy
    """
    z = np.array(z)
    flux = np.array(flux)
    
    dL = cosmo.luminosity_distance(z).to(u.meter).value # Luminosity distance
    luminosity = 4.0*np.pi*dL*dL/(1.0+z) * flux * 1e-26
    
    return z, luminosity

"""
Explanation: How to calculate luminosity:
$$L_{\nu} (\nu_{e}) = \frac{4 \pi D_{L}^2}{1+z} \cdot S_{\nu} (\nu_{o})$$
Notes: 
- Calculate Luminosity or Power in a specific wavelength (without k-correction e.g. using spectral index)
- $L_{\nu}$ in watt/Hz, in emited freq
- $S_{\nu}$ in watt/m$^2$/Hz, in observed freq
- $D_L$ is luminosity distance, calculated using astropy.cosmology function
- need to calculate distance in meter
- need to convert Jy to watt/m$^2$/Hz ----- $\times 10^{-26}$
End of explanation
"""


z3, l3 = calc_power(z3, f3)
z6, l6 = calc_power(z6, f6)
z7, l7 = calc_power(z7, f7)

plt.figure(figsize=(15,10))

plt.subplot(221)
plt.plot(z3, np.log10(l3), 'r*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B3")

plt.subplot(222)
plt.plot(z6, np.log10(l6), 'g*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B6")

plt.subplot(223)
plt.plot(z7, np.log(l7), 'b*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B7")

plt.subplot(224)
plt.plot(z3, np.log10(l3), 'r*', z6, np.log10(l6), 'g*', z7, np.log10(l7), 'b*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B3, B6, B7")

"""
Explanation: Plot $\log_{10}(L)$ vs $z$
End of explanation
"""


z3, l3 = calc_power(z3, f3)
z6, l6 = calc_power(z6, f6)
z7, l7 = calc_power(z7, f7)

plt.figure(figsize=(15,10))

plt.subplot(221)
plt.plot(z3, l3, 'r*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B3")

plt.subplot(222)
plt.plot(z6, l6, 'g*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B6")

plt.subplot(223)
plt.plot(z7, l7, 'b*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B7")

plt.subplot(224)
plt.plot(z3, l3, 'r*', z6, l6, 'g*', z7, l7, 'b*', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B3, B6, B7")

"""
Explanation: Without log10
End of explanation
"""


