<jupyter_start><jupyter_text>Need to Edit next panel for your own use  Look for Maps below in www.astro.caltech.edu/viero/Jason/cosmos/<jupyter_code>## Map Directories 

dir_spitzer_maps = '/home/agn/viero/Jason/cosmos/'
dir_spire_maps = '/home/agn/viero/Jason/cosmos/'
dir_pacs_maps = '/home/agn/viero/Jason/cosmos/'
dir_scuba_maps = '/home/agn/viero/Jason/cosmos/'
dir_aztec_maps = '/home/agn/viero/Jason/cosmos/'
dir_aztec_maps = '/home/agn/viero/Jason/cosmos/'

## Dictionary Names
library_keys =['mips24'
               ,'pacs_green'
               ,'pacs_red'
               ,'spire_PSW'
               ,'spire_PMW'
               ,'scuba_450'
               ,'spire_PLW'
               ,'scuba_850'
               ,'aztec'
              ]

wavelength=[24.,100,160,250,350,450,500,850,1100]
nwv = np.sum(wv0) 
fwhm =[6.32,7.4, 11.3,18.1, 25.2,7., 36.6,15.,18.]
efwhm=[6.32,6.7, 11.2,17.6, 23.9,7.8, 35.2,14.5,18.] # want to the measured effective FWHM later
color_correction=[1.25,23.58,23.82,1.018,0.9914,1e-3,0.95615,1e-3,1.0]
beam_area = [1.55e-09,1.,1.,1.,1.,1.62e-09,1.,5.6e-09,1.] #sr
beam_area = [1.55e-09,1.,1.,1.,1.,1.,1.,1.,1.] #sr

pixsize_suffix = '6.0_arcsec_pixels'

maps = [dir_spitzer_maps+'mips_24_GO3_sci_10.cutout.fits' 
        ,dir_pacs_maps+'pep_COSMOS_green_Map.DR1.sci.cutout.fits'
        ,dir_pacs_maps+'pep_COSMOS_red_Map.DR1.sci.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PSW.signal.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PMW.signal.cutout.fits'
        ,dir_scuba_maps+'map450_new_header.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PLW.signal.cutout.fits'
        ,dir_scuba_maps+'map850_new_header.cutout.fits'
        ,dir_aztec_maps+'cosmos_jcmt_kscott20100925_new_header_6.0_arcsec_pixels_map.fits'
       ]
noises = [dir_spitzer_maps+'mips_24_GO3_unc_10.cutout.fits' 
        ,dir_pacs_maps+'pep_COSMOS_green_Map.DR1.err.cutout.fits'
        ,dir_pacs_maps+'pep_COSMOS_red_Map.DR1.err.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PSW.noise.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PMW.noise.cutout.fits'
        ,dir_scuba_maps+'map450_new_header_rms.cutout.fits'
        ,dir_spire_maps+'cosmos-uvista-hipe12_itermap_10_iterations_'+pixsize_suffix+'_PLW.noise.cutout.fits'
        ,dir_scuba_maps+'map850_new_header_rms.cutout.fits'
        ,dir_aztec_maps+'cosmos_jcmt_kscott20100925_new_header_6.0_arcsec_pixels_noise.fits'
       ]<jupyter_output><empty_output><jupyter_text>Reading Maps into Objects, and Objects into Dictionaries<jupyter_code>sky_library={}

for t in indstack[0]:
    sky_library[library_keys[t]] = Skymaps(maps[t],noises[t],efwhm[t],color_correction=color_correction[t])
    sky_library[library_keys[t]].add_wavelength(wavelength[t])
    sky_library[library_keys[t]].add_fwhm(efwhm[t]) 
    if beam_area[t] != 1: sky_library[library_keys[t]].beam_area_correction(beam_area[t])<jupyter_output><empty_output><jupyter_text>Catalog Directories and Names Look for Maps below in www.astro.caltech.edu/viero/Jason/cosmos/<jupyter_code>path_catalog = '/data/maps_cats_models/catalogs/UVISTA/'
#Catalog must have the following (case sensitive) Columns:
# - RA
# - DEC
# - z_peak
# - LMASS
# - rf_U_V
# - rf_V_J
#ALPHA_J2000,DELTA_J2000,ZPDF,ZPDF_L68,ZPDF_H68,CLASS,MASS_MED,MASS_MED_MIN68,MASS_MED_MAX68
file_catalog = 'COSMOS2015_Laigle+_v1.1.csv'
tbl = pd.read_table(path_catalog+file_catalog,sep=',')
tbl['ID']=range(len(tbl['ALPHA_J2000']))
tbl['ra']=tbl['ALPHA_J2000']
tbl['dec']=tbl['DELTA_J2000']
tbl['z_peak']=tbl['ZPDF']
tbl['sfg']=tbl['CLASS']
tbl['LMASS']=tbl['MASS_MED']
tbl.keys()
uVista = Field_catalogs(tbl)
#Stack in redshift bins, with Layers divided by Stellar Mass and Star-Forming/Quiescent Galaxies
pop = ['sf','qt']
npop=2
all_stacked_fluxes = np.zeros([nwv,nz,nm,npop])
all_luminosity_temp = np.zeros([nz,nm,npop,2])

for iz in range(nz):
    zn = z_nodes[iz:iz+2]
    uVista.get_sf_qt_mass_redshift_bins(zn,m_nodes)
    radec_m_z_p = uVista.subset_positions(uVista.id_z_ms)
    stacked_fluxes_  =  None
    n_sources_max = None
    
    stacked_fluxes = stack_libraries_in_redshift_slices(
        sky_library,
        radec_m_z_p)
    
    #pdb.set_trace()
    args = radec_m_z_p.keys()
    for iwv in range(nwv):
        stacked_fluxes_wv = stacked_fluxes[str(wavelength[indstack[0][iwv]])]
        for j in range(nm):
            for k in range(npop):
                arg = 'z_'+str(zn[0])+'-'+str(zn[1])+'__m_'+str(m_nodes[j])+'-'+str(m_nodes[j+1])+'_'+pop[k]
                all_stacked_fluxes[iwv,iz,j,k] = stacked_fluxes_wv[arg.replace('.','p').replace('-','_')].value
    
    #pdb.set_trace()
    #PLOT
    plt.figure()
    plt.ylim([-1e-1,3e0])
    plt.ylim([1e-2,9e1])
    plt.xlim([20,510])
    plt.xlim([20,1800])
    plt.xlim([20,2500])
    plt.yscale('log')
    plt.ylabel('Flux Density [mJy]',fontsize=14)
    plt.xlabel('Wavelength [microns]',fontsize=14)
    ln=['','--']
    for k in range(npop):
        for j in range(nm)[::-1]: 
            arg = 'z_'+str(round(zn[0],3))+'-'+str(round(zn[1],3))+'__m_'+str(round(m_nodes[j],3))+'-'+str(round(m_nodes[j+1],3))+'_'+pop[k]
            ng = len(radec_m_z_p[arg.replace('.','p').replace('-','_')][0])
            plt.plot(np.array(wavelength)[indstack],1e3*all_stacked_fluxes[:,iz,j,k],ln[k],color=popcolor[npop*j+k],label=[str(ng)+' '+pop[k],str(m_mid[j])])
            plt.legend()
    plt.show()
    #pdb.set_trace()
def nghack(a):
    b = a
    for i in range(len(a)):
        if a[i] < 0:
            b[i] = 1e-5
    return b
for iwv in range(nwv):
    plt.figure()
    for j in range(nm):
        plt.yscale('log')
        plt.ylabel('Flux Density [mJy]',fontsize=14)
        plt.xlabel('Redshift',fontsize=14)
        plt.title('Star-Forming Galaxies, '+str(wavelength[indstack[0][iwv]]))
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]),'o')
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]))
    plt.show()
for iwv in range(nwv):
    plt.figure()
    for j in range(nm):
        plt.yscale('log')
        plt.ylabel('Flux Density [mJy]',fontsize=14)
        plt.xlabel('Redshift',fontsize=14)
        plt.title('Star-Forming Galaxies, '+str(wavelength[indstack[0][iwv]]))
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]),'o')
        plt.plot(z_mid,1e3*nghack(all_stacked_fluxes[iwv,:,j,0]))
    plt.show()
for iwv in range(nwv):
    plt.figure()
    for j in range(nm):
        plt.yscale('log')
        plt.ylim([1e-2,9e1])
        plt.ylabel('Flux Density [mJy]',fontsize=14)
        plt.xlabel('Redshift',fontsize=14)
        plt.title('Quiescent Galaxies, '+str(wavelength[indstack[0][iwv]]))
        plt.plot(z_mid,1e3*(all_stacked_fluxes[iwv,:,j,1]),'o')
        plt.plot(z_mid,1e3*(all_stacked_fluxes[iwv,:,j,1]))
    plt.show()
wvs = np.array(wavelength)[indstack]
luminosities = np.zeros([nz,nm])
temperatures = np.zeros([nz,nm])
luminosity_dictionary = {}
lumfit_dictionary = {}

for iz in range(nz):
    zn = z_nodes[iz:iz+2] 
    plt.ylim([9,14])
    for j in range(nm)[::-1]:
        flux_wv = (all_stacked_fluxes[:,iz,j,0] )   
        m = fast_sed_fitter(wvs,flux_wv,covar=np.sqrt(abs(flux_wv)))
        #sed = fast_sed(m, wvs)
        L_mod = fast_Lir(m, np.mean(zn))
        luminosities[iz,j] = L_mod.value
        temperatures[iz,j] = m['T_observed'].value * np.mean(zn)
        #print temperatures[iz,k] 
        
for j in range(nm)[::-1]:
    luminosity_dictionary['L_sed'+str(j)] = luminosities[:,j]
    luminosity_dictionary['T_sed'+str(j)] = temperatures[:,j]
        #print LIR.value
        #plt.plot(np.mean(zn),np.log10(LIR.value),'o')
    plt.plot(z_mid,np.log10(luminosity_dictionary['L_sed'+str(j)]))
        #plt.yscale('log')
        
    plt.ylabel('Luminosity [L_star]',fontsize=14)
    plt.xlabel('Redshift',fontsize=14)
#plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1730745456793176, 'fraction_numerical': 0.04413400914822599, 'mean_word_length': 4.0, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '166375', 'license': 'MIT', 'n_tokens_mistral': 3759, 'n_tokens_neox': 3500, 'n_words': 366, 'path': 'notebooks/.ipynb_checkpoints/simstack_with_libraries_example_Laigle-checkpoint.ipynb', 'repo_name': 'marcoviero/simstack', 'size': '8089'}