import numpy as np
from matplotlib import pyplot as plt
from scipy import integrate

tidVol = 500 #ml
#conversion is   1 millilitre = 0.001225 grams
tidVolkg = tidVol * 0.001225/1000

respTime = 4 #seconds
freq = respTime / (2*np.pi)
time = np.arange(0,2,0.005)

amplitude = 0.0004812
flowLPM = amplitude*np.sin(time/freq)

flow = lambda time: (amplitude*np.sin(time/freq)) #kg/s
a,b=integrate.quad(flow, 0, 2) #take half cycle as it's symmetrical
#inhaled tidal volume is 500ml per breath
print("inhaled tidal volume is", a, "kg", "which is", a*1000, "g")
print("The average flow rate is", a/2, "kg/s")
plt.plot(time,flowLPM)

timeA = time[0:int(len(time)/2)]
constLPM = np.ones(len(timeA))*a/2
plt.plot(timeA,constLPM)
plt.ylabel("Flowrate [L/s]")
plt.xlabel('time [s]')
print("max inhalation is", amplitude*60, "LPM")

# write transient profile -mass FR
fileName = "massFlowRate.prof"
nTerms = len(time)
fid = open(fileName,"w")      # Open the file for writing
fid.write(f"( (sin_2secBreath transient {nTerms} 1)\n") #1 =periodic 0=not periodic
fid.write("(time ")
fid.write(" ".join(str(format(i, ".5f")) for i in time))
fid.write(")\n")
fid.write("(mass-flow ")
fid.write(" ".join(str(format(i, ".12f")) for i in flowLPM))
fid.write(")\n")
fid.write(")\n")
fid.close()