Fix LTO engine calculations

src/IO/read_input.jl

@@ -5,7 +5,7 @@ export read_aircraft_model, load_default_model
     read_input(k::String, dict::AbstractDict=data, 
     default_dict::AbstractDict = default)
 
-Reads the input from a given dictonary (typically parsed from a TOML file).
+Reads the input from a given dictionary (typically parsed from a TOML file).
 If requested input does not exist in dictionary, looks for value in default input
 and stores default value into the given dictionary (primarily for later output/
 saving as an aircraft model file)

src/TASOPT.jl

@@ -17,10 +17,11 @@ using PyPlot
 using Dates
 using ForwardDiff
 using CSV, Tables
-using DocStringExtensions
+using DocStringExtensions, TOML
 
 #convenient directories
-const __TASOPTroot__ = @__DIR__
+const __version__ = TOML.parsefile(joinpath(pkgdir(TASOPT), "Project.toml"))["version"]
+const __TASOPTroot__ = @__DIR__ #points to TASOPT/src for relative imports
 const __TASOPTindices__ = joinpath(__TASOPTroot__,"misc/index.inc") #include(__TASOPTindices__) in REPL
 export __TASOPTroot__, __TASOPTindices__
 

src/mission/LTO.jl

@@ -1,45 +1,120 @@
-function LTO(name, NPSS::Base.Process, ifirst, pare; fileout = stdout)
-nTshaft = parpt[ipt_nTshaft]
-#LTO values
-LTOpoints = [1.0, 0.85, 0.3, 0.07]
-LTOEINOx = zero(4)
-LTOmdotf = zero(4)
-LHV = parg[igLHVfuel]
-# Get Foo based on Tt4 or T/O thrust
-Tt4 = pare[ieTt4, ipstatic]
-Foo = pare[ieFe , ipstatic]
-ip = iptest
-EIs = zeros(4)
-mfs = zeros(4)
-NPSS_success, Foo, η, P, Hrej, heatexcess, 
-mdotf, deNOx, EINOx1, EINOx2, FAR, Tt3, OPR,
-Wc3, Tt41, EGT = NPSS_TEsysOD(NPSS, 0.0, 0.0,
-Foo, 0.0, 0.0, 0.0, 0.0, 0.0, ifirst, parg, parpt, pare, ip)
-println(fileout, "Aircraft LTO characteristics \t\t"*Dates.format(now(), DateFormat("u dd yyyy")))
-println(fileout, name)
-println(fileout, @sprintf("Foo = %.3f [kN]", Foo/1000))
-println(fileout, @sprintf("%2s) %4s %15s %15s %15s %15s %15s %15s %8s %16s", 
-"#", "TS", "Fn[kN]", "̇mf_total[kg/s]", "̇mf[kg/s]", "EI(NOₓ)[g/kg]", "Total_NOₓ[g/s]", "NOₓ[g/s]", "deNOₓ[%]", "JetAeqEI(NOₓ)" ))
-
-    for (i,TS) in enumerate(LTOpoints)
-        println(@sprintf("Running LTO point %.1f %%", TS*100))
-        NPSS_success, Ftotal, η, P, Hrej, heatexcess, 
-        mdotf, deNOx, EINOx1, EINOx2, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TEsysOD(NPSS, 0.0, 0.0, Foo*TS, 0.0, 0.0, 0.0, 0.0, 0.0, ifirst, parg, parpt, pare, ip)
-        ifirst  = false
-
-        println(fileout, 
-        @sprintf("%2d) %4.2f %15.5f %15.5f %15.5f %15.5f %15.5f %15.5f %8.2f %16.5f",
-        i, TS, Ftotal/1000, mdotf, mdotf/nTshaft, EINOx2, EINOx2*mdotf, EINOx2*mdotf/nTshaft, deNOx*100, EINOx2*43/LHV))
-        EIs[i] = EINOx2
-        mfs[i] = mdotf/nTshaft
-    end
-    # for TS in LTOpoints[end:-1:1]
-    #     NPSS_success, Ftotal, η, P, Hrej, heatexcess, 
-    #     mdotf, deNOx, EINOx1, EINOx2, FAR, Tt3, OPR, Wc3, Tt41, EGT = NPSS_TEsysOD(NPSS, 0.0, 0.0, Foo*TS, 0.0, 0.0, 0.0, 0.0, 0.0, ifirst, parg, parpt, pare, ip)
-    # end
-    println(fileout, "------------------------- AEIC ENG_EI output --------------------------")
-    for (i,j) in zip([3,2,1,4], [1,2,3,4])
-        println(fileout,@sprintf("ZIAENG, %d, 0.00, 0.00, %10.5f, 0.00, 1.00,  %10.5f, ZIAENG, 1.0", j, EIs[i], mfs[i]))
-    end
+"""
+    LTO(name, ac; fileout = stdout)
 
+Prints out LTO EI(NOₓ) values
+"""
+function LTO(name, ac; fileout = stdout, method = "cubic")
+
+    #LTO values
+    LTOpoints = [1.0, 0.85, 0.3, 0.07]
+
+    LHV = ac.parg[igLHVfuel]
+    # Get Foo based on Tt4 or T/O thrust
+    Tt4 = ac.pared[ieTt4, ipstatic]
+    Foo = ac.pared[ieFe , ipstatic]
+
+    ip = iptest
+
+    ac.pared[:, ip] = ac.pared[:, ipstatic]
+    ac.parad[:, ip] = ac.parad[:, ipstatic]
+
+    EIs = zeros(4)
+    mfs = zeros(4)
+
+    icall = 2
+    icool = 1
+    initeng = 0
+
+    tfcalc!(ac.pari, ac.parg, view(ac.parad, :, ip), 
+                                view(ac.pared, :, ip), ip, icall, icool, initeng)
+
+
+    println(fileout, "Aircraft LTO characteristics \t\t"*Dates.format(now(), DateFormat("u dd yyyy")))
+    println(fileout, name)
+    println(fileout, @sprintf("Foo = %.3f [kN]", Foo/1000))
+    println(fileout, @sprintf("%2s) %4s %15s %15s %15s %15s", 
+    "#", "TS", "Fn[kN]", "̇mf[kg/s]", "EI(NOₓ)[g/kg]",  "NOₓ[g/s]" ))
+
+        for (i,TS) in enumerate(LTOpoints)
+            # println(@sprintf("Running LTO point %.1f %%", TS*100))
+            ac.pared[ieFe, iptest] = Foo*TS
+            tfcalc!(ac.pari, ac.parg, view(ac.parad, :, iptest), 
+                                    view(ac.pared, :, iptest), 
+                                    ip, icall, icool, initeng)
+
+            F_total = ac.pared[ieFe, iptest]
+            mdotf = ac.pared[ieff, iptest] * ac.pared[iemcore]
+            P3_kPa = ac.pared[iept3, iptest]/1000.0
+            T3_K   = ac.pared[ieTt3, iptest]
+            EI = EINOx(ac, iptest; method=method)
+            EIs[i] = EI
+            mfs[i] = mdotf
+
+            println(fileout, 
+            @sprintf("%2d) %4.2f %15.5f %15.5f %15.5f %15.5f",
+            i, TS, F_total/1000, mdotf, EI, EI*mdotf))
+
+        end
+
+        println(fileout, "------------------------- AEIC ENG_EI output --------------------------")
+        for (i,j) in zip([3,2,1,4], [1,2,3,4])
+            println(fileout,@sprintf("%s, %d, 0.00, 0.00, %10.5f, 0.00, 1.00,  %10.5f, %s, 1.0",
+                                        name, j, EIs[i], mfs[i], name))
+        end
+
+        return EIs, mfs
+end
+
+"""
+    EINOx3(P3_kPa,T3_K, sp_humidity = 0.00634)
+
+
+Returns the EI(NOₓ) for a CFM56 level engine.
+Uses a third order polynomial in T₃ 
+Assumes a default specific humidity of 0.00634 kg water/kg dry air per 
+ICAO Annex 16 Vol. II (part 2.1.4.1)
+"""
+function EINOx3(P3_kPa, T3_K, sp_humidity = 0.00634)
+    # Constants derived using a CRN model for a CFM56 tech level engine
+    a = 6.25528852e-08
+    b = -1.17064467e-04
+    c = 7.36953400e-02
+    d = -1.50392850e+01
+
+    H = -19.0*(sp_humidity - 0.00634)
+
+    EI_NOx = exp(H)*P3_kPa^0.4 * (a*T3_K^3 + b*T3_K^2 + c*T3_K + d)
+
+    return EI_NOx
+end
+"""
+    EINOx4(P3_kPa, T3_K, sp_humidity = 0.00634)
+
+Returns the EI(NOₓ) for a CFM56 level engine.
+Uses a fourth order polynomial which behaves a little better at low T₃ than the cubic
+Assumes a default specific humidity of 0.00634 kg water/kg dry air per 
+ICAO Annex 16 Vol. II (part 2.1.4.1)
+"""
+function EINOx4(P3_kPa, T3_K, sp_humidity = 0.00634)
+    a = 4.85354237e-11
+    b = -6.51089333e-08
+    c =  7.19366066e-06
+    d = 2.06850617e-02
+    e =  -6.69412110e+00
+
+    H = -19.0*(sp_humidity - 0.00634)
+
+    EI_NOx = exp(H)*P3_kPa^0.4 * (a*T3_K^4 + b*T3_K^3 + c*T3_K^2 + d*T3_K + e)
+    return EI_NOx
+end  # function EINOx4
+
+function EINOx(ac::aircraft, ip::Int; sp_humidity = 0.00634, method="cubic")
+       P3_kPa = ac.pared[iept3, ip]/1000.0
+       T3_K   = ac.pared[ieTt3, ip]
+       if lowercase(method) == "cubic"
+            EI = EINOx3(P3_kPa, T3_K, sp_humidity)
+       elseif lowercase(method) == "quartic"
+            EI = EINOx4(P3_kPa, T3_K, sp_humidity)
+       end
+       return EI
 end

test/regression_test_wsize.jl

@@ -30,6 +30,23 @@
 
     @test ac.parm[imPFEI] ≈  0.919121257897844
 
+    @testset "LTO" begin
+        EIs, mfs = TASOPT.LTO("Default A/C + Engine", ac)
+        test_EIs = [35.76509639768726, 26.70084308861548, 10.22948188721507, 6.480767538587702]
+        test_mfs = [1.401100172004722, 1.2535426491385873, 0.7015473738995113, 0.4123593208289815]
+        for i in eachindex(test_EIs)
+            @test test_EIs[i] ≈ EIs[i]
+            @test test_mfs[i] ≈ mfs[i]
+        end
+
+        EIs, mfs = TASOPT.LTO("Default A/C + Engine", ac, method = "quartic")
+        test_EIs = [35.72641687463056, 26.27891743413049, 10.095612208446196, 6.180472562714481]
+        for i in eachindex(test_EIs)
+            @test test_EIs[i] ≈ EIs[i]
+        end
+
+    end
+
 end
 
 @testset "Wide sizing" verbose=true begin

utils/modify_index.py

@@ -0,0 +1,60 @@
+import re
+import click
+
+@click.command()
+@click.option('--input_file', '-i', required=True, 
+               help="Path to index file that you want to de-fortranify")
+@click.option('--variable_names', '-vars', required=True, type=str,
+              help="Enter a list of index name to be removed separated by commas.")
+@click.option('--output_file', '-o', default=None, 
+              help="Give output file name")
+def remove_variable_and_renumber_indices(input_file, output_file,
+                                          variable_names):
+    print("Input file:",input_file)
+    if output_file == None:
+        output_file = input_file + "_modified"
+    print("Output file:",output_file)
+    # print(variable_names)
+    variable_names = variable_names.split(',')
+    print("Indices to remove:",variable_names,"\n")
+
+    # Read the contents of the file
+    with open(input_file, 'r') as file:
+        lines = file.readlines()
+
+    for variable_name in variable_names:
+        var_start = variable_name[0:2] #Store this so can only act on vars of this category
+
+        # Find the index of the variable to be removed
+        removed_index = None
+        removed_line_no = None
+        n_lines = len(lines)
+        for i, line in zip(range(0,n_lines), lines):
+            match = re.search(rf'\b{variable_name}\b\s*=\s*(\d+)\s*', line)
+            if match:
+                # print(i,line)
+                removed_index = int(match.group(1))
+                removed_line_no = i
+                print(f"Removed Ind {variable_name} [{removed_index}] from '{var_start}' category")
+                lines[i] = "### Good Riddance ###\n"
+                break
+
+        # Renumber the indices for the remaining variables in that category
+        for i, line in zip(range(removed_line_no,n_lines), lines[removed_line_no:-1]):
+            # Only match variables of the same category
+            match = re.search(rf'\b{var_start}\w+\s*=\s*(\d+)\s*', line)
+            if match:
+                index = int(match.group(1))
+                if index > removed_index:
+                    # print("Index",index)
+                    new_index = index - 1
+                    lines[i] = re.sub(rf'\b{index}\b', str(new_index), line)
+
+    # Write the modified list to a new text file
+    print("\nWriting to file...")
+    with open( output_file, 'w') as file:
+        file.write(''.join(lines))
+    print("Done!")
+
+if __name__ == '__main__':
+    remove_variable_and_renumber_indices()