<jupyter_start><jupyter_text>Gawa Simulations This jn is intended to create simulations of dwarf galaxies and globular clusters using as field stars the catalog of DES. These simulations will be later copied to gawa jn, a pipeline to detect stellar systems with field's stars. In principle this pipeline read a table in data base with g and r magnitudes, subtract the extinction in each band, and randomize the positions in RA and DEC in order to avoid stellar systems in the FoV. The star clusters are inserted later, centered in each HP pixel with specific nside.To complete all the steps you just have to run all the cells below in sequence. Firstly, install the packages not available in the image via terminal. Restart the kernel and so you can run the cell bellow.<jupyter_code># Rodar no terminal e depois reiniciar o kernel:
# python -m pip install -U --user pip
# python -m pip install -U --user scikit-image
# pip install --user parsl

import os
from os import system
import healpy as hp
import numpy as np
import astropy.io.fits as fits
from astropy import units as u
import astropy.coordinates as coord
from astropy.coordinates import SkyCoord
from astropy.coordinates import match_coordinates_sky
from astropy.io.fits import getdata
from astropy.table import Table
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.path as mpath
from matplotlib.ticker import ScalarFormatter
import collections
from matplotlib import rc
import warnings
import sqlalchemy
warnings.simplefilter('ignore', UserWarning)
%matplotlib inline
from scipy.stats import expon
from itertools import compress
import parsl
from parsl.app.app import python_app, bash_app
from parsl.configs.local_threads import config
parsl.load(config)
mpl.rcParams['legend.numpoints'] = 1
cmap = mpl.cm.get_cmap("inferno_r")
cmap.set_under('dimgray')<jupyter_output><empty_output><jupyter_text>Below are the items of the configuration for field stars and simulations. A small description follows as a comment.<jupyter_code># Main settings:
nside_ini = 64 # the nside in which pixel star clusters will be inserted in the center
nside_ftp = 4096 # nside of the footprint map
nside3 = 2 ** 18 # pixelization of the star's position
ra_min, ra_max = 45, 55 # 5., 85. # limits in ra of the simulated map
dec_min, dec_max = -30, -20 # -60., -20. # limits in dec of the simulated map
mmin, mmax, cmin, cmax = 17.00, 24.00, -0.4, 1.6 # magnitude and color limits of stars
ell_min, ell_max = 0.00, 0.20 # limits in ellipticity of star clusters simulated
pa_min, pa_max = 0., 180. # limits in position angle
log10_mass_min, log10_mass_max = 3.00, 3.10 # limits for the visible mass of the clusters,
# avoiding simulate many stars unseen in the catalog given the limiting magnitude
mM_min, mM_max = 20.001, 20.002 # limits in modulus distance
log10_rexp_min, log10_rexp_max = 0.60, 1.00 # limits in exponential radius
survey = 'des' # survey name
frac_bin = 0.5 # fraction of binaries (amount of stars in binaries = Nbinaries / Ntotal)
IMF_author = 'Kroupa' # selection of initial mass function
border_extract = 1. # remove star clusters in the edge of the map, in degrees
file_iso = 'iso_input2.dat' # file with the isochrone
iso_info = np.loadtxt(file_iso, usecols=(1, 2), unpack=True) 
# just to print the metalicity and age of the input file:
FeH_iso = iso_info[0][0]
logAge_iso = iso_info[1][0]
hdu_ngp = fits.open("/archive/external_catalogs/SDSS3/schlegel/schlegel/SFD_dust_4096_ngp.fits", memmap=True)
ngp = hdu_ngp[0].data
hdu_sgp = fits.open("/archive/external_catalogs/SDSS3/schlegel/schlegel/SFD_dust_4096_sgp.fits", memmap=True)
sgp = hdu_sgp[0].data<jupyter_output><empty_output><jupyter_text>Now, printing the age and metalicity of the model:<jupyter_code>print('[M/H]={:.2f}, Age={:.2f} Gyr'.format(FeH_iso, 10**(logAge_iso-9)))<jupyter_output>[M/H]=-2.00, Age=13.49 Gyr<jupyter_text>Bellow there are many functions that will be used later. A small description follows after the definition of each function.<jupyter_code>def dist_ang(ra1, dec1, ra_ref, dec_ref):
    """ Calculated the angular distance between (ra1, dec1) and (ra_ref, dec_ref)
    ra-dec in degrees
    ra1-dec1 can be arrays
    ra_ref-dec_ref are scalars
    output is in degrees
    Parameters
    ----------
    ra1, dec1 : lists or arrays
        Positions in the sky
    ra_ref, dec_ref : floats
        Reference position in the sky
    
    Returns
    -------
    list
        a list of floats with the angular distance of the first couple of 
        inputs to the position of reference.

    """
    costheta = np.sin(np.radians(dec_ref)) * np.sin(np.radians(dec1)) + \
               np.cos(np.radians(dec_ref)) * np.cos(np.radians(dec1)) * np.cos(np.radians(ra1-ra_ref))
    dist_ang = np.arccos(costheta)
    return np.rad2deg(dist_ang) # degrees<jupyter_output><empty_output><jupyter_text>Add a function to return the extinction in the V band since the catalog used for detection is extinction free. In the future, may be read the magnitude already corrected by extinction in the VAC or catalog.<jupyter_code>def get_av(gal_l, gal_b):
    """Return extinction (A) in V band based on l, b position. Rv is defined as a constant and equal to 3.1.

    Parameters
    ----------
    gal_l : list
        Galactic longitude of the objects (degrees)
    gal_b : list
        Galactic latitude of the objects (degrees)

    Returns
    -------
    av : list
        a list of Galactic extinction in the V band to each position
    """
    lt = np.radians(gal_l)
    bt = np.radians(gal_b)

    n = np.array(bt > 0).astype(int)
    n[n == 0] = -1
    x = 2048 * ((1 - (n * np.sin(bt))) ** (0.5)) * np.cos(lt) + 2047.5
    y = -2048 * n * ((1 - (n * np.sin(bt))) ** (0.5)) * np.sin(lt) + 2047.5

    x, y = np.round(x, 0).astype(int), np.round(y, 0).astype(int)
    av = np.empty(len(gal_l))

    for idx, num in enumerate(bt):
        if num > 0:
            av[idx] = 3.1 * ngp[y[idx], x[idx]]
        else:
            av[idx] = 3.1 * sgp[y[idx], x[idx]]
    '''
    for i in range(len(bt)):
        if bt[i] > 0:
            av[i] = 3.1 * ngp[y[i], x[i]]
        else:
            av[i] = 3.1 * sgp[y[i], x[i]]
    '''
    return av
def make_footprint():
    """Creates a partial HealPix map based on the area selected,
    with nside = nside_ftp

    Parameters
    ----------
    ra_min, ra_max, dec_min, dec_max : float (global)
        Limits in ra and dec (degrees)
    nside_ftp : int (global)
        Nside of the footprint map.

    Returns
    -------
    list
        a list of ipixels in the area selected (inclusive=False)
    """

    global ra_min, ra_max, dec_min, dec_max, nside_ftp

    vertices = hp.ang2vec([ra_min, ra_max, ra_max, ra_min], [dec_min, dec_min, dec_max, dec_max], lonlat=True)

    hp_sample = hp.query_polygon(nside_ftp, vertices, inclusive=False, fact=64, nest=True, buff=None)

    # m = np.bincount(hp_sample, minlength=hp.nside2npix(nside_ftp))

    # hp.mollview(m, nest=True, flip='astro')
    SIGNAL = np.ones(len(hp_sample))

    col0 = fits.Column(name='HP_PIXEL_NEST_4096', format='J', array=hp_sample)
    col1 = fits.Column(name='SIGNAL', format='E', array=SIGNAL)

    cols = fits.ColDefs([col0, col1])
    tbhdu = fits.BinTableHDU.from_columns(cols)
    tbhdu.writeto('ftp_4096_nest.fits', overwrite=True)
    
    return hp_sample
    
hpx_ftp = make_footprint()<jupyter_output><empty_output><jupyter_text>Reading the catalog and writing as a fits file (to avoid read from the DB many times in the case the same catalog will be used multiple times).<jupyter_code>tablename = 'vac_ga_y6.catalog_6048'
engine = sqlalchemy.create_engine('postgresql://untrustedprod:<EMAIL>:5432/prod_gavo')
conn = engine.connect()

query = 'select ra, dec, mag_g, magerr_g, mag_r, magerr_r from %s where (ra > %s) and (ra <%s) and (dec > %s) and (dec < %s)' % (tablename, ra_min, ra_max, dec_min, dec_max)
stm = sqlalchemy.sql.text(query)
stm_get = conn.execute(stm)   
stm_result = stm_get.fetchall()  
table = Table(rows=stm_result, names=('ra', 'dec', 'mag_g', 'magerr_g', 'mag_r', 'magerr_r'))

RA = np.array(table['ra'])
DEC = np.array(table['dec'])
MAG_G = np.array(table['mag_g'])
MAGERR_G = np.array(table['magerr_g'])
MAG_R = np.array(table['mag_r'])
MAGERR_R = np.array(table['magerr_r'])
print(len(RA))
'''
dist = np.zeros(len(RA))
for i in range(len(RA)):
    cond = (RA > RA[i] - 0.1)&(RA < RA[i] + 0.1)&(DEC > DEC[i] - 0.1)&(DEC < DEC[i] + 0.1)
    xxx = np.sort(dist_ang(RA[cond], DEC[cond], RA[i], DEC[i]))
    dist[i] = xxx[1]
plt.hist(dist, bins=100)
plt.show()
'''<jupyter_output>542317<jupyter_text>The cells below reads the position, calculates the extinction using the previous function and correct the aparent magnitude (top of the Galaxy), filter the stars for magnitude and color ranges,  and writes a file with the original position of the stars and corrected magnitudes.<jupyter_code>c = SkyCoord(ra=RA * u.degree, dec=DEC * u.degree, frame='icrs')
L = c.galactic.l.degree
B = c.galactic.b.degree

MAG_G -= 1.19863 * get_av(L, B)
MAG_R -= 0.83734* get_av(L, B)

cond = (MAG_G < mmax)&(MAG_G > mmin)&(MAG_G-MAG_R > cmin)&(MAG_G-MAG_R < cmax)

RA = RA[cond]
DEC = DEC[cond]
MAG_G = MAG_G[cond]
MAG_R = MAG_R[cond]
MAGERR_G = MAGERR_G[cond]
MAGERR_R = MAGERR_R[cond]

col1 = fits.Column(name='RA', format='D', array=RA)
col2 = fits.Column(name='DEC', format='D', array=DEC)
col3 = fits.Column(name='MAG_G', format='E', array=MAG_G)
col4 = fits.Column(name='MAG_R', format='E', array=MAG_R)
col5 = fits.Column(name='MAGERR_G', format='E', array=MAGERR_G)
col6 = fits.Column(name='MAGERR_R', format='E', array=MAGERR_R)

cols = fits.ColDefs([col1, col2, col3, col4, col5, col6])
tbhdu = fits.BinTableHDU.from_columns(cols)
tbhdu.writeto('DES_Y6_Gold_v1.fits', overwrite=True)
'''
hdu = fits.open("DES_Y6_Gold_v1.fits", memmap=True)
RA = hdu[1].data.field('RA')
DEC = hdu[1].data.field('DEC')
MAG_G = hdu[1].data.field('MAG_G')
MAG_R = hdu[1].data.field('MAG_R')
MAGERR_G = hdu[1].data.field('MAGERR_G')
MAGERR_R = hdu[1].data.field('MAGERR_R')
'''<jupyter_output><empty_output><jupyter_text>This is the function to distribute stars in the real catalog (stars not belonging to the simulated clusters).<jupyter_code>def d_star_real_cat(hpx_ftp, len_):
    """ This function distributes the set of stars with lenght = len_
    in a catalog based on the set of hpx_ftp (healpix footprint).
    
    Parameters
    ----------
    hpx_ftp : list
        The set of ipixels in the footprint
    len_ : int
        Total amount of stars in the real catalog

    Returns
    -------
    ra_mw_stars, dec_mw_stars : lists
        The position of the stars in the catalog (degrees)
    """

    global nside3, nside_ftp
    
    f2 = nside3 / nside_ftp
    A = np.repeat(hpx_ftp, f2 ** 2)
    a = np.random.choice(np.arange(f2 ** 2), len(hpx_ftp), replace=True)
    set_pixels_nside3 = A * (f2 ** 2) + a
    hpx_star = np.random.choice(set_pixels_nside3, len_, replace=False)
    np.random.shuffle(hpx_star)
    ra_mw_stars, dec_mw_stars = hp.pix2ang(nside3, hpx_star, nest=True, lonlat=True)
    return ra_mw_stars, dec_mw_stars<jupyter_output><empty_output><jupyter_text>Simulation of dwarf galaxies and globular clustersIn fact, the dwarf galaxies and globular clusters are very similar in terms of stellar populations. Dwarf galaxieshave a half-light radius larger than globular clusters (given the amount of dark matter) with the same absolute magnitude. The code below simulates stars using a Kroupa or Salpeter IMF, and an exponential radius for the 2D distribution of stars.<jupyter_code>def IMF_(author):
    """ Defines dictionary for Kroupa and Salpeter initial mass functions.

    Parameters
    ----------
    author : str
        The name of the initial mass function (IMF)

    Returns
    -------
    dictionary
        a dict with the alpha values and mass breaks
    """
    if author == 'Kroupa':
        return {
               'IMF_alpha_1' : -1.3,
               'IMF_alpha_2'   : -2.3,
               'IMF_mass_break': 0.5
               }
    if author == 'Salpeter':
        return {
               'IMF_alpha_1' : -2.3,
               'IMF_alpha_2'   : -2.3,
               'IMF_mass_break': 0.5
               }


def apply_err(mag, mag_table, err_table):
    """ This function returns magnitude errors for the 'mag' variable
    based on mag_table and err_table.
    
    Parameters
    ----------
    mag : list
        The list of magnitudes to be calculated
    mag_table : list
        List of magnitudes
    err_table : List
        List of magnitude errors (1-sigma) respective to mag_table

    Returns
    -------
    list
        a list of magnitude errors following Normal distribution with
        1-sigma error as informed
    """
    err_interp = np.interp(mag, mag_table, err_table)
    return np.abs(err_interp * np.random.randn(len(err_interp)))


def faker_bin(total_bin, IMF_author, file_in, dist):
    """Calculates the fraction of binaries in the simulated clusters.

    Parameters
    ----------
    total_bin : float
        The amount of binaries. Definition: N is the total amount of
        stars (take care to count a system of a binary as two stars),
        and B is the amount of stars in binary systems.
        so bin_frac = B / N
    IMF_author : str
        Name of the IMF (see function above)
    file_in : str
        The name of the file with star's masses and magnitudes
    dist : float
        Distance cluster-observer in parsecs

    Returns
    -------
    binaries[:,0]
        a list of magnitudes of the binaries in the first band
    binaries[:,1]
        a list of magnitudes of the binaries in the second band
    """
    global cmin, cmax, mmin, mmax

    mass, mag1, mag2 = np.loadtxt(file_in, usecols=(3, 29, 30), unpack=True)

    # bin in mass (solar masses)
    binmass = 5.e-4

    mag1 += 5 * np.log10(dist) - 5
    mag2 += 5 * np.log10(dist) - 5
    
    IMF = IMF_(IMF_author)

    # amostra is an array with the amount of stars in each bin of mass. ex.: [2,3,4,1,2]
    massmin = np.min(mass)
    massmax = np.max(mass)
    bins_mass = int((massmax-massmin)/binmass)
    amostra = np.zeros(bins_mass)

    for i in range(bins_mass):
        if ((i*binmass) + massmin <= IMF['IMF_mass_break']):
            amostra[i] = round((massmin + i * binmass)**(IMF['IMF_alpha_1']))
        else:
            amostra[i] = round((massmin + i * binmass)**(IMF['IMF_alpha_2']))
    # Soma is the total amount of stars (float), the sum of amostra
    soma = np.sum(amostra)
    # Now normalizing the array amostra
    #for idx, num in enumerate(amostra):
    amostra = np.multiply(amostra, total_bin / soma)

    massa_calculada = np.zeros(int(total_bin))

    count = 0

    for j in range(bins_mass): # todos os intervalos primarios de massa
        for k in range(int(amostra[j])): # amostra() eh a amostra de estrelas dentro do intervalo de massa
            massa_calculada[count] = massmin + (j * binmass) + (k * binmass/ amostra[j])
            # massa calculada eh a massa de cada estrela
            count += 1
            
    # mag1 mag1err unc1 mag2 mag2err unc2
    binaries = np.zeros((total_bin, 2))

    for i in range(total_bin):
        for k in range(len(mass)-1): # abre as linhas do arquivo em massa
            # se a massa estiver no intervalo das linhas
            if (mass[k] < massa_calculada[i])&(mass[k+1] > massa_calculada[i]):
            # vai abrir tantas vezes quantas forem as estrelas representadas
                intervalo = (massa_calculada[i] - mass[k]) / (mass[k+1] - mass[k]) # intervalo entre zero e um
                binaries[i,0] = mag1[k] - (mag1[k] - mag1[k+1]) * intervalo
                binaries[i,1] = mag2[k] - (mag2[k] - mag2[k+1]) * intervalo
    return binaries[:,0], binaries[:,1]

def unc(mag, mag_table, err_table):
    """Interpolates the uncertainty in magnitude for a specific magnitude
    using magnitude and error from table.

    Parameters
    ----------
    mag : float or list
        The magnitude to be interpolated
    mag_table : list
        List of magnitudes in table
    err_table : list
        List of magnitude errors in table
    
    Returns
    -------
    err_interp : float or list
        Magnitudes interpolated
    """
    err_interp = np.interp(mag, mag_table, err_table)
    return err_interp


@python_app
def faker(N_stars_cmd, frac_bin, IMF_author, x0, y0, rexp, ell_, pa, dist, hpx):
    """Creates an array with positions, magnitudes, magnitude errors and magnitude
    uncertainties for the simulated stars in two bands. The stars belong to a simple
    stellar population and they are spatially distributed following an exponential profile.
    The code firstly simulates the stars in the CMDs and finally simulates only the
    companions of the binaries (it does not matter how massive the companions are) to
    join to the number of points in the CMD.
    Bear in mind these quantities (definitions):
    N_stars_cmd = number of stars seen in the CMD. The binaries are seen as a single star.
    N_stars_single = amount of stars that are single stars.
    N_stars_bin = amount of stars that are binaries in the CMD. For each of these kind of
    stars, a companion should be calculated later.

    Parameters
    ----------
    N_stars_cmd : int
        Points in simulated cmd given the limiting magnitude. Some of this stars are single,
        some of the are binaries. This amount obeys the following relation:
        N_stars_cmd = N_stars_single + N_stars_bin, where N_stars_single are the single stars
        in the cmd and N_stars_bin are the points in CMD that are binaries. A single star in
        each system are accounted for. In this case, the total amount of stars simulated is
        N_stars_single + 2 * N_stars_bin
    frac_bin : float (0-1)
        Fraction of binaries. This is the total amount of stars in the CMD that belongs to a
        binary system (= 2 * N_stars_bin / total amount of stars).
    IMF_author : str
        Name of the IMF (see function above)
    x0 : float (degrees)
        RA position of the center of cluster
    y0 : float(degrees)
        DEC position of the center of cluster
    rexp : float (degrees)
        Exponential radii of the cluster following the exponential law of density:
        N = A * exp(-r/rexp)
    ell_ : float
        Ellipticity of the cluster (ell_=sqrt((a^2-b^2)/(a^2)))
    pa : float
        Positional angle (from North to East), in degrees
    dist : float
        Distance to the cluster in parsecs
    hpx : int
        Pixel where the cluster resides (nested)
    
    """

    global cmin, cmax, mmin, mmax, mag_, err1_, err2_, file_iso

    mass, mag1, mag2 = np.loadtxt(file_iso, usecols=(3, 29, 30), unpack=True)

    # bin in mass (solar masses)
    binmass = 5.e-4

    mag1 += 5 * np.log10(dist) - 5
    mag2 += 5 * np.log10(dist) - 5

    # Warning: cut in mass to avoid faint stars with high errors showing up in the
    # bright part of magnitude. The mass is not the total mass of the cluster,
    # only a lower limit for the total mass.
    cond = (mag1 <= mmax + 0.5)
    mass, mag1, mag2 = mass[cond], mag1[cond], mag2[cond]

    # amostra is an array with the amount of stars in each bin of mass. ex.: [2,3,4,1,2]
    massmin = np.min(mass)
    massmax = np.max(mass)
    bins_mass = int((massmax-massmin)/binmass)
    amostra = np.zeros(bins_mass)

    IMF = IMF_(IMF_author)

    for i in range(bins_mass):
        if ((i*binmass) + massmin <= IMF['IMF_mass_break']):
            amostra[i] = round((massmin + i * binmass)**(IMF['IMF_alpha_1']))
        else:
            amostra[i] = round((massmin + i * binmass)**(IMF['IMF_alpha_2']))
    # Soma is the total amount of stars (float), the sum of amostra
    soma = np.sum(amostra)
    # Now normalizing the array amostra
    for i in range(len(amostra)):
        amostra[i] = N_stars_cmd * amostra[i] / soma
    
    massa_calculada = np.zeros(int(N_stars_cmd))
    
    count = 0
  
    for j in range(bins_mass): # todos os intervalos primarios de massa
        for k in range(int(amostra[j])): # amostra() eh a amostra de estrelas dentro do intervalo de massa
            massa_calculada[count] = massmin + (j * binmass) + (k * binmass/ amostra[j])
            # massa calculada eh a massa de cada estrela
            count += 1
                
    # 0-RA, 1-DEC, 2-mag1, 3-mag1err, 4-unc1, 5-mag2, 6-mag2err, 7-unc2
    star = np.zeros((N_stars_cmd, 8))
    
    for i in range(N_stars_cmd):
        for k in range(len(mass)-1): # abre as linhas do arquivo em massa
            # se a massa estiver no intervalo das linhas
            if (mass[k] < massa_calculada[i])&(mass[k+1] > massa_calculada[i]):
            # vai abrir tantas vezes quantas forem as estrelas representadas
                intervalo = (massa_calculada[i] - mass[k]) / (mass[k+1] - mass[k]) # intervalo entre zero e um
                star[i,2] = mag1[k] - (mag1[k] - mag1[k+1]) * intervalo
                star[i,5] = mag2[k] - (mag2[k] - mag2[k+1]) * intervalo

    # apply binarity
    # definition of binarity: fb = N_stars_in_binaries / N_total
    N_stars_bin = int(N_stars_cmd / ((2. / frac_bin) - 1))
    mag1_bin, mag2_bin = faker_bin(N_stars_bin, 'Kroupa', file_iso, dist)
        
    j = np.random.randint(N_stars_cmd, size=N_stars_bin)
    k = np.random.randint(N_stars_bin, size=N_stars_bin)

    for j, k in zip(j,k):
        star[j,2] = -2.5 * np.log10(10.**(-0.4 * star[j,2]) + 10.**(-0.4 * mag1_bin[k]))
        star[j,5] = -2.5 * np.log10(10.**(-0.4 * star[j,5]) + 10.**(-0.4 * mag2_bin[k]))                
        
    star[:,3] = apply_err(star[:,2], mag1_, err1_)
    star[:,6] = apply_err(star[:,5], mag1_, err2_)
    
    star[:,4] = unc(star[:,2], mag1_, err1_)
    star[:,7] = unc(star[:,5], mag1_, err2_)
        
    mag_ref_comp = 22.5
    comp_mag_ref = 1.
    comp_mag_max = 0.10
    dy_dx = (comp_mag_max-comp_mag_ref)/(mmax-mag_ref_comp)
    p_values = np.zeros(len(star[:,0]))
    cond = (star[:,2]+star[:,3] > mag_ref_comp)
    p_values[cond] = np.abs((comp_mag_ref - dy_dx * mag_ref_comp) + dy_dx * (star[:,2][cond] + star[:,3][cond]))
    p_values[star[:,2] > mmax] = 1.e-9
    p_values[star[:,2] < mag_ref_comp] = 1.

    aaaa = np.random.choice(len(star[:,0]), N_stars_cmd, replace=False, p=p_values/np.sum(p_values))

    r = expon.rvs(size=N_stars_cmd)
    r *= rexp
                
    rexp = (180/np.pi) * np.arctan(rexp/dist) # in deg
    r = (180/np.pi) * np.arctan(r/dist) # in deg
        
    phi = 2 * np.pi * np.random.rand(N_stars_cmd)
    
    X = r * np.sin(phi)
    Y = r * np.cos(phi)
    y_ell = Y
    x_ell = np.multiply(X, (1.-ell_))
    r_ell = np.sqrt(x_ell ** 2 + y_ell ** 2)
    phi_ell = np.arctan(-x_ell / y_ell)
    phi_ell[x_ell < 0.] += np.pi
    phi_ell += np.deg2rad(pa)
    star[:,0] = x0 + (r_ell * np.sin(phi_ell)) / np.cos(np.deg2rad(y0))
    star[:,1] = y0 + r_ell * np.cos(phi_ell)
    
    aaa = open(str(hpx) + '_clus.dat', 'w')

    for ii in aaaa: #range(len(star[:,0])):
        cor = star[ii,2]+star[ii,3]-(star[ii,5]+star[ii,6])
        mmag = star[ii,2]+star[ii,3]
        if (mmag < mmax)&(mmag > mmin)&(cor >= cmin)&(cor <= cmax):
            print(star[ii,0], star[ii,1], star[ii,2]+star[ii,3], star[ii,4], star[ii,5]+star[ii,6], \
                  star[ii,7], star[ii,3], star[ii,4], file=aaa)
    aaa.close()<jupyter_output><empty_output><jupyter_text>Calculating the properties of the sample of simulated clusters. Write the properties to file 'objects.dat'.<jupyter_code>cell_area = hp.nside2pixarea(nside_ini, degrees=True)

area = (dec_max - dec_min) * np.cos(np.deg2rad((ra_max + ra_min)/2.)) * (ra_max - ra_min)
    
vertices = hp.ang2vec([ra_min + border_extract, ra_max - border_extract, ra_max - border_extract, \
                       ra_min + border_extract], [dec_min + border_extract, dec_min + border_extract, \
                                                  dec_max - border_extract, dec_max - border_extract], lonlat=True)

hp_sample_un = hp.query_polygon(nside_ini, vertices, inclusive=False, nest=True, buff=None)

RA_pix, DEC_pix = hp.pix2ang(nside_ini, hp_sample_un, nest=True, lonlat=True)

c = SkyCoord(ra=RA_pix * u.degree, dec=DEC_pix * u.degree, frame='icrs')
L = c.galactic.l.degree
B = c.galactic.b.degree

obj_file = open('objects.dat', 'w')

mM = mM_min + np.random.rand(len(hp_sample_un)) * (mM_max - mM_min)
r_exp = 10**(log10_rexp_min * (log10_rexp_max / log10_rexp_min) ** np.random.rand(len(hp_sample_un)))
mass = 10**(log10_mass_min * (log10_mass_max / log10_mass_min) ** np.random.rand(len(hp_sample_un)))
dist = 10 ** ((mM/5) + 1)

ell = ell_min + np.random.rand(len(hp_sample_un)) * (ell_max - ell_min)
pa = pa_min + np.random.rand(len(hp_sample_un)) * (pa_max - pa_min)

for i in range(len(hp_sample_un)):
    print('{:d} {:.4f} {:.4f} {:.4f} {:.4f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}'.format(
        hp_sample_un[i], L[i], B[i], RA_pix[i], DEC_pix[i], r_exp[i], ell[i], pa[i], \
        mass[i], dist[i]), file=obj_file)

obj_file.close()<jupyter_output><empty_output><jupyter_text>Dist starsReading data from magnitude and errors.<jupyter_code>mag1_, err1_, err2_ = np.loadtxt('errors_Y6.dat', usecols=(0, 1, 2), unpack=True)
err1_ += 0.015; err2_ += 0.015<jupyter_output><empty_output><jupyter_text>Now simulating the clusters using 'faker' function.<jupyter_code>for jjj in range(len(hp_sample_un)):
    N_stars_cmd = int(mass[jjj] / 0.54)
    # os.register_at_fork(after_in_child=lambda: _get_font.cache_clear())
    faker(N_stars_cmd, frac_bin, 'Kroupa', RA_pix[jjj], DEC_pix[jjj], r_exp[jjj], ell[jjj], pa[jjj], dist[jjj], hp_sample_un[jjj])<jupyter_output><empty_output><jupyter_text>Now functions to join catalogs of simulated clusters and field stars, and to estimate signal-to-noise ratio.<jupyter_code>def JoinCat():
    """
    This is a function created only to join the catalog of real stars
    with random position and
    all the small catalogs of simulated clusters into a single one.
    There is a cut in the coord limits implemented below.
    Only global parameters are needed (see its description above)
    and the code is intended to write a single fits file.

    """
    global ra_min, ra_max, dec_min, dec_max, hp_sample_un, survey, RA, DEC, MAG_G, MAG_R, MAGERR_G, MAGERR_R, nside_ini
    
    GC = np.zeros(len(RA), dtype=int)

    for j in range(len(hp_sample_un)):
        try:
            RA_clus, DEC_clus, MAG1_clus, MAGERR1_clus, MAG2_clus, MAGERR2_clus = \
            np.loadtxt(str(hp_sample_un[j]) + '_clus.dat', usecols=(0,1,2,3,4,5), unpack=True)

            pr_limit = (RA_clus >= ra_min)&(RA_clus <= ra_max)&(DEC_clus >= dec_min)&(DEC_clus <= dec_max)&\
            (MAG1_clus <= mmax)&(MAG1_clus >= mmin)&(MAG1_clus - MAG2_clus >= cmin)&(MAG1_clus - MAG2_clus <= cmax)
        
            RA_clus, DEC_clus, MAG1_clus, MAG2_clus, MAGERR1_clus, MAGERR2_clus = \
            RA_clus[pr_limit], DEC_clus[pr_limit], MAG1_clus[pr_limit], MAG2_clus[pr_limit], \
            MAGERR1_clus[pr_limit], MAGERR2_clus[pr_limit]
          
            GC_clus = np.ones(len(RA_clus), dtype=int)
            GC = np.concatenate((GC, GC_clus), axis=0)
            RA = np.concatenate((RA, RA_clus), axis=0)
            DEC = np.concatenate((DEC, DEC_clus), axis=0)
            MAG_G = np.concatenate((MAG_G, MAG1_clus), axis=0)
            MAG_R = np.concatenate((MAG_R, MAG2_clus), axis=0)
            MAGERR_G = np.concatenate((MAGERR_G, MAGERR1_clus), axis=0)
            MAGERR_R = np.concatenate((MAGERR_R, MAGERR2_clus), axis=0)
        except:
            print('zero stars in ', hp_sample_un[j])
    
    HPX64 = hp.ang2pix(nside_ini, RA, DEC, nest=True, lonlat=True)
    col0 = fits.Column(name='GC', format='I', array=GC)
    col1 = fits.Column(name='ra', format='D', array=RA)
    col2 = fits.Column(name='dec', format='D', array=DEC)
    col3 = fits.Column(name='mag_g_with_err', format='E', array=MAG_G)
    col4 = fits.Column(name='mag_r_with_err', format='E', array=MAG_R)
    col5 = fits.Column(name='magerr_g', format='E', array=MAGERR_G)
    col6 = fits.Column(name='magerr_r', format='E', array=MAGERR_R)
    col7 = fits.Column(name='HPX64', format='K', array=HPX64)
    cols = fits.ColDefs([col0, col1, col2, col3, col4, col5, col6, col7])
    tbhdu = fits.BinTableHDU.from_columns(cols)
    tbhdu.writeto(survey + '_mockcat_for_detection.fits', overwrite=True)


def SNR_estimate(RA__, DEC__, G__, GR__, PIX_sim, nside1, mM_, inner_circle, rin_annulus, rout_annulus):
    """
    Estimate the SNR (Signal to Noise Ratio) of the simulated cluster in density-number.


    Parameters
    ----------
    RA__, DEC__ : lists
        The coordinates of the stars in the catalog
    G__, GR__ : lists
        Apparent magnitude in g band and g-r color of the stars.
    PIX_sim : int
        The pixel where the simulated cluster resided (nside=nside1)
        Nest=True.
    nside1 : int
        Nside of the pixelization of the sky reagarding the distribution of simulated clusters.
    mM_ : float
        Modulus distance of the simulated cluster
    inner_circle : float
        Radius of the inner circle (in deg) where the signal is estimated
    rin_annulus : float
        Inner radius where the noise is estimated (in deg)
    rout_annulus : float
        Outer radius where the noise is estimated (in deg)
    
    Returns
    -------
    float
        The SNR estimated following the SNR = N_inner / sqrt(N_outter_circle)
        where N_inner is the amount of stars in the inner_circle
        and N_outter_circle is the star counts in the outter circle
        normalized by the area of the inner_circle.
    """
    ra_center, dec_center = hp.pix2ang(nside1, PIX_sim, nest=True, lonlat=True)
    # loading data from isochronal mask
    gr_mask, g_mask, kind_mask = np.loadtxt('gr_g_model_D0.asc', usecols=(0, 1, 2), unpack=True)
    g_mask += mM_
    
    _g, _gerr = np.loadtxt('des_y6_g_gerr.asc', usecols=(0,1), unpack=True)
    _r, _rerr = np.loadtxt('des_y6_r_rerr.asc', usecols=(0,1), unpack=True)

    for i in range(len(gr_mask)):
        err_ = np.sqrt(np.interp(g_mask[i], _g, _gerr) ** 2. + np.interp(g_mask[i] - gr_mask[i], _r, _rerr) ** 2.)
        if err_ <= 1.00:
            if kind_mask[i] == 0.:
                gr_mask[i] -= err_
            elif kind_mask[i] == 1.:
                gr_mask[i] += err_

    points = np.column_stack([GR__, G__])
    verts = np.array([gr_mask, g_mask]).T
    path = mpath.Path(verts)
    t = path.contains_points(points)
    RA_ = list(compress(RA__, t))
    DEC_ = list(compress(DEC__, t))
    
    area_inner_circle = np.pi * (inner_circle ** 2)
    area_annulus = np.pi * ((rout_annulus ** 2) - (rin_annulus ** 2))
    r_star = np.sqrt(((RA_ - ra_center)*np.cos(np.deg2rad(dec_center)))**2 + (DEC_ - dec_center)**2)
    N_bg_equal_area = float(len(r_star[(r_star > rin_annulus)&(r_star < rout_annulus)]) * area_inner_circle / area_annulus)

    return len(r_star[r_star < inner_circle]) / np.sqrt(N_bg_equal_area)


def WriteSimClusFeatures():
    """
    Write a few features of the clusters in a file called 'N_stars,dat'.
    The columns of the file are:
    - the ipix of the cluster that serves as an ID (hp_sample_un[j]),
    - star counts (len(RA[cond_clus])),
    - absolute magnitude in V band (M_abs_V), and
    - signal-to-noise ratio (SNR).
    Only global parameters are needed (see its description in the first cells).
    The file must be written after the simulation because the absolute
    magnitude of the cluster and the number of stars are estimated after
    the simulation. The absolute magnitude depends strongly on the brightest
    stars and the star counts may vary in a few counts for instance
    for two clusters with the same mass (stars are a numerical realization
    within an IMF).
    """
    global hp_sample_un, nside_ini, survey, mM

    hdu = fits.open(survey + '_mockcat_for_detection.fits', memmap=True)
    GC = hdu[1].data.field('GC')
    RA = hdu[1].data.field('ra')
    DEC = hdu[1].data.field('dec')
    MAG_G = hdu[1].data.field('mag_g_with_err')
    MAG_R = hdu[1].data.field('mag_r_with_err')
    MAGERR_G = hdu[1].data.field('magerr_g')
    MAGERR_R = hdu[1].data.field('magerr_r')
    HPX64 = hdu[1].data.field('HPX64')
    
    ccc = open('N_stars.dat', 'w')
    
    for j in range(len(hp_sample_un)):
        # try:
        cond = (HPX64 == hp_sample_un[j])
        RA__, DEC__, MAGG__, MAGR__ = RA[cond], DEC[cond], MAG_G[cond], MAG_R[cond]
        # plt.scatter(RA__, DEC__)
        # plt.show()
        SNR = SNR_estimate(RA__, DEC__, MAGG__, MAGG__-MAGR__, hp_sample_un[j], nside_ini, mM[j], \
                               2./60, 10./60, 25./60)
            
        cond_clus = (cond)&(GC == 1)
        RA_clus, DEC_clus, MAGG_clus, MAGR_clus = RA[cond_clus], DEC[cond_clus], \
        MAG_G[cond_clus], MAG_R[cond_clus]
        flux_g = 10 ** (-0.4 * MAGG_clus)
        flux_r = 10 ** (-0.4 * MAGR_clus)
        M_abs_g = -2.5 * np.log10(np.sum(flux_g)) - mM[j]
        M_abs_r = -2.5 * np.log10(np.sum(flux_r)) - mM[j]
        M_abs_V = M_abs_g - 0.58*(M_abs_g-M_abs_r) - 0.01 # in V band following Jester 2005

        print('{:d} {:d} {:.2f} {:.2f}'.format(hp_sample_un[j], len(RA[cond_clus]), M_abs_V, SNR), file=ccc)
        #except:
        #    print(hp_sample_un[j], 0.000, 99.999, 0.000, file=ccc)
    
    ccc.close()
    
def SplitOutputHPX(file_in, out_dir):
    
    try:
        os.mkdir(out_dir)
    except:
        print('Folder already exists or is an invalid name.')

    data = getdata(file_in)
    GC = data['GC']
    ra = data['ra']
    dec = data['dec']
    mag_g_with_err = data['mag_g_with_err']
    mag_r_with_err = data['mag_r_with_err']
    magerr_g = data['magerr_g']
    magerr_r = data['magerr_r']
    HPX64 = data['HPX64']

    HPX_un = set(HPX64)

    for i in HPX_un:

        cond = (HPX64 == i)

        col0 = fits.Column(name='GC', format='I', array=GC[cond])
        col1 = fits.Column(name='ra', format='D', array=ra[cond])
        col2 = fits.Column(name='dec', format='D', array=dec[cond])
        col3 = fits.Column(name='mag_g_with_err', format='E', array=mag_g_with_err[cond])
        col4 = fits.Column(name='mag_r_with_err', format='E', array=mag_r_with_err[cond])
        col5 = fits.Column(name='magerr_g', format='E', array=magerr_g[cond])
        col6 = fits.Column(name='magerr_r', format='E', array=magerr_r[cond])
        col7 = fits.Column(name='HPX64', format='K', array=np.repeat(int(i), len(GC[cond])))
        cols = fits.ColDefs([col0, col1, col2, col3, col4, col5, col6, col7])
        tbhdu = fits.BinTableHDU.from_columns(cols)
        tbhdu.writeto(out_dir + '/' + str(int(i)) + '.fits', overwrite=True)

def SplitFtpHPX(file_in, out_dir, nside_in=4096, nest_in=True, nside_out=64, nest_out=True):
    
    try:
        os.mkdir(out_dir)
    except:
        print('Folder already exists or is an invalid name.')

    data = getdata(file_in)
    HPX4096 = data['HP_PIXEL_NEST_4096']
    SIGNAL = data['SIGNAL']

    ra, dec = hp.pix2ang(nside_in, HPX4096, nest=nest_in, lonlat=True)

    HPX64 = hp.ang2pix(nside_out, ra, dec, nest=nest_out, lonlat=True)

    HPX_un = set(HPX64)

    for i in HPX_un:

        cond = (HPX64 == i)

        col0 = fits.Column(name='HP_PIXEL_NEST_4096', format='K', array=HPX4096[cond])
        col1 = fits.Column(name='SIGNAL', format='E', array=SIGNAL[cond])
        cols = fits.ColDefs([col0, col1])
        tbhdu = fits.BinTableHDU.from_columns(cols)
        tbhdu.writeto(out_dir + '/' + str(int(i)) + '.fits', overwrite=True)<jupyter_output><empty_output><jupyter_text>Now applying the two functions.<jupyter_code>JoinCat()
SplitOutputHPX(survey + '_mockcat_for_detection.fits', 'HPX_cats')
# SplitFtpHPX('ftp_4096_nest.fits', 'HPX_ftp')
WriteSimClusFeatures()<jupyter_output>Folder already exists or is an invalid name.<jupyter_text>Merge both files in a single file using join command.<jupyter_code>os.system('join --nocheck-order N_stars.dat objects.dat > star_clusters_simulated.dat')<jupyter_output><empty_output><jupyter_text>A few plots to characterize the simulated clusters.<jupyter_code>PIX_sim, NSTARS, MAG_ABS_V, RA, DEC, R_EXP, ELL, PA, MASS, DIST = np.loadtxt('star_clusters_simulated.dat', \
                                                              usecols=(0, 1, 2, 6, 7, 8, 9, 10, 11, 12), \
                                                              unpack=True)
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5), sharey=True)
ax1.scatter(np.log10(1.7*R_EXP[MAG_ABS_V < 0.]), MAG_ABS_V[MAG_ABS_V < 0.])
ax1.set_ylabel('M(V)')
ax1.set_xlabel('log10(h-l radii(pc))')

ax2.scatter(MASS, MAG_ABS_V)
ax2.set_xlabel('mass(Msun)')
plt.show()
def plot_ftp():
    """Plot footprint map to check area.
    """
    nside = 4096
    npix = hp.nside2npix(nside)
    
    data = getdata('ftp_4096_nest.fits')
    pix_ftp = data['HP_PIXEL_NEST_4096']
    
    ra_pix_ftp, dec_pix_ftp = hp.pix2ang(nside, pix_ftp, nest=True, lonlat=True)
    map_ftp = np.zeros(hp.nside2npix(nside))
    map_ftp[pix_ftp] = 1

    test = hp.cartview(map_ftp, nest=True, lonra=[np.min(ra_pix_ftp), np.max(ra_pix_ftp)], \
               latra=[np.min(dec_pix_ftp), np.max(dec_pix_ftp)], hold=True, cbar=False, \
               title='', return_projected_map=True)
    # plt.clf()
    
    PIX_sim, NSTARS, MAG_ABS_V, RA, DEC, R_EXP, ELL, PA, MASS, DIST = np.loadtxt('star_clusters_simulated.dat', \
                                                              usecols=(0, 1, 2, 6, 7, 8, 9, 10, 11, 12), \
                                                              unpack=True)
    for i in range(len(RA)):
        hp.projtext(RA[i], DEC[i], str(PIX_sim[i]), lonlat=True, fontsize=10, c='k',\
                    horizontalalignment='center')
        hp.projscatter(RA[i], DEC[i], lonlat=True, coord='C', s=1., color='k', lw=0.1)

    data = getdata(survey + '_mockcat_for_detection.fits')
    RA_star, DEC_star = data['RA'], data['DEC']
    fig, axs = plt.subplots(1, 1, figsize=(10,10))
    axs.imshow(test, origin='lower', extent=(ra_max, ra_min, dec_min, dec_max), interpolation = 'none')
    axs.scatter(RA, DEC, s=20, c='k', marker='s', label='Simulated clusters')
    axs.scatter(RA_star, DEC_star, s=0.1, c='k', marker='o', label='Simulated stars')
    axs.set_xlim([ra_min, ra_max])
    axs.set_ylim([dec_min, dec_max])
    axs.set_xlabel('RA (deg)')
    axs.set_ylabel('DEC (deg)')
    axs.set_title('Distribution of stars on Footprint Map')
    axs.grid()
    plt.legend()
    plt.show()
    print(len(PIX_sim))
    
plot_ftp()
def plots_ang_size():
    """Plots to analyze the simulated clusters.
    """
    
    global mmin, mmax, cmin, cmax
    
    cmap = mpl.cm.get_cmap("inferno")
    cmap.set_under('dimgray')
    cmap.set_bad('black')
    
    hp_sample_un, NSTARS, MAG_ABS_V, RA_pix, DEC_pix, r_exp, ell, pa, mass, \
    dist = np.loadtxt('star_clusters_simulated.dat',
                      usecols=(0, 1, 2, 6, 7, 8, 9, 10,11,12), 
                      unpack=True)
    
    for ii in hp_sample_un:
        star = np.loadtxt(str(int(ii)) + '_clus.dat')
        '''
        plt.scatter(star[:,2]-star[:,4], star[:,2], color='b')
        plt.title('HPX ' + str(int(ii)) + ', N=' + str(len(star[:,2])))
        plt.ylim([mmax, mmin])
        plt.xlim([cmin, cmax])
        plt.xlabel('mag1-mag2')
        plt.ylabel('mag1')
        plt.savefig(str(int(ii)) + '_cmd.png')
        plt.close()
        '''
        h1, xedges, yedges, im1 = plt.hist2d(star[:,2]-star[:,4], star[:,2], bins=50,
                                             range=[[cmin, cmax], [mmin, mmax]],
                                             # norm=mpl.colors.LogNorm(),
                                             cmap=cmap)
        plt.clf()
        plt.title('HPX ' + str(int(ii)) + ', N=' + str(len(star[:,2])))
        im1 = plt.imshow(h1.T, interpolation='None', origin='lower', vmin=0.1, vmax=np.max(h1),
                         extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]], aspect='auto', cmap=cmap)
        plt.ylim([mmax, mmin])
        plt.xlim([cmin, cmax])
        plt.xlabel('mag1-mag2')
        plt.ylabel('mag1')
        plt.colorbar(im1, cmap=cmap, orientation='vertical', label='stars per bin') 
        plt.savefig(str(int(ii)) + '_cmd.png')
        plt.close()


    ang_size = 60 * np.rad2deg(np.arctan(1.7 * r_exp / dist))
    plt.hist(dist / 1000, bins=50)
    plt.xlabel('Distance (kpc)')
    plt.ylabel('N objects')
    plt.show()
    plt.hist(ang_size, bins=50)
    plt.xlabel('Half-light radii (arcmin)')
    plt.ylabel('N objects')
    plt.show()
    plt.scatter(dist/1000, ang_size)
    plt.xlabel('Distance (kpc)')
    plt.ylabel('Half-light radii (arcmin)')
    plt.show()
    
    plt.scatter(mass, NSTARS)
    plt.xlabel('MASS(MSun)')
    plt.ylabel('N stars')
    plt.show()
    
    plt.scatter(mass, MAG_ABS_V)
    plt.xlabel('MASS(MSun)')
    plt.ylabel('MAG_ABS_V')
    plt.show()
    
plots_ang_size()
def radec2GCdist(ra, dec, dist_kpc):
    """
    Return Galactocentric distance from ra, dec, D_sun_kpc.

    Parameters
    ----------
    ra, dec : float or list
        Coordinates of the objects (in deg)
    dist_kpc : float or list
        Distance in kpc of the objects

    Returns
    -------
    float of list
        the Galactocentric distance to the object[s]
    """
    c1 = coord.SkyCoord(ra = ra*u.degree,
                        dec = dec*u.degree,
                        distance = dist_kpc*u.kpc,
                        frame='icrs')
    x, y, z = c1.transform_to(coord.Galactocentric).x.value, \
              c1.transform_to(coord.Galactocentric).y.value, \
              c1.transform_to(coord.Galactocentric).z.value

    return np.sqrt(x*x + y*y + z*z)


def plots_ref():
    """Make a few plots about the simulated clusters

    """
    
    ra_DG, dec_DG, dist_kpc_obj, Mv_obj, rhl_pc_obj, FeH_DG = np.loadtxt('catalogos/objects_in_ref.dat',
                                                                           usecols=(0, 1, 4, 8, 10, 11), unpack=True)
    name_obj = np.loadtxt('catalogos/objects_in_ref.dat', dtype=str, usecols=(2), unpack=True)
    # 0-Name 1-L 2-B 3-R_gc	4-Fe/H 5-M-M 6-Mv 7-rhl arcmin
    R_MW_GC, FeH_GC, mM_GC, Mv_GC, rhl_arcmin_GC = np.loadtxt('catalogos/Harris_updated.dat',
                                                           usecols=(3, 4, 5, 6, 7), unpack=True)
    dist_kpc_GC = 10 ** (mM_GC/5 - 2)
    PIX_sim, NSTARS, MAG_ABS_V, RA, DEC, R_EXP, ELL, PA, MASS, DIST = np.loadtxt('star_clusters_simulated.dat', \
                                                                      usecols=(0, 1, 2, 6, 7, 8, 9, 10, 11, 12), \
                                                                      unpack=True)
    LOG10_RHL_PC_SIM = np.log10(1.7 * R_EXP)
    rhl_pc_GC = 1000 * dist_kpc_GC * (rhl_arcmin_GC/(57.3*60))
    
    MW_center_distance_DG_kpc = radec2GCdist(ra_DG, dec_DG, dist_kpc_obj)
    
    fig, axs = plt.subplots(2, 2, figsize=(15, 8))
    axs[0, 0].hist(MAG_ABS_V, bins=20, range=(-16, 0.), histtype='stepfilled', label='Sim', color='grey', ls='--')
    axs[0, 0].hist(Mv_obj, bins=20, range=(-16, 0.), histtype='step', label='DG', color='r')
    axs[0, 0].hist(Mv_GC, bins=20, range=(-16, 0.), histtype='step', label='GC', color='b')
    axs[0, 0].set_xlabel(r'$M_V$')
    axs[0, 0].set_ylabel('N')
    axs[0, 0].legend(loc=2)
    
    axs[0, 1].hist(LOG10_RHL_PC_SIM, bins=20, histtype='stepfilled', range=(0, 4.0), label='Sim', color='grey', ls='--')
    axs[0, 1].hist(np.log10(rhl_pc_obj), bins=20, histtype='step', range=(0, 4.0), label='DG', color='r')
    axs[0, 1].hist(np.log10(rhl_pc_GC), bins=20, histtype='step', range=(0, 4.0), label='GC', color='b')
    axs[0, 1].set_xlabel('log10(rhl[pc])')
    axs[0, 1].legend(loc=1)
    
    axs[1, 0].hist(DIST/1000, bins=20, range=(0, 400.), histtype='stepfilled', label='Sim', color='grey', ls='--')
    axs[1, 0].hist(dist_kpc_obj, bins=20, range=(0, 400.), histtype='step', label='DG', color='r')
    axs[1, 0].hist(dist_kpc_GC, bins=20, range=(0, 400.), histtype='step', label='GC', color='b')
    axs[1, 0].set_xlabel('Distance (kpc)')
    axs[1, 0].legend(loc=1)

    axs[1, 1].hist(np.repeat(FeH_iso, len(MAG_ABS_V)), bins=20, range=(-3, 1.), histtype='stepfilled',
                   label='Sim', color='grey', ls='--')
    axs[1, 1].hist(FeH_GC, bins=20, range=(-3, 1.), histtype='step', label='GC', color='r')
    axs[1, 1].hist(FeH_DG, bins=20, range=(-3, 1.), histtype='step', label='DG', color='b')
    axs[1, 1].set_xlabel('[Fe/H]')
    axs[1, 1].legend(loc=1)

    plt.suptitle('Physical features of 58 Dwarf Gal + 152 GC + ' + str(len(PIX_sim)) + ' Simulations', fontsize=16)
    fig.tight_layout()
    plt.subplots_adjust(top=0.92)
    plt.savefig('_01_real_objects.png')
    plt.close()
    
    rhl = np.logspace(np.log10(1.8), np.log10(1800), 10, endpoint=True)
    m_v = np.linspace(1, -14, 10, endpoint=True)
    
    plt.scatter(1.7*R_EXP, MAG_ABS_V, marker='s', color='grey', label='Sim')
    plt.scatter(rhl_pc_obj, Mv_obj, marker='^', color='r', label='DG')
    plt.scatter(rhl_pc_GC, Mv_GC, marker='x', color='b', label='GC')
    plt.plot(np.logspace(np.log10(1.8), np.log10(1800), 10, endpoint=True), \
             np.linspace(1, -14, 10, endpoint=True), color='b', ls=':')
    plt.plot(np.logspace(np.log10(4.2), np.log10(4200), 10, endpoint=True), \
             np.linspace(1, -14, 10, endpoint=True), color='b', ls=':')
    plt.plot(np.logspace(np.log10(11), np.log10(11000), 10, endpoint=True), \
             np.linspace(1, -14, 10, endpoint=True), color='b', ls=':')
    plt.plot(np.logspace(np.log10(28), np.log10(28000), 10, endpoint=True), \
             np.linspace(1, -14, 10, endpoint=True), color='b', ls=':')
    plt.text(300, -7.9, r'$\mu_V=27\ mag/arcsec$', rotation=40)
    plt.text(400, -4.2, r'$\mu_V=31\ mag/arcsec$', rotation=40)
    plt.xscale('log')
    plt.xlim([0.4,4000])
    plt.ylim([1, -14])
    plt.ylabel(r'$M_V$')
    plt.xlabel(r'$r_h\ (pc)$')
    plt.legend(loc=2, frameon=True)
    plt.savefig('_03_mv_rh.png')
    plt.close()
    
    plt.scatter(DIST/1000, np.repeat(FeH_iso, len(DIST)), label='Sim', color='grey', marker='x', lw=1.)
    plt.scatter(MW_center_distance_DG_kpc, FeH_DG, label='DG', color='r')
    plt.scatter(R_MW_GC, FeH_GC, label='GC', color='b')
    plt.xlabel('Distance to the Galactic center (kpc)')
    plt.ylabel('[Fe/H]')
    plt.ylim([-3.5, 0])
    plt.legend()
    plt.grid()
    plt.savefig('_02_feh_rgc.png')
    plt.close()
    
plots_ref()<jupyter_output><empty_output><jupyter_text>Copying the main files to the folder where the Gawa code will be able to detect the simulated clusters.<jupyter_code># os.system('cp ftp_4096_nest.fits ../ga-wazpy_detect_sim/raw_ftp/')
# os.system('cp des*_detection.fits ../ga-wazpy_detect_sim/raw_cats_corr/')
# os.system('cp star_clusters_simulated.dat ../ga-wazpy_detect_sim/')
# os.system('cp _*.png ../ga-wazpy_detect_sim/')
# os.system('cp *_cmd.png ../ga-wazpy_detect_sim/')
# os.system('ls *_cmd.png > list_cmds.dat')
# os.system('python format_html.py > sim_main.html')
# os.system('cp sim_main.html ../ga-wazpy_detect_sim/')
def plot_err():
    """Plot the magnitude and error of the simulated clusters compared to the
    real stars, in log scale.

    """
    hdu = fits.open('des_mockcat_for_detection.fits', memmap=True)
    GC = hdu[1].data.field('GC')
    ra = hdu[1].data.field('ra')
    dec = hdu[1].data.field('dec')
    mag_g_with_err = hdu[1].data.field('mag_g_with_err')
    mag_r_with_err = hdu[1].data.field('mag_r_with_err')
    magerr_g = hdu[1].data.field('magerr_g')
    magerr_r = hdu[1].data.field('magerr_r')
    hdu.close()
    
    plt.scatter(mag_g_with_err[GC == 0], magerr_g[GC == 0], label='Field stars', c='k')
    plt.scatter(mag_g_with_err[GC == 1], magerr_g[GC == 1], label='Simulated stars', c='r', zorder=10)
    plt.yscale('log')
    plt.xlabel('mag_g_with_err')
    plt.ylabel('magerr_g')
    plt.legend()
    plt.show()

plot_err()
def RemoveCloseStars(name_cat, output_cat, PSF_factor):
    """ This function removes the stars closer than PSF_factor * PSF_size
    This is an observational bias of the DES since the photometric pipeline
    is set to join regions closer than an specific distance.
    In this sense, if many objects are crowded in a small area,
    the pipeline joins the detections into a single object.
    
    This function may be improved to run in parallel using many ipixels
    since it takes too long to check and remove stars.

    Parameters
    ----------
    name_cat : str
        The file with the information of the objects
    output_cat : str
        The output file to be written with the stars that survived the test.
    PSF_factor : float
        A factor to multiply the PSF size. The result is the minimal distance
        between two objects to be detected as two objects. If the distance
        is less than PSF_factor * PSF_size, none of the objects survives.

    """
    global nside_ini

    PSF_size = 0.8 # arcsec
    hdu = fits.open(name_cat, memmap=True)
    GC = hdu[1].data.field('GC')
    ra = hdu[1].data.field('ra')
    dec = hdu[1].data.field('dec')
    mag_g_with_err = hdu[1].data.field('mag_g_with_err')
    mag_r_with_err = hdu[1].data.field('mag_r_with_err')
    magerr_g = hdu[1].data.field('magerr_g')
    magerr_r = hdu[1].data.field('magerr_r')
    HPX64 = hdu[1].data.field('HPX64')
    hdu.close()

    idx = []
    seplim=(1./3600) * u.degree
    for i in range(len(ra)):
        cond = (ra > ra[i] - 0.05)&(ra < ra[i] + 0.05)&(dec > dec[i] - 0.05)&(dec < dec[i] + 0.05)
        
        c = SkyCoord(ra = ra[cond] * u.degree, dec = dec[cond] * u.degree)
        idx_, sep2d, dist3d = match_coordinates_sky(c, c, nthneighbor=2, storekdtree='kdtree_sky')
        cond2 = [sep2d < seplim]
        if len(idx_[cond2]) == 0:
            idx.append(i)
        #dists = dist_ang(ra[cond], dec[cond], ra[i], dec[i])
        #dist = np.sort(list(dists))[1]
        #if dist > (PSF_factor * PSF_size/3600.):
        #    idx.append(i)
            
    print(len(ra), len(idx))
    HPX64 = hp.ang2pix(nside_ini, ra, dec, nest=True, lonlat=True)
    col0 = fits.Column(name='GC', format='I', array=np.asarray([GC[i] for i in idx]))
    col1 = fits.Column(name='ra', format='D', array=np.asarray([ra[i] for i in idx]))
    col2 = fits.Column(name='dec', format='D', array=np.asarray([dec[i] for i in idx]))
    col3 = fits.Column(name='mag_g_with_err', format='E', array=np.asarray([mag_g_with_err[i] for i in idx]))
    col4 = fits.Column(name='mag_r_with_err', format='E', array=np.asarray([mag_r_with_err[i] for i in idx]))
    col5 = fits.Column(name='magerr_g', format='E', array=np.asarray([magerr_g[i] for i in idx]))
    col6 = fits.Column(name='magerr_r', format='E', array=np.asarray([magerr_r[i] for i in idx]))
    col7 = fits.Column(name='HPX64', format='K', array=np.asarray([HPX64[i] for i in idx]))
    cols = fits.ColDefs([col0, col1, col2, col3, col4, col5, col6, col7])
    tbhdu = fits.BinTableHDU.from_columns(cols)
    tbhdu.writeto(output_cat, overwrite=True)

# RemoveCloseStars('des_mockcat_for_detection.fits', 'des_mockcat_for_detection_wo_closest.fits', 1.0)<jupyter_output>/opt/conda/lib/python3.6/site-packages/ipykernel_launcher.py:45: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
{'fraction_non_alphanumeric': 0.14038922850731253, 'fraction_numerical': 0.025709974464133714, 'mean_word_length': 3.0131200993711666, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 92, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '625613', 'license': 'MIT', 'n_tokens_mistral': 20146, 'n_tokens_neox': 18855, 'n_words': 4709, 'path': 'gawa_simulations.ipynb', 'repo_name': 'linea-it/ga_sim', 'size': '51714'}