# Lets define all imports and  matplotlib figure defaults

import numpy as np
import matplotlib.pyplot as plt
from __future__ import division
import math
import matplotlib
%matplotlib inline

matplotlib.rc('figure', figsize=(14, 7))
matplotlib.rc('font', weight='bold', size=12)

# Thermal cyclic analysis of turbojet engine
# define constants
prgcon={"llp":10,"llt":4,"cp":1005,"r":287,"g":9.81,"tlr":6.5,
        "cfuel":45868.0,"prc":2.0,"tet":800,"pa":1.1013,"ta":291.15,
        "masc":0.1289,"effc":0.95,"effcom":0.99,"plcom":0.05,
        "bleed":0.01,"efmech":0.98,"effct":0.84,"effn":0.95,"powreq":1.0
        }


def toyturbojet(prgcon):
    result={"temps":[],"prs":[],"thrust":0.0,"spthrust":0.0,"sfc":0.0}
    prc=prgcon["prc"]
    t03=prgcon["tet"]
    gam=(prgcon["cp"]/prgcon["r"])/(prgcon["cp"]/prgcon["r"]-1)
    gamx=prgcon["cp"]/prgcon["r"]
    gamy=(gam-1.0)/2.0

    macfl=0.0
    r1=1.0
    pa=prgcon["pa"]
    ta=prgcon["ta"]

    p01=pa*r1*((1.0+gamy*macfl**2)**gamx)
    t01=ta*(1.0+gamy*macfl**2)
    den1=(p01*10000.0*prgcon["g"])/(prgcon["r"]*t01)

#print p01,t01,den1
    trais=prc**(1.0/gamx)
    tdrpis=t01*(trais-1.0)
    tdrp=tdrpis/prgcon["effc"]
    powerc=prgcon["masc"]*prgcon["cp"]*tdrp/1000.0
    p02=p01*prc
    t02=t01+tdrp

    mascom=prgcon["masc"]-prgcon["bleed"]*prgcon["masc"]
    tracom=t03-t02

    fara1=1.0
    if(tracom > 000  and tracom < 400):
        fara1=0.99 *((t03-t02-10)*(1.+t02/3250.))/(prgcon["cfuel"]*prgcon["effcom"])
    elif(tracom > 400 and tracom < 9000):
       fara1=1.10*((t03-t02-50.)*(1.+t02/3250.))/(prgcon["cfuel"]*prgcon["effcom"])

        
    masfuel=fara1*mascom
    p03=p02*(1.0-prgcon["plcom"])

    powert=powerc/prgcon["efmech"]
    masct=mascom+masfuel
    cp34=950.0+0.21*t03
    tdrpct=powert/(masct*cp34/1000.0)
    gamxc=cp34/prgcon["r"]
    tdctis=tdrpct/prgcon["effct"]
    t04is=t03-tdctis
    pract=(t03/t04is)**gamxc
    p04=p03/pract
    t04=t03-tdrpct

#print p01,p02,p03,p04,t01,t02,t03,t04

    rthst=prgcon["powreq"]
    p05i=p04
    t05=t04
    cp45=950.0+0.21*t04
    gamxp=cp45/prgcon["r"]
    gamn=gamxp/(gamxp-1.0)
    popan=p05i/pa
    if p05i < 1.1*pa:
        result["temps"]=([t01,t02,t03,t04,t05])
        result["prs"]=([p01,p02,p03,p04,0.0])
        result["thrust"]=0.0
        result["sfc"]=0.0
        result["spthrust"]=0.0
        return result
    emach=math.sqrt(2*(popan**(1/gamxp)-1)/(gamn-1))
    t5i=t05/(1+0.5*(gamn-1)*emach*emach)
    t5=t05-prgcon["effn"]*(t05-t5i)
    vjet=math.sqrt(2*cp45*(t05-t5))
    p05=pa*((t05/t5)**(gamn/(gamn-1)))
    massj=masct+prgcon["bleed"]+prgcon["masc"]
    thst=massj*vjet/prgcon["g"]
    ssp=thst/prgcon["masc"]
    sfc=masfuel/thst*3600.

    result["temps"]=([t01,t02,t03,t04,t05])
    result["prs"]=([p01,p02,p03,p04,p05])
    result["thrust"]=thst
    result["sfc"]=sfc
    result["spthrust"]=ssp

    return result

#print p01,p02,p03,p04,p05,t01,t02,t03,t04,t05,thst,ssp,sfc


res=toyturbojet(prgcon)
print res
plt.plot(res["prs"])
plt.xticks([0,1,2,3,4])
plt.title("Pressure at each stage ")
plt.ylabel("Pressure (bar)")




plt.plot(res["temps"])
plt.xticks([0,1,2,3,4])
plt.title("Temperature at each stage ")
plt.ylabel("Teperature (K)")


pr1 = np.linspace(2.0,24.0,20)
tet1=np.linspace(800.0,1400.0,6)
pr= np.zeros((pr1.shape[0],tet1.shape[0]))


#print pr.shape[0]
sp=np.zeros((pr.shape[0],tet1.shape[0]))
sfc=np.zeros((pr.shape[0],tet1.shape[0]))
thr=np.zeros((pr.shape[0],tet1.shape[0]))


for j, tet in enumerate(tet1): 
    prgcon["tet"]=tet
    pr[:,j]=pr1
    for i,pp in enumerate(pr1):
        prgcon["prc"]=pp
        res =toyturbojet(prgcon)
        #print i,pp
        sp[i,j]=res["spthrust"]
        sfc[i,j]=res["sfc"]
        thr[i,j]=res["thrust"]

plt.subplot(3,1,1)
plt.plot(pr,sfc)
plt.title("sfc vs pr ")
plt.ylabel("sfc")
plt.xlabel("pr")

plt.legend="best"
plt.subplot(3,1,2)
plt.plot(pr,sp)
plt.title("Sp. Thrust vs pr ")
plt.ylabel("Sp. Thrust")
plt.xlabel("pr")

plt.subplot(3,1,3)
plt.legend="best"
plt.plot(pr,thr)
plt.title("Thrust vs pr ")
plt.ylabel("Thrust")
plt.xlabel("pr")