#------------------------------------------------------------
#Description:
#
#This script is a modified version of the calculation in
#PureCO2OLR.py, designed as a crude approximation to
#Venusian conditions.  The OLR calculation is set up
#in the same way as in PureCO2OLR.py, but the shortwave
#bands (where HITRAN isn't accurate for Venus) are eliminated
#from the band table and replace by an equivalent blackbody
#emitter at pressure pRad.  The atmospheric temperature is
#assumed to be on the ideal-gas dry adiabat, as before.
#------------------------------------------------------------

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

from ClimateUtilities import *
import phys,planets
import math

import miniClimtFancy as radmodel

#Path to the workbook datasets
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
#Path to the exponential sum tables
EsumPath = datapath + 'Chapter4Data/ExpSumTables/'


radmodel.Tstar = 900. #Temperature scaling parameter
#
#Initialize the radiation model parameters
#
radmodel.BackgroundGas = phys.air #The transparent background gas (not used)
radmodel.GHG =phys.CO2 #The greenhouse gas
radmodel.g = planets.Venus.g
radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableAllWave.260K.100mb.self.data')
#Eliminate transparent shortwave bands, as kludge for high-temperature Venus
radmodel.bandData = radmodel.bandData[:-9] #Uncomment this for "Venus"
radmodel.SelfContinuum = radmodel.CO2SelfContinuum
radmodel.ForeignContinuum = radmodel.CO2ForeignContinuum

#Define adiabat functions
#Set needed thermodynamic constants
stuff = radmodel.GHG
L = stuff.L_sublimation # Latent heat
cp = stuff.cp #Specific heat
R = stuff.R #Gas constant for the atmosphere

def OLR(ps,Tg,pRad,n=20):
    p = radmodel.setpLog(ps,ptop,n)
    q = Numeric.ones(n,Numeric.Float)#1. #Pure greenhouse gas
    T = Tg*(p/ps)**(R/cp)
    TRad = Tg*(pRad/ps)**(R/cp)
    #Handle condensation
    flux,heat = radmodel.radCompLW(p,T,Tg,q)
    #Add in contribution of part of spectrum outside the band table
    nu1 = radmodel.bandData[0].nu1
    nu2 = radmodel.bandData[-1].nu2
    m = romberg(radmodel.Planck) #Evaluates definite integral
    #Instead of adding in the shortwave emission for temperature
    #Tg, we add in the emission for TRad
    flux1 = m([.001,nu1],Tg) + m([nu2,50000.],TRad)
    return flux[0]+flux1

ptop = 1. #Top of atmosphere pressure,Pa.  A global default

#Call the OLR function until we match the absorbed solar for Venus
#You can do this interactively, or use a small Newton's method
#calculation to find the answer.

S = (1.-planets.Venus.albedo)*planets.Venus.L/4.