#Script for plotting the comparison of the spectra
#of absorbed Solar energy vs. the spectra of emitted
#energy for the planets of the Solar system
#

#Data on section of text which this script is associated with
Chapter = '3.**'
Figure = '**'
#

import phys,math
from ClimateUtilities import *
from planets import *

#Under construction. After trying a number of things,
#the main point is illustrated best by plotting on a
#log-wavenumber axis, and transforming the density
#accordingly.  It would also be possible to use
#the cumulative emission function, but that might
#be harder to read.  For the publication version,
#only selected planets are shown, to make the graph
#easier to read.

#ToDo: Consider using observed planetary albedos rather
#than a hypothetical planet with zero albedo

#Planet data
Tstar = 5800.
SolarRadius = 695.e6 #Solar radius in m
PlanetList = [Mercury,Venus,Earth,Mars,Jupiter,Saturn,Uranus,Neptune,Pluto]

#Specify list of temperatures.

Tlist = [Tstar*(1.-p.albedo)**.25*math.sqrt(1./(2.*(p.rsm/SolarRadius)))
         for p in PlanetList]
Tlist.reverse() #Put temperature in ascending order
#Append Solar temperature
Tlist.append(Tstar)
#List of names
names = [p.name for p in PlanetList]
names.reverse()
names.append('Sun')

#Generate a good range of frequencies for plotting
#This assumes that the first temperature in the list
#is the minimum, and the last is the maximum
nuA = phys.k*Tlist[0]/phys.h
dnu = nuA/2.
nuB = phys.k*Tlist[-1]/phys.h
N = int(20.*nuB/dnu)

nuList = [i*dnu for i in range(1,N)]
c = Curve()
wavenum = [.01*nu/phys.c for nu in nuList] #Wavenumber in 1/cm
c.addCurve(wavenum,'n')

#Note: All spectra are normalized to have unit total emission
#Note also that we transform the density so that it is
#properly defined when plotted on a log axis.
names.reverse() #To use as pop-off list for labeling
for T in Tlist:
  Blist = [nu*math.pi*phys.B(nu,T)/(phys.sigma*T**4) for nu in nuList]
  label = names.pop()
  print label,T
  c.addCurve(Blist,label,label)
c.XlogAxis = True
w= plot(c)