<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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.
Step1:  Example, retrieve all the calibrator with a flux > 0.1 Jy
Step2:  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
Step3:  We can write a "report file" or only use the "resume data", some will have redshift data retrieved from NED.
Step4:  Sometimes there is no redshift information found in NED
Retrieve object nonALMACAL and ALMACAL and combining data in resume and listcal
Step5:  Select objects which has redshift
collect the flux, band, freq, and obsdate
plot based on the Band
Step6:  Plot Flux vs Redshift, NonALMACAL vs ALMACAL
same object will located in the same z
some of them will not have flux in all 3 bands.
nonALMACAL = red
ALMACAL = blue
Step7:  Distribution of redshift
only B3
Step8:  Plot log(Luminosity) vs redshift, nonALMACAL vs ALMACAL
Step10:  How to calculate luminosity
Step11:  Plot $\log_{10}(L)$ vs $z$
Step12:  Black-dashed line are for 0.1 Jy flux.
Distribution for Luminosity
Only Band 3
<ASSISTANT_TASK:>
Python Code:
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:60., 6:60., 7:60.}, \
                                    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
report_non, resume_non = q.select_object_from_sqldb("calibrators_brighterthan_0.1Jy_20180419.db", \
                                    maxFreqRes=999999999, array='12m', \
                                    excludeCycle0=True, \
                                    selectPol=False, \
                                    minTimeBand={3:60., 6:60., 7:60.}, \
                                    nonALMACAL = True, \
                                    silent=True)
report_almacal, resume_almacal = q.select_object_from_sqldb("calibrators_brighterthan_0.1Jy_20180419.db", \
                                    maxFreqRes=999999999, array='12m', \
                                    excludeCycle0=True, \
                                    selectPol=False, \
                                    minTimeBand={3:60., 6:60., 7:60.}, \
                                    onlyALMACAL = True, \
                                    silent=True)
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"
for i, obj in enumerate(resume_non):
    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"
for i, obj in enumerate(resume_almacal):
    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
Retrieve object nonALMACAL and ALMACAL and combining data in resume and listcal
End of explanation
def collect_z_and_flux(Band, resume_array):
    z = []
    flux = []
    for idata in resume_array:
        if idata[6] is not None: # select object which has redshift information
            fluxnya = idata[-1][0]
            bandnya = idata[-1][1]
            freqnya = idata[-1][2]
            datenya  = idata[-1][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(3, resume)
print("Number of seleted source in B3: ", len(z3))
z6, f6 = collect_z_and_flux(6, resume)
print("Number of seleted source in B6: ", len(z6))
z7, f7 = collect_z_and_flux(7, resume)
print("Number of seleted source in B7: ", len(z7))
znon3, fnon3 = collect_z_and_flux(3, resume_non)
print("Number of seleted source in B3: ", len(znon3))
znon6, fnon6 = collect_z_and_flux(6, resume_non)
print("Number of seleted source in B6: ", len(znon6))
znon7, fnon7 = collect_z_and_flux(7, resume_non)
print("Number of seleted source in B7: ", len(znon7))
zalm3, falm3 = collect_z_and_flux(3, resume_almacal)
print("Number of seleted source in B3: ", len(zalm3))
zalm6, falm6 = collect_z_and_flux(6, resume_almacal)
print("Number of seleted source in B6: ", len(zalm6))
zalm7, falm7 = collect_z_and_flux(7, resume_almacal)
print("Number of seleted source in B7: ", len(zalm7))
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(znon3, fnon3, 'r.', zalm3, falm3, 'b.')
plt.xlabel("z")
plt.ylabel("Flux density (Jy)")
plt.title("B3")
plt.subplot(222)
plt.plot(znon6, fnon6, 'r.', zalm6, falm6, 'b.')
plt.xlabel("z")
plt.ylabel("Flux density (Jy)")
plt.title("B6")
plt.subplot(223)
plt.plot(znon7, fnon7, 'r.', zalm7, falm7, 'b.')
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 Flux vs Redshift, NonALMACAL vs ALMACAL
same object will located in the same z
some of them will not have flux in all 3 bands.
nonALMACAL = red
ALMACAL = blue
End of explanation
# from matplotlib import rcParams
# rcParams['font.family'] = ['Humor Sans', 'Comic Sans MS']
# rcParams['font.size'] = 14.0
# rcParams['axes.linewidth'] = 1.5
# rcParams['lines.linewidth'] = 2.0
# rcParams['figure.facecolor'] = 'white'
# rcParams['grid.linewidth'] = 0.0
# rcParams['axes.unicode_minus'] = False
# rcParams['xtick.major.size'] = 8
# rcParams['xtick.major.width'] = 3
# rcParams['ytick.major.size'] = 8
# rcParams['ytick.major.width'] = 3
plt.xkcd()
plt.figure(figsize=(10,6))
#plt.subplot(121)
plt.hist(z3, bins=20, label="all sample", histtype="step", linewidth=4)
plt.hist(znon3, bins=20, label="our final target", histtype="step", linewidth=4)
plt.hist(zalm3, bins=20, label="ALMACAL I", histtype="step", linewidth=4)
plt.xlabel(r"Redshift $(z)$")
plt.ylabel("Number of Calibrator")
plt.title("Fig 1. Redshift distribution of sample")
plt.legend()
# plt.subplot(122)
# plt.hist(z3, bins=15, label="all", normed=True, cumulative=True, histtype='step')
# plt.hist(znon3, bins=15, label="nonALMACAL", normed=True, cumulative=True, histtype='step')
# plt.hist(zalm3, bins=15, label="ALMACAL", normed=True, cumulative=True, histtype='step')
# plt.legend()
# plt.xlabel(r"$z$")
# plt.ylabel("Normalized-Cumulative N")
# plt.show()
plt.savefig("sample_redshift.png", transparent=True, bbox_tight=True)
plt.title?
Explanation: Distribution of redshift
only B3
End of explanation
from astropy.cosmology import FlatLambdaCDM
cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)
Explanation: Plot log(Luminosity) vs redshift, nonALMACAL vs ALMACAL
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)
znon3, lnon3 = calc_power(znon3, fnon3)
znon6, lnon6 = calc_power(znon6, fnon6)
znon7, lnon7 = calc_power(znon7, fnon7)
zalm3, lalm3 = calc_power(zalm3, falm3)
zalm6, lalm6 = calc_power(zalm6, falm6)
zalm7, lalm7 = calc_power(zalm7, falm7)
zdummy = np.linspace(0.001, 2.5, 100)
fdummy = 0.1 # Jy, our cut threshold
zdummy, Ldummy0 = calc_power(zdummy, fdummy)
zdummy, Ldummy3 = calc_power(zdummy, np.max(f3))
zdummy, Ldummy6 = calc_power(zdummy, np.max(f6))
zdummy, Ldummy7 = calc_power(zdummy, np.max(f7))
plt.figure(figsize=(15,10))
plt.subplot(221)
plt.plot(znon3, np.log10(lnon3), 'ro', zalm3, np.log10(lalm3), 'bo', \
         zdummy, np.log10(Ldummy0), 'k--', zdummy, np.log10(Ldummy3), 'r--', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B3")
plt.subplot(222)
plt.plot(znon6, np.log10(lnon6), 'ro', zalm6, np.log10(lalm6), 'bo', \
         zdummy, np.log10(Ldummy0), 'k--', zdummy, np.log10(Ldummy6), 'g--', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B6")
plt.subplot(223)
plt.plot(znon7, np.log10(lnon7), 'ro', zalm7, np.log10(lalm7), 'bo', \
         zdummy, np.log10(Ldummy0), 'k--', zdummy, np.log10(Ldummy7), 'b--', alpha=0.5)
plt.xlabel(r"$z$"); plt.ylabel(r"$\log_{10}(L_{\nu_e})$"); plt.title("B7")
# plt.subplot(224)
# plt.plot(znon3, np.log10(lnon3), 'r*', znon6, np.log10(lnon6), 'g*', znon7, np.log10(lnon7), 'b*', \
#          zdummy, np.log10(Ldummy0), 'k--', \
#          zdummy, np.log10(Ldummy3), 'r--', \
#          zdummy, np.log10(Ldummy6), 'g--', \
#          zdummy, np.log10(Ldummy7), '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
plt.figure(figsize=(15,5))
plt.subplot(121)
plt.hist(l3, bins=20, label="all", histtype="step")
plt.hist(lnon3, bins=20, label="nonALMACAL", histtype="step")
plt.hist(lalm3, bins=20, label="ALMACAL", histtype="step")
plt.legend()
plt.xlabel(r"$L_{\nu}$")
plt.ylabel("N")
plt.subplot(122)
plt.hist(l3, bins=15, label="all", normed=True, cumulative=True, histtype='step')
plt.hist(lnon3, bins=15, label="nonALMACAL", normed=True, cumulative=True, histtype='step')
plt.hist(lalm3, bins=15, label="ALMACAL", normed=True, cumulative=True, histtype='step')
plt.legend()
plt.xlabel(r"$L_{\nu}$")
plt.ylabel("Normalized-Cumulative N")
plt.show()
from scipy.stats import ks_2samp
print(ks_2samp(z3, znon3))
print(ks_2samp(z3, zalm3))
print(ks_2samp(znon3, zalm3))
print(ks_2samp(l3, lnon3))
print(ks_2samp(l3, lalm3))
print(ks_2samp(lnon3, lalm3))
Explanation: Black-dashed line are for 0.1 Jy flux.
Distribution for Luminosity
Only Band 3
End of explanation
<END_TASK>