From 6e9bfca961dfa497742ee9fdd433e5a55609171c Mon Sep 17 00:00:00 2001
From: Tatiana Acero Cuellar <cuellart@jupyterMiLAB>
Date: Sun, 2 May 2021 16:49:32 -0500
Subject: [PATCH] add some functions in .py files
---
Functions/constants.py | 19 ++++
Functions/functions.py | 226 +++++++++++++++++++++++++++++++++++++++++
Functions/imports.py | 6 ++
3 files changed, 251 insertions(+)
create mode 100644 Functions/constants.py
create mode 100644 Functions/functions.py
create mode 100644 Functions/imports.py
diff --git a/Functions/constants.py b/Functions/constants.py
new file mode 100644
index 0000000..fb58e3b
--- /dev/null
+++ b/Functions/constants.py
@@ -0,0 +1,19 @@
+# -- traer los imports del archivo "imports.py"
+from imports import *
+
+
+# -- constante H0, ecuación 7 del paper y pie de nota en página 4 del paper
+H0 = 70
+# -- H0 esta en km/(s Mpc) convertirlo a km/(s pc)
+H0 = ((H0 * units.km / (units.s * units.Mpc)).to('km / s pc')).value
+G = ((constants.G * constants.M_sun / constants.pc).to('km2 /s2') * units.pc /units.M_sun).value
+
+# -- parametros iniciales para el ajuste (leastsquare) para el modelo NFW
+# -- p0: son los valores de la galaxia simulada Dwarf1 del paper
+M200 = (3e10 * units.M_sun).value
+c = 11
+
+# -- parametros iniciales para el ajuste (leastsquare) para el modelo ISO
+# -- p0: tienen "inspiracion" de la gráfica derecha de la tabla 14 (pag 20) del paper
+Rc = ((10.78 * units.kpc).to('pc')).value
+rho0 = (np.exp(15.278)*1e-1 * units.M_sun / units.pc **3 ).value
diff --git a/Functions/functions.py b/Functions/functions.py
new file mode 100644
index 0000000..271950e
--- /dev/null
+++ b/Functions/functions.py
@@ -0,0 +1,226 @@
+# -- trae las constantes del archivo "constants.py"
+from constants import *
+
+def read_col_names(path):
+ """
+ Recibe el path donde esta la galaxia, ej: "Data/rcugc7876.txt".
+ El formato del archivo que se lee es:
+ primera fila: nombre de las columnas. ej: r(kpc), v(km/s), etc
+ segunda fila y hasta abajo: datos.
+
+ La función lee desde la segunda fila (los datos numéricos), y aparte lee la primera fila,
+ donde estan los nombres de las columnas, le quita los espacios a estos nombres.
+ ej: nombre original "r (kpc)" --> nombre posterior "r(kpc)".
+
+
+ Retorna los datos como un dataframe.
+ """
+
+ # -- read the first row, that contains the name of each column, but with spaces
+ header_raw = pd.read_csv(path,nrows=1,header=None)
+ # -- convert the list with all names, into a numpy array
+ header_raw = header_raw.iloc[0].to_numpy()
+ # -- split the str array, if a double space appears
+ h = np.array(header_raw[0].split(" "))
+ # -- extract unique names and the index of these uniques in the original array
+ header, index = np.unique(h,return_index=True)
+ # -- extract the unique header, in the correct order
+ header = h[sorted(index)]
+ # -- not consider space a name column
+ header = header[1:]
+ # -- convert the names into dataframe
+ dfh = pd.DataFrame(list(header), columns = ["header"])
+ # -- delete space inside each name
+ dfh["header"] = dfh["header"].str.replace(" ","")
+ columns = dfh["header"].to_numpy()
+
+ return columns
+
+def vnfw (MM, c, x, G = G, H0 = H0):
+ '''
+ Velocity for NFW model.
+ Recibe como parametros el M_200 y la concentracion c. G y H0 son valores dados y fijos.
+ Recibe la variables x como x = r/r200, el r200 depende de M_200 y se define abajo como
+ otra función.
+
+ La M200 debe estar en unidades de masa solar.
+ x = r/r200 donde r y r200 deben estar en unidades parsec (pc)
+ La concentración no tiene unidades
+
+ Retorna el valor de la velocidad en km/s para los parámetros de M200 y c dados.
+ '''
+ v = (10*G*H0*MM)**(1/3) * ((np.log(1+c*x) - c*x/(1+c*x)) /(x*(np.log(1+c) - c/(1+c))) )**(1/2)
+ return v
+
+def viso (rr, x, G = G):
+ '''
+ Velocity for ISO model.
+ Recibe como parámetros el rho0 y Rc. El "rr" dentro de "v = ..." es rr = (rho0)(Rc^2).
+ Recibe la variable x como x = r/Rc.
+
+ rho0 esta dado en unidades de Msun / (pc^3)
+ x = r/Rc donde r y Rc esta dado en parsec.
+
+
+ Retorna el valor de la velocidad en km/s para los parámetros de rho0 y Rc dados.
+ '''
+ v = ( 4*np.pi*G*rr*x**(-1) *(x - np.arctan(x) ) )**(1/2)
+ return v
+
+def funcline(params,x,i = 0):
+ """
+ Recibe los parametros para cada modelo (NFW o ISO), y valores del eje x
+ i = 0: si es modelo NFW
+ i = 1: si es modelo ISO
+
+ Por default esta NFW.
+
+ Retorna la función del modelo (NFW o ISO) al cual se le quiere hacer ajuste por
+ least square.
+ """
+ if i == 0:
+ r200 = ((G * params[0] / (10*G*H0*params[0])**(2/3))*units.pc).value
+ function = vnfw(params[0], params[1], x/r200, G = G, H0 = H0)
+ else:
+ function = viso (params[0]*params[1]**2, x/params[1], G = G)
+ return function
+
+def Errorlineal(tpl,x,y,incer=1,i=0):
+ """
+ Recibe los parámetros (tpl), los valores de x, para calcular la velocidad "teórica" de acuerdo a los parámetros
+ y los valores de la velocidad observada (y).
+ Recibe el valor de la incertidumbre (en este caso dado por la columna "err_v(km/s)"),
+ por defaul esta en 1 por si no queremos considerar o no hay errores.
+ Recibe i = 0 (NFW) o i = 1 (ISO)
+ """
+ return (funcline(tpl,x,i)-y)/ np.sqrt(abs(incer))
+
+def calculator(galaxy,p0,Errorlineal,incer=1,i=0):
+ """
+ Recibe la dataframe de los datos de la galaxia (galaxy)
+ Recibe los parametros (p0) de partida para least square
+ Recibe la función de Errorlineal.
+ Recibe el valor de la incertidumbre (en este caso dado por la columna "err_v(km/s)"),
+ por defaul esta en 1 por si no queremos considerar o no hay errores.
+ Recibe i = 0 (NFW) o i = 1 (ISO)
+
+ Retorna el mejor fit a los datos observados mediante la función de least square.
+ Retorna la "matriz de covarianza" **no es exactamente esa matriz
+ """
+
+ x = galaxy["r(kpc)"].to_numpy()*1000
+ best,covx, infodict, errmsg, success= leastsq(Errorlineal, p0,
+ args=(x,galaxy["v(km/s)"].to_numpy(),incer,i),
+ full_output=1, epsfcn=0.0001)
+ return best,x,covx
+
+
+def uncertainty(best,x,covx, galaxy, p0,i=0):
+
+ """
+ [1](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html)
+ To obtain the covariance matrix of the parameters x, cov_x must be multiplied by the
+ variance of the residuals .
+
+
+ [2](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)
+ The estimated covariance of popt. The diagonals provide the variance of the parameter
+ estimate. To compute one standard deviation errors on the parameters use
+ perr = np.sqrt(np.diag(pcov)).
+
+ Recibe el mejor fit a los datos, la matriz de covarianza, el dataframe de la galaxia a trabajar,
+ los parámetros inciales y i = 0 (NFW) o i = 1 (ISO)
+
+ Retorna una lista del mejor parámetro y su respectiva incertidumbre formal
+ ej: M200 = (10e10, 1.0e9)
+
+ """
+
+ if i == 0:
+ residual_err = (Errorlineal(best,x,galaxy["v(km/s)"].to_numpy(),
+ incer=galaxy["err_v(km/s)"].to_numpy(),i=0)**2).sum() / (len(galaxy["v(km/s)"]) - len(p0))
+ else:
+ residual_err = (Errorlineal(best,x,galaxy["v(km/s)"].to_numpy(),
+ incer=galaxy["err_v(km/s)"].to_numpy(),i=1)**2).sum() / (len(galaxy["v(km/s)"]) - len(p0))
+
+ covm = covx*residual_err
+ err = np.sqrt(np.diag(abs(covm)))
+
+ return list(zip(best,err))
+
+def determine(gal,label):
+ """
+ Recibe el dataframe de la galaxia a trabajar, puede ser considerando todos
+ los puntos RA (receiding - approaching) o solo R o solo A y
+ el label = ["ra","r","a"].
+
+ Hace el fit para modelo NFW ( i = 0) y para ISO ( i = 1) y para cada uno de
+ estos, considera es "ra", "r" o "a"
+
+ Retorna una lista indicando si se hizo para NFW (0) o ISO (1), el mejor parámetro,
+ su respectiva incertidumbre formal y si es "ra", "r" o "a"
+
+ ej: Modelo NFW, parametro M200
+ = ("0", 10e10, 1.0e9, "ra",),("0", 11e10, 2.0e9,"r"),("0", 10.5e10, 1.3e9,"a")
+ """
+ AA = []
+ for ii in [0,1]:
+ if ii == 0:
+ p0 = [M200,c]
+ else:
+ p0 = [rho0,Rc]
+
+ best,x,covx = calculator(gal, p0, Errorlineal, gal["err_v(km/s)"].to_numpy(),ii)
+ try:
+ unce = uncertainty(best,x,covx,gal,p0,ii)
+ except:
+ unce = list(zip(best,np.zeros(len(best))))
+ pass
+
+
+ AA.append((ii,unce,label))
+ return AA
+
+def return_values(AA):
+ """
+ Recibe la lista AA dada por la función determine(gal,label) y la desgloza para construir un data frame
+ donde se reporte el parámetro del mejor fit con su respectivo error formal y para cada modelo (NFW o ISO)
+
+ Retorna dataframe con el
+ * nombre de la galaxia,
+ * side: "ra", "r" o "a"
+ * M200_NFW(Msun)
+ * err_M200_NFW(Msun)
+ * c_NFW
+ * err_c_NFW
+ * rho0_ISO(Msun/pc3)
+ * err_rho0_ISO(Msun/pc3)
+ * Rc_ISO(pc)
+ * err_Rc_ISO(pc)
+ """
+
+ A = np.array(AA)
+ ID = A[:,0]
+ RA =[i[0][-1] for i in A[:,1]]
+ M200_NFW = [i[0][1][0][0] for i in A[:,1]]
+ err_M200_NFW = [i[0][1][0][1] for i in A[:,1]]
+ c_NFW = [i[0][1][1][0] for i in A[:,1]]
+ err_c_NFW = [i[0][1][1][1] for i in A[:,1]]
+ rho0_ISO = [i[1][1][0][0] for i in A[:,1]]
+ err_rho0_ISO = [i[1][1][0][1] for i in A[:,1]]
+ Rc_ISO = [i[1][1][1][0] for i in A[:,1]]
+ err_Rc_ISO = [i[1][1][1][1] for i in A[:,1]]
+
+ datas = pd.DataFrame({"ID":ID,"Side":RA,
+ "M200_NFW(Msun)":M200_NFW,"err_M200_NFW(Msun)":err_M200_NFW,
+ "c_NFW":c_NFW,"err_c_NFW":err_c_NFW,
+ "rho0_ISO(Msun/pc3)":rho0_ISO,"err_rho0_ISO(Msun/pc3)":err_rho0_ISO,
+ "Rc_ISO(pc)":Rc_ISO,"err_Rc_ISO(pc)":err_Rc_ISO})
+
+ return datas
+
+# -- esta funcioncita cálcula el r200 de acuerdo a la ecuación 1.14 de la tesis, donde M200 (pm) viene
+# -- dado por la tabla de datos del paper (M_halo) y V200 viene dados por la ecuación 7 del paper
+# -- y pie de nota en página 4 del paper
+r200 = lambda pm : ((G * pm / (10*G*H0*pm)**(2/3))*units.pc).value
+
diff --git a/Functions/imports.py b/Functions/imports.py
new file mode 100644
index 0000000..6a98d7c
--- /dev/null
+++ b/Functions/imports.py
@@ -0,0 +1,6 @@
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+from scipy.optimize import leastsq
+from astropy import units
+from astropy import constants
--
GitLab