#!/usr/bin/python
# Filename: calculate_normalised_sigmaM.py

from __future__ import division

import numpy as np
import pylab
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
#import sim_cosmology as cosmo
import sys

# Basic parameters

mean_density = 27.755e10*0.272#cosmo.rho_mean() # Mean density of the Universe
#h = 0.704
h = 0.697
massUnit = 1.e10
mass_in_8 = 4 / 3 * np.pi * 8**3 * mean_density * massUnit # Mass in 8 Mpc [physical] sphere
#sigma8 = 0.81 # For WMAP7 cosmology
sigma8 = 0.817
mmin = 1.e6
mmax = 10**16

if len(sys.argv) > 1:

    # Load sigmaM files: assuming first row has been removed

    sigma = np.genfromtxt(sys.argv[1], usecols = 1)
    mass = np.genfromtxt(sys.argv[1], usecols=0)
    f = interp1d(mass, sigma)
    norm = f(mass_in_8)
    sigma = sigma / norm * sigma8

    # Plotting tools
    
    
    plt.plot(mass, sigma, linewidth = 2, label = sys.argv[2])
    plt.legend(loc='lower left')
    plt.axvline(x = 10**13.3, color = 'k', linestyle = '--')
    plt.xscale('log')
    pylab.xlim([mmin, mmax])
    pylab.ylim([0,9])

    plt.xlabel(r'$M\,[M_\odot]$')
    plt.ylabel(r'$\sigma\,(M)$')
    


    # Optional: Find the closest match to d_c = 1.686
    
    x = np.linspace(10**11,10**14,100)
    g = interp1d(mass, sigma)
    test_sigmas =  g(x)
    mask = abs(test_sigmas - 1.686).argmin() # Spherical top-hat overdensity
    print(x[mask])

def normed_sigma(data_file):
    
    "Calcualate the sigma(M) normalised to sigma_8 for this cosmology."
    
    sigma = np.genfromtxt(data_file, usecols = 1)
    mass = np.genfromtxt(data_file, usecols=0)
    f = interp1d(mass, sigma)
    norm = f(mass_in_8)
    sigma = sigma / norm * sigma8
    print(mass_in_8)
    return mass, sigma