#---------------------Description--------------------------------------
#
#**DESCRIPTION 
#
#
#-----------------------------------------------------------------------
#ToDo:
#       *Figure out some way to deal with the small cross section
#        of atomic H, when it's mixed with another gas.
#

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

from ClimateUtilities import *

import planets
import phys
import math

#Define some atomic gases not in phys
Hydrogen = Dummy()
Hydrogen.R = phys.Rstar
Hydrogen.MolecularWeight = 1.
Hydrogen.name = 'Atomic H'
Oxygen = Dummy()
Oxygen.MolecularWeight = 16.
Oxygen.R = phys.Rstar/16.
Oxygen.name = 'Atomic O'
Nitrogen = Dummy()
Nitrogen.name = 'Atomic N'
Nitrogen.R = phys.Rstar/14.
Nitrogen.MolecularWeight = 14.
#

#Data
chi = math.pi*(125.e-12)**2 #Molecular cross-section
mu = 1./(1000.*phys.N_avogadro) #Atomic mass unit (kg.)

ps = 1.e5 #Total surface pressure in Pa
planet = planets.Earth
gas1 =  Nitrogen
gas2 = Hydrogen
eta = 0. #Molar concentration of second gas
Tlower = 300. #Mean temperature of lower atmosphere (homosphere)


#Functions for iteration to find exobase height

def f(r,T):
    #ns = (ps/(gas.R*Tlower))/(mu*gas.MolecularWeight)
    #Ho = gas.R*Tlower/planet.g
    eta = n(r)[1]
    Mbar = ((1.-eta)*gas1.MolecularWeight + eta*gas2.MolecularWeight)
    R = phys.Rstar/Mbar
    g =  planet.g*(planet.a/r)**2
    H = R*T/g
    nex = 1./(math.sqrt(2.)*chi*H)
    #Find the radius at which the exospheric density is attained
    #We can do this analytically using the following statement, for
    #a one-component atmosphere:
    #r1 = planet.a/(1. - (Ho/planet.a)*math.log(ns/nex))
    #but instead we do the following iteration, to allow for more complex
    #multicomponent atmospheres.
    ans1.setParams(nex)
    r1 = ans1(planet.a)
    return r1-r

#Particle density model (use this to generalize code to multicomponent case)
#Re-write this in terms of homopause density and homopause gravity
nH = 1.2e19 #Homopause density (Eventually put in correct scaling)
rH = planet.a+100.e3 #Temporary. Put in right value for smaller planets
gH = planet.g*(planet.a/rH)**2 

#In this version r is distance relative to homopause
def n(r):
    Mbar = ((1.-eta)*gas1.MolecularWeight + eta*gas2.MolecularWeight)
    R = phys.Rstar/Mbar
    H1 = gas1.R*Tlower/gH
    H2 = gas2.R*Tlower/gH
    nn = nH*(1.-eta)*math.exp((rH/H1)*(rH/r - 1.)) \
             + nH*eta*math.exp((rH/H2)*(rH/r - 1.))
    if nn > 0.:
        eta2 = nH*eta*math.exp((rH/H2)*(rH/r - 1.))/nn
    else:
        eta2 = -1.
    return nn,eta2


def fn(r,nex):
    return n(r)[0]- nex

ans = newtSolve(f)
ans1 = newtSolve(fn)

rexGuess = planet.a
#ToDo: Fix this to work properly with multiple constituents
def tLoss(T):
    global rexGuess
    m1 = mu*gas1.MolecularWeight
    m2 = mu*gas2.MolecularWeight
    ans.setParams(T)
    rex = ans(rexGuess)
    if (type(rex) == type(1.)):
        if rex > 0.:
            rexGuess = rex
    g = planet.g*(planet.a/rex)**2
    nex,eta = n(rex)
    lam1 = g*rex*m1/(phys.k*T)
    lam2 = g*rex*m2/(phys.k*T)
    wex1 = math.sqrt(phys.k*T/m1)*.5*math.sqrt(2./math.pi)*(1.+lam1)*math.exp(-lam1)
    wex2 = math.sqrt(phys.k*T/m2)*.5*math.sqrt(2./math.pi)*(1.+lam2)*math.exp(-lam2)
    flux1 = nex*(1.-eta)*wex1
    flux2 = nex*eta*wex2
    Mbar = ((1.-eta)*gas1.MolecularWeight + eta*gas2.MolecularWeight)
    Natmos = (ps/planet.g)/(Mbar*mu)
    if flux1 > 0.:
        t1 = (1.-eta)*Natmos/flux1
        t1 = t1/(24.*3600.*365.*1.e9)
    else:
        t1 = 1.e500

    #Compute escape flux per surface area of second constituent.
    #We do it this way to allow for cases where the atmosphere
    #dissociates above the homopause
    Phi = flux2*(rex/planet.a)**2
    
    return wex1,t1,wex2,'%e'%Phi, (rex-planet.a)/1000.

#Computes the mean free path at infinity (needs to be modified
#to take into account the small collision diameter of atomic hydrogen
def mfpInf():
    ninf = n(1.e30*planet.a)[0]
    return 1./(2.**.5 * chi*ninf)


#Function for computing the sum of the collision radii from
#data on the binary diffusion parameter D*n = b = A*T**s .
#M1 and M2 are the molecular weights of the diffusing substance
#and the background gas. The collision radii are determined by
#fitting b to the formula for a hard-sphere gas. The function
#returns the average of the effective radii of the two species.
def molrad(M1,M2,A,s,T):
	mu = 1./(1000.*phys.N_avogadro)
	b = A*T**s
	f = (3./64.)*(1.*math.pi*phys.k*T*(M1+M2)/(M1*M2*mu))**.5
	d = ((4./math.pi)*(f/b))**.5
	return .5*d/1.e-12


#Functions to compute photon flux (useful for photolysis rate)
def fPhoton(nu):
	return math.pi*phys.B(nu,6000.)/(phys.h*nu)

#lam1 and lam2 are wavelength in nanometers
def photons(planet,lam1,lam2,tolerance = 1.e6):
    rSun = 6.955e8
    nu1 = phys.c/(lam1*1.e-9)
    nu2 = phys.c/(lam2*1.e-9)
    mPhoton = romberg(fPhoton)
    flux = mPhoton([nu1,nu2],None,tolerance)
    return flux*(rSun/planet.rsm)**2

cases = [(planets.Earth,phys.N2,1.e5,300.,300.),(planets.Earth,Oxygen,.2e5,300.,1000.),
         (planets.Venus,phys.CO2,90.e5,500.,300.),(planets.Venus,phys.water,1.e5,400.,300.),
         (planets.Mars,phys.CO2,2.e5,250.,300.),(planets.Mars,Oxygen,1.e5,250.,1000.),
         (planets.Titan,phys.N2,1.5e5,100.,200.),(planets.Titan,Nitrogen,1.5e5,100.,200.),
         (planets.Moon,phys.N2,1.e5,260.,300.)]

#Loop to make table of exosphere properties
##gas2 = Hydrogen
##for case in cases:
##	planet,gas1,ps,Tlower,Tex = case
##	result = tLoss(Tex)
##	wJ = result[0]
##	t = result[1]
##	zex = result[-1]
##	wJH = result[2]
##	print planet.name,gas1.name,ps/1.e5,Tlower,Tex,'%6.0f'%zex,'%2.0e'%t,wJH