from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

import sys, math
import scipy.integrate as integrate
#from scipy import inf
import numpy as num

__author__ = 'becker@astro.washington.edu (Andrew Becker)'

# becker@astro.washington.edu

# Many of these adapted from
# astro-ph/9905116

# README :
# Version 1.1 fixes bug in self.Tfunc (lookback time); thanks to Neil Crighton for catching it!
# Version 1.1.1 added unit conversions and changing cosmology at init (Conor Mancone, cmancone@ufl.edu)
# Version 1.1.2 added a couple extra functions of various utility (Conor Mancone)

VERSION = '1.1.2'


class Cosmology:
    lookup_tol = 1e-8

    # Note, all the distance units come from Dh, so this returns everything in meters

    def __init__(self, Om=0.272, Ol=0.728, h=0.704, w=-1):
        self.Om = Om
        self.Ol = Ol
        self.w = w
        self.h = h

        # CONSTANTS
        self.c = 2.9979E5  # km/s
        self.G = 6.67259E-11  # m^3 / kg / s^2
        self.Msun = 1.98892E30  # kg
        self.pc = 3.085677E16  # m

        # Functions for integration
        # Eqn 14 from Hogg, adding in 'w' to Ol.
        self.Efunc = lambda z, self=self : num.sqrt( (self.Om   * (1. + z)**3 +                 \
                                                       self.Ok() * (1. + z)**2 +                 \
                                                       self.Ol   * (1. + z)**(3 * (1. + self.w)) \
                                                       )**-1 )
        # Eqn 30
        self.Tfunc = lambda z, self=self: self.Efunc(z) / (1. + z)

        # Omega total
    def Otot(self):
        return self.Om + self.Ol

    # Curvature
    def Ok(self):
        return 1. - self.Om - self.Ol

    # Hubble constant, km / s / Mpc
    def H0(self):
        return 100 * self.h

    # Hubble constant at a particular epoch
    # Not sure if this is correct
    #def Hz(self, z):
    #    return self.H0() * (1. + z) * num.sqrt(1 + self.Otot() * z)

    # Got this from Jake
    def Hz(self, z):
        return self.H0() / self.Efunc(z)

    # Scale factor
    def a(self, z):
        return 1. / (1. + z)

    # Hubble distance, c / H0
    # Returns meters
    def Dh(self, meter=False, pc=False, kpc=False, mpc=False, cm=False):
        d = self.c / self.H0()  # km / s / (km / s / Mpc) = Mpc
        d *= self.pc * 1e6  # m
        return d * self.lengthConversion(cm=cm,
                                         meter=meter,
                                         pc=pc,
                                         kpc=kpc,
                                         mpc=mpc)

    # Hubble time, 1 / H0
    # Returns seconds
    def Th(self, s=False, yr=False, myr=False, gyr=False):
        t = 1. / self.H0()  # Mpc s / km
        t *= self.pc * 1e3
        return t * self.timeConversion(s=s, yr=yr, myr=myr, gyr=gyr)

    # Lookback time
    # Difference between the age of the Universe now and the age at z
    def Tl(self, z, s=False, yr=False, myr=False, gyr=False):
        #integ = integrate.romberg(self.Tfunc, 0, z)
        integ, _ = integrate.quad(self.Tfunc, 0, z)
        return self.Th() * integ * self.timeConversion(s=s, yr=yr, myr=myr, gyr=gyr)

    # Line of sight comoving distance
    # Remains constant with epoch if objects are in the Hubble flow
    def Dc(self, z, cm=False, meter=False, pc=False, kpc=False, mpc=False):
        #integ = integrate.romberg(self.Efunc, 0, z)
        integ, _ = integrate.quad(self.Efunc, 0, z)
        return self.Dh() * integ * self.lengthConversion(cm=cm, meter=meter, pc=pc, kpc=kpc, mpc=mpc)

    # Transverse comoving distance
    # At same redshift but separated by angle dtheta; Dm * dtheta is transverse comoving distance
    def Dm(self, z, cm=False, meter=False, pc=False, kpc=False, mpc=False):
        Ok = self.Ok()
        sOk = num.sqrt(num.abs(Ok))
        Dc = self.Dc(z)
        Dh = self.Dh()

        conversion = self.lengthConversion(cm=cm,
                                           meter=meter,
                                           pc=pc,
                                           kpc=kpc,
                                           mpc=mpc)

        if Ok > 0:
            return Dh / sOk * num.sinh(sOk * Dc / Dh) * conversion
        elif Ok == 0:
            return Dc * conversion
        else:
            return Dh / sOk * num.sin(sOk * Dc / Dh) * conversion

    # Angular diameter distance
    # Ratio of an objects physical transvserse size to its angular size in radians
    def Da(self, z, cm=False, meter=False, pc=False, kpc=False, mpc=False):
        return self.Dm(z) / (1. + z) * self.lengthConversion(cm=cm,
                                                             meter=meter,
                                                             pc=pc,
                                                             kpc=kpc,
                                                             mpc=mpc)

    # Angular diameter distance between objects at 2 redshifts
    # Useful for gravitational lensing
    def Da2(self,
            z1,
            z2,
            cm=False,
            meter=False,
            pc=False,
            kpc=False,
            mpc=False):
        # does not work for negative curvature
        assert (self.Ok()) >= 0

        # z1 < z2
        if (z2 < z1):
            foo = z1
            z1 = z2
            z2 = foo
        assert (z1 <= z2)

        Dm1 = self.Dm(z1)
        Dm2 = self.Dm(z2)
        Ok = self.Ok()
        Dh = self.Dh()

        return 1. / (1 + z2) * (
            Dm2 * num.sqrt(1. + Ok * Dm1**2 / Dh**2) - Dm1 *
            num.sqrt(1. + Ok * Dm2**2 / Dh**2)) * self.lengthConversion(
                cm=cm, meter=meter, pc=pc,
                kpc=kpc, mpc=mpc)

    # Luminosity distance
    # Relationship between bolometric flux and bolometric luminosity
    def Dl(self, z, cm=False, meter=False, pc=False, kpc=False, mpc=False):
        return (1. + z) * self.Dm(z) * self.lengthConversion(cm=cm,
                                                             meter=meter,
                                                             pc=pc,
                                                             kpc=kpc,
                                                             mpc=mpc)

    # Distance modulus
    # Recall that Dl is in m
    def DistMod(self, z):
        return 5. * num.log10(self.Dl(z) / self.pc / 10)

    # added 12/20/11 - Conor Mancone
    def Tu(self, s=False, yr=False, myr=False, gyr=False):
        return self.Th() * integrate.quad(
            self.Tfunc, 0, num.inf)[0] * self.timeConversion(s=s,
                                                         yr=yr,
                                                         myr=myr,
                                                         gyr=gyr)

    # added 12/16/11 - Conor Mancone
    # returns conversion from arcseconds to physical angular size
    def scale(self, z, cm=False, meter=False, pc=False, kpc=False, mpc=False):
        return math.tan(1.0 / 3600 * math.pi / 180) * self.Da(z,
                                                              cm=cm,
                                                              meter=meter,
                                                              pc=pc,
                                                              kpc=kpc,
                                                              mpc=mpc)

    # added 12/05/11 - Conor Mancone  Ages should be in years
    def GetZ(self, incoming_ages, zf):

        # deal with incoming scalar values
        is_scalar = False
        if len(num.asarray(incoming_ages).shape) == 0:
            incoming_ages = num.asarray([incoming_ages])
        else:
            incoming_ages = num.asarray(incoming_ages)

        # negative ages are not allowed
        if incoming_ages.min() < 0:
            raise ValueError('Negative ages are not allowed in call to GetZ')

        # also nothing older than the lookback time
        if incoming_ages.max() > self.Tl(zf, yr=True):
            raise ValueError(
                'Ages older than the lookback time of zf are not allowed in call to GetZ')

        # there is no cosmology routine to calculate redshift given formation redshift and age, so we must work the problem backwards.
        # Make a semi-regular grid of zs, calculate age given zf, and then interpolate on that.

        # Keep track of zfs we've looked up and their stored interpolation data
        if not hasattr(self, 'lookup_zfs'):
            self.lookups = []
            self.lookup_zfs = num.array([])

        # if we have looked up this formation redshift before then figure out where it is in the lookups array
        fetched = False
        if len(self.lookup_zfs):
            ind = num.abs(self.lookup_zfs - zf).argmin()
            if num.abs(self.lookup_zfs - zf)[ind] < self.lookup_tol:
                ages = self.lookups[ind]['ages']
                zs = self.lookups[ind]['zs']
                fetched = True

        # if we haven't found the stored index then we have to look this one up now
        if not fetched:

            # make a list of redshifts/ages out to zf.
            # use finer sampling at low redshift.  The details aren't physically motivated - so sue me.
            zs = num.arange(0.0, 0.01, 0.001)
            if zf >= 0.01:
                zs = num.append(zs, num.arange(0.01, num.min([0.1, zf + 1e-5]),
                                               0.005))
            if zf >= 0.1:
                zs = num.append(zs, num.arange(0.1, num.min([2.0, zf + 1e-5]),
                                               0.025))
            if zf >= 2.0: zs = num.append(zs, num.arange(2.0, zf + 1e-5, 0.1))

            # reverse so that age is monotonically increasing
            zs = zs[::-1]

            # now calculate ages for all those redshifts given formation redshift
            ages = num.array(
                [self.Tl(zf, yr=True) - self.Tl(z, yr=True) for z in zs])

            # store reversed so age is monotonically increasing
            self.lookup_zfs = num.append(self.lookup_zfs, zf)
            self.lookups.append({'zs': zs, 'ages': ages})

        # now use the data in self.lookup_zfs[stored_ind] to do the lookup
        # return nan if the age is out of bounds (i.e. z < 0 or z > zf)
        outgoing_zs = num.zeros(incoming_ages.size)
        outgoing_zs[:] = num.nan
        m = (incoming_ages >= ages.min()) & (incoming_ages <= ages.max())
        if m.sum(): outgoing_zs[m] = num.interp(incoming_ages, ages, zs)

        # all done!
        return outgoing_zs

    # added 12/11/09 - Conor Mancone
    def timeConversion(self, s=False, yr=False, myr=False, gyr=False):
        # all times come in as seconds.  Convert accordingly
        if yr: return 1.0 / (3600.0 * 24 * 365.25)
        elif myr: return 1.0 / (3600.0 * 24 * 365.25 * 1e6)
        elif gyr: return 1.0 / (3600.0 * 24 * 365.25 * 1e9)
        else: return 1.0

    # added 12/11/09 - Conor Mancone
    def lengthConversion(self,
                         cm=False,
                         meter=False,
                         pc=False,
                         kpc=False,
                         mpc=False):
        # all lengths come in as meters.  Convert accordingly
        if cm: return 100.0
        elif pc: return 1.0 / self.pc
        elif kpc: return 1.0 / (self.pc * 1e3)
        elif mpc: return 1.0 / (self.pc * 1e6)
        else: return 1.0


if __name__ == '__main__':

    c = Cosmology()

    if len(sys.argv) < 2:
        print('Usage : cosmology.py z1 [z2]')
    z1 = float(sys.argv[1])

    print('Cosmology : H0           =', c.H0())
    print('Cosmology : Omega Matter =', c.Om)
    print('Cosmology : Omega Lambda =', c.Ol)
    print()

    print('Hubble distance                %.2f Mpc' % (c.Dh() / c.pc / 1e6))
    print('Hubble time                    %.2f Gyr' %
          (c.Th() / 3600 / 24 / 365.25 / 1e9))
    print()

    print('For z = %.2f:' % (z1))
    print('Lookback time                  %.2f Gyr' %
          (c.Tl(z1) / 3600 / 24 / 365.25 / 1e9))
    print('Scale Factor a                 %.2f' % (c.a(z1)))
    print('Comoving L.O.S. Distance (w)   %.2f Mpc' % (c.Dc(z1) / c.pc / 1e6))
    print('Angular diameter distance      %.2f Mpc' % (c.Da(z1) / c.pc / 1e6))
    print('Luminosity distance            %.2f Mpc' % (c.Dl(z1) / c.pc / 1e6))
    print('Distance modulus               %.2f mag' % (c.DistMod(z1)))
