add some functions in .py files

Functions/constants.py

+# -- 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

Functions/functions.py

+# -- 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
+

Functions/imports.py

+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