#Utility script to make table of critical impact
#erosion parameters



from ClimateUtilities import *
import planets,phys
import math

mEarth = 5.9636e24
mMoon =  7.3477e22

#Impactor mass spectrum function
mp = .1*mMoon #Maximum impactor mass, kg
q = 1.5 #Power law coeffficient
def mTotFact(mc):
    return mp*((1-q)/(2-q))/(1.-(mc/mp)**(1.-q))


def impact(planet,T,ps,gas):
    ell2 = gas.R*T/(planet.g/planet.a)
    mc = ell2* (ps/planet.g)
    #For satellite case, multiply by ratio (v_esc/v_impact)**2
    #(e.g. .1 for Titan orbiting Saturn)
    rSilica = ((.75/math.pi)*mc/3000.)**(1./3.)
    rIce = ((.75/math.pi)*mc/960.)**(1./3.)
    N = 1./(.25*ell2/planet.a**2) #Number of impacts for erosion
    mTot = N*mTotFact(mc)
    #For satellite case multiply mTot by ratio of area of primary
    #to area of satellite (e.g. 500 for Saturn/Titan
    return '%2.1e'%mc,rSilica/1000.,rIce/1000.,N,'%2.1e'%(mTot/mEarth)

#Define a hypothetical Super Earth (a bit like Gliese)
SuperEarth = Dummy()
SuperEarth.g = 16.75
SuperEarth.a = planets.Earth.a*(5.**(1./3.))

cases = [(planets.Earth,280.,1.e5,phys.N2), (planets.Mars,280.,2.e5,phys.CO2),\
         (planets.Mars,220.,1.e4,phys.CO2),(planets.Venus,700.,90.e5,phys.CO2),\
         (planets.Venus,280.,1.e5,phys.N2),(planets.Titan,80.,1.5e5,phys.N2),\
         (SuperEarth,280.,1.e5,phys.N2)]

for case in cases:
    planet,T,ps,gas = case
    print impact(planet,T,ps,gas)