||
||	a simple test function and its period,
||	for testing periodic function integration.
||
testf :: num -> num
  testf(x) = cos (cos x) ^ 5.78653
             ||where
             ||tan x = sin x / cos x

period_testf :: num
  period_testf = 2*pi



||
||	fundamental physical constants, and popular light frequencies.
||	Note: the colour names are somewhat arbitrary!
||
c = 2.99792458e8
h = 6.626e-34
k = 1.38e-23

redfreq    = c / 700e-9
yellowfreq = (2*redfreq + bluefreq) / 3
greenfreq  = (redfreq + 2*bluefreq) / 3
bluefreq   = c / 400e-9



||
||	the black-body photon count distribution against frequency.
||	(Input is frequency in units of kT/h; arbitrary output units.)
||
bbphotons :: num -> num
  bbphotons(f) = f^2 / (exp f - 1),                                  if f>0
               = 0,                                                  if f=0
               = error "bbphotons: frequency must be non-negative!", if f<0



||
||	the black-body energy distribution against frequency.
||	(Input is frequency in units of kT/h; arbitrary output units.)
||
bbenergy :: num -> num
  bbenergy(f)  = f^3 / (exp f - 1),                                  if f>0
               = 0,                                                  if f=0
               = error "bbenergy: frequency must be non-negative!",  if f<0



||
||	fraction of energy delivered as visible light,
||	for a given temperature
||
frac_visible :: num -> num
  frac_visible temp = integrate4 10000 redfreqratio bbenergy bluefreqratio /
                      whole_spectrum

                      where

                      tempfreq = k*temp / h
                      redfreqratio = redfreq / tempfreq
                      bluefreqratio = bluefreq / tempfreq

                      whole_spectrum = integrate4 10000 0 bbenergy 40



||
||	colour ratio for a given temperature, defined somewhat arbitrarily
||	as the ratio of delivered energy in the "top third" to "bottom third"
||	of the visible spectrum
||
colour_ratio :: num -> num
  colour_ratio temp = integrate4 10000 greenfreqratio bbenergy bluefreqratio /
                      integrate4 10000 redfreqratio bbenergy yellowfreqratio

                      where

                      tempfreq = k*temp / h
                      redfreqratio = redfreq / tempfreq
                      yellowfreqratio = yellowfreq / tempfreq
                      greenfreqratio = greenfreq / tempfreq
                      bluefreqratio = bluefreq / tempfreq



||
||	front end to frac_visible and colour_ratio above,
||	for running them on a list of temperatures
||
run_frac_visible :: [num] -> [(num,num)]
  run_frac_visible temps = zip2 temps (map frac_visible temps)

run_colour_ratio :: [num] -> [(num,num)]
  run_colour_ratio temps = zip2 temps (map colour_ratio temps)



||
||	some test lists of temperatures
||	(the smaller ones are just in case the big one runs out of heap space!)
||
test_temps :: [num]
  test_temps = [1000,2000,2500,3000,4000,5000,5780,6000,6500,7000,7500,8000,8500,9000,10000,15000,20000,50000,100000,1e6,1e7]

test_temps_1 :: [num]
  test_temps_1 = [1000,2000,2500,3000,4000,5000,5780,6000]

test_temps_2 :: [num]
  test_temps_2 = [6000,6500,7000,7500,8000,8500,9000,10000]

test_temps_3 :: [num]
  test_temps_3 = [10000,15000,20000,50000,100000,1e6,1e7]



||
||	average takes a function and computes the weighted average value
||	of that function for a blackbody photon.
||	e.g. the average frequency of a blackbody photon (in units of kT/h)
||	     is average id ;
||	     the average wavelength (in units of hc/kT) is average (1/).
||
average :: (num -> num) -> num
  average f = integrate4 10000 0.0001 weighted_f 40.0001 /
              integrate4 10000 0.0001 bbphotons 40.0001
              where
              weighted_f(freq) = bbphotons(freq) * f(freq)



exact_integral :: num -> num
||
||	well, OK, not *exact* as such since we chop at 10000, but:
||	the integral (from freq = zero to infinity) of
||	freq^d / (exp freq - 1), by using 1 + r + r^2 + ... = 1/(1-r),
||	comes to d! * [1^-(d+1) + 2^-(d+1) + 3^-(d+1) + ...].
||

exact_integral d = product [1..d] * sum (map (^-(d+1)) [1..10000])



||
||	OK, now the numerical integration stuff itself!
||



||
||	non-periodic: order 2, 3, 4 respectively.
||



integrate2 :: num -> num -> (num -> num) -> (num -> num)
||
||	no. of points -> start -> FUNCTION -> (finish -> value_of_INTEGRAL)
||

integrate2 n a f b = (sum (map f (samples n a b))
                      - (
                         1/2*f(a) + 1/2*f(b)
                        )
                     ) * h
                     where h = (b-a)/(n-1)



integrate3 :: num -> num -> (num -> num) -> (num -> num)
||
||	no. of points -> start -> FUNCTION -> (finish -> value_of_INTEGRAL)
||

integrate3 n a f b = (sum (map f (samples n a b))
                      - (
                         7/12*f(a) - 1/12*f(a+h) - 1/12*f(b-h) + 7/12*f(b)
                        )
                     ) * h
                     where h = (b-a)/(n-1)



integrate4 :: num -> num -> (num -> num) -> (num -> num)
||
||	no. of points -> start -> FUNCTION -> (finish -> value_of_INTEGRAL)
||

integrate4 n a f b
  = (sum (map f (samples n a b))
     - (
          31/48*f(a) - 11/48*f(a+h) + 5/48*f(a+2*h) - 1/48*f(a+3*h)
        + 31/48*f(b) - 11/48*f(b-h) + 5/48*f(b-2*h) - 1/48*f(b-3*h)
       )
    ) * h
    where h = (b-a)/(n-1)



||
||	periodic: order here is meaningless (technically infinite!).
||



integrateP :: (num -> num) -> num -> (num -> num)
||
||	FUNCTION -> period -> no. of points -> value_of_INTEGRAL
||

integrateP f p n = sum (map f (take n [0,p/n..])) * p/n



errorP :: (num -> num) -> num -> [num]
||
||	for smoothly periodic functions only.
||	[not functions made artificially periodic by repetition!]
||
|| FUNCTION -> period -> [error_values]
||

errorP f p = map error_for trials
             where
             exact = integrateP f p 1024	|| OK, fairly exact! :-)
             trials = map (2^) [0..9]
             error_for n = integrateP f p n - exact



chain :: [num] -> [num]
||
||	the successive ratios; for estimating an "order" of convergence
||	(which doesn't always exist: sometimes convergence is *faster* than
||	any given order!!)
||
chain (a:b:rest) = (a/b) : chain (b:rest)
chain [x] = []
chain []  = []



||
||	front end to errorP and chain above.
||
go :: ([num], [num])
  go = (results, chain results)
       where
       results = errorP testf period_testf



||
||	support functions.
||



samples :: num -> num -> num -> [num]
||
||	no. of points -> start -> finish -> [sample_points]
||

samples n a b = take n [a,a+h..]
                where
                h = (b-a)/(n-1)
                ||
                ||	we use "take n" rather than an explicit finish to avoid
                ||	round-off errors losing us an entire point!!
                ||
