#Computes condensible equivalent path
#on the moist adiabat

from ClimateUtilities import *
import phys
import planets

condensible = phys.CH4
noncon = phys.N2
g = planets.Titan.g
selfrat = 1. #6.

#-----------
ps = 1.2e5
Ts = 100.
m = phys.MoistAdiabat(condensible,noncon)
p,T,molarCon,q = m(ps,Ts)

#Plot temperature
c = Curve()
c.addCurve(p,'p','Pressure (Pa)')
c.addCurve(T,'T','Temperature (K)')
c.switchXY = c.reverseY = c.YlogAxis = True
plot(c)

#Plot mixing ratio
c1 = Curve()
c1.switchXY = c1.reverseY = c1.YlogAxis = c1.XlogAxis = True
c1.addCurve(p,'p','Pressure (Pa)')
c1.addCurve(molarCon+1.e-30,'molarCon','Molar concentration of condensible')
plot(c1)

#Path computation
pref = 1.e4
pathSelf = [0.]
pathForeign = [0.]
pb = [p[0]]
for i in range(len(p)-1):
	pbar = .5*(p[i]+p[i+1])
	rbar = .5*(molarCon[i]+molarCon[i+1])
	qbar = .5*(q[i]+q[i+1])
	dp = p[i]-p[i+1]
	#Factor of 2 in following is from slant angle term
	pathForeign.append(pathForeign[-1] + 2.*qbar*pbar*(1-rbar)*dp/(pref*g))
	pathSelf.append(pathSelf[-1] + 2.*qbar*pbar*rbar*dp/(pref*g))
	pq = (1.- rbar)*p + selfrat*rbar*p
	pb.append(pbar)
pathSelf = Numeric.array(pathSelf)
pathForeign = Numeric.array(pathForeign)
pathForeign = pathForeign[-1] - pathForeign
pathSelf = pathSelf[-1] - pathSelf 

c2 = Curve()
c2.addCurve(pb[:-1])
c2.addCurve(pathSelf[:-1]+1.e-10,'Self')
c2.addCurve(pathForeign[:-1]+1.e-10,'Foreign')
c2.XlogAxis = c2.YlogAxis = c2.reverseY = c2.switchXY = True
plot(c2)