Page MenuHomec4science
Files F88375846

No OneTemporary
Actions

File Metadata

Created
Fri, Oct 18, 11:40

View Options

diff --git a/ParticlesToSurfaceBrightness/ARRAKIHS.py b/ParticlesToSurfaceBrightness/ARRAKIHS.py
new file mode 100644
index 0000000..026b907
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/ARRAKIHS.py
@@ -0,0 +1,621 @@
+#!/usr/bin/env python3
+###########################################################################################
+# package: pNbody
+# file: __init__.py
+# brief: init file
+# copyright: GPLv3
+# Copyright (C) 2019 EPFL (Ecole Polytechnique Federale de Lausanne)
+# LASTRO - Laboratory of Astrophysics of EPFL
+# author: Yves Revaz <yves.revaz@epfl.ch>
+#
+# This file is part of pNbody.
+###########################################################################################
+
+
+import numpy as np
+
+from astropy import constants as cte
+from astropy import units as u
+
+from pNbody import *
+from pNbody import ic
+from pNbody import mapping
+from pNbody import geometry
+
+import argparse
+
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+
+from scipy import signal
+
+import stars_class
+
+
+def SetTelecope(opt):
+
+
+ ##################
+ # set telescope
+
+ if opt.telescope=='arrakihs_vis':
+
+ opt.ccd_shape=[4096,4096]
+ opt.ccd_size=[1.5,1.5]
+
+ if opt.filter is None:
+ opt.filter = "VIS"
+
+ elif opt.telescope=='arrakihs_nir':
+
+ opt.ccd_shape=[2048,2048]
+ opt.ccd_size=[1.3,1.3]
+
+ if opt.filter is None:
+ opt.filter = "UKIRT_J"
+
+
+ elif opt.telescope=='euclid_vis':
+
+ opt.ccd_shape=[24576,24576]
+ opt.ccd_size=[0.787,0.787]
+
+
+ ##################
+ # set filters
+
+ if opt.filter=='F475X':
+ opt.mlim = 30.12
+ opt.mlmin = 29.50
+ opt.mlmax = opt.mlim
+
+ elif opt.filter=='VISeuclid':
+ opt.mlim = 33 #30.23
+ opt.mlmin = 29.5 #29.50 # Euclid preparation: XVI. Borlaff 2021
+ opt.mlmax = opt.mlim
+
+ elif opt.filter=='Yeuclid':
+ opt.mlim = 29.59
+ opt.mlmin = 29.50
+ opt.mlmax = opt.mlim
+
+ elif opt.filter=='Jeuclid':
+ opt.mlim = 31.37
+ opt.mlmin = 29.50
+ opt.mlmax = opt.mlim
+
+
+
+
+
+
+ ##################################
+ # ccd size
+ ##################################
+
+ # degree to arcsec
+ opt.ccd_size = [opt.ccd_size[0]*3600,opt.ccd_size[1]*3600]
+ opt.ccd_pixel_size = [opt.ccd_size[0]/opt.ccd_shape[0],opt.ccd_size[1]/opt.ccd_shape[1]]
+ opt.ccd_pixel_area = opt.ccd_pixel_size[0]*opt.ccd_pixel_size[1]
+
+
+
+
+def GenerateDwarf(opt):
+ '''
+ Generate a dwarf model based on a plummer sphere.
+ Note that the plummer scaling length is equivalent to
+ the projected half light radius.
+ '''
+
+ # set maximal radius
+ rmax = opt.e * opt.rmax_e_ratio
+
+ # convert to int
+ opt.N = int(opt.N)
+
+ # generate the model
+ nb = ic.plummer(opt.N,1,1,1,opt.e,rmax,M=opt.M,irand=1,vel='no',name=opt.outputfilename,ftype='swift')
+
+ if opt.outputfilename is not None:
+ nb.write()
+
+ return nb
+
+
+
+def gaussian_filter(kernel_size, sigma=1, muu=0):
+
+ # Initializing value of x,y as grid of kernel size
+ # in the range of kernel size
+
+ x, y = np.meshgrid(np.linspace(-kernel_size, kernel_size, kernel_size),np.linspace(-kernel_size, kernel_size, kernel_size))
+ dst = np.sqrt(x**2+y**2)
+
+
+ # lower normal part of gaussian
+ #normal = 1/np.sqrt(2 * np.pi * sigma**2)
+
+ # Calculating Gaussian filter
+ gauss = np.exp(-((dst-muu)**2 / (2.0 * sigma**2)))
+ gauss = gauss/gauss.sum()
+
+ return gauss
+
+
+
+def Project(opt):
+ '''
+ Project a stellar model
+ '''
+
+
+ ##################################
+ # factor conversion kpc to arcsec
+ ##################################
+
+ # Mpc to kpc
+ distance = opt.distance * 1000
+
+ # define a conversion function : kpc -> arcsec
+ fct_kpc2arcsec = np.vectorize(lambda x: 3600*180/np.pi* np.arctan(x/distance) )
+ opt.fct_kpc2arcsec = lambda x: 3600*180/np.pi* np.arctan(x/distance)
+
+ # define a conversion function : arcsec -> kpc
+ fct_arcsec2kpc = np.vectorize(lambda x: distance*np.tan(x*np.pi/3600/180) )
+ opt.fct_arcsec2kpc = lambda x: distance*np.tan(x*np.pi/3600/180)
+
+
+ # image properties
+ if opt.fov is not None:
+ xmin = -opt.fov/2.
+ xmax = +opt.fov/2.
+ ymin = -opt.fov/2.
+ ymax = +opt.fov/2.
+
+ nx = int(opt.fov/opt.ccd_pixel_size[0])+1
+ ny = int(opt.fov/opt.ccd_pixel_size[1])+1
+
+
+ elif opt.size is not None:
+ nx = opt.size[0] + 1 # +1 as the resulting matrix will be -1
+ ny = opt.size[1] + 1 # +1 as the resulting matrix will be -1
+
+ xmin = -nx*opt.ccd_pixel_size[0] /2.
+ xmax = +nx*opt.ccd_pixel_size[0] /2.
+ ymin = -ny*opt.ccd_pixel_size[1] /2.
+ ymax = +ny*opt.ccd_pixel_size[1] /2.
+
+ else:
+ nx = int(opt.ccd_shape[0] *opt.ccd_field_fraction)
+ ny = int(opt.ccd_shape[1] *opt.ccd_field_fraction)
+
+ xmin = -opt.ccd_size[0]/2 *opt.ccd_field_fraction
+ xmax = +opt.ccd_size[0]/2 *opt.ccd_field_fraction
+ ymin = -opt.ccd_size[1]/2 *opt.ccd_field_fraction
+ ymax = +opt.ccd_size[1]/2 *opt.ccd_field_fraction
+
+
+
+ # extension in kpc
+ opt.xmin_kpc = fct_arcsec2kpc(xmin)
+ opt.xmax_kpc = fct_arcsec2kpc(xmax)
+ opt.ymin_kpc = fct_arcsec2kpc(ymin)
+ opt.ymax_kpc = fct_arcsec2kpc(ymax)
+
+
+
+ ##################################
+ # open the model
+ ##################################
+
+ if opt.file is None:
+ nb = opt.nb
+ else:
+ nb = Nbody(opt.file)
+ nb.mass = nb.mass*opt.toMsol
+
+ # remove doublets
+ u,idx = np.unique(nb.rxyz(),return_index=True)
+ nb = nb.selectp(lst=nb.num[idx])
+
+ if not nb.has_array("rsp"):
+ # compute Hsml
+ #nb.set_tpe(0)
+ #nb.InitSphParameters(DesNumNgb=32, MaxNumNgbDeviation=2)
+ #nb.getTree()
+ #nb.rsp = nb.get_rsp_approximation()
+ nb.ComputeRsp(5)
+
+ if not nb.has_array("age"):
+ # compute Age [in Gyr]
+ print("Compute ages...")
+ nb.age = nb.StellarAge(units="Gyr")
+ print("done.")
+
+ # rotate
+ if opt.los is not None:
+ nb.align(opt.los)
+
+
+ # compute mass
+ mass = nb.mass
+
+
+
+
+ # scale the model : kpc -> arcsec
+ x = fct_kpc2arcsec(nb.x())
+ y = fct_kpc2arcsec(nb.y())
+ z = fct_kpc2arcsec(nb.z())
+
+ nb.pos = np.transpose((x,y,z))
+ nb.pos = nb.pos.astype(np.float32)
+
+
+ binsx = np.linspace(xmin,xmax,nx)
+ binsy = np.linspace(ymin,ymax,ny)
+
+
+
+
+ # get the luminosity:
+ if opt.filter=="V" or opt.filter=="VIS" or opt.filter=="UKIRT_J" :
+
+ # here we assume the mass of the particle to contain its luminosity
+
+ # luminosity (L_sun)
+ #L = nb.mass.sum()
+
+ # map from numpy
+ #L_map,xe,ye = np.histogram2d(x, y,bins=[binsx,binsy],weights=nb.mass)
+ #L_map = np.rot90(L_map)
+
+ # map from pnbody
+ params = {}
+ params['size'] = (xmax,ymax)
+ params['shape'] = (nx,ny)
+ params['rendering'] = 'map'
+ params['obs'] = None
+ params['xp'] = None
+ params['view'] = opt.view
+ params['mode'] = 'm'
+ params['exec'] = None
+ params['macro'] = None
+ params['frsp'] = opt.frsp
+ params['filter_name'] = None
+
+ L_map = nb.CombiMap(params)
+ L_map = np.rot90(L_map)
+ L_map = np.flipud(L_map)
+
+ # psf convolution
+ if opt.convolve:
+ psf = opt.psf/opt.ccd_pixel_size[0]
+ # gaussian filter
+ psf_map = gaussian_filter(L_map.shape[0],psf)
+ # convolve
+ L_map = signal.convolve2d(L_map, psf_map, mode='same', boundary='fill', fillvalue=0)
+
+
+ # compute the absolute magnitude in each pixel (vband)
+ Mvega = 4.81
+ M_map = Mvega - 2.5*np.log10(L_map)
+
+ # filter stuffs
+ M_V = M_map
+ L_V = L_map
+
+ M_R = M_V -0.45
+ M_I = M_R -0.5
+ M_J = M_I -0.6
+ L_R = np.log10( (M_R-4.43)/(-2.5) )
+ L_I = np.log10( (M_R-4.10)/(-2.5) )
+ L_J = np.log10( (M_J-3.65)/(-2.5) )
+
+ L_VIS = L_V + L_R + L_I
+ M_VIS = 3.216 - 2.5*np.log10(L_VIS)
+
+ if opt.filter=="V":
+ M_map = M_V
+ elif opt.filter=="VIS":
+ M_map = M_VIS
+ elif opt.filter=="UKIRT_J":
+ M_map = M_J
+
+
+ # compute the apparent magnitude in each pixel
+ d = distance*100 # kpc -> 10pc
+ m_map = M_map + 5*np.log10(d)
+
+ # compute the surface brightness in each pixel
+ S_map = m_map + 2.5*np.log10(opt.ccd_pixel_area)
+
+
+ #S_map = L_map
+
+
+
+ elif opt.filter=="F475X" or opt.filter=="VISeuclid" or opt.filter=="Yeuclid" or opt.filter=="Jeuclid":
+
+ if not nb.has_array("MagF475X") or not nb.has_array("MagVISeuclid") or not nb.has_array("MagYeuclid") or not nb.has_array("MagJeuclid"):
+
+
+ # compute a magnitude for the filter in each mass bin
+ print("Compute magnitudes...")
+
+ if opt.filter=="F475X":
+ M = stars_class.HST475X_fun(None,nb.age,nb.MH())
+
+ elif opt.filter=="VISeuclid":
+ M = stars_class.VISeuclid_fun(None,nb.age,nb.MH())
+
+ elif opt.filter=="Yeuclid":
+ M = stars_class.Yeuclid_fun(None,nb.age,nb.MH())
+
+ elif opt.filter=="Jeuclid":
+ M = stars_class.Jeuclid_fun(None,nb.age,nb.MH())
+
+ # convert to flux (ignore the zero point)
+ # to get the correct flux, see notes in Garrotxa
+ F = 10**(-M/2.5)
+
+ # get the number of stars in each mass bin
+ Nstars = stars_class.Stars_fun(mass,None,None, 'normed_3slope')
+
+ # sum the contribution of the mass bins
+ F = np.sum(F*Nstars, axis=0)
+
+ # compute the absolute magnitude in each pixel (as before we ignore the zero point)
+ M = - 2.5*np.log10(F)
+
+ print("done.")
+
+ else:
+
+ if opt.filter=="F475X":
+ M = nb.MagF475X
+
+ elif opt.filter=="VISeuclid":
+ M = nb.MagVISeuclid
+
+ elif opt.filter=="Yeuclid":
+ M = nb.MagYeuclid
+
+ elif opt.filter=="Jeuclid":
+ M = nb.MagJeuclid
+
+
+
+
+ # convert to flux (ignore the zero point)
+ F = 10**(-M/2.5)
+
+ # store in the mass field
+ nb.mass = F.astype(np.float32)
+
+ # map from numpy
+ #L_map,xe,ye = np.histogram2d(x, y,bins=[binsx,binsy],weights=nb.mass)
+ #L_map = np.rot90(L_map)
+
+ # map from pnbody
+ params = {}
+ params['size'] = (xmax,ymax)
+ params['shape'] = (nx,ny)
+ params['rendering'] = 'map'
+ params['obs'] = None
+ params['xp'] = None
+ params['view'] = opt.view
+ params['mode'] = 'm'
+ params['exec'] = None
+ params['macro'] = None
+ params['frsp'] = opt.frsp
+ params['filter_name'] = None
+
+ L_map = nb.CombiMap(params)
+ L_map = np.rot90(L_map)
+ L_map = np.flipud(L_map)
+
+ # psf convolution
+ if opt.convolve:
+ psf = opt.psf/opt.ccd_pixel_size[0]
+ # gaussian filter
+ psf_map = gaussian_filter(L_map.shape[0],psf)
+ # convolve
+ L_map = signal.convolve2d(L_map, psf_map, mode='same', boundary='fill', fillvalue=0)
+
+
+ # compute the absolute magnitude in each pixel (as before we ignore the zero point)
+ M_map = - 2.5*np.log10(L_map)
+
+ # compute the apparent magnitude in each pixel
+ d = distance*100 # kpc -> 10pc
+ m_map = M_map + 5*np.log10(d)
+
+ # compute the surface brightness in each pixel
+ S_map = m_map + 2.5*np.log10(opt.ccd_pixel_area)
+
+
+
+
+
+
+
+
+
+
+ return S_map,xmin,xmax,ymin,ymax
+
+
+
+
+def SkySurfaceBrightness(opt,shape):
+ '''
+ Generate a sky surface brightness
+
+ sky background on earth : 21.8 mag/arcsec2 in V
+ '''
+
+ mean = opt.sky_mean
+ std = (mean - opt.mlim)/3.
+ sky_image = np.random.normal(mean,std,shape)
+
+ print(" mean :",sky_image.mean())
+ print(" 3*std :",3*sky_image.std())
+ print(" mean-3*std :",sky_image.mean()-3*sky_image.std())
+
+ return sky_image
+
+
+
+
+
+def AddSkyBackground(opt,image):
+ """
+ !!! this must be corrected !!! the filter is bad
+ """
+
+ if opt.sky_background:
+ sky = SkySurfaceBrightness(opt,image.shape)
+
+
+ M0 = 4.81
+ Fg = 10**((-image+M0)/2.5)
+ Fs = 10**((-sky +M0)/2.5)
+
+ image = -2.5*np.log10(Fg + Fs) + M0
+
+ plt.imshow(image)
+ plt.show()
+
+
+
+ return image
+
+
+def Display(opt,ax,image,xmin,xmax,ymin,ymax):
+
+
+ mmin = 25 # min magnitude
+ mmax = 36 # max magnitude
+
+
+ if opt.split_colorbar:
+ mlim = opt.mlim # limit magnitude
+
+ n1 = int(256* (mlim-mmin)/(mmax-mmin))
+ n2 = 256-n1
+
+ colors1 = plt.cm.gist_heat_r(np.linspace(0.75, 0.25, n1))
+ #colors2 = plt.cm.binary(np.linspace(0.25, 0.0, n2))
+ colors2 = plt.cm.binary(np.linspace(0.25, 0.0, n2))
+
+ # combine them and build a new colormap
+ colors = np.vstack((colors1, colors2))
+ mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)
+
+ else:
+ mymap = plt.cm.Greys_r
+
+ # use kpc units
+ xmin=opt.xmin_kpc
+ xmax=opt.xmax_kpc
+ ymin=opt.ymin_kpc
+ ymax=opt.ymax_kpc
+
+ im = ax.imshow(image,aspect='equal',extent=(xmin,xmax,ymin,ymax),cmap=mymap,interpolation=None,vmin=mmin,vmax=mmax)
+
+
+ return im
+
+
+
+
+def DisplayDiff(opt,ax,image,xmin,xmax,ymin,ymax):
+
+
+ mmin = 25 # min magnitude
+ mmax = 36 # max magnitude
+
+
+ if opt.split_colorbar:
+ mlim = opt.mlim # limit magnitude
+
+ mlmin = opt.mlmin
+ mlmax = opt.mlmax
+
+ n1 = int(256* (mlmin-mmin)/(mmax-mmin))
+ n2 = int(256* (mlmax-mmin)/(mmax-mmin)) - n1
+ n3 = 256 - int(256* (mlmax-mmin)/(mmax-mmin))
+
+
+ colors1 = plt.cm.binary(np.linspace(0.0, 0.0, n1))
+ colors2 = plt.cm.gist_heat_r(np.linspace(0.7, 0.3, n2))
+ colors3 = plt.cm.binary(np.linspace(0.0, 0.0, n3))
+
+
+ # combine them and build a new colormap
+ colors = np.vstack((colors1,colors2,colors3))
+ mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)
+
+ else:
+ mymap = plt.cm.Greys_r
+
+ # use kpc units
+ xmin=opt.xmin_kpc
+ xmax=opt.xmax_kpc
+ ymin=opt.ymin_kpc
+ ymax=opt.ymax_kpc
+
+ im = ax.imshow(image,aspect='equal',extent=(xmin,xmax,ymin,ymax),cmap=mymap,interpolation=None,vmin=mmin,vmax=mmax)
+
+
+ return im
+
+
+
+
+def DisplaySingle(opt,ax,image,xmin,xmax,ymin,ymax):
+
+
+
+ mmin = opt.mlmin
+ mmax = opt.mlmax
+
+ opt.split_colorbar=False
+
+ if opt.split_colorbar:
+ #mlim = opt.mlim # limit magnitude
+
+ #mlmin = opt.mlmin
+ #mlmax = opt.mlmax
+
+ #n1 = int(256* (mlmin-mmin)/(mmax-mmin))
+ #n2 = int(256* (mlmax-mmin)/(mmax-mmin)) - n1
+ #n3 = 256 - int(256* (mlmax-mmin)/(mmax-mmin))
+
+
+ #colors1 = plt.cm.binary(np.linspace(0.0, 0.0, n1))
+ colors2 = plt.cm.gist_heat(np.linspace(0.3, 1.0, 255))
+ #colors3 = plt.cm.binary(np.linspace(0.0, 0.0, n3))
+
+
+ # combine them and build a new colormap
+ #colors = np.vstack((colors1,colors2,colors3))
+ mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors2)
+
+ #mymap = plt.cm.gist_heat
+
+ else:
+ mymap = plt.cm.gist_heat
+
+ # use kpc units
+ xmin=opt.xmin_kpc
+ xmax=opt.xmax_kpc
+ ymin=opt.ymin_kpc
+ ymax=opt.ymax_kpc
+
+ im = ax.imshow(image,aspect='equal',extent=(xmin,xmax,ymin,ymax),cmap=mymap,interpolation=None,vmin=mmin,vmax=mmax)
+
+
+ return im
diff --git a/ParticlesToSurfaceBrightness/IMF.py b/ParticlesToSurfaceBrightness/IMF.py
new file mode 100644
index 0000000..f82e34d
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/IMF.py
@@ -0,0 +1,232 @@
+import numpy as np
+
+####################
+ ## IMF subroutines
+####################
+
+def slope_imf(x,p1,p2,p3,kn1,kn2):
+
+#Is calculating a three slope IMF
+#INPUT:
+# x = An array of masses for which the IMF should be calculated
+# p1..p3 = the slopes of the power law
+# kn1, kn2 = Where the breaks of the power law are
+#OUTPUT:
+# An array of frequencies matching the mass base array x
+ if(x > kn2):
+ t = (pow(kn2,p2)/pow(kn2,p3))*pow(x,p3+1)
+ elif (x < kn1):
+ t = (pow(kn1,p2)/pow(kn1,p1))*pow(x,p1+1)
+ else:
+ t = pow(x,p2+1)
+ return t
+
+
+def lifetime(m,Z):
+
+#here we will calculate the MS lifetime of the star after Argast et al., 2000, A&A, 356, 873
+#INPUT:¡
+# m = mass in Msun
+# Z = metallicity in Zsun
+
+#OUTPUT:
+# returns the lifetime of the star in Gyrs
+
+ lm = np.log10(m)
+ a0 = 3.79 + 0.24*Z
+ a1 = -3.10 - 0.35*Z
+ a2 = 0.74 + 0.11*Z
+ tmp = a0 + a1*lm + a2*lm*lm
+ return np.divide(np.power(10,tmp),1000)
+
+
+class IMF(object):
+
+#This class represents the IMF normed to 1 in units of M_sun.
+#Input for initialisation:
+
+# mmin = minimal mass of the IMF
+
+# mmax = maximal mass of the IMF
+
+# intervals = how many steps inbetween mmin and mmax should be given
+
+#Then one of the IMF functions can be used
+
+# self.x = mass base
+
+# self.dn = the number of stars at x
+
+# self.dm = the masses for each mass interval x
+
+ def __init__(self, mmin = 0.08 , mmax = 100., intervals = 5000):
+ self.mmin = mmin
+ self.mmax = mmax
+ self.intervals = intervals
+ self.x = np.linspace(mmin,mmax,intervals)
+ self.dx = self.x[1]-self.x[0]
+
+ def normed_3slope(self,paramet = (-1.3,-2.2,-2.7,0.5,1.0)):
+
+# Three slope IMF, Kroupa 1993 as a default
+
+ s1,s2,s3,k1,k2 = paramet
+ u = np.zeros_like(self.x)
+ v = np.zeros_like(self.x)
+ for i in range(len(self.x)):
+ u[i] = slope_imf(self.x[i],s1,s2,s3,k1,k2)
+ v = np.divide(u,self.x)
+ self.dm = np.divide(u,sum(u))
+ self.dn = np.divide(self.dm,self.x)
+ return(self.dm,self.dn)
+
+
+ def Chabrier_1(self, paramet = (0.69, 0.079, -2.3)):
+
+# Chabrier IMF from Chabrier 2003 equation 17 field IMF with variable high mass slope and automatic normalisation
+
+ sigma, m_c, expo = paramet
+ dn = np.zeros_like(self.x)
+ for i in range(len(self.x)):
+ if self.x[i] <= 1:
+ index_with_mass_1 = i
+ dn[i] = (1. / float(self.x[i])) * np.exp(-1*(((np.log10(self.x[i] / m_c))**2)/(2*sigma**2)))
+ else:
+ dn[i] = (pow(self.x[i],expo))
+ # Need to 'attach' the upper to the lower branch
+ derivative_at_1 = dn[index_with_mass_1] - dn[index_with_mass_1 - 1]
+ target_y_for_m_plus_1 = dn[index_with_mass_1] + derivative_at_1
+ rescale = target_y_for_m_plus_1 / dn[index_with_mass_1 + 1]
+ dn[np.where(self.x>1.)] *= rescale
+ # Normalizing to 1 in mass space
+ self.dn = np.divide(dn,sum(dn))
+ dm = dn*self.x
+ self.dm = np.divide(dm,sum(dm))
+ self.dn = np.divide(self.dm,self.x)
+ return(self.dm,self.dn)
+
+
+ def Chabrier_2(self,paramet = (22.8978, 716.4, 0.25,-2.3)):
+
+# Chabrier IMF from Chabrier 2001, IMF 3 = equation 8 parameters from table 1
+
+
+ A,B,sigma,expo = paramet
+ expo -= 1. ## in order to get an alpha index normalisation
+ dn = np.zeros_like(self.x)
+ for i in range(len(self.x)):
+ dn[i] = A*(np.exp(-pow((B/self.x[i]),sigma)))*pow(self.x[i],expo)
+ self.dn = np.divide(dn,sum(dn))
+ dm = dn*self.x
+ self.dm = np.divide(dm,sum(dm))
+ self.dn = np.divide(self.dm,self.x)
+ return(self.dm,self.dn)
+
+
+ def salpeter(self, alpha = (2.35)):
+
+# Salpeter IMF
+
+# Input the slope of the IMF
+
+ self.alpha = alpha
+ temp = np.power(self.x,-self.alpha)
+ norm = sum(temp)
+ self.dn = np.divide(temp,norm)
+ u = self.dn*self.x
+ self.dm = np.divide(u,sum(u))
+ self.dn = np.divide(self.dm,self.x)
+ return (self.dm,self.dn)
+
+
+ def BrokenPowerLaw(self, paramet):
+ breaks,slopes = paramet
+ if len(breaks) != len(slopes)-1:
+ print("error in the precription of the power law. It needs one more slope than break value")
+ else:
+ dn = np.zeros_like(self.x)
+ self.breaks = breaks
+ self.slopes = slopes
+ self.mass_range = np.hstack((self.mmin,breaks,self.mmax))
+ for i,item in enumerate(self.slopes):
+ cut = np.where(np.logical_and(self.x>=self.mass_range[i],self.x<self.mass_range[i+1]))
+ dn[cut] = np.power(self.x[cut],item)
+ if i != 0:
+ renorm = np.divide(last_dn,dn[cut][0])
+ dn[cut] = dn[cut]*renorm
+ last_dn = dn[cut][-1]
+ last_x = self.x[cut][-1]
+ self.dn = np.divide(dn,sum(dn))
+ u = self.dn*self.x
+ self.dm = np.divide(u,sum(u))
+ self.dn = np.divide(self.dm,self.x)
+ return (self.dm,self.dn)
+
+
+ def imf_mass_fraction(self,mlow,mup):
+
+# Calculates the mass fraction of the IMF sitting between mlow and mup
+
+ norm = sum(self.dm)
+ cut = np.where(np.logical_and(self.x>=mlow,self.x<mup))
+ fraction = np.divide(sum(self.dm[cut]),norm)
+ return(fraction)
+
+ def imf_number_fraction(self,mlow,mup):
+
+# Calculating the number fraction of stars of the IMF sitting between mlow and mup
+
+ norm = sum(self.dn)
+ cut = np.where(np.logical_and(self.x>=mlow,self.x<mup))
+ fraction = np.divide(sum(self.dn[cut]),norm)
+ return(fraction)
+
+ def imf_number_stars(self,mlow,mup):
+ cut = np.where(np.logical_and(self.x>=mlow,self.x<mup))
+ number = sum(self.dn[cut])
+ return(number)
+
+ def stochastic_sampling(self, mass):
+
+# The analytic IMF will be resampled according to the mass of the SSP.
+# The IMF will still be normalised to 1
+
+# Stochastic sampling is realised by fixing the number of expected stars and then drawing from the probability distribution of the number density
+# Statistical properties are tested for this sampling and are safe: number of stars and masses converge.
+
+ number_of_stars = int(round(sum(self.dn) * mass))
+ self.dm_copy = np.copy(self.dm)
+ self.dn_copy = np.copy(self.dn)
+
+ #self.dn_copy = np.divide(self.dn_copy,sum(self.dn_copy))
+ random_number = np.random.uniform(low = 0.0, high = sum(self.dn_copy), size = number_of_stars)
+ self.dm = np.zeros_like(self.dm)
+ self.dn = np.zeros_like(self.dn)
+
+ ### This could be favourable if the number of stars drawn is low compared to the imf resolution
+# for i in range(number_of_stars):
+ ### the next line randomly draws a mass according to the number distribution of stars
+# cut = np.where(np.abs(np.cumsum(self.dn_copy)-random_number[i])== np.min(np.abs(np.cumsum(self.dn_copy)-random_number[i])))
+# x1 = self.x[cut][0]
+ #self.dn[cut] += 0.5
+# self.dn[cut[0]] += 1
+# self.dm[cut[0]] += x1 + self.dx/2.
+# t.append(x1 + self.dx/2.)
+
+ counting = np.cumsum(self.dn_copy)
+ for i in range(len(counting)-1):
+ if i == 0:
+ cut = np.where(np.logical_and(random_number>0.,random_number<=counting[i]))
+ else:
+ cut = np.where(np.logical_and(random_number>counting[i-1],random_number<=counting[i]))
+ number_of_stars_in_mass_bin = len(random_number[cut])
+ self.dm[i] = number_of_stars_in_mass_bin * self.x[i]
+ if number_of_stars:
+ self.dm = np.divide(self.dm,sum(self.dm))
+ else:
+ self.dm = np.divide(self.dm, 1)
+ self.dn = np.divide(self.dm,self.x)
+ return (self.dm,self.dn)
+########
+ ## End IMF subroutines
+########
diff --git a/ParticlesToSurfaceBrightness/Readme b/ParticlesToSurfaceBrightness/Readme
new file mode 100644
index 0000000..223cea5
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/Readme
@@ -0,0 +1,48 @@
+# requirement : pNbody https://lastro.epfl.ch/projects/pNbody/
+
+
+#################################################
+# download snapshots containing only stars
+# we need to find another repo for big data
+
+stars_1keV_1000kpc.hdf5
+stars_3keV_1000kpc.hdf5
+stars_10keV_1000kpc.hdf5
+stars_30keV_1000kpc.hdf5
+
+#################################################
+# add additional information (age, magnitudes, html)
+
+./arrakihs_addfields.py -o stars_1keV_1000kpc.hdf5 stars_1keV_1000kpc.hdf5
+./arrakihs_addfields.py -o stars_3keV_1000kpc.hdf5 stars_3keV_1000kpc.hdf5
+./arrakihs_addfields.py -o stars_10keV_1000kpc.hdf5 stars_10keV_1000kpc.hdf5
+./arrakihs_addfields.py -o stars_30keV_1000kpc.hdf5 stars_30keV_1000kpc.hdf5
+
+
+################################################
+# compute a set of images
+
+./arrakihs_compute_images.py --telescope arrakihs_vis --filter F475X --fov 5000 --distance 25 --frsp 0.25 --nlos 16 -o F475X_1keV.pkl stars_1keV_1000kpc.hdf5
+./arrakihs_compute_images.py --telescope arrakihs_vis --filter VISeuclid --fov 5000 --distance 25 --frsp 0.25 --nlos 16 -o VISeuclid_3keV.pkl stars_3keV_1000kpc.hdf5
+./arrakihs_compute_images.py --telescope arrakihs_vis --filter Yeuclid --fov 5000 --distance 25 --frsp 0.25 --nlos 16 -o Yeuclid_1keV.pkl stars_1keV_1000kpc.hdf5
+./arrakihs_compute_images.py --telescope arrakihs_vis --filter Jeuclid --fov 5000 --distance 25 --frsp 0.25 --nlos 16 -o Jeuclid_1keV.pkl stars_1keV_1000kpc.hdf5
+
+################################################
+# display or save images in the pkl file
+
+./arrakihs_show_images.py F475X_1keV.pkl -o F475X_1keV.png
+./arrakihs_show_images.py VISeuclid_3keV.pkl -o VISeuclid_3keV.png
+./arrakihs_show_images.py Yeuclid_1keV.pkl -o Yeuclid_1keV.png
+./arrakihs_show_images.py Jeuclid_1keV.pkl -o Jeuclid_1keV.png
+
+###############################################
+# compute and compare the fraction of illuminated pixels
+
+./arrakihs_compute_illuminated.py -o F475X_coverage_fraction.png F475X_1keV.pkl F475X_3keV.pkl F475X_10keV.pkl F475X_30keV.pkl
+./arrakihs_compute_illuminated.py -o VISeuclid_coverage_fraction.png VISeuclid_1keV.pkl VISeuclid_3keV.pkl VISeuclid_10keV.pkl VISeuclid_30keV.pkl
+./arrakihs_compute_illuminated.py -o Yeuclid_coverage_fraction.png Yeuclid_1keV.pkl Yeuclid_3keV.pkl Yeuclid_10keV.pkl Yeuclid_30keV.pkl
+./arrakihs_compute_illuminated.py -o Jeuclid_coverage_fraction.png Jeuclid_1keV.pkl Jeuclid_3keV.pkl Jeuclid_10keV.pkl Jeuclid_30keV.pkl
+
+
+
+
diff --git a/ParticlesToSurfaceBrightness/arrakihs_addfields.py b/ParticlesToSurfaceBrightness/arrakihs_addfields.py
new file mode 100755
index 0000000..102dfe5
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/arrakihs_addfields.py
@@ -0,0 +1,154 @@
+#!/usr/bin/env python3
+
+import argparse
+import numpy as np
+from pNbody import *
+
+import stars_class
+
+
+####################################################################
+# option parser
+####################################################################
+
+description=""
+epilog =""""""
+
+parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
+
+
+parser.add_argument(action="store",
+ dest="files",
+ metavar='FILE',
+ type=str,
+ default=None,
+ nargs='*',
+ help='a file name')
+
+parser.add_argument("-o",
+ action="store",
+ type=str,
+ dest="outputfilename",
+ default=None,
+ help="Name of the output file")
+
+####################################################################
+# main
+####################################################################
+
+
+if __name__ == '__main__':
+
+ opt = parser.parse_args()
+
+ for f in opt.files:
+ nb = Nbody(f,ftype="arepo")
+
+
+ # remove doublets
+ u,idx = np.unique(nb.rxyz(),return_index=True)
+ nb = nb.selectp(lst=nb.num[idx])
+
+
+ ###################################
+ # Stellar ages
+ ###################################
+ print("Compute ages...")
+ nb.age = nb.StellarAge(units="Gyr")
+ print("done.")
+
+ ###################################
+ # compute Hsml
+ ###################################
+ print("Compute Rsp...")
+ #nb.set_tpe(0)
+ #nb.InitSphParameters(DesNumNgb=32, MaxNumNgbDeviation=2)
+ #nb.getTree()
+ #nb.rsp = nb.get_rsp_approximation()
+ #nb.set_tpe(4)
+ nb.ComputeRsp(5)
+ print("done.")
+
+ ###################################
+ # compute Magnitudes
+ ###################################
+ print("Compute magnitudes...")
+ MH = nb.MH()
+
+ # mass in Msol
+ mass = nb.mass * 1e10
+
+ # get the number of stars in each mass bin
+ Nstars = stars_class.Stars_fun(mass,None,None, 'normed_3slope')
+
+ ##############################
+ # F475X magnitude
+
+ M = stars_class.HST475X_fun(None,nb.age,MH)
+ # convert to flux (ignore the zero point)
+ F = 10**(-M/2.5)
+ # sum the contribution of the mass bins
+ F = np.sum(F*Nstars, axis=0)
+ # compute the absolute magnitude in each pixel (as before we ignore the zero point)
+ M = - 2.5*np.log10(F)
+ nb.MagF475X = M
+
+ ##############################
+ # VIS Euclid magnitude
+
+ M = stars_class.VISeuclid_fun(None,nb.age,MH)
+ # convert to flux (ignore the zero point)
+ F = 10**(-M/2.5)
+ # sum the contribution of the mass bins
+ F = np.sum(F*Nstars, axis=0)
+ # compute the absolute magnitude in each pixel (as before we ignore the zero point)
+ M = - 2.5*np.log10(F)
+ nb.MagVISeuclid = M
+
+ ##############################
+ # Y Euclid magnitude
+
+ M = stars_class.Yeuclid_fun(None,nb.age,MH)
+ # convert to flux (ignore the zero point)
+ F = 10**(-M/2.5)
+ # sum the contribution of the mass bins
+ F = np.sum(F*Nstars, axis=0)
+ # compute the absolute magnitude in each pixel (as before we ignore the zero point)
+ M = - 2.5*np.log10(F)
+ nb.MagYeuclid = M
+
+ ##############################
+ # J Euclid magnitude
+
+ M = stars_class.Jeuclid_fun(None,nb.age,MH)
+ # convert to flux (ignore the zero point)
+ F = 10**(-M/2.5)
+ # sum the contribution of the mass bins
+ F = np.sum(F*Nstars, axis=0)
+ # compute the absolute magnitude in each pixel (as before we ignore the zero point)
+ M = - 2.5*np.log10(F)
+ nb.MagJeuclid = M
+
+
+ print("done.")
+
+
+ if opt.outputfilename:
+ nb.rename(opt.outputfilename)
+ nb.write()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ParticlesToSurfaceBrightness/arrakihs_compute_illuminated.py b/ParticlesToSurfaceBrightness/arrakihs_compute_illuminated.py
new file mode 100755
index 0000000..45c31d6
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/arrakihs_compute_illuminated.py
@@ -0,0 +1,180 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+from astropy import constants as cte
+from astropy import units as u
+
+import argparse
+
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+
+from pNbody import *
+from pNbody import ic
+import ARRAKIHS
+
+import pickle
+
+
+####################################################################
+# option parser
+####################################################################
+
+description=""
+epilog =""""""
+
+parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
+
+
+
+
+
+parser.add_argument(action="store",
+ dest="files",
+ metavar='FILE',
+ type=str,
+ default=None,
+ nargs='*',
+ help='a file name')
+
+parser.add_argument("-o",
+ action="store",
+ type=str,
+ dest="outputfilename",
+ default=None,
+ help="Name of the output file")
+
+
+
+parser.add_argument('--mlim',
+ action="store",
+ dest="mlim",
+ metavar='FLOAT',
+ type=float,
+ default=31,
+ help='magnitude limit')
+
+
+
+
+
+parser.add_argument("--mlmin",
+ action="store",
+ dest="mlmin",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='mag lim min')
+
+parser.add_argument("--mlmax",
+ action="store",
+ dest="mlmax",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='mag lim max')
+
+parser.add_argument("--nlos",
+ action="store",
+ dest="nlos",
+ metavar='INT',
+ type=int,
+ default=1,
+ help='number of los')
+
+
+
+
+####################################################################
+# main
+####################################################################
+
+cs = ['r','g','b','y']
+nmagbins = 50
+mag0 = 33
+magmax = mag0
+
+if __name__ == '__main__':
+
+ opt = parser.parse_args()
+
+
+ for j,f in enumerate(opt.files):
+
+ print(f)
+
+ fd = open(f,"rb")
+ images = pickle.load(fd)
+ fd.close()
+
+ # loop over images of the same model
+
+ Npixs = np.zeros((len(images),nmagbins))
+
+
+ for k in range(len(images)):
+
+ image,xmin,xmax,ymin,ymax = images[k]
+
+
+ Npixtot = image.shape[0]*image.shape[1]
+
+ # clean
+ image = np.where(np.isinf(image),100,image)
+
+
+
+ magbins = np.linspace(26,magmax,nmagbins)
+ Npix = np.zeros(nmagbins)
+
+ for i in tqdm(range(nmagbins)):
+ img1 = np.where(image<magbins[i],1,0)
+ Npix[i] = np.ravel(img1).sum()
+
+ # normalize
+ #Npix = Npix/Npixtot
+ #idx = np.argmin(np.fabs(magbins - mag0))
+ #Npix = Npix/Npix[idx]
+ Npix = Npix/Npixtot
+
+
+
+ #plt.plot(magbins,Npix,color=cs[j])
+
+
+ # sum
+ Npixs[k] = Npix
+
+
+
+ mean = Npixs.mean(axis=0)
+ std = Npixs.std(axis=0)
+
+ label = os.path.basename(f)
+ label = label[label.find("_")+1:]
+ label = label[:-4]
+
+
+ plt.plot(magbins,mean,color=cs[j],label=r"$%s$"%label,lw=2)
+ #plt.plot(magbins,mean+1*std,color=cs[j],lw=5)
+ #plt.plot(magbins,mean-1*std,color=cs[j],lw=5)
+
+ ax = plt.gca()
+ ax.fill_between(magbins,mean+1*std,mean-1*std,facecolor=cs[j],alpha=0.1)
+
+
+
+ plt.legend()
+ plt.xlabel(r"$\rm{VIS}\,\,\rm{magnitude}$")
+ plt.ylabel(r"$\rm{coverage\,\,fraction}$")
+
+ if opt.outputfilename is None:
+ plt.show()
+ else:
+ plt.savefig(opt.outputfilename)
+
+
+
+
+
diff --git a/ParticlesToSurfaceBrightness/arrakihs_compute_images.py b/ParticlesToSurfaceBrightness/arrakihs_compute_images.py
new file mode 100755
index 0000000..e8a3055
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/arrakihs_compute_images.py
@@ -0,0 +1,294 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+from astropy import constants as cte
+from astropy import units as u
+
+import argparse
+
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+
+from pNbody import *
+from pNbody import ic
+import ARRAKIHS
+
+import pickle
+
+
+####################################################################
+# option parser
+####################################################################
+
+description=""
+epilog =""""""
+
+parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
+
+
+
+
+
+parser.add_argument(action="store",
+ dest="file",
+ metavar='FILE',
+ type=str,
+ default=None,
+ nargs='*',
+ help='a file name')
+
+parser.add_argument("-o",
+ action="store",
+ type=str,
+ dest="outputfilename",
+ default=None,
+ help="Name of the output file")
+
+
+parser.add_argument("--ccd_shape",
+ action="store",
+ dest="ccd_shape",
+ metavar='INT INT',
+ type=int,
+ default=[4096,4096],
+ nargs=2,
+ help='ccd shape')
+
+parser.add_argument("--size",
+ action="store",
+ dest="size",
+ metavar='INT INT',
+ type=int,
+ default=None,
+ nargs=2,
+ help='image size')
+
+
+
+parser.add_argument("--ccd_size",
+ action="store",
+ dest="ccd_size",
+ metavar='FLOAT FLOAT',
+ type=float,
+ default=[1.5,1.5],
+ nargs=2,
+ help='ccd size in degree')
+
+
+parser.add_argument("--ccd_field_fraction",
+ action="store",
+ dest="ccd_field_fraction",
+ metavar='FLOAT',
+ type=float,
+ default=1,
+ help='fraction of the field displayed')
+
+
+parser.add_argument("--fov",
+ action="store",
+ dest="fov",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='field of view in arcsec')
+
+
+
+parser.add_argument("--psf",
+ action="store",
+ dest="psf",
+ metavar='FLOAT',
+ type=float,
+ default=0.8,
+ help='psf in arcsec')
+
+
+parser.add_argument("--convolve",
+ action="store_true",
+ dest="convolve",
+ default=False,
+ help='convolve with the pdf')
+
+
+parser.add_argument('--distance',
+ action="store",
+ dest="distance",
+ metavar='FLOAT',
+ type=float,
+ default=20,
+ help='distance of the object in Mpc')
+
+
+
+parser.add_argument('--mlim',
+ action="store",
+ dest="mlim",
+ metavar='FLOAT',
+ type=float,
+ default=31,
+ help='magnitude limit')
+
+
+parser.add_argument("--sky_background",
+ action="store_true",
+ dest="sky_background",
+ default=False,
+ help='add sky background')
+
+
+parser.add_argument("--sky_N",
+ action="store",
+ dest="sky_N",
+ metavar='FLOAT',
+ type=float,
+ default=1e6,
+ help='number of particles in the sky')
+
+parser.add_argument("--sky_M",
+ action="store",
+ dest="sky_M",
+ metavar='FLOAT',
+ type=float,
+ default=3.4e5,
+ help='total sky mass [Msol]')
+
+parser.add_argument("--sky_e",
+ action="store",
+ dest="sky_e",
+ metavar='FLOAT',
+ type=float,
+ default=3,
+ help='sky size [kpc]')
+
+parser.add_argument("--sky_mean",
+ action="store",
+ dest="sky_mean",
+ metavar='FLOAT',
+ type=float,
+ default=34,
+ help='sky mean surface brighness')
+
+
+
+parser.add_argument("--telescope",
+ action="store",
+ dest="telescope",
+ metavar='STR',
+ type=str,
+ default='arrakihs_vis',
+ help='telescope name')
+
+parser.add_argument("--filter",
+ action="store",
+ dest="filter",
+ metavar='STR',
+ type=str,
+ default=None,
+ help='filter name')
+
+
+parser.add_argument("--toMsol",
+ action="store",
+ dest="toMsol",
+ metavar='FLOAT',
+ type=float,
+ default=1e10,
+ help='mass convertion to solar mass')
+
+parser.add_argument("--frsp",
+ action="store",
+ dest="frsp",
+ metavar='FLOAT',
+ type=float,
+ default=0,
+ help='factor for the image smoothing')
+
+
+parser.add_argument("--view",
+ action="store",
+ dest="view",
+ metavar='STR',
+ type=str,
+ default='xy',
+ help='view xy, yz, zx')
+
+parser.add_argument("--mlmin",
+ action="store",
+ dest="mlmin",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='mag lim min')
+
+parser.add_argument("--mlmax",
+ action="store",
+ dest="mlmax",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='mag lim max')
+
+parser.add_argument("--nlos",
+ action="store",
+ dest="nlos",
+ metavar='INT',
+ type=int,
+ default=1,
+ help='number of los')
+
+
+
+
+####################################################################
+# main
+####################################################################
+
+def get_axes(n=7,irand=0):
+
+ np.random.seed(irand)
+ random2 = np.random.random(n)
+ random3 = np.random.random(n)
+
+ phi = random2 * np.pi * 2.
+ costh = 1. - 2. * random3
+
+ sinth = np.sqrt(1. - costh**2)
+
+ x = sinth * np.cos(phi)
+ y = sinth * np.sin(phi)
+ z = costh
+
+ pos = np.transpose(np.array([x, y, z]))
+
+ return pos
+
+
+
+
+if __name__ == '__main__':
+
+ opt = parser.parse_args()
+
+ naxes = opt.nlos
+ axes = get_axes(naxes)
+
+
+ ARRAKIHS.SetTelecope(opt)
+
+ images = []
+
+ # loop over axes
+ for i in range(naxes):
+ opt.los=axes[i]
+ image = ARRAKIHS.Project(opt) # here we get: image,xmin,xmax,ymin,ymax
+ images.append(image)
+
+
+ if opt.outputfilename is not None:
+
+ f = open(opt.outputfilename,'wb')
+ pickle.dump(images,f)
+ f.close()
+
+
diff --git a/ParticlesToSurfaceBrightness/arrakihs_show_images.py b/ParticlesToSurfaceBrightness/arrakihs_show_images.py
new file mode 100755
index 0000000..1ef266a
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/arrakihs_show_images.py
@@ -0,0 +1,203 @@
+#!/usr/bin/env python3
+
+import numpy as np
+
+from astropy import constants as cte
+from astropy import units as u
+
+import argparse
+
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+
+from pNbody import *
+from pNbody import ic
+import ARRAKIHS
+
+import pickle
+import os
+
+####################################################################
+# option parser
+####################################################################
+
+description=""
+epilog =""""""
+
+parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
+
+
+
+
+
+parser.add_argument(action="store",
+ dest="file",
+ metavar='FILE',
+ type=str,
+ default=None,
+ nargs='*',
+ help='a file name')
+
+parser.add_argument("-o",
+ action="store",
+ type=str,
+ dest="outputfilename",
+ default=None,
+ help="Name of the output file")
+
+
+
+
+
+
+parser.add_argument('--mlim',
+ action="store",
+ dest="mlim",
+ metavar='FLOAT',
+ type=float,
+ default=31,
+ help='magnitude limit')
+
+
+
+
+
+parser.add_argument("--mlmin",
+ action="store",
+ dest="mlmin",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='mag lim min')
+
+parser.add_argument("--mlmax",
+ action="store",
+ dest="mlmax",
+ metavar='FLOAT',
+ type=float,
+ default=None,
+ help='mag lim max')
+
+parser.add_argument("--nlos",
+ action="store",
+ dest="nlos",
+ metavar='INT',
+ type=int,
+ default=1,
+ help='number of los')
+
+
+parser.add_argument("--filter",
+ action="store",
+ dest="filter",
+ metavar='STR',
+ type=str,
+ default=None,
+ help='filter name')
+
+
+parser.add_argument("--MDM",
+ action="store",
+ dest="MDM",
+ metavar='FLOAT',
+ type=float,
+ default=1,
+ help='dark mattermass')
+
+
+
+####################################################################
+# main
+####################################################################
+
+
+if __name__ == '__main__':
+
+ opt = parser.parse_args()
+
+
+ if len(opt.file)==0:
+ print("need to provie one file")
+ exit()
+
+
+ f = open(opt.file[0],"rb")
+ images = pickle.load(f)
+ f.close()
+
+ nx = int(np.sqrt(len(images)))
+ ny = nx
+
+
+ fig, ax = plt.subplots(nx,ny)
+
+ # plot
+ fig = plt.gcf()
+ fig.set_size_inches(12*1.15, 12)
+
+ fig.subplots_adjust(left=0.1)
+ fig.subplots_adjust(right=1)
+ fig.subplots_adjust(bottom=0.05)
+ fig.subplots_adjust(top=0.95)
+ fig.subplots_adjust(wspace=0.02)
+ fig.subplots_adjust(hspace=0.02)
+
+
+ colors2 = plt.cm.gist_heat(np.linspace(0.3, 1.0, 255))
+ mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors2)
+
+ for i in range(len(images)):
+
+ image,xmin,xmax,ymin,ymax = images[i]
+
+ #image = np.where(image>26,34,image)
+
+
+ ix = i//nx
+ iy = i- (ix*nx)
+
+
+ #ax[i] = plt.gca()
+ ax[ix,iy].set_aspect('equal')
+
+ mmin = 24
+ mmax = 33
+
+ im = ax[ix,iy].imshow(image,aspect='equal',extent=(xmin,xmax,ymin,ymax),interpolation=None,vmin=mmin,vmax=mmax,cmap=mymap)
+ #ax[ix,iy].contour(image,levels=[26],colors='r')
+
+
+ ax[ix,iy].set_xlabel(r"$x\,[\rm{kpc}]$")
+ ax[ix,iy].set_ylabel(r"$y\,[\rm{kpc}]$")
+
+
+ if i != nx*(ny-1):
+ ax[ix,iy].get_xaxis().set_visible(False)
+ ax[ix,iy].get_yaxis().set_visible(False)
+ ax[ix,iy].get_xaxis().set_ticks([])
+ ax[ix,iy].get_yaxis().set_ticks([])
+
+
+ plt.colorbar(im,label="surface brightness [mag/arcsec]",ax=ax)
+
+ basename = os.path.basename(opt.file[0])
+ opt.filter = basename.split("_")[0]
+ opt.MDM = basename.split("_")[1].split("keV")[0]
+
+ txt = r"$\rm{%s}\,\,:\,\,%s\,\rm{keV}$"%(opt.filter,opt.MDM)
+ fig.suptitle(txt,fontsize=20)
+
+
+
+
+ if opt.outputfilename:
+ plt.savefig(opt.outputfilename)
+ else:
+ plt.show()
+
+
+
+
+
+
+
diff --git a/ParticlesToSurfaceBrightness/stars_class.py b/ParticlesToSurfaceBrightness/stars_class.py
new file mode 100644
index 0000000..811eda3
--- /dev/null
+++ b/ParticlesToSurfaceBrightness/stars_class.py
@@ -0,0 +1,1449 @@
+import numpy as np
+from IMF import *
+
+
+
+def Stars_fun(mass, age, metals, imf = 'normed_3slope'):
+ """
+ This function calculates the number of stars for a given stellar mass particle, for different IMFs:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ Metallicity of the particle [M/H]
+ imf: str.
+ Initial Mass Function. Eligible: 'normed_3slope', 'Chabier_1', 'Chabier_2', 'sapeter' or 'BrokenPowerLaw'. By default: 'normed_3slope'.
+
+ """
+ init_imf=IMF(mmin=0.05, mmax=100.0, intervals=5000)
+ if imf=='Chabrier_1':
+ IMF_type=init_imf.Chabrier_1() # Chabrier_1() # Chabrier_2(), salpeter(), BrokenPowerLaw(),
+ # normed_3slope() that is Kroupa.
+ if imf=='Chabrier_2':
+ IMF_type=init_imf.Chabrier_2()
+
+ if imf=='salpeter':
+ IMF_type=init_imf.salpeter()
+
+ if imf=='normed_3slope':
+ IMF_type=init_imf.normed_3slope()
+
+
+ M1=init_imf.imf_mass_fraction(0.0,0.1)
+ M2=init_imf.imf_mass_fraction(0.1,0.2)
+ M3=init_imf.imf_mass_fraction(0.2,0.3)
+ M4=init_imf.imf_mass_fraction(0.3,0.4)
+ M5=init_imf.imf_mass_fraction(0.4,0.5)
+ M6=init_imf.imf_mass_fraction(0.5,0.6)
+ M7=init_imf.imf_mass_fraction(0.6,0.7)
+ M8=init_imf.imf_mass_fraction(0.7,0.8)
+ M9=init_imf.imf_mass_fraction(0.8,0.9)
+ M10=init_imf.imf_mass_fraction(0.9,1)
+ M11=init_imf.imf_mass_fraction(1,2)
+ M12=init_imf.imf_mass_fraction(2,3)
+ M13=init_imf.imf_mass_fraction(3,4)
+ M14=init_imf.imf_mass_fraction(4,5)
+ M15=init_imf.imf_mass_fraction(5,6)
+ M16=init_imf.imf_mass_fraction(6,7)
+ M17=init_imf.imf_mass_fraction(7,8)
+ M18=init_imf.imf_mass_fraction(8,101.)
+
+
+
+ n_1=[M1*i/0.1 for i in mass]
+ n_2=[M2*i/0.2 for i in mass]
+ n_3=[M3*i/0.3 for i in mass]
+ n_4=[M4*i/0.4 for i in mass]
+ n_5=[M5*i/0.5 for i in mass]
+ n_6=[M6*i/0.6 for i in mass]
+ n_7=[M7*i/0.7 for i in mass]
+ n_8=[M8*i/0.8 for i in mass]
+ n_9=[M9*i/0.9 for i in mass]
+ n_10=[M10*i/1 for i in mass]
+ n_11=[M11*i/2 for i in mass]
+ n_12=[M12*i/3 for i in mass]
+ n_13=[M13*i/4 for i in mass]
+ n_14=[M14*i/5 for i in mass]
+ n_15=[M15*i/6 for i in mass]
+ n_16=[M16*i/7 for i in mass]
+ n_17=[M17*i/8 for i in mass]
+ n_18=[M18*i/100 for i in mass]
+
+
+ return np.array([n_1,n_2, n_3,n_4,n_5,n_6, n_7,n_8,n_9,n_10, n_11, n_12,n_13, n_14, n_15,n_16,n_17,n_18])
+
+
+def Temp_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the temperature of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_euclid.txt',comments='#')
+ MH,logAge,Mini,Te,logL=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4]
+ # Difererentes edades de las estrellas del cmd
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+ t1=Te[k+1]
+ t2=Te[k+2]
+ t3=Te[k+3]
+ t4=Te[k+4]
+ t5=Te[k+5]
+ t6=Te[k+6]
+ t7=Te[k+7]
+ t8=Te[k+8]
+ t9=Te[k+9]
+ t10=Te[k+10]
+ t11=Te[k+11]
+ t12=Te[k+12]
+ t13=Te[k+13]
+ t14=Te[k+14]
+ t15=Te[k+15]
+ t16=Te[k+16]
+ t17=Te[k+17]
+ t18=Te[k+18]
+
+
+ return t1,t2,t3,t4,t5,t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18
+
+
+
+def Luminosity_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the intrinsic luminosity of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_euclid.txt',comments='#')
+ MH,logAge,Mini,Te,logL=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1=10**logL[k+1]
+ l2=10**logL[k+2]
+ l3=10**logL[k+3]
+ l4=10**logL[k+4]
+ l5=10**logL[k+5]
+ l6=10**logL[k+6]
+ l7=10**logL[k+7]
+ l8=10**logL[k+8]
+ l9=10**logL[k+9]
+ l10=10**logL[k+10]
+ l11=10**logL[k+11]
+ l12=10**logL[k+12]
+ l13=10**logL[k+13]
+ l14=10**logL[k+14]
+ l15=10**logL[k+15]
+ l16=10**logL[k+16]
+ l17=10**logL[k+17]
+ l18=10**logL[k+18]
+
+
+ return l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18
+
+
+#####################################
+########## HST FUNCTIONS
+#####################################
+
+
+
+
+def HST475X_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on HST-F475X from WFC3 of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_hst_wfc3.txt',comments='#')
+ MH,logAge,Mini,Te,logL, HST475Xmag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,5]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+ # yr: the following lines gives the wrong results !!!
+ '''
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+ '''
+
+ # yr
+ get_age_idx = np.vectorize(lambda x:np.argmin( np.fabs( x*1e9 -10**age_l )))
+ get_met_idx = np.vectorize(lambda x:np.argmin( np.fabs( x -met_l )))
+
+ idx_age = get_age_idx(age)
+ idx_met = get_met_idx(metals)
+
+ k=(idx_met*(age_len+1)*(mass_len+1)+idx_age*(mass_len+1))
+ k = k-1 # really ???
+
+
+
+ l1 = HST475Xmag[k+1]
+ l2 = HST475Xmag[k+2]
+ l3 = HST475Xmag[k+3]
+ l4 = HST475Xmag[k+4]
+ l5 = HST475Xmag[k+5]
+ l6 = HST475Xmag[k+6]
+ l7 = HST475Xmag[k+7]
+ l8 = HST475Xmag[k+8]
+ l9 = HST475Xmag[k+9]
+ l10 = HST475Xmag[k+10]
+ l11 = HST475Xmag[k+11]
+ l12 = HST475Xmag[k+12]
+ l13 = HST475Xmag[k+13]
+ l14 = HST475Xmag[k+14]
+ l15 = HST475Xmag[k+15]
+ l16 = HST475Xmag[k+16]
+ l17 = HST475Xmag[k+17]
+ l18 = HST475Xmag[k+18]
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+#####################################
+########## EUCLID FUNCTIONS
+#####################################
+
+def VISeuclid_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on Euclid-VIS of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_euclid.txt',comments='#')
+ MH,logAge,Mini,Te,logL, VISmag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,5]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+ # yr: the following lines gives the wrong results !!!
+ '''
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+ '''
+
+ # yr
+ get_age_idx = np.vectorize(lambda x:np.argmin( np.fabs( x*1e9 -10**age_l )))
+ get_met_idx = np.vectorize(lambda x:np.argmin( np.fabs( x -met_l )))
+
+ idx_age = get_age_idx(age)
+ idx_met = get_met_idx(metals)
+
+ k=(idx_met*(age_len+1)*(mass_len+1)+idx_age*(mass_len+1))
+ k = k-1 # really ???
+
+
+
+ l1 = VISmag[k+1]
+ l2 = VISmag[k+2]
+ l3 = VISmag[k+3]
+ l4 = VISmag[k+4]
+ l5 = VISmag[k+5]
+ l6 = VISmag[k+6]
+ l7 = VISmag[k+7]
+ l8 = VISmag[k+8]
+ l9 = VISmag[k+9]
+ l10 = VISmag[k+10]
+ l11 = VISmag[k+11]
+ l12 = VISmag[k+12]
+ l13 = VISmag[k+13]
+ l14 = VISmag[k+14]
+ l15 = VISmag[k+15]
+ l16 = VISmag[k+16]
+ l17 = VISmag[k+17]
+ l18 = VISmag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+def Yeuclid_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on Euclid-Y of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_euclid.txt',comments='#')
+ MH,logAge,Mini,Te,logL, Ymag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,7]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+ # yr: the following lines gives the wrong results !!!
+ '''
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+ '''
+
+ # yr
+ get_age_idx = np.vectorize(lambda x:np.argmin( np.fabs( x*1e9 -10**age_l )))
+ get_met_idx = np.vectorize(lambda x:np.argmin( np.fabs( x -met_l )))
+
+ idx_age = get_age_idx(age)
+ idx_met = get_met_idx(metals)
+
+ k=(idx_met*(age_len+1)*(mass_len+1)+idx_age*(mass_len+1))
+ k = k-1 # really ???
+
+
+ l1 = Ymag[k+1]
+ l2 = Ymag[k+2]
+ l3 = Ymag[k+3]
+ l4 = Ymag[k+4]
+ l5 = Ymag[k+5]
+ l6 = Ymag[k+6]
+ l7 = Ymag[k+7]
+ l8 = Ymag[k+8]
+ l9 = Ymag[k+9]
+ l10 = Ymag[k+10]
+ l11 = Ymag[k+11]
+ l12 = Ymag[k+12]
+ l13 = Ymag[k+13]
+ l14 = Ymag[k+14]
+ l15 = Ymag[k+15]
+ l16 = Ymag[k+16]
+ l17 = Ymag[k+17]
+ l18 = Ymag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+def Jeuclid_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on Euclid-J of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_euclid.txt',comments='#')
+ MH,logAge,Mini,Te,logL, Jmag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,6]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+ # yr: the following lines gives the wrong results !!!
+ '''
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+ '''
+
+ # yr
+ get_age_idx = np.vectorize(lambda x:np.argmin( np.fabs( x*1e9 -10**age_l )))
+ get_met_idx = np.vectorize(lambda x:np.argmin( np.fabs( x -met_l )))
+
+ idx_age = get_age_idx(age)
+ idx_met = get_met_idx(metals)
+
+ k=(idx_met*(age_len+1)*(mass_len+1)+idx_age*(mass_len+1))
+ k = k-1 # really ???
+
+
+ l1 = Jmag[k+1]
+ l2 = Jmag[k+2]
+ l3 = Jmag[k+3]
+ l4 = Jmag[k+4]
+ l5 = Jmag[k+5]
+ l6 = Jmag[k+6]
+ l7 = Jmag[k+7]
+ l8 = Jmag[k+8]
+ l9 = Jmag[k+9]
+ l10 = Jmag[k+10]
+ l11 = Jmag[k+11]
+ l12 = Jmag[k+12]
+ l13 = Jmag[k+13]
+ l14 = Jmag[k+14]
+ l15 = Jmag[k+15]
+ l16 = Jmag[k+16]
+ l17 = Jmag[k+17]
+ l18 = Jmag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+
+
+
+#####################################
+########## SLOAN FUNCTIONS
+#####################################
+
+
+def rsloan_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on sloan r of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_sloan.txt',comments='#')
+ MH,logAge,Mini,Te,logL, rmag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,5]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = rmag[k+1]
+ l2 = rmag[k+2]
+ l3 = rmag[k+3]
+ l4 = rmag[k+4]
+ l5 = rmag[k+5]
+ l6 = rmag[k+6]
+ l7 = rmag[k+7]
+ l8 = rmag[k+8]
+ l9 = rmag[k+9]
+ l10 = rmag[k+10]
+ l11 = rmag[k+11]
+ l12 = rmag[k+12]
+ l13 = rmag[k+13]
+ l14 = rmag[k+14]
+ l15 = rmag[k+15]
+ l16 = rmag[k+16]
+ l17 = rmag[k+17]
+ l18 = rmag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+def gsloan_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on sloan g of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_sloan.txt',comments='#')
+ MH,logAge,Mini,Te,logL, gmag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,6]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = gmag[k+1]
+ l2 = gmag[k+2]
+ l3 = gmag[k+3]
+ l4 = gmag[k+4]
+ l5 = gmag[k+5]
+ l6 = gmag[k+6]
+ l7 = gmag[k+7]
+ l8 = gmag[k+8]
+ l9 = gmag[k+9]
+ l10 = gmag[k+10]
+ l11 = gmag[k+11]
+ l12 = gmag[k+12]
+ l13 = gmag[k+13]
+ l14 = gmag[k+14]
+ l15 = gmag[k+15]
+ l16 = gmag[k+16]
+ l17 = gmag[k+17]
+ l18 = gmag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+def isloan_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on sloan g of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_sloan.txt',comments='#')
+ MH,logAge,Mini,Te,logL, imag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,7]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = imag[k+1]
+ l2 = imag[k+2]
+ l3 = imag[k+3]
+ l4 = imag[k+4]
+ l5 = imag[k+5]
+ l6 = imag[k+6]
+ l7 = imag[k+7]
+ l8 = imag[k+8]
+ l9 = imag[k+9]
+ l10 = imag[k+10]
+ l11 = imag[k+11]
+ l12 = imag[k+12]
+ l13 = imag[k+13]
+ l14 = imag[k+14]
+ l15 = imag[k+15]
+ l16 = imag[k+16]
+ l17 = imag[k+17]
+ l18 = imag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+#####################################
+########## GALEX FUNCTIONS
+#####################################
+
+def NUVmags_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on GALEX NUV of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_galex.txt',comments='#')
+ MH,logAge,Mini,Te,logL, NUVmag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,7]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+ l1 = NUVmag[k+1]
+ l2 = NUVmag[k+2]
+ l3 = NUVmag[k+3]
+ l4 = NUVmag[k+4]
+ l5 = NUVmag[k+5]
+ l6 = NUVmag[k+6]
+ l7 = NUVmag[k+7]
+ l8 = NUVmag[k+8]
+ l9 = NUVmag[k+9]
+ l10 = NUVmag[k+10]
+ l11 = NUVmag[k+11]
+ l12 = NUVmag[k+12]
+ l13 = NUVmag[k+13]
+ l14 = NUVmag[k+14]
+ l15 = NUVmag[k+15]
+ l16 = NUVmag[k+16]
+ l17 = NUVmag[k+17]
+ l18 = NUVmag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+def FUVmags_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on GALEX FUV of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_galex.txt',comments='#')
+ MH,logAge,Mini,Te,logL, FUVmag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,6]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+ l1 = FUVmag[k+1]
+ l2 = FUVmag[k+2]
+ l3 = FUVmag[k+3]
+ l4 = FUVmag[k+4]
+ l5 = FUVmag[k+5]
+ l6 = FUVmag[k+6]
+ l7 = FUVmag[k+7]
+ l8 = FUVmag[k+8]
+ l9 = FUVmag[k+9]
+ l10 = FUVmag[k+10]
+ l11 = FUVmag[k+11]
+ l12 = FUVmag[k+12]
+ l13 = FUVmag[k+13]
+ l14 = FUVmag[k+14]
+ l15 = FUVmag[k+15]
+ l16 = FUVmag[k+16]
+ l17 = FUVmag[k+17]
+ l18 = FUVmag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+#####################################
+########## LSST FUNCTIONS
+#####################################
+
+def ulsst_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on LSST-U of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_lsst.txt',comments='#')
+ MH,logAge,Mini,Te,logL, umag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,5]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+ l1 = umag[k+1]
+ l2 = umag[k+2]
+ l3 = umag[k+3]
+ l4 = umag[k+4]
+ l5 = umag[k+5]
+ l6 = umag[k+6]
+ l7 = umag[k+7]
+ l8 = umag[k+8]
+ l9 = umag[k+9]
+ l10 = umag[k+10]
+ l11 = umag[k+11]
+ l12 = umag[k+12]
+ l13 = umag[k+13]
+ l14 = umag[k+14]
+ l15 = umag[k+15]
+ l16 = umag[k+16]
+ l17 = umag[k+17]
+ l18 = umag[k+18]
+
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+def glsst_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on LSST-G of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_lsst.txt',comments='#')
+ MH,logAge,Mini,Te,logL, gmag=Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,6]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+ l1 = gmag[k+1]
+ l2 = gmag[k+2]
+ l3 = gmag[k+3]
+ l4 = gmag[k+4]
+ l5 = gmag[k+5]
+ l6 = gmag[k+6]
+ l7 = gmag[k+7]
+ l8 = gmag[k+8]
+ l9 = gmag[k+9]
+ l10 = gmag[k+10]
+ l11 = gmag[k+11]
+ l12 = gmag[k+12]
+ l13 = gmag[k+13]
+ l14 = gmag[k+14]
+ l15 = gmag[k+15]
+ l16 = gmag[k+16]
+ l17 = gmag[k+17]
+ l18 = gmag[k+18]
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+
+
+
+#####################################
+########## DECAM FUNCTIONS
+#####################################
+
+
+
+
+def g_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on DECAM-g of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_decam.txt',comments='#')
+ MH,logAge,Mini,Te,logL, gmag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,5]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = gmag[k+1]
+ l2 = gmag[k+2]
+ l3 = gmag[k+3]
+ l4 = gmag[k+4]
+ l5 = gmag[k+5]
+ l6 = gmag[k+6]
+ l7 = gmag[k+7]
+ l8 = gmag[k+8]
+ l9 = gmag[k+9]
+ l10 = gmag[k+10]
+ l11 = gmag[k+11]
+ l12 = gmag[k+12]
+ l13 = gmag[k+13]
+ l14 = gmag[k+14]
+ l15 = gmag[k+15]
+ l16 = gmag[k+16]
+ l17 = gmag[k+17]
+ l18 = gmag[k+18]
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+def r_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on DECAM-r of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_decam.txt',comments='#')
+ MH,logAge,Mini,Te,logL, rmags = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,7]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = rmags[k+1]
+ l2 = rmags[k+2]
+ l3 = rmags[k+3]
+ l4 = rmags[k+4]
+ l5 = rmags[k+5]
+ l6 = rmags[k+6]
+ l7 = rmags[k+7]
+ l8 = rmags[k+8]
+ l9 = rmags[k+9]
+ l10 = rmags[k+10]
+ l11 = rmags[k+11]
+ l12 = rmags[k+12]
+ l13 = rmags[k+13]
+ l14 = rmags[k+14]
+ l15 = rmags[k+15]
+ l16 = rmags[k+16]
+ l17 = rmags[k+17]
+ l18 = rmags[k+18]
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+
+def z_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on DECAM-z of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_decam.txt',comments='#')
+ MH,logAge,Mini,Te,logL, zmag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,6]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = zmag[k+1]
+ l2 = zmag[k+2]
+ l3 = zmag[k+3]
+ l4 = zmag[k+4]
+ l5 = zmag[k+5]
+ l6 = zmag[k+6]
+ l7 = zmag[k+7]
+ l8 = zmag[k+8]
+ l9 = zmag[k+9]
+ l10 = zmag[k+10]
+ l11 = zmag[k+11]
+ l12 = zmag[k+12]
+ l13 = zmag[k+13]
+ l14 = zmag[k+14]
+ l15 = zmag[k+15]
+ l16 = zmag[k+16]
+ l17 = zmag[k+17]
+ l18 = zmag[k+18]
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+def i_fun(mass, age, metals):
+ '''''''''''''''
+ This function calculates the magnitude on DECAM-i of stars from the synthetic CMD:
+ Parameters
+ -------
+ mass: numpy array.
+ Mass of the particle.
+ age: numpy array.
+ Age of the particle
+ metals: numpy array.
+ '''''''''''''''
+ Isoc=np.genfromtxt('./Isochrones/new_CMDs/Isoc_new_decam.txt',comments='#')
+ MH,logAge,Mini,Te,logL, imag = Isoc[:,0],Isoc[:,1],Isoc[:,2],Isoc[:,3],Isoc[:,4], Isoc[:,8]
+ # Difererentes edades de las estrellas del cmd
+
+ logAge_iter=[0]
+ logAge_iter[0]=logAge[0]
+ j=0
+ for i in range(0,len(logAge)-1):
+ if (logAge_iter[j] != logAge[i]) & (logAge_iter[j] <= logAge[i]):
+ j=j+1
+ logAge_iter.append(logAge[i])
+
+ # Difererentes metalicidades de las estrellas del cmd
+ MH_iter=[0]
+ MH_iter[0]=MH[0]
+ j=0
+ for i in range(0,len(MH)-1):
+ if MH_iter[j] != MH[i]:
+ j=j+1
+ MH_iter.append(MH[i])
+
+ met_l = np.unique(MH)
+ age_l = np.unique(logAge)
+ Mass_l = np.unique(Mini)
+ met_len=len(met_l)-1
+ age_len=len(age_l)-1
+ mass_len=len(Mass_l)-1
+
+
+ m=((age_len)*((age*1e9-10**np.min(age_l))/(10**np.max(age_l)-10**np.min(age_l)))).round(0).astype(int)
+ m[m<0]=0
+ m[m>=age_len]=age_len
+
+
+ j=((met_len)*(metals-np.min(met_l))/(np.max(met_l)-np.min(met_l))).astype(int)
+
+ j[j<0]=0
+ j[j>=met_len]=met_len
+
+ k=(j*(age_len+1)*(mass_len+1)+m*(mass_len+1))
+
+
+ l1 = imag[k+1]
+ l2 = imag[k+2]
+ l3 = imag[k+3]
+ l4 = imag[k+4]
+ l5 = imag[k+5]
+ l6 = imag[k+6]
+ l7 = imag[k+7]
+ l8 = imag[k+8]
+ l9 = imag[k+9]
+ l10 = imag[k+10]
+ l11 = imag[k+11]
+ l12 = imag[k+12]
+ l13 = imag[k+13]
+ l14 = imag[k+14]
+ l15 = imag[k+15]
+ l16 = imag[k+16]
+ l17 = imag[k+17]
+ l18 = imag[k+18]
+
+ return np.array([l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,l11,l12, l13,l14,l15,l16,l17,l18])
+
+
+

Event Timeline

c4science · Help