#!/usr/bin/python
#McDermott
#2-13-2016

from __future__ import division # make floating point division default as in Matlab, e.g., 1/2=0.5
import math
import numpy as np
import scipy.special as sp
import matplotlib.pyplot as plt
from matplotlib import rc
rc('font',**{'family':'sans-serif','sans-serif':['Helvetica'],'size':16})

def tsqramp(t0,t1,dt):
    t = t0
    while t <= t1:
        T = t-t0
        F = ((t-t0)/(t1-t0))**2
        t = t+dt
        f.write("&RAMP ID='fuel ramp', T="+str(T)+", F="+str(F)+"/\n")
    return

HRR_PEAK = 516. # kW
BURNER_AREA = math.pi*(36.*0.0254)**2/4. # circular sand burner area [m]
BL = math.sqrt(BURNER_AREA)
HRRPUA = HRR_PEAK/BURNER_AREA
L = 18.28
W = 12.19
H = 6.10

f = open('FM_SNL_4.fds','w')

f.write("&HEAD CHID='FM_SNL_4', TITLE='ENFP 426 HW2 FDS input file powered by Python'/\n")
f.write("\n")
f.write("&MESH IJK=48,32,16, XB="+format(-L/2,'.3f')+","+format(L/2,'.3f')+
    ","+format(-W/2,'.3f')+","+format(W/2,'.3f')+
    ","+format(0,'.3f')+","+format(H,'.3f')+"/\n")
f.write("\n")
f.write("&TIME T_END=0/\n") # total time is 1210 s
f.write("\n")
f.write("&MISC TMPA=15/\n")
f.write("\n")
f.write("&REAC FUEL='PROPYLENE', SOOT_YIELD=0.02/\n")
f.write("\n")
f.write("&SURF ID='burner', HRRPUA="+format(HRRPUA,'.2f')+", RAMP_Q='fuel ramp', COLOR='RED'/\n")
f.write("\n")
tsqramp(0,240,60)
f.write("&RAMP ID='fuel ramp', T=600, F=1/\n")
f.write("&RAMP ID='fuel ramp', T=601, F=0/\n")
f.write("\n")

CFM=800
VF=CFM*0.0283168/60./6 # volume flow per vent [m3/s]
f.write("&SURF ID='port', VOLUME_FLOW="+format(-VF,'.3e')+", COLOR='LIGHT BLUE'/\n")
f.write("\n")

x0 = 10*12*0.0254
f.write("&OBST XB="+format(x0-BL/2,'.3f')+","+format(x0+BL/2,'.3f')+
    ","+format(-BL/2,'.3f')+","+format(BL/2,'.3f')+
    ",0,.2, THICKEN=.TRUE., SURF_IDS='burner','INERT','INERT'/\n")
f.write("\n")

# ventilation ports (assume 0.3 m x 0.3 m x 1.2 m)
f.write("&MULT ID='m1', DX=6.096,DY=6.096, I_UPPER=2,J_UPPER=1/\n")
f.write("&OBST XB=-6.396,-5.796,-3.3,-2.7,4.9,6.1, MULT_ID='m1', SURF_IDS='INERT','INERT','port'/\n")
f.write("\n")

# exhaust port
f.write("&VENT XB=-9,-8.4,-0.9,0.9,6.1,6.1, SURF_ID='OPEN'/\n")
f.write("\n")

# slice files
f.write("&SLCF PBY=0, QUANTITY='VELOCITY', VECTOR=.TRUE./\n")
f.write("&SLCF PBY=0, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./\n")
f.write("\n")

f.write("\n")
f.write("&TAIL/\n")
f.close()

