import Data.Time
getCurrentTime

import IHaskell.Display
import Graphics.Rendering.Chart.Backend.Diagrams

import Graphics.Rendering.Chart.Easy

-- For the purposes of this notebook you do not need to understand this
-- block. It may be useful to point out that in Haskell, comments are
-- introduced by "--" (or bracketed by "{-" and "-}", but I don't plan
-- to use that style here).
--
instance IHaskellDisplay (Renderable a) where
  display renderable = do
    (svg, _) <- renderableToSVG renderable 450 300
    -- now let blaze worry about displaying the SVG
    display svg

-- Calculate f(z; omega_m). The first argument is the matter density
-- and the second the redshift.
f :: Double -> Double -> Double
f om z = let ol = 1 - om
             t = (1 + z)^2 * (1 + om*z) - (2 + z) * ol * z
         in 1 / sqrt t

g = f 0.3

:type g
:type f

g 0

g 1.2

:type map

:type map g

map g [0, 0.4, 3.2]

zs = [0::Double, 0.1 .. 10]
gs = map g zs

toRenderable (do
    layout_title .= "f"
    layout_x_axis . laxis_title .= "Redshift"
    layout_y_axis . laxis_title .= "f(z)"
    plot (line "O_matter=0.3 O_lambda=0.7" [zip zs gs])
    )

import Numeric.Integration.TanhSinh

absolute 1.0e-6 (simpson g 0 2)

absolute 1.0e-6 (trap g 0 2)

result (absolute 1.0e-6 (trap g 0 2))

:type simpson
:type trap

:type absolute
:type relative

take 3 (trap g 0 2)

take 3 (simpson g 0 2)

daH om z = let fs = result (absolute 1.0e-6 (trap (f om) 0 z))
           in fs / (1+z)

:type daH

-- This uses the list of redshifts defined earlier.
-- 
-- (daH om) has type Double -> Double
-- zs has type [Double]
-- so map (daH om) zs has type [Double]
--
-- zip's signature is [a] -> [b] -> [(a,b)]
-- that is, it pairs up corresponding elements of
-- the two input arrays. This means that
-- the output of zip zs (map (daH om) zs)
-- is [(Double, Double)], where the first element
-- of eachpair is the redshift, and the second element
-- the anular-diameter distance for this redshift
-- (in units of the Hubble length).
--
calcDAH om = zip zs (map (daH om) zs)

toRenderable (do
    layout_title .= "Angular-diameter distance"
    layout_x_axis . laxis_title .= "Redshift"
    layout_y_axis . laxis_title .= "dA(z)"
    plot (line "O_m=0.1 O_l=0.9" [calcDAH 0.1])
    plot (line "O_m=0.3 O_l=0.7" [calcDAH 0.3])
    plot (line "O_m=0.9 O_l=0.1" [calcDAH 0.9])
    )

-- Calculate the angular-diameter distance for a flat Cosmology
-- with a matter density of om, Hubble constant of H0 (in km/s/Mpc),
-- and redshift. The result is in Mpc.
--
da om h0 z = daH om z * c / h0
  where
    -- speed of light in km/s
    c = 299792.458

map (da 0.3 70) [0, 0.5, 1.0, 1.2, 2.0]

nedWrightDA = [1259.0, 1651.8, 1710.5, 1726.3]

-- zipWith has a type of (a -> b -> c) -> [a] -> [b] -> [c], and
-- applies the function (first argument) to the elements in the two input
-- lists, creating the output list. Here I supply a list of redshifts
-- and a list of the values from CosmoCalc for those redshifts, and
-- then return the angular-diameter distance calculated for the
-- redshift divided by the CosmoCalc value.
--
-- The (\lambda z nw -> da 0.3 70 z / nw) fragment is Haskell syntax
-- for creating an anonymous function that takes two arguments (z and nw
-- here) and then applying the function body (in this case divide the
-- value of "da 0.3 70 z" by nw).
--
-- I could have written this in several lines, and avoided the need for
-- all these comments!
rel = zipWith (\z nw -> da 0.3 70 z / nw) [0.5, 1.0, 1.2, 2.0] nedWrightDA

rel

calcDA om = zip zs (map (da om 70) zs)

toRenderable (do
    layout_title .= "Angular-diameter distance; H_0 = 70 km/s/Mpc"
    layout_x_axis . laxis_title .= "Redshift"
    layout_y_axis . laxis_title .= "dA(z) [Mpc]"
    plot (line "O_m=0.1 O_l=0.9" [calcDA 0.1])
    plot (line "O_m=0.3 O_l=0.7" [calcDA 0.3])
    plot (line "O_m=0.9 O_l=0.1" [calcDA 0.9])
    )

-- Given Omega_m, H_0 in km/s/Mpc, z, and an angle in arcminutes,
-- return the linear scale in kpc.
adist om h0 z angle = 1000 * d * a
    where
      -- convert angle, in arcminutes, into radians
      a = pi * angle / (180 * 60)
      
      -- d is in Mpc
      d = da om h0 z

adist 0.3 70 2 (1/60)

import qualified Numeric.Units.Dimensional.TF.Prelude as P

c = 299792458.0 P.*~ (P.meter P./ P.second)

:type c

c :: P.Velocity Double
c = 299792458.0 P.*~ (P.meter P./ P.second)

:type c

c

c P./~ (P.kilo P.meter P./ P.second)

c P./ P._2

c P./ (2 P.*~ P.one)

arcmin = (1::Double) P.*~ P.minuteOfArc

:type arcmin

arcmin + c

arcmin

pi / (60 * 180)

arcmin P./~ P.radian == (pi / (60*180))

arcmin == (pi / (60*180))

parsec :: P.Unit P.DLength Double
parsec = P.prefix 3.085678e16 P.metre

3.2 P.*~ P.mega parsec

hubbleConstant :: P.Unit P.DFrequency Double
hubbleConstant = P.kilo P.meter P./ P.second P./ P.mega parsec

wrongConstant :: P.Unit P.DFrequency Double
wrongConstant = P.kilo P.meter P./ P.mega parsec

70 P.*~ hubbleConstant

:type 70 P.*~ hubbleConstant

type HConstant = P.Quantity P.DFrequency Double
type Length = P.Quantity P.DLength Double

da2 :: Double -> HConstant -> Double -> Length 
da2 om h0 z = (daH om z P.*~ P.one) P.* c P./ h0

adist2 om h0 z angle = da2 om h0 z P.* angle

:type adist2

adist3 :: Double -> HConstant -> Double -> P.Quantity P.DPlaneAngle Double -> Length
adist3 = adist2

adist2 0.3 (70 P.*~ hubbleConstant) 2 (1 P.*~ P.meter)

adist3 0.3 (70 P.*~ hubbleConstant) 2 (1 P.*~ P.meter)

adist3 0.3 (70 P.*~ hubbleConstant) 2 (1 P.*~ P.secondOfArc)

adist3 0.3 (70 P.*~ hubbleConstant) 2 (1 P.*~ P.secondOfArc) P./~ P.kilo parsec

adist3 0.3 (70 P.*~ hubbleConstant) 2 (2.4 P.*~ P.minuteOfArc) P./~ P.mega parsec

dlH om z = let fs = result (absolute 1.0e-6 (trap (f om) 0 z))
           in fs * (1+z)
           
dl :: Double -> HConstant -> Double -> Length
dl om h0 z = (dlH om z P.*~ P.one) P.* c P./ h0

dl 0.3 (70 P.*~ hubbleConstant) 2

dl 0.3 (70 P.*~ hubbleConstant) 2 P./~ P.mega parsec

import Numeric.Units.Dimensional.TF.NonSI (year)

lightYear = let d = c P.* (1 P.*~ year) P./~ P.meter
            in P.prefix d P.meter

dl 0.3 (70 P.*~ hubbleConstant) 2 P./~ P.giga lightYear

luminosity :: Double -> HConstant -> Double -> P.Irradiance Double -> P.RadiantIntensity Double
luminosity om h0 z flux = P._4 P.* P.pi P.* d P.* d P.* flux
  where
    d = dl om h0 z

erg = P.prefix (1.0e-7::Double) P.joule

ergLum = erg P./ P.second

ergFlux = ergLum P./ P.square (P.centi P.meter)

lum = luminosity 0.3 (70 P.*~ hubbleConstant) 2 (2e-13 P.*~ ergFlux)

lum

lum P./~ P.watt

lum P./~ ergLum

bigG = 6.67428e-11 P.*~ (P.metre P.^ P.pos3 P./ P.kilo P.gram P./ P.second P.^ P.pos2)

criticalDensity h = P._3 P.* h P.^ P.pos2 P./ (P._8 P.* P.pi P.* bigG)

c70 = criticalDensity (70 P.*~ hubbleConstant)
c70

c70 P./~ (P.kilo P.metricTon P./ P.astronomicalUnit P.^ P.pos3)