text:
@doc """
    COV_var(x::Array{Float64})::Array{Float64}

Determine Coefficient of Variation, COV
    """ ->
function COV_var(x::Array{Float64})::Array{Float64}
    @assert size(x,1) >= 15 "Minimum of 15 readings required to confirm metal loss is general - provide more readings'"
        index = collect(1:1:size(x,1))
        tam = sum(x) * 1/size(x,1)
        # Calculate (Trdi - Tam) and (Trdi - Tam)^2
trdi_tam = zeros(size(x,1))
trdi_tam_2power = zeros(size(x,1))
for i = 1:size(x,1)
  trdi_tam[i] = x[i] - tam
  trdi_tam_2power[i] = (x[i] - tam)^2
end

S = sum(trdi_tam_2power) # sum (Trdi - Tam)^2 sqaured differences
N = size(x,1)

# COV
COV = (1/tam)*((S/(N-1))^0.5)

    return out = [tam, COV]
end # function

@doc """
    part4_PTR_Level_1_Assessment(x::Array{Float64}; annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0)

Average Measured Thickness from Point Thickness Readings (PTR)

Variables\n
equipment_group = "piping" # "vessel", "tank"\n
annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
# "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
# "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]\n
flaw_location = "external" # "External","Internal"\n
units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"\n
tnom = .3 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.\n
FCAml = .05 # Future Corrosion Allowance applied to the region of metal loss.\n
FCA = .05 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).\n
LOSS = 0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.\n
Do = 3.5 # Outside Diameter\n
D = Do - 2*(tnom)  # inside diameter of the shell corrected for FCAml , as applicable\n
P = 1480.0 # internal design pressure.\n
S = 20000.0 # allowable stress.\n
E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.\n
MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.\n
Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .\n
t = tnom - LOSS - FCA  # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.\n
tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).\n

    """ ->
function part4_PTR_Level_1_Assessment(x::Array{Float64}; annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0)
    @assert any(annex2c_tmin_category .== ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
    "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
    "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]) "Invalid entry must be any of the following tmin groups: 'Cylindrical Shell','Spherical Shell','Hemispherical Head','Elliptical Head','Torispherical Head','Conical Shell','Toriconical Head','Conical Transition','Nozzles Connections in Shells',
    'Junction Reinforcement Requirements at Conical Transitions','Tubesheets','Flat head to cylinder connections','Bolted Flanges','Straight Pipes Subject To Internal Pressure','Boiler Tubes','Pipe Bends Subject To Internal Pressure',
    'MAWP for External Pressure','Branch Connections','API 650 Storage Tanks' "
    @assert any(equipment_group .== ["piping", "vessel", "tank"]) "Invalid input: please enter either: 'piping','vessel','tank' "
    @assert any(flaw_location .== ["external", "internal"]) "Invalid input: please enter either: 'external shells', 'internal shells' "
    @assert any(units .== ["lbs-in-psi", "nmm-mm-mpa"]) "Invalid input: please enter either: 'lbs-in-psi', 'nmm-mm-mpa' "
    @assert size(x,1) >= 15 "Minimum of 15 readings required to confirm metal loss is general - provide more readings'"

    print("Begin -- Part 4 LTA Level 1 Assessment (Point Thickness Readings) - API 579 June 2016 Edition\n")

    # STEP 1 – Take the point thickness reading data in accordance with paragraph 4.3.3 Tam , and the
    # Coefficient Of Variation (COV). A template for computing the COV is provided in Table 4.3.
    # Tam = average measured wall thickness of the component based on the point thickness readings (PTR) measured at the time of the inspection.
    # COV = Coefficient Of Variation.
    # Trd = Thickness reading within corroded
    step_1 = COV_var(x)
    tmm = minimum(x)
    tam = step_1[1]
    c_o_v = step_1[2]

    # STEP 2 – If the COV from STEP 1 is less than or equal to 0.1, then proceed to STEP 3 to complete the
    # assessment using the average thickness, tam . If the COV is greater than 0.1 then the use of thickness
    # profiles should be considered for the assessment (see paragraph 4.4.2.2).
    if c_o_v <= 0.1
        print("COV is <= 0.1 - Procced to Step 3 to complete the assessment using the average thickness tam :: COV = ", round(step_1[2],digits =2),"\n")
        step3_satisfied = 1
    else
        print("COV is > 0.1 - Thickness Profiles should be conisdered for assessment :: (see paragraph 4.4.2.2) :: COV = ", round(step_1[2],digits =2),"\n")
        step3_satisfied = 0
    end

    # STEP 3 – The acceptability of the component for continued operation can be established using the Level 1
    # criteria in Table 4.4, Table 4.5, Table 4.6, and Table 4.7. The averaged measured thickness or
    # MAWP acceptance criterion may be used. In either case, the minimum thickness criterion shall be
    # satisfied. For MAWP acceptance criterion (see Part 2, paragraph 2.4.2.2.e) to determine the
    # acceptability of the equipment for continued operation
    if (step3_satisfied == 1) # begin
        if (annex2c_tmin_category == "Cylindrical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Spherical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Hemispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Elliptical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Torispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Toriconical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Transition")
            #tmin here
        elseif (annex2c_tmin_category == "Nozzles Connections in Shells")
            #tmin here
        elseif (annex2c_tmin_category == "Junction Reinforcement Requirements at Conical Transitions")
            #tmin here
        elseif (annex2c_tmin_category == "Tubesheets")
            #tmin here
        elseif (annex2c_tmin_category == "Flat head to cylinder connections")
            #tmin here
        elseif (annex2c_tmin_category == "Bolted Flanges")
            #tmin here
        elseif (annex2c_tmin_category == "Straight Pipes Subject To Internal Pressure")
            # calculate MAWP and Tmin
            # Tmin
            tcmin = Pipingtcmin(P, Do=Do, S=S, E=E, Yb31=Yb31, MA=MA)
            tlmin = Pipingtlmin(P, Do=Do, S=S, E=E, Yb31=Yb31, tsl=tsl, MA=MA)
            tmin = maximum([tcmin,tlmin])
            # MAWP
            MAWPc = PipingMAWPc(S, E=E, t=t, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Piping MAWPc = ",round(MAWPc, digits=3),"psi\n")
            MAWPl = PipingMAWPl(S; E=E, t=t, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Piping MAWPl = ",round(MAWPl, digits=3),"psi\n")
            MAWP = minimum([MAWPc,MAWPl])
            print("Final MAWP = ",round(MAWP, digits=3),"psi\n")

            # acceptance crtiera
            # tmin level 1 acceptance for PTR
            # The averaged measured thickness or MAWP acceptance criterion may be used. In either case, the minimum thickness criterion shall be satisfied.
            if (tam - FCAml) >= tcmin
                print("Average Measured Thickness from Point Thickness Readings (PTR) - (tam - FCAml) >= tcmin == True\n")
                print("Check the minimum measured thickness criterion.\n")
                tmin_accept=1
            else
                print("Average Measured Thickness from Point Thickness Readings (PTR) - Level 1 assessment has not been satisfied :: (tam - FCAml) >= tcmin == False\n")
                tmin_accept=0
                tmm_accept = 0
            end

            #++ Alternatively ++#
            # mawp level 1 acceptance for PTR
            # adjust tam by: (tam - FCAml)
            print("Alternatively, determine MAWPcr using (tam - FCAml)\n")
            t_adj = (tam - FCAml)
            MAWPcr = PipingMAWPc(S, E=E, t=t_adj, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Reduced MAWPcr in circumferential or hoop direction = ",round(MAWPcr, digits=3),"psi\n")
            MAWPlr = PipingMAWPl(S; E=E, t=t_adj, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Reduced MAWPlr in longitudinal direction = ",round(MAWPlr, digits=3),"psi\n")
            if (MAWPcr >= P)
                print("MAWP from Point Thickness Readings (PTR) - (MAWPcr >= P) == True\n")
                print("Check the minimum measured thickness criterion.\n")
                mawp_accept=1
            else
                print("MAWP from Point Thickness Readings (PTR) - Level 1 assessment has not been satisfied :: (MAWPcr >= P) == False\n")
                mawp_accept=0
                tmm_accept = 0
            end

            # Minimum Measured Thickness level 1 for PTR
            tmin = maximum([tcmin,tlmin])
            if units == "lbs-in-psi"
            tlim = maximum([.2*tnom,0.05])
            elseif units ==  "nmm-mm-mpa"
            tlim = maximum([.2*tnom,1.3])
            end
            if tmin_accept == 1 || mawp_accept == 1
                    if (tmm - FCAml) >= maximum([0.5*tmin,tlim])
                        print("Minimum Measured Thickness - (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == True\n")
                        tmm_accept = 1
                    else
                        print("Minimum Measured Thickness - Level 1 assessment has not been satisfied :: (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == False\n")
                        tmm_accept = 0
                    end
                end

        if tmin_accept + tmm_accept == 2
            print("Part 4 level 1 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 1 assessment (Averaged Measured Thickness)\n")
        elseif tmin_accept + tmm_accept != 2
            print("Part 4 level 1 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 1 assessment\n")
        end
        if mawp_accept + tmm_accept == 2
                print("Part 4 level 1 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 1 assessment (MAWP)\n")
        elseif mawp_accept + tmm_accept != 2
                print("Part 4 level 1 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 1 assessment\n")
        end

        # end piping level 1

        elseif (annex2c_tmin_category == "Boiler Tubes")
            #tmin here
        elseif (annex2c_tmin_category == "Pipe Bends Subject To Internal Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "MAWP for External Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "Branch Connections")
            #tmin here
        elseif (annex2c_tmin_category == "API 650 Storage Tanks")
            #tmin here
        end # end equations
        labels = ["tnom","Do","D","FCA","FCAml","tmm", "tam", "COV", "tcmin", "tlmin", "MAWPc", "MAWPl", "MAWPcr","MAWPlr", "P", "S"]
        global part_4_level_1_ctp_calculated_parameters = [tnom,Do,D,FCA,FCAml,tmm, tam, c_o_v, tcmin, tlmin, MAWPc, MAWPl, MAWPcr, MAWPlr, P, S]
        out = hcat(labels,part_4_level_1_ctp_calculated_parameters)
        return out
    elseif (step3_satisfied == 0)
        print("Use CTPs")
    end # PTR logical check
end # end function

# Level 2 Assessment
@doc """
    part4_PTR_Level_2_Assessment(x::Array{Float64}; annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0, RSFa::Float64=0.9)

Average Measured Thickness from Point Thickness Readings (PTR)

Variables\n
equipment_group = "piping" # "vessel", "tank"\n
annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
# "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
# "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]\n
flaw_location = "external" # "External","Internal"\n
units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"\n
tnom = .3 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.\n
FCAml = .05 # Future Corrosion Allowance applied to the region of metal loss.\n
FCA = .05 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).\n
LOSS = 0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.\n
Do = 3.5 # Outside Diameter\n
D = Do - 2*(tnom)  # inside diameter of the shell corrected for FCAml , as applicable\n
P = 1480.0 # internal design pressure.\n
S = 20000.0 # allowable stress.\n
E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.\n
MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.\n
Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .\n
t = tnom - LOSS - FCA  # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.\n
tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).\n

    """ ->
function part4_PTR_Level_2_Assessment(x::Array{Float64}; annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0, RSFa::Float64=0.9)
    @assert any(annex2c_tmin_category .== ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
    "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
    "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]) "Invalid entry must be any of the following tmin groups: 'Cylindrical Shell','Spherical Shell','Hemispherical Head','Elliptical Head','Torispherical Head','Conical Shell','Toriconical Head','Conical Transition','Nozzles Connections in Shells',
    'Junction Reinforcement Requirements at Conical Transitions','Tubesheets','Flat head to cylinder connections','Bolted Flanges','Straight Pipes Subject To Internal Pressure','Boiler Tubes','Pipe Bends Subject To Internal Pressure',
    'MAWP for External Pressure','Branch Connections','API 650 Storage Tanks' "
    @assert any(equipment_group .== ["piping", "vessel", "tank"]) "Invalid input: please enter either: 'piping','vessel','tank' "
    @assert any(flaw_location .== ["external", "internal"]) "Invalid input: please enter either: 'external shells', 'internal shells' "
    @assert any(units .== ["lbs-in-psi", "nmm-mm-mpa"]) "Invalid input: please enter either: 'lbs-in-psi', 'nmm-mm-mpa' "
    @assert size(x,1) >= 15 "Minimum of 15 readings required to confirm metal loss is general - provide more readings'"

    print("Begin -- Part 4 LTA Level 2 Assessment (Point Thickness Readings) - API 579 June 2016 Edition\n")

    # STEP 1 – Take the point thickness reading data in accordance with paragraph 4.3.3 Tam , and the
    # Coefficient Of Variation (COV). A template for computing the COV is provided in Table 4.3.
    # Tam = average measured wall thickness of the component based on the point thickness readings (PTR) measured at the time of the inspection.
    # COV = Coefficient Of Variation.
    # Trd = Thickness reading within corroded
    step_1 = COV_var(x)
    tmm = minimum(x)
    tam = step_1[1]
    c_o_v = step_1[2]

    # STEP 2 – If the COV from STEP 1 is less than or equal to 0.1, then proceed to STEP 3 to complete the
    # assessment using the average thickness, tam . If the COV is greater than 0.1 then the use of thickness
    # profiles should be considered for the assessment (see paragraph 4.4.2.2).
    if c_o_v <= 0.1
        print("COV is <= 0.1 - Procced to Step 3 to complete the assessment using the average thickness tam :: COV = ", round(step_1[2],digits =2),"\n")
        step3_satisfied = 1
    else
        print("COV is > 0.1 - Thickness Profiles should be conisdered for assessment :: (see paragraph 4.4.2.2) :: COV = ", round(step_1[2],digits =2),"\n")
        step3_satisfied = 0
    end

    # STEP 3 – The acceptability of the component for continued operation can be established using the Level 1
    # criteria in Table 4.4, Table 4.5, Table 4.6, and Table 4.7. The averaged measured thickness or
    # MAWP acceptance criterion may be used. In either case, the minimum thickness criterion shall be
    # satisfied. For MAWP acceptance criterion (see Part 2, paragraph 2.4.2.2.e) to determine the
    # acceptability of the equipment for continued operation
    if (step3_satisfied == 1) # begin
        if (annex2c_tmin_category == "Cylindrical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Spherical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Hemispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Elliptical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Torispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Toriconical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Transition")
            #tmin here
        elseif (annex2c_tmin_category == "Nozzles Connections in Shells")
            #tmin here
        elseif (annex2c_tmin_category == "Junction Reinforcement Requirements at Conical Transitions")
            #tmin here
        elseif (annex2c_tmin_category == "Tubesheets")
            #tmin here
        elseif (annex2c_tmin_category == "Flat head to cylinder connections")
            #tmin here
        elseif (annex2c_tmin_category == "Bolted Flanges")
            #tmin here
        elseif (annex2c_tmin_category == "Straight Pipes Subject To Internal Pressure")
            # calculate MAWP and Tmin
            P = P * RSFa # adjust P for tmin calculation
            # Tmin
            tcmin = Pipingtcmin(P, Do=Do, S=S, E=E, Yb31=Yb31, MA=MA)
            tlmin = Pipingtlmin(P, Do=Do, S=S, E=E, Yb31=Yb31, tsl=tsl, MA=MA)
            tmin = maximum([tcmin,tlmin])
            # MAWP
            MAWPc = PipingMAWPc(S, E=E, t=t, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Piping MAWPc = ",round(MAWPc, digits=3),"psi\n")
            MAWPl = PipingMAWPl(S; E=E, t=t, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Piping MAWPl = ",round(MAWPl, digits=3),"psi\n")
            MAWP = minimum([MAWPc,MAWPl])
            print("Final MAWP = ",round(MAWP, digits=3),"psi\n")

            #########++++++ acceptance crtiera ++++++#########
            # tmin level 1 acceptance for PTR
            # The averaged measured thickness or MAWP acceptance criterion may be used. In either case, the minimum thickness criterion shall be satisfied.
            if (tam - FCAml) >= maximum([tcmin,tlmin])
                print("Average Measured Thickness from Point Thickness Readings (PTR) - (tam - FCAml) >= max([tcmin,tlin]) == True\n")
                print("Check the minimum measured thickness criterion.\n")
                tmin_accept=1
            else
                print("Average Measured Thickness from Point Thickness Readings (PTR) - Level 1 assessment has not been satisfied :: (tam - FCAml) >= max([tcmin,tlin]) == False\n")
                tmin_accept=0
                tmm_accept = 0
            end

            #++ Alternatively ++#
            # mawp level 1 acceptance for PTR
            # adjust tam by: (tam - FCAml)
            print("Alternatively, determine MAWPcr using (tam - FCAml)\n")
            t_adj = (tam - tsl - FCAml)
            MAWPcr = PipingMAWPc(S, E=E, t=t_adj, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Reduced MAWPcr in circumferential or hoop direction = ",round(MAWPcr, digits=3),"psi\n")
            MAWPlr = PipingMAWPl(S; E=E, t=t_adj, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Reduced MAWPlr in longitudinal direction = ",round(MAWPlr, digits=3),"psi\n")
            if ((minimum([MAWPcr,MAWPlr]) / RSFa) >= P)
                print("MAWP from Point Thickness Readings (PTR) - ((min([MAWPcr,MAWPlr]) / RSFa) >= P) == True\n")
                print("Check the minimum measured thickness criterion.\n")
                mawp_accept=1
            else
                print("MAWP from Point Thickness Readings (PTR) - Level 1 assessment has not been satisfied :: ((min([MAWPcr,MAWPlr]) / RSFa) >= P) == False\n")
                mawp_accept=0
                tmm_accept = 0
            end

            # Minimum Measured Thickness level 1 for PTR
            tmin = maximum([tcmin,tlmin])
            if units == "lbs-in-psi"
            tlim = maximum([.2*tnom,0.05])
            elseif units ==  "nmm-mm-mpa"
            tlim = maximum([.2*tnom,1.3])
            end
            if tmin_accept == 1 || mawp_accept == 1
                    if (tmm - FCAml) >= maximum([0.5*tmin,tlim])
                        print("Minimum Measured Thickness - (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == True\n")
                        tmm_accept = 1
                    else
                        print("Minimum Measured Thickness - Level 1 assessment has not been satisfied :: (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == False\n")
                        tmm_accept = 0
                    end
                end

        if tmin_accept + tmm_accept == 2
            print("Part 4 level 2 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 2 assessment (Averaged Measured Thickness)\n")
        elseif tmin_accept + tmm_accept != 2
            print("Part 4 level 2 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 2 assessment\n")
        end
        if mawp_accept + tmm_accept == 2
                print("Part 4 level 2 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 2 assessment (MAWP)\n")
        elseif mawp_accept + tmm_accept != 2
                print("Part 4 level 2 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 2 assessment\n")
        end

        # end piping level 1

        elseif (annex2c_tmin_category == "Boiler Tubes")
            #tmin here
        elseif (annex2c_tmin_category == "Pipe Bends Subject To Internal Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "MAWP for External Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "Branch Connections")
            #tmin here
        elseif (annex2c_tmin_category == "API 650 Storage Tanks")
            #tmin here
        end # end equations
        labels = ["tnom","Do","D","FCA","FCAml","tmm", "tam", "COV", "tcmin", "tlmin", "MAWPc", "MAWPl", "MAWPcr","MAWPlr", "P", "S"]
        global part_4_level_1_ctp_calculated_parameters = [tnom,Do,D,FCA,FCAml,tmm, tam, c_o_v, tcmin, tlmin, MAWPc, MAWPl, MAWPcr, MAWPlr, P, S]
        out = hcat(labels,part_4_level_1_ctp_calculated_parameters)
        return out
    elseif (step3_satisfied == 0)
        print("Use CTPs")
    end # PTR logical check
end # end function

# CTPs

@doc """
part4_CTP_Level_1_Assessment(CTPGrid::Array{Float64,2}; FCA_string::String="external",annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0, s_spacings::Float64=0.5, c_spacings::Float64=0.5, RSFa::Float64=0.9)

Average CTP Thickness from Critial Thickness Profile (CTP)

Variables\n
FCA_string = "external" # Application of FCA external / internal\n
equipment_group = "piping" # "vessel", "tank"\n
annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
# "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
# "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]\n
flaw_location = "external" # "External","Internal"\n
units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"\n
tnom = .3 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.\n
FCAml = .05 # Future Corrosion Allowance applied to the region of metal loss.\n
FCA = .05 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).\n
LOSS = 0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.\n
Do = 3.5 # Outside Diameter\n
D = Do - 2*(tnom)  # inside diameter of the shell corrected for FCAml , as applicable\n
P = 1480.0 # internal design pressure.\n
S = 20000.0 # allowable stress.\n
E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.\n
MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.\n
Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .\n
t = tnom - LOSS - FCA  # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.\n
tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).\n
s_spacings = 0.5 # grid spacing in longitudinal direction
c_spacings = 0.5 # grid circumferential in longitudinal direction
RSFa = 0.9 # remaining strength factor - consult API 579 is go lower than 0.9\n
    """ ->
function part4_CTP_Level_1_Assessment(CTPGrid::Array{Float64,2}; FCA_string::String="external",annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0, s_spacings::Float64=0.5, c_spacings::Float64=0.5, RSFa::Float64=0.9)
    @assert any(annex2c_tmin_category .== ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
    "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
    "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]) "Invalid entry must be any of the following tmin groups: 'Cylindrical Shell','Spherical Shell','Hemispherical Head','Elliptical Head','Torispherical Head','Conical Shell','Toriconical Head','Conical Transition','Nozzles Connections in Shells',
    'Junction Reinforcement Requirements at Conical Transitions','Tubesheets','Flat head to cylinder connections','Bolted Flanges','Straight Pipes Subject To Internal Pressure','Boiler Tubes','Pipe Bends Subject To Internal Pressure',
    'MAWP for External Pressure','Branch Connections','API 650 Storage Tanks' "
    @assert any(equipment_group .== ["piping", "vessel", "tank"]) "Invalid input: please enter either: 'piping','vessel','tank' "
    @assert any(FCA_string .== ["external", "internal"]) "Invalid input: please enter either: 'external', 'internal' "
    @assert any(flaw_location .== ["external", "internal"]) "Invalid input: please enter either: 'external shells', 'internal shells' "

print("Begin -- Part 4 LTA Level 1 Assessment (Critical Thickness Profile) - API 579 June 2016 Edition\n")

    # STEP 1 – Determine the thickness profile data in accordance with paragraph 4.3.3.3 and determine the minimum measured thickness, tmm.
    tmm, long_CTP, circ_CTP = CTP_Grid_Part4(CTPGrid) # tmm = minimum measured thickness determined at the time of the inspection.

    # STEP 2 – Determine the wall thickness and diameter to be used in the assessment using Equation (4.2)
    # and Equation (4.3) or Equation (4.4).
    tml = tnom - FCAml

    # Adjust the FCA by internal and external as below
    if (FCA_string == "internal")
    Dml = D + (2*FCAml) # inside diameter of the shell corrected for FCAml, as applicable.
    elseif (FCA_string == "external")
        Dml = D # inside diameter of the shell corrected for FCAml , as applicable.
    end

    # STEP 3 – Compute the remaining thickness ratio, Rt .
    Rt = (tmm[1]-FCAml) / tml # remaining thickness ratio. # (5.5)

    # STEP 4 – Compute the length for thickness averaging, L where the parameter Q is evaluated using Table 4.8
    # handle cases for Q = 50
    if (Rt < RSFa)
        Q = 1.123*( ( ( ((1-Rt) / (1-(Rt/RSFa))) ^2) -1)^0.5)
    elseif Rt >= RSFa
        Q = 50.0
    end
    L = Q*sqrt(Dml*tml)

    # STEP 5 – Establish the Critical Thickness Profiles (CTP’s) from the thickness profile data (see paragraph 4.3.3.3).
    # Determine the average measured thickness tsam based on the longitudinal CTP and the average measured thickness tcam based on the circumferential CTP.
    # The average measured thicknesses tsam and tcam shall be based on the length L determined in STEP 4 (see Figure 4.19).
    # The length L shall be located on the respective CTP such that the resulting average thickness is a minimum.

    # create index position of grid
    long_index = collect(0:s_spacings:length(long_CTP))
    circ_index = collect(0:c_spacings:length(circ_CTP))
    bounds_grid_error = 0

    # Longitudinal CTP (tsam)
    # find index position of minimum
    long_min = minimum(long_CTP)

    min_index = Int64.(zeros(size(long_CTP,1)))
    for i = 1:length(long_CTP)
        if long_CTP[i] == long_min
             min_index[i] = i
        end
    end

    # reduce
    min_index = min_index[min_index .!= 0][1]

# find thickness to use in averaging
    if L/2 <= s_spacings
        if min_index > 1
        t1 = long_min + (long_CTP[min_index-1] - long_min) * ((L / 2) / s_spacings)
        bounds_grid_error = 0
elseif min_index == 1
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
    @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
    bounds_grid_error = 1
        t1 = long_min + (long_CTP[min_index-1] - long_min) * ((L / 2) / s_spacings)
end

    if min_index < length(long_CTP)
        t2 = long_min + (long_CTP[min_index+1] - long_min) * ((L / 2) / s_spacings)
            bounds_grid_error = 0
    elseif min_index == length(long_CTP)
        @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
        @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
            bounds_grid_error = 1
            t2 = long_min + (long_CTP[min_index+1] - long_min) * ((L / 2) / s_spacings)
        end

elseif L/2 > s_spacings # find the next thickness over the thickness averaging length/2
    # initialize array
    out = zeros(size(long_index,1))
    for i = 1:length(long_index)
        if L/2 > long_index[i]
            out[i] = true
        end
    end
s_index_pos_for_l = Int64.(sum(out))
# error check for bounds errors
if min_index - s_index_pos_for_l > 0
    t1 = long_min + (long_CTP[min_index-s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
        bounds_grid_error = 0
elseif min_index - s_index_pos_for_l <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t1 exceeds grid input size"
    @info "Please increase data grid in longitudinal direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t1 = long_min + (long_CTP[min_index-s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
end

if min_index+s_index_pos_for_l > 0
    t2 = long_min + (long_CTP[min_index+s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
        bounds_grid_error = 0
elseif min_index+s_index_pos_for_l <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t2 exceeds grid input size"
    @info "Please increase data grid in longitudinal direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t2 = long_min + (long_CTP[min_index+s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
end
end

A1 = ((t1 + long_min) / 2) * (L / 2)
A2 = ((t2 + long_min) / 2) * (L / 2)
A = A1 + A2
tsam = A / L

# Circumferential CTP (tcam)
# find index position of minimum
circ_min = minimum(circ_CTP)

min_index = Int64.(zeros(size(circ_CTP,1)))
for i = 1:length(circ_CTP)
    if circ_CTP[i] == circ_min
         min_index[i] = i
    end
end

# reduce
min_index = min_index[min_index .!= 0][1]

# find thickness to use in averaging
    if L/2 <= c_spacings
        if min_index > 1
    t3 = circ_min + (circ_CTP[min_index-1] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
elseif min_index == 1
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
    @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t3 = circ_min + (circ_CTP[min_index-1] - circ_min) * ((L / 2) / c_spacings)
end

    if min_index < length(circ_CTP)
    t4 = circ_min + (circ_CTP[min_index+1] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
    elseif min_index == length(circ_CTP)
        @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
        @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
            bounds_grid_error = 1
            t4 = circ_min + (circ_CTP[min_index+1] - circ_min) * ((L / 2) / c_spacings)
        end

    elseif L/2 > c_spacings
    # initialize array
    out = zeros(size(circ_index,1))
    for i = 1:length(circ_index)
    if L/2 > circ_index[i]
        out[i] = true
        end
    end

s_index_pos_for_c = Int64.(sum(out))
# error check for bounds errors
if min_index - s_index_pos_for_c > 0
    t3 = long_min + (circ_CTP[min_index-s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
elseif min_index - s_index_pos_for_c <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
    @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t3 = long_min + (circ_CTP[min_index-s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
end

if min_index+s_index_pos_for_c > 0
    t4 = long_min + (circ_CTP[min_index+s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
elseif min_index+s_index_pos_for_c <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t4 exceeds grid input size"
    @info "Please increase data grid in circumferential direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t4 = long_min + (circ_CTP[min_index+s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
end
end

A3 = ((t3 + circ_min) / 2) * (L / 2)
A4 = ((t4 + circ_min) / 2) * (L / 2)
A = A3 + A4
tcam = A / L

# STEP 6 – Based on the values of tsam  and ctascm  from STEP 5, determine the acceptability of the component for continued
# operation using the Level 1 criteria in Table 4.4, Table 4.5, Table 4.6, and Table 4.7, as applicable. The averaged measured
# thickness or MAWP acceptance criterion may be used. In either case, the minimum measured thickness, tmm , shall satisfy the criterion in Table 4.4, Table 4.5,
# Table 4.6 and Table 4.7. For MAWP acceptance criterion (see Part 2, paragraph 2.4.2.2.e) to
# determine the acceptability of the equipment for continued operation.
        if (annex2c_tmin_category == "Cylindrical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Spherical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Hemispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Elliptical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Torispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Toriconical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Transition")
            #tmin here
        elseif (annex2c_tmin_category == "Nozzles Connections in Shells")
            #tmin here
        elseif (annex2c_tmin_category == "Junction Reinforcement Requirements at Conical Transitions")
            #tmin here
        elseif (annex2c_tmin_category == "Tubesheets")
            #tmin here
        elseif (annex2c_tmin_category == "Flat head to cylinder connections")
            #tmin here
        elseif (annex2c_tmin_category == "Bolted Flanges")
            #tmin here
        elseif (annex2c_tmin_category == "Straight Pipes Subject To Internal Pressure")
            # calculate MAWP and Tmin
            # Tmin
            tcmin = Pipingtcmin(P, Do=Do, S=S, E=E, Yb31=Yb31, MA=MA)
            tlmin = Pipingtlmin(P, Do=Do, S=S, E=E, Yb31=Yb31, tsl=tsl, MA=MA)
            tmin = maximum([tcmin,tlmin])
            # MAWP
            MAWPc = PipingMAWPc(S, E=E, t=t, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Piping MAWPc = ",round(MAWPc, digits=3),"psi\n")
            MAWPl = PipingMAWPl(S; E=E, t=t, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Piping MAWPl = ",round(MAWPl, digits=3),"psi\n")
            MAWP = minimum([MAWPc,MAWPl])
            print("Final MAWP = ",round(MAWP, digits=3),"psi\n")

            # acceptance crtiera
            # tmin level 1 acceptance for CTP
            # The averaged measured thickness or MAWP acceptance criterion may be used. In either case, the minimum thickness criterion shall be satisfied.
            if (tsam - FCAml) >= tcmin
                print("Average Measured Thickness Longitudinal CTP (Circumferential Plane) :: (tsam - FCAml) >= tcmin == True\n")
                print("Check the minimum measured thickness criterion.\n")
                tcmin_accept=1
            else
                print("Average Measured Thickness from Longitudinal CTP (Circumferential Plane) - Level 1 assessment has not been satisfied :: (tsam - FCAml) >= tcmin == False\n")
                tcmin_accept=0
                tmm_accept = 0
            end

            if (tcam - FCAml) >= tlmin
                print("Average Measured Thickness from Circumferential CTP (Longitudinal Plane) :: (tcam - FCAml) >= tlmin == True\n")
                print("Check the minimum measured thickness criterion.\n")
                tlmin_accept=1
            else
                print("Average Measured Thickness from Circumferential CTP (Longitudinal Plane) - Level 1 assessment has not been satisfied :: (tcam - FCAml) >= tcmin == False\n")
                tlmin_accept=0
                tmm_accept = 0
            end

            #++ Alternatively ++#
            # mawp level 1 acceptance for PTR
            # adjust tam by: (tam - FCAml)
            print("Alternatively, determine MAWPcr using (tsam - FCAml)\n")
            t_adj = (tsam - FCAml)
            MAWPcr = PipingMAWPc(S, E=E, t=t_adj, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Reduced MAWPcr in circumferential or hoop direction = ",round(MAWPcr, digits=3),"psi\n")
            t_adj = (tcam - FCAml)
            MAWPlr = PipingMAWPl(S; E=E, t=t_adj, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Reduced MAWPlr in longitudinal direction = ",round(MAWPlr, digits=3),"psi\n")
            final_mawp = minimum([MAWPcr,MAWPlr])
            if (final_mawp >= P)
                print("MAWP from CTP :: (min[MAWPcr,MAWPlr] >= P) == True\n")
                print("Check the minimum measured thickness criterion.\n")
                mawp_accept=1
            else
                print("MAWP from CTP - Level 1 assessment has not been satisfied :: (MAWPcr >= P) == False\n")
                mawp_accept=0
                tmm_accept = 0
            end

            # Minimum Measured Thickness level 1 for CTP
            tmin = maximum([tcmin,tlmin])
            if units == "lbs-in-psi"
            tlim = maximum([.2*tnom,0.05])
            elseif units ==  "nmm-mm-mpa"
            tlim = maximum([.2*tnom,1.3])
            end
            if (tcmin_accept == 1 && tlmin_accept == 1) || mawp_accept == 1
                    if (tmm - FCAml) >= maximum([0.5*tmin,tlim])
                        print("Minimum Measured Thickness - (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == True\n")
                        tmm_accept = 1
                    else
                        print("Minimum Measured Thickness - Level 1 assessment has not been satisfied :: (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == False\n")
                        tmm_accept = 0
                    end
                end

        if tcmin_accept + tlmin_accept + tmm_accept == 3
            print("Part 4 level 1 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 1 assessment (Averaged CTP (over L) Thickness)\n")
        elseif tcmin_accept + tlmin_accept + tmm_accept != 3
            print("Part 4 level 1 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 1 assessment (Averaged CTP (over L) Thickness)\n")
        end
        if mawp_accept + tmm_accept == 2
                print("Part 4 level 1 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 1 assessment (Averaged CTP (over L) MAWP)\n")
        elseif mawp_accept + tmm_accept != 2
                print("Part 4 level 1 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 1 assessment (Averaged CTP (over L) MAWP)\n")
        end

        # end piping level 1

        elseif (annex2c_tmin_category == "Boiler Tubes")
            #tmin here
        elseif (annex2c_tmin_category == "Pipe Bends Subject To Internal Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "MAWP for External Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "Branch Connections")
            #tmin here
        elseif (annex2c_tmin_category == "API 650 Storage Tanks")
            #tmin here
        end # end equations

        labels = ["tnom","Do","D","Dml","FCA","FCAml","tmm", "tml", "Rt", "L","tsam","tcam", "tcmin", "tlmin", "MAWPc", "MAWPl", "MAWPcr","MAWPlr", "P", "S", "RSFa"]
        global part_4_level_1_ctp_calculated_parameters = [tnom,Do,D,Dml,FCA,FCAml,tmm, tml, Rt, L ,tsam,tcam, tcmin, tlmin, MAWPc, MAWPl, MAWPcr, MAWPlr, P, S, RSFa]
        out = hcat(labels,part_4_level_1_ctp_calculated_parameters)
        return out
end # function

# Level 2 CTP

@doc """
part4_CTP_Level_2_Assessment(CTPGrid::Array{Float64,2}; FCA_string::String="external",annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0, s_spacings::Float64=0.5, c_spacings::Float64=0.5, RSFa::Float64=0.9)

Average CTP Thickness from Critial Thickness Profile (CTP)

Variables\n
FCA_string = "external" # Application of FCA external / internal\n
equipment_group = "piping" # "vessel", "tank"\n
annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
# "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
# "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]\n
flaw_location = "external" # "External","Internal"\n
units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"\n
tnom = .3 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.\n
FCAml = .05 # Future Corrosion Allowance applied to the region of metal loss.\n
FCA = .05 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).\n
LOSS = 0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.\n
Do = 3.5 # Outside Diameter\n
D = Do - 2*(tnom)  # inside diameter of the shell corrected for FCAml , as applicable\n
P = 1480.0 # internal design pressure.\n
S = 20000.0 # allowable stress.\n
E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.\n
MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.\n
Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .\n
t = tnom - LOSS - FCA  # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.\n
tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).\n
s_spacings = 0.5 # grid spacing in longitudinal direction
c_spacings = 0.5 # grid circumferential in longitudinal direction
RSFa = 0.9 # remaining strength factor - consult API 579 is go lower than 0.9\n
    """ ->
function part4_CTP_Level_2_Assessment(CTPGrid::Array{Float64,2}; FCA_string::String="external",annex2c_tmin_category::String, equipment_group::String="piping",flaw_location::String="external",units::String="lbs-in-psi",tnom::Float64=0.0,
    FCAml::Float64=0.0,Do::Float64=0.0,D::Float64=0.0,P::Float64=0.0,S::Float64=0.0,E::Float64=0.0,MA::Float64=0.0,Yb31::Float64=0.0, tsl::Float64=0.0, t::Float64=0.0, s_spacings::Float64=0.5, c_spacings::Float64=0.5, RSFa::Float64=0.9)
    @assert any(annex2c_tmin_category .== ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
    "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
    "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]) "Invalid entry must be any of the following tmin groups: 'Cylindrical Shell','Spherical Shell','Hemispherical Head','Elliptical Head','Torispherical Head','Conical Shell','Toriconical Head','Conical Transition','Nozzles Connections in Shells',
    'Junction Reinforcement Requirements at Conical Transitions','Tubesheets','Flat head to cylinder connections','Bolted Flanges','Straight Pipes Subject To Internal Pressure','Boiler Tubes','Pipe Bends Subject To Internal Pressure',
    'MAWP for External Pressure','Branch Connections','API 650 Storage Tanks' "
    @assert any(equipment_group .== ["piping", "vessel", "tank"]) "Invalid input: please enter either: 'piping','vessel','tank' "
    @assert any(FCA_string .== ["external", "internal"]) "Invalid input: please enter either: 'external', 'internal' "
    @assert any(flaw_location .== ["external", "internal"]) "Invalid input: please enter either: 'external shells', 'internal shells' "
    @assert any(units .== ["lbs-in-psi", "nmm-mm-mpa"]) "Invalid input: please enter either: 'lbs-in-psi', 'nmm-mm-mpa' "

print("Begin -- Part 4 LTA Level 2 Assessment (Critical Thickness Profile) - API 579 June 2016 Edition\n")

    # STEP 1 – Determine the thickness profile data in accordance with paragraph 4.3.3.3 and determine the minimum measured thickness, tmm.
    tmm, long_CTP, circ_CTP = CTP_Grid_Part4(CTPGrid) # tmm = minimum measured thickness determined at the time of the inspection.

    # STEP 2 – Determine the wall thickness and diameter to be used in the assessment using Equation (4.2)
    # and Equation (4.3) or Equation (4.4).
    tml = tnom - FCAml

    # Adjust the FCA by internal and external as belows
    if (FCA_string == "internal")
    Dml = D + (2*FCAml) # inside diameter of the shell corrected for FCAml, as applicable.
    elseif (FCA_string == "external")
        Dml = D # inside diameter of the shell corrected for FCAml , as applicable.
    end

    # STEP 3 – Compute the remaining thickness ratio, Rt .
    Rt = (tmm[1]-FCAml) / tml # remaining thickness ratio. # (5.5)

    # STEP 4 – Compute the length for thickness averaging, L where the parameter Q is evaluated using Table 4.8
    # handle cases for Q = 50
    if (Rt < RSFa)
        Q = 1.123*( ( ( ((1-Rt) / (1-(Rt/RSFa))) ^2) -1)^0.5)
    elseif Rt >= RSFa
        Q = 50.0
    end
    L = Q*sqrt(Dml*tml)

    # STEP 5 – Establish the Critical Thickness Profiles (CTP’s) from the thickness profile data (see paragraph 4.3.3.3).
    # Determine the average measured thickness tsam based on the longitudinal CTP and the average measured thickness tcam based on the circumferential CTP.
    # The average measured thicknesses tsam and tcam shall be based on the length L determined in STEP 4 (see Figure 4.19).
    # The length L shall be located on the respective CTP such that the resulting average thickness is a minimum.

    # create index position of grid
    long_index = collect(0:s_spacings:length(long_CTP))
    circ_index = collect(0:c_spacings:length(circ_CTP))
    bounds_grid_error = 0

    # Longitudinal CTP (tsam)
    # find index position of minimum
    long_min = minimum(long_CTP)

    min_index = Int64.(zeros(size(long_CTP,1)))
    for i = 1:length(long_CTP)
        if long_CTP[i] == long_min
             min_index[i] = i
        end
    end

    # reduce
    min_index = min_index[min_index .!= 0][1]

# find thickness to use in averaging
    if L/2 <= s_spacings
        if min_index > 1
        t1 = long_min + (long_CTP[min_index-1] - long_min) * ((L / 2) / s_spacings)
        bounds_grid_error = 0
elseif min_index == 1
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
    @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
    bounds_grid_error = 1
        t1 = long_min + (long_CTP[min_index-1] - long_min) * ((L / 2) / s_spacings)
end

    if min_index < length(long_CTP)
        t2 = long_min + (long_CTP[min_index+1] - long_min) * ((L / 2) / s_spacings)
            bounds_grid_error = 0
    elseif min_index == length(long_CTP)
        @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
        @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
            bounds_grid_error = 1
            t2 = long_min + (long_CTP[min_index+1] - long_min) * ((L / 2) / s_spacings)
        end

elseif L/2 > s_spacings # find the next thickness over the thickness averaging length/2
    # initialize array
    out = zeros(size(long_index,1))
    for i = 1:length(long_index)
        if L/2 > long_index[i]
            out[i] = true
        end
    end
s_index_pos_for_l = Int64.(sum(out))
# error check for bounds errors
if min_index - s_index_pos_for_l > 0
    t1 = long_min + (long_CTP[min_index-s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
        bounds_grid_error = 0
elseif min_index - s_index_pos_for_l <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t1 exceeds grid input size"
    @info "Please increase data grid in longitudinal direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t1 = long_min + (long_CTP[min_index-s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
end

if min_index+s_index_pos_for_l > 0
    t2 = long_min + (long_CTP[min_index+s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
        bounds_grid_error = 0
elseif min_index+s_index_pos_for_l <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t2 exceeds grid input size"
    @info "Please increase data grid in longitudinal direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t2 = long_min + (long_CTP[min_index+s_index_pos_for_l] - long_min) * ((L / 2) / s_spacings)
end
end

A1 = ((t1 + long_min) / 2) * (L / 2)
A2 = ((t2 + long_min) / 2) * (L / 2)
A = A1 + A2
tsam = A / L

# Circumferential CTP (tcam)
# find index position of minimum
circ_min = minimum(circ_CTP)

min_index = Int64.(zeros(size(circ_CTP,1)))
for i = 1:length(circ_CTP)
    if circ_CTP[i] == circ_min
         min_index[i] = i
    end
end

# reduce
min_index = min_index[min_index .!= 0][1]

# find thickness to use in averaging
    if L/2 <= c_spacings
        if min_index > 1
    t3 = circ_min + (circ_CTP[min_index-1] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
elseif min_index == 1
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
    @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t3 = circ_min + (circ_CTP[min_index-1] - circ_min) * ((L / 2) / c_spacings)
end

    if min_index < length(circ_CTP)
    t4 = circ_min + (circ_CTP[min_index+1] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
    elseif min_index == length(circ_CTP)
        @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
        @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
            bounds_grid_error = 1
            t4 = circ_min + (circ_CTP[min_index+1] - circ_min) * ((L / 2) / c_spacings)
        end

    elseif L/2 > c_spacings
    # initialize array
    out = zeros(size(circ_index,1))
    for i = 1:length(circ_index)
    if L/2 > circ_index[i]
        out[i] = true
        end
    end

s_index_pos_for_c = Int64.(sum(out))
# error check for bounds errors
if min_index - s_index_pos_for_c > 0
    t3 = long_min + (circ_CTP[min_index-s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
elseif min_index - s_index_pos_for_c <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t3 exceeds grid input size"
    @info "Please increase data grid in the circumferential direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t3 = long_min + (circ_CTP[min_index-s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
end

if min_index+s_index_pos_for_c > 0
    t4 = long_min + (circ_CTP[min_index+s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
        bounds_grid_error = 0
elseif min_index+s_index_pos_for_c <= 0
    @error "part4_general_metal_loss_functions.jl: Thickness needed for t4 exceeds grid input size"
    @info "Please increase data grid in circumferential direction so that the correct averaging thickness's are used in the assessment"
        bounds_grid_error = 1
    t4 = long_min + (circ_CTP[min_index+s_index_pos_for_c] - circ_min) * ((L / 2) / c_spacings)
end
end

A3 = ((t3 + circ_min) / 2) * (L / 2)
A4 = ((t4 + circ_min) / 2) * (L / 2)
A = A3 + A4
tcam = A / L

# STEP 6 – Based on the values of tsam  and ctascm  from STEP 5, determine the acceptability of the component for continued
# operation using the Level 1 criteria in Table 4.4, Table 4.5, Table 4.6, and Table 4.7, as applicable. The averaged measured
# thickness or MAWP acceptance criterion may be used. In either case, the minimum measured thickness, tmm , shall satisfy the criterion in Table 4.4, Table 4.5,
# Table 4.6 and Table 4.7. For MAWP acceptance criterion (see Part 2, paragraph 2.4.2.2.e) to
# determine the acceptability of the equipment for continued operation.
        if (annex2c_tmin_category == "Cylindrical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Spherical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Hemispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Elliptical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Torispherical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Shell")
            #tmin here
        elseif (annex2c_tmin_category == "Toriconical Head")
            #tmin here
        elseif (annex2c_tmin_category == "Conical Transition")
            #tmin here
        elseif (annex2c_tmin_category == "Nozzles Connections in Shells")
            #tmin here
        elseif (annex2c_tmin_category == "Junction Reinforcement Requirements at Conical Transitions")
            #tmin here
        elseif (annex2c_tmin_category == "Tubesheets")
            #tmin here
        elseif (annex2c_tmin_category == "Flat head to cylinder connections")
            #tmin here
        elseif (annex2c_tmin_category == "Bolted Flanges")
            #tmin here
        elseif (annex2c_tmin_category == "Straight Pipes Subject To Internal Pressure")
            # calculate MAWP and Tmin
            P = P * RSFa
            # Tmin
            tcmin = Pipingtcmin(P, Do=Do, S=S, E=E, Yb31=Yb31, MA=MA)
            tlmin = Pipingtlmin(P, Do=Do, S=S, E=E, Yb31=Yb31, tsl=tsl, MA=MA)
            tmin = maximum([tcmin,tlmin])
            # MAWP
            MAWPc = PipingMAWPc(S, E=E, t=t, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Piping MAWPc = ",round(MAWPc, digits=3),"psi\n")
            MAWPl = PipingMAWPl(S; E=E, t=t, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Piping MAWPl = ",round(MAWPl, digits=3),"psi\n")
            MAWP = minimum([MAWPc,MAWPl])
            print("Final MAWP = ",round(MAWP, digits=3),"psi\n")

            # acceptance crtiera
            # tmin level 1 acceptance for CTP
            # The averaged measured thickness or MAWP acceptance criterion may be used. In either case, the minimum thickness criterion shall be satisfied.
            if (tsam - FCAml) >= tcmin
                print("Average Measured Thickness Longitudinal CTP (Circumferential Plane) :: (tsam - FCAml) >= tcmin == True\n")
                print("Check the minimum measured thickness criterion.\n")
                tcmin_accept=1
            else
                print("Average Measured Thickness from Longitudinal CTP (Circumferential Plane) - Level 2 assessment has not been satisfied :: (tsam - FCAml) >= tcmin == False\n")
                tcmin_accept=0
                tmm_accept = 0
            end

            if (tcam - FCAml) >= tlmin
                print("Average Measured Thickness from Circumferential CTP (Longitudinal Plane) :: (tcam - FCAml) >= tlmin == True\n")
                print("Check the minimum measured thickness criterion.\n")
                tlmin_accept=1
            else
                print("Average Measured Thickness from Circumferential CTP (Longitudinal Plane) - Level 2 assessment has not been satisfied :: (tcam - FCAml) >= tcmin == False\n")
                tlmin_accept=0
                tmm_accept = 0
            end

            #++ Alternatively ++#
            # mawp level 1 acceptance for PTR
            # adjust tam by: (tam - FCAml)
            print("Alternatively, determine MAWPcr using (tsam - FCAml)\n")
            t_adj = (tsam - FCAml)
            MAWPcr = PipingMAWPc(S, E=E, t=t_adj, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.147)
            print("Reduced MAWPcr in circumferential or hoop direction = ",round(MAWPcr, digits=3),"psi\n")
            t_adj = (tcam - tsl - FCAml)
            MAWPlr = PipingMAWPl(S; E=E, t=t_adj, tsl=tsl, MA=MA, Do=Do, Yb31=Yb31) # eq (2C.150)
            print("Reduced MAWPlr in longitudinal direction = ",round(MAWPlr, digits=3),"psi\n")
            final_mawp = minimum([MAWPcr,MAWPlr])
            if ((final_mawp / RSFa) >= P)
                print("MAWP from CTP :: ((min[MAWPcr,MAWPlr] / RSFa) >= P) == True\n")
                print("Check the minimum measured thickness criterion.\n")
                mawp_accept=1
            else
                print("MAWP from CTP - Level 2 assessment has not been satisfied :: (MAWPcr >= P) == False\n")
                mawp_accept=0
                tmm_accept = 0
            end

            # Minimum Measured Thickness level 1 for CTP
            tmin = maximum([tcmin,tlmin])
            if units == "lbs-in-psi"
            tlim = maximum([.2*tnom,0.05])
            elseif units ==  "nmm-mm-mpa"
            tlim = maximum([.2*tnom,1.3])
            end
            if (tcmin_accept == 1 && tlmin_accept == 1) || mawp_accept == 1
                    if (tmm - FCAml) >= maximum([0.5*tmin,tlim])
                        print("Minimum Measured Thickness - (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == True\n")
                        tmm_accept = 1
                    else
                        print("Minimum Measured Thickness - Level 2 assessment has not been satisfied :: (tmm - FCAml) >= maximum([0.5*tmin,tlim]) == False\n")
                        tmm_accept = 0
                    end
                end

        if tcmin_accept + tlmin_accept + tmm_accept == 3
            print("Part 4 level 2 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 2 assessment (Averaged CTP (over L) Thickness)\n")
        elseif tcmin_accept + tlmin_accept + tmm_accept != 3
            print("Part 4 level 2 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 2 assessment (Averaged CTP (over L) Thickness)\n")
        end
        if mawp_accept + tmm_accept == 2
                print("Part 4 level 2 - Satisfies June 2016 API 579-1/ASME FFS-1 Level 2 assessment (Averaged CTP (over L) MAWP)\n")
        elseif mawp_accept + tmm_accept != 2
                print("Part 4 level 2 - Does not satisfy June 2016 API 579-1/ASME FFS-1 Level 2 assessment (Averaged CTP (over L) MAWP)\n")
        end

        # end piping level 1

        elseif (annex2c_tmin_category == "Boiler Tubes")
            #tmin here
        elseif (annex2c_tmin_category == "Pipe Bends Subject To Internal Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "MAWP for External Pressure")
            #tmin here
        elseif (annex2c_tmin_category == "Branch Connections")
            #tmin here
        elseif (annex2c_tmin_category == "API 650 Storage Tanks")
            #tmin here
        end # end equations

        labels = ["tnom","Do","D","Dml","FCA","FCAml","tmm", "tml", "Rt", "L","tsam","tcam", "tcmin", "tlmin", "MAWPc", "MAWPl", "MAWPcr","MAWPlr", "P", "S", "RSFa"]
        global part_4_level_1_ctp_calculated_parameters = [tnom,Do,D,Dml,FCA,FCAml,tmm, tml, Rt, L ,tsam,tcam, tcmin, tlmin, MAWPc, MAWPl, MAWPcr, MAWPlr, P, S, RSFa]
        out = hcat(labels,part_4_level_1_ctp_calculated_parameters)
        return out
end # function


text:
include("./utils.jl")
include("neural_net.jl")
using Flux.Tracker
using Flux: params, mse, crossentropy
using ProgressMeter
using CuArrays

function total_loss(y, ŷ, ps)
    (π, z), (p, v)  = y, ŷ
    crossentropy(p, π) + mse(z, v) + 0.0001f0*sum([sum(i.^2) for i in ps])
end

@inline function all_equal(x)
    e1 = x[1]
    @inbounds for i=2:30
        x[i] == e1 || return false
    end
    return true
end

# function retrain(curr_nn::NeuralNet, epoch::Int=1000)
#     curr_nn = mapleaves(cu, curr_nn)
#     total = 0.
#     losses = []
#     opt = ADAM()
#     anim, plt = get_animPlot()
#     @showprogress "training network: " for i=1:epoch
#         data = load_data()
#         x = data[1] |> gpu
#         π = data[2] |> gpu
#         z = data[3] |> gpu
#         y = (π, z)
#         ps, ŷ = params(curr_nn), curr_nn(x)
#         @sync begin
#             @async begin
#                 loss = total_loss(y, ŷ, ps).data
#                 total += loss
#                 push!(plt, [loss; total/i])
#                 push!(losses, [loss; total/i])
#                 frame(anim, plt)
#             end
#             @async begin gs = gradient(ps) do
#                 total_loss(y, ŷ, ps)
#                 end
#                 Tracker.update!(opt, ps, gs)
#             end
#         end
#         i>100 && all_equal(losses[end-30:end]) && break
#     end
#     save_retrain(anim, plt, losses)
#     save_model(curr_nn, true)
#     return curr_nn
# end

ext: jl


text:
module Fahy_examples
using FinEtools
using Gaston
import Arpack: eigs

function fahy_H20_example()
    println("""
    Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
    Frank J. Fahy, Paolo Gardonio, page 483.

    Hexahedral mesh.
    """)

    t0 = time()

    rho=1.21*1e-9;# mass density
    c =343.0*1000;# millimeters per second
    bulk= c^2*rho;
    L=500.0;# length of the box, millimeters
    A=200.0; # cross-sectional area of the box
    n=40;#
    neigvs=8;
    OmegaShift=10.0;

    fens,fes = H8block(L,sqrt(A),sqrt(A),n,1,1); # Mesh
    fens,fes = H8toH20(fens,fes)

    geom = NodalField(fens.xyz)
    P = NodalField(zeros(size(fens.xyz,1),1))

    numberdofs!(P)

    femm = FEMMAcoust(IntegData(fes, GaussRule(3, 3)), MatAcoustFluid(bulk, rho))

    S = acousticstiffness(femm, geom, P);
    C = acousticmass(femm, geom, P);

    d,v,nev,nconv = eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    println("Total time elapsed = ",time() - t0,"s")

    File =  "fahy_H20.vtk"
    en = 5;
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.H20;
    scalars=[("Pressure_mode_$en", v[:,en])])
    @async run(`"paraview.exe" $File`)
    println("Done")
    true

end # fahy_H20_example

function fahy_H27_example()
    println("""
    Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
    Frank J. Fahy, Paolo Gardonio, page 483.

    Hexahedral mesh. 27-node elements.
    """)

    t0 = time()

    rho=1.21*1e-9;# mass density
    c =343.0*1000;# millimeters per second
    bulk= c^2*rho;
    L=500.0;# length of the box, millimeters
    A=200.0; # cross-sectional area of the box
    n=40;#
    neigvs=8;
    OmegaShift=10.0;

    fens,fes = H27block(L,sqrt(A),sqrt(A),n,1,1); # Mesh

    geom = NodalField(fens.xyz)
    P = NodalField(zeros(size(fens.xyz,1),1))

    numberdofs!(P)

    femm = FEMMAcoust(IntegData(fes, GaussRule(3, 3)), MatAcoustFluid(bulk, rho))

    S = acousticstiffness(femm, geom, P);
    C = acousticmass(femm, geom, P);

    d,v,nev,nconv = eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    println("Total time elapsed = ",time() - t0,"s")

    println("Done")
    true

end # fahy_H27_example

function fahy_H8_example()
    println("""
    Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
    Frank J. Fahy, Paolo Gardonio, page 483.

    Hexahedral mesh.
    """)

    t0 = time()

    rho=1.21*1e-9;# mass density
    c =343.0*1000;# millimeters per second
    bulk= c^2*rho;
    L=500.0;# length of the box, millimeters
    A=200.0; # cross-sectional area of the box
    n=40;#
    neigvs=8;
    OmegaShift=10.0;

    fens,fes = H8block(L,sqrt(A),sqrt(A),n,1,1); # Mesh

    geom = NodalField(fens.xyz)
    P = NodalField(zeros(size(fens.xyz,1),1))

    numberdofs!(P)

    femm = FEMMAcoust(IntegData(fes, GaussRule(3, 2)), MatAcoustFluid(bulk, rho))

    S = acousticstiffness(femm, geom, P);
    C = acousticmass(femm, geom, P);

    d,v,nev,nconv = eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    println("Total time elapsed = ",time() - t0,"s")

    File =  "fahy_H8.vtk"
    en = 5;
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.H8;
    scalars=[("Pressure_mode_$en", v[:,en])])
    println("Done")
    true

end # fahy_H8_example

function fahy_L2_example()
    println("""
    Example from Sound and Structural Vibration, Second Edition: Radiation, Transmission and Response [Paperback]
    Frank J. Fahy, Paolo Gardonio, page 483.

    1D mesh.
    """)

    t0  =  time()

    rho = 1.21*1e-9;# mass density
    c  = 343.0*1000;# millimeters per second
    bulk =  c^2*rho;
    L = 500.0;# length of the box, millimeters
    A = 200.0; # cross-sectional area of the box
    graphics =  true;# plot the solution as it is computed?
    n = 40;#
    neigvs = 8;
    OmegaShift = 10.0;

    fens,fes  =  L2block(L,n); # Mesh

    geom = NodalField(fens.xyz)
    P = NodalField(zeros(size(fens.xyz,1),1))

    numberdofs!(P)

    femm = FEMMAcoust(IntegData(fes, GaussRule(1, 2)), MatAcoustFluid(bulk, rho))

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    d,v,nev,nconv  = eigs(C+OmegaShift*S, S; nev = neigvs, which = :SM)
    d  =  d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    println("Total time elapsed  =  ",time() - t0,"s")

    # @pyimport matplotlib.pyplot as plt
    # plt.style[:use]("seaborn-whitegrid")
    # fig = plt.figure()
    # ax = plt.axes()
    # en = 2
    # ix = sortperm(geom.values[:])
    # ax[:plot](geom.values[:][ix], v[:,en][ix], color = "blue")
    # ax[:set_xlabel]("x")
    # ax[:set_ylabel]("P")
    # plt.show()

    set(axis="normal", plotstyle="linespoints", linewidth=2, pointsize = 2, color = "black", xlabel = "x", ylabel = "P", grid="on", title = "")
    ix = sortperm(geom.values[:])
    plot(geom.values[:][ix], v[:,2][ix], legend = "Pressure mode", marker = "edmd")

    true

end # fahy_L2_example

function allrun()
    println("#####################################################")
    println("# fahy_H20_example ")
    fahy_H20_example()
    println("#####################################################")
    println("# fahy_H27_example ")
    fahy_H27_example()
    println("#####################################################")
    println("# fahy_H8_example ")
    fahy_H8_example()
    println("#####################################################")
    println("# fahy_L2_example ")
    fahy_L2_example()
end # function allrun

end # module Fahy_examples


text:
function init_all_OHT(n_vars::Int64,func_evals::Int64,nof_funcs::Int64,nof_optimizers::Int64)
    println(Dates.format(now(), "HH:MM:SS")," >>>>>>>>>>>>>START----------------------------------------------->")
    bb_methods = copy(BlackBoxOptim.MethodNames)
    history_all=zeros(Float64,func_evals,17+3,nof_funcs).+1e10
    domain_width=10
    lb=-domain_width*ones(Float64,n_vars)
    ub=-2lb
    all_bounds=Tuple{Float64,Float64}[(lb[z], ub[z]) for z in 1:n_vars]
    optis=zeros(Float64,nof_optimizers+1+3)
    times_all=zeros(nof_optimizers+1+3)
    x_evals=zeros(func_evals,n_vars)
    opti_evals=zeros(func_evals).+1e100
    return history_all, lb, ub, all_bounds, optis, times_all, bb_methods, x_evals, opti_evals
end
function init_vars_OHT(func::Symbol,lb::Array{Float64,1},ub::Array{Float64,1},history_all::Array{Float64,3},n_vars::Int64,func_evals::Int64,fff::Int64,x_evals::Array{Float64,2},opti_evals::Array{Float64,1})
    x=lb+(ub-lb).*rand(n_vars)
    x_evals[1,:]=lb
    x_evals[2,:]=x
    x_evals[3,:]=ub
    opti_evals[1]=eval(func)(lb)
    opti_evals[2]=eval(func)(x)
    opti_evals[3]=eval(func)(ub)
    iopt=sortperm(opti_evals[1:3])
    opti_f=opti_evals[iopt[1]]
    opti_x=x_evals[iopt[1],:]
    x=copy(opti_x)
    iter=3
    history_all[1:3,17,fff].=opti_evals[1:3]
    return x, func, x_evals, opti_evals, opti_f, opti_x, iter, history_all
end
function run_OHT(n_vars::Int64,func_evals::Int64,ibb::Array{Int64,2},funcs::Array{Symbol,1})
    history_all, lb, ub, all_bounds, optis, times_all, bb_methods, x_evals, opti_evals = init_all_OHT(n_vars,func_evals,length(funcs),length(ibb))
    times_all_step=Array{Float64,1}(undef,0);iter=0;min_iter=Inf
    SR_all=zeros(length(funcs))
    alpha_plus_x=0; alpha_plus_f=0; alpha_minus_f=0; i_true=trues(func_evals)
    cdf_y=[];inp_yy=[];inp_y=[];inp_x=[];pdf_y=[];opti_f=Inf;indsBEST=[];var_xx=[];alpha1=[];var_yy=[];
    delta_fu_all=zeros(Float64,func_evals,n_vars)
    # delta_fu_all[1:3,:].=1.0
    for fff=1:length(funcs)
        nof_a1=[]
        alpha_plus_x=0; alpha_plus_f=0; alpha_minus_f=0; i_true=trues(func_evals)
        cdf_y=[];inp_yy=[];inp_y=[];inp_x=[];pdf_y=[];opti_f=Inf;indsBEST=[];var_xx=[];alpha1=[];var_yy=[];
        # try
            t1=now()
            x, func, x_evals, opti_evals, opti_f, opti_x, iter, history_all = init_vars_OHT(funcs[fff],lb,ub,history_all,n_vars,func_evals,fff,x_evals,opti_evals)
            alpha_plus_x=0; alpha_plus_f=0; alpha_minus_f=0; i_true=trues(func_evals,n_vars)
            alpha1s=[]
            while iter<func_evals
                for i=1:n_vars
                    for fght=1:1
                        alpha1=copy(iter)
                        if iter>func_evals/5 alpha1=Int64(floor((func_evals/5+10)*(1-iter/func_evals)^1)) +12 end
                        if alpha1>iter alpha1=copy(iter) end; push!(alpha1s,alpha1)
                        indsBEST=sortperm(opti_evals[1:iter])[1:alpha1]

                        if alpha1>50
                            curr_dists=x_evals[indsBEST,:].-repeat(opti_x',alpha1,1)
                            i_others=deleteat!(collect(1:n_vars),i)
                            curr_dists=sum(abs.(curr_dists[:,i_others]),dims=2)[:,1]
                            i_near=sortperm(curr_dists)[1:Int64(floor(alpha1/10))]
                            indsBEST=indsBEST[i_near]
                        end

                        inp_x=x_evals[indsBEST,i]; inp_y=opti_evals[indsBEST]; inp_d=delta_fu_all[indsBEST,i]
                        ind_sort_x=sortperm(inp_x); inp_x=inp_x[ind_sort_x]; inp_y=inp_y[ind_sort_x];

                        # if iter in [40;100;500;900] display(scatter(inp_x,inp_y,size=(320,300),legend=false,title=string(iter," ",i))) end

                        rr=rand()
                        if rr>1-iter/func_evals+0.1  # 0.5
                            x[i]=inp_x[1] .+ rand().*(inp_x[end].-inp_x[1])
                        else
                            x[i]=lb[i] .+ rand().*(ub[i].-lb[i])
                        end

                        iter+=1; x_evals[iter,:].=copy(x); fu=eval(func)(x);opti_evals[iter]=fu

                        if fu<opti_f opti_x.=copy(x); opti_f=copy(fu);  else x.=copy(vec(opti_x)) end
                        history_all[iter,17,fff]=opti_f;if iter==func_evals @goto  escape_label   end

                    end
                end
            end
            @label escape_label
            push!(times_all_step,convert(Int64,Dates.value(now()-t1))/1000)
            # println("unique f=",length(unique(opti_evals[1:iter]))," iter=",iter," opti_f=",opti_f," opti_x=",opti_x)
            if iter<func_evals history_all[iter:func_evals,17,fff].=opti_f end
            if iter<min_iter min_iter=copy(iter) end
            optis[end-3]=copy(opti_f)
            times_all[end-3]=convert(Int64,Dates.value(now()-t1))/1000
            # display(scatter(alpha1s,size=(320,300)))
            indt=Int64(0)
            for j=ibb
                indt+=1
                t1=now()
                matt,opti1,sol1 = call_bbox(bb_methods,all_bounds,iter,j,func,times_all[end])
                # println(sol1)
                history_all[1:iter,j,fff]=matt[1:iter]'
                optis[indt]=opti1
                times_all[indt]=convert(Int64,Dates.value(now()-t1))/1000
            end
            t1=now()
            res = Optim.optimize(getfield(Main, func), fill(lb[1], n_vars), fill(ub[1], n_vars), x_evals[2,:], NelderMead(), Optim.Options(f_calls_limit=func_evals,store_trace =true))
            # println(res.minimizer)
            hh=Optim.f_trace(res);hh[findfirst(hh.==minimum(hh)):end].=minimum(hh);matt=repeat(hh', trunc(Int, func_evals/length(hh))+1, 1)[:]
            history_all[1:iter,18,fff]=matt[1:iter]'
            optis[end-2]=res.minimum
            times_all[end-2]=convert(Int64,Dates.value(now()-t1))/1000

            res = Optim.optimize(getfield(Main, func), fill(lb[1], n_vars), fill(ub[1], n_vars), x_evals[2,:], ParticleSwarm(), Optim.Options(f_calls_limit=func_evals,store_trace =true))
            # println(res.minimizer)
            hh=Optim.f_trace(res);hh[findfirst(hh.==minimum(hh)):end].=minimum(hh);matt=repeat(hh', trunc(Int, func_evals/length(hh))+1, 1)[:];
            history_all[1:iter,19,fff]=matt[1:iter]'
            optis[end-1]=res.minimum
            times_all[end-1]=convert(Int64,Dates.value(now()-t1))/1000

            res = Optim.optimize(getfield(Main, func), fill(lb[1], n_vars), fill(ub[1], n_vars), x_evals[2,:], SAMIN(verbosity=0), Optim.Options(f_calls_limit=func_evals,store_trace =true));
            # println(res.minimizer)
			matt=[]
            try hh=Optim.f_trace(res);hh[findfirst(hh.==minimum(hh)):end].=minimum(hh);matt=repeat(hh', trunc(Int, func_evals/length(hh))+1, 1)[:]; catch
                matt=collect(range(opti_evals[2],length=func_evals,stop=res.minimum)); end;
            history_all[1:iter,20,fff]=matt[1:iter]';
            optis[end]=res.minimum;
            times_all[end]=convert(Int64,Dates.value(now()-t1))/1000

            SR_all[fff]=(opti_f-minimum(optis))/maximum(opti_evals)
            println(fff," of ",length(funcs)," ",Dates.format(now(), "HH:MM:SS")," func=", string(funcs[fff]), " ",optis," first_fu= ",convert(Int64,floor(opti_evals[1]))," times= ",times_all," SR_all[fff]= ",SR_all[fff],"\n")

            # catch ex
        #     if isa(ex, InterruptException) @warn("InterruptException..."); break
        #     else println(iter," ",ex," ",inp_y," ",pdf_y)
        #         if iter<min_iter min_iter=copy(iter) end end
        # end
    end
    # display(bar(SR_all,size=(320,300)))
    # println(delta_fu_all)
    return history_all, times_all_step, x_evals, opti_evals, cdf_y,inp_yy,inp_x,pdf_y,min_iter,i_true
end

# inp_yy[inp_yy.<1e-5].=0.0
# pdf_y=exp.(-(20-16iter/func_evals)inp_yy.^2)
# pdf_y=exp.(-4inp_yy.^2)
# pdf_y=1 ./(.1 .+inp_yy.^2)
# pdf_y=1 ./(inp_yy .+1e-100 .+(1e-10)*(1-iter/func_evals))
# pdf_y=1 ./(inp_yy .+1e-2)
# pdf_y=1 ./(inp_yy .+1e-2*(1-iter/func_evals))


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666],
     y=[0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512],
     color=[0.2774, 0.800248, 0.934636, 0.527809, 0.866783, 0.803046, 0.278974, 0.92097, 0.292635, 0.743778],
     Scale.color_sqrt)


text:
export af_open, FilePlayer, rewind, samplerate

const SFM_READ = int32(0x10)
const SFM_WRITE = int32(0x20)

const SF_FORMAT_WAV =  0x010000
const SF_FORMAT_FLAC = 0x170000
const SF_FORMAT_OGG =  0x200060

const SF_FORMAT_PCM_S8 = 0x0001 # Signed 8  bit data
const SF_FORMAT_PCM_16 = 0x0002 # Signed 16 bit data
const SF_FORMAT_PCM_24 = 0x0003 # Signed 24 bit data
const SF_FORMAT_PCM_32 = 0x0004 # Signed 32 bit data

const SF_SEEK_SET = 0
const SF_SEEK_CUR = 1
const SF_SEEK_END = 2

const EXT_TO_FORMAT = [
    ".wav" => SF_FORMAT_WAV,
    ".flac" => SF_FORMAT_FLAC
]

type SF_INFO
    frames::Int64
    samplerate::Int32
    channels::Int32
    format::Int32
    sections::Int32
    seekable::Int32

    function SF_INFO(frames::Integer, samplerate::Integer, channels::Integer,
                     format::Integer, sections::Integer, seekable::Integer)
        new(int64(frames), int32(samplerate), int32(channels), int32(format),
            int32(sections), int32(seekable))
    end
end

type AudioFile
    filePtr::Ptr{Void}
    sfinfo::SF_INFO
end

samplerate(f::AudioFile) = f.sfinfo.samplerate

# AudioIO.open is part of the public API, but is not exported so that it
# doesn't conflict with Base.open
function open(path::String, mode::String = "r",
            sampleRate::Integer = 44100, channels::Integer = 1,
            format::Integer = 0)
    @assert channels <= 2

    sfinfo = SF_INFO(0, 0, 0, 0, 0, 0)
    file_mode = SFM_READ

    if mode == "w"
        file_mode = SFM_WRITE
        sfinfo.samplerate = sampleRate
        sfinfo.channels = channels
        if format == 0
            _, ext = splitext(path)
            sfinfo.format = EXT_TO_FORMAT[ext] | SF_FORMAT_PCM_16
        else
            sfinfo.format = format
        end
    end

    filePtr = ccall((:sf_open, libsndfile), Ptr{Void},
                    (Ptr{Uint8}, Int32, Ptr{SF_INFO}),
                    path, file_mode, &sfinfo)

    if filePtr == C_NULL
        errmsg = ccall((:sf_strerror, libsndfile), Ptr{Uint8}, (Ptr{Void},), filePtr)
        error(bytestring(errmsg))
    end

    return AudioFile(filePtr, sfinfo)
end

function Base.close(file::AudioFile)
    err = ccall((:sf_close, libsndfile), Int32, (Ptr{Void},), file.filePtr)
    if err != 0
        error("Failed to close file")
    end
end

function open(f::Function, args...)
    file = AudioIO.open(args...)
    try
        f(file)
    finally
        close(file)
    end
end

function af_open(args...)
    warn("af_open is deprecated, please use AudioIO.open instead")
    AudioIO.open(args...)
end

# TODO: we should implement a general read(node::AudioNode) that pulls data
# through an arbitrary render chain and returns the result as a vector
function Base.read(file::AudioFile, nframes::Integer, dtype::Type)
    @assert file.sfinfo.channels <= 2
    # the data comes in interleaved
    arr = zeros(dtype, file.sfinfo.channels, nframes)

    if dtype == Int16
        nread = ccall((:sf_readf_short, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Int16}, Int64),
                        file.filePtr, arr, nframes)
    elseif dtype == Int32
        nread = ccall((:sf_readf_int, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Int32}, Int64),
                        file.filePtr, arr, nframes)
    elseif dtype == Float32
        nread = ccall((:sf_readf_float, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Float32}, Int64),
                        file.filePtr, arr, nframes)
    elseif dtype == Float64
        nread = ccall((:sf_readf_double, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Float64}, Int64),
                        file.filePtr, arr, nframes)
    end

    return arr[:, 1:nread]'
end

Base.read(file::AudioFile, dtype::Type) = Base.read(file, file.sfinfo.frames, dtype)
Base.read(file::AudioFile, nframes::Integer) = Base.read(file, nframes, Int16)
Base.read(file::AudioFile) = Base.read(file, Int16)

function Base.write{T}(file::AudioFile, frames::Array{T})
    @assert file.sfinfo.channels <= 2
    nframes = int(length(frames) / file.sfinfo.channels)

    if T == Int16
        return ccall((:sf_writef_short, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Int16}, Int64),
                        file.filePtr, frames, nframes)
    elseif T == Int32
        return ccall((:sf_writef_int, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Int32}, Int64),
                        file.filePtr, frames, nframes)
    elseif T == Float32
        return ccall((:sf_writef_float, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Float32}, Int64),
                        file.filePtr, frames, nframes)
    elseif T == Float64
        return ccall((:sf_writef_double, libsndfile), Int64,
                        (Ptr{Void}, Ptr{Float64}, Int64),
                        file.filePtr, frames, nframes)
    end
end

function Base.seek(file::AudioFile, offset::Integer, whence::Integer)
    new_offset = ccall((:sf_seek, libsndfile), Int64,
        (Ptr{Void}, Int64, Int32), file.filePtr, offset, whence)

    if new_offset < 0
        error("Could not seek to $(offset) in file")
    end

    new_offset
end

# Some convenience methods for easily navigating through a sound file
Base.seek(file::AudioFile, offset::Integer) = seek(file, offset, SF_SEEK_SET)
rewind(file::AudioFile) = seek(file, 0, SF_SEEK_SET)

type FileRenderer <: AudioRenderer
    file::AudioFile

    function FileRenderer(file::AudioFile)
        node = new(file)
        finalizer(node, n -> close(n.file))
        return node
    end
end

typealias FilePlayer AudioNode{FileRenderer}
FilePlayer(file::AudioFile) = FilePlayer(FileRenderer(file))
FilePlayer(path::String) = FilePlayer(AudioIO.open(path))

function render(node::FileRenderer, device_input::AudioBuf, info::DeviceInfo)
    @assert node.file.sfinfo.samplerate == info.sample_rate

    # Keep reading data from the file until the output buffer is full, but stop
    # as soon as no more data can be read from the file
    audio = Array(AudioSample, 0, node.file.sfinfo.channels)
    while true
        read_audio = read(node.file, info.buf_size-size(audio, 1), AudioSample)
        audio = vcat(audio, read_audio)
        if size(audio, 1) >= info.buf_size || size(read_audio, 1) <= 0
            break
        end
    end

    # if the file is stereo, mix the two channels together
    if node.file.sfinfo.channels == 2
        return (audio[:, 1] / 2) + (audio[:, 2] / 2)
    else
        return audio
    end
end

function play(filename::String, args...)
    player = FilePlayer(filename)
    play(player, args...)
end

function play(file::AudioFile, args...)
    player = FilePlayer(file)
    play(player, args...)
end


text:
using OptiMo, Bazinga, ScSTO
using Printf, PyPlot, PyCall
using DelimitedFiles

foldername = "/home/albertodm/Documents/ScSTO.jl/demo/"

"""
Maintenance optimization problem

    [SAL14]     Sun, Aw, Loxton, Teo, "An optimal machine maintenance problem with
    probabilistic state constraints" (2014)

    state : machine state `x` and cumulated profit `y`

 3-dynamics
 (i) production : full output, no maintenance
 (ii) minor maintenance : reduced output, slow, cheap maintenance, minor improvement
 (iii) major maintenance : no output, fast maintenance, major improvement
"""
# time interval
t0 = 0.0
tf = 1.0

# initial state
x0 = [1.0; 0.0]
nx = length(x0)

# control sequence
uvec = Matrix(repeat([0.0; 1.0; 0.0; 1.0; 0.0; 2.0], 3, 1)')
uvec = uvec[:, 1:end-1]
nu = size(uvec, 1)
N = size(uvec, 2)

# cost matrices
Q = Matrix([0.0 0.0; 0.0 0.0])
E = Matrix([0.0 0.0; 0.0 -1.0])

# system dynamics
function dynam(x::Vector{Float64}, u::Vector{Float64})
    k = [2.0; 25.0; 1.0; 1.0; 8.0; 2.0; 50.0; 40.0]
    if u[1] < 0.5 # no maintenance
        f = [-k[1] * x[1]; k[2] * x[1]^2 - k[3]]
    elseif u[1] < 1.5 # minor maintenance
        f = [k[4] * (1 - x[1]); k[5] * x[1]^2 - k[6]]
    else # major maintenance
        f = [k[7] * (1 - x[1]^2); -k[8]]
    end
    return f
end
function d_dynam(x::Vector{Float64}, u::Vector{Float64})
    k = [2.0; 25.0; 1.0; 1.0; 8.0; 2.0; 50.0; 40.0]
    if u[1] < 0.5 # no maintenance
        df = [
            -k[1] 0
            2 * k[2] * x[1] 0
        ]
    elseif u[1] < 1.5 # minor maintenance
        df = [
            -k[4] 0
            2 * k[5] * x[1] 0
        ]
    else # major maintenance
        df = [
            -2 * k[7] * x[1] 0
            0 0
        ]
    end
    return df
end

# fixed grid points
ngrid = 400

# STO problem with switching costs
prob = scstoproblem(
    x0,
    dynam,
    d_dynam,
    uvec,
    ngrid = ngrid,
    t0 = t0,
    tf = tf,
    Q = Q,
    E = E,
    swc = 0.0,
);

# solve
solver = Bazinga.ZEROFPR(
    max_iter = 100,
    verbose = true,
    gamma_min = 1e-16,
)
out = solver(prob)
print(out)

# switching intervals
delta0 = copy(prob.meta.x0)
swdelta = out.x
# switching times
tau0, _ = gettau(prob, delta0)
swtau, swtf = gettau(prob, swdelta)
# simulation
nt = 1000
xsim0, _, objective0, tsim = simulate(prob, tau0, length = nt)
usim0, _ = simulateinput(prob, tau0, tsim)
xsim, _, objective, _ = simulate(prob, swtau, tsim)
usim, _ = simulateinput(prob, swtau, tsim)
objf_0 = prob.repo.objf

tmpmat = Array{Float64}(undef, nt, 1+2*3)
tmpmat[:,1] = tsim
tmpmat[:,2] = xsim0[1, :]
tmpmat[:,3] = xsim0[2, :]
tmpmat[:,4] = usim0[1, :]
tmpmat[:,5] = xsim[1, :]
tmpmat[:,6] = xsim[2, :]
tmpmat[:,7] = usim[1, :]

writedlm( foldername * "omm_init_traj.csv",  tmpmat, ',')

#=
figure()
subplot(3, 1, 1)
plot(tsim, xsim0[1, :], linestyle = "--")
plot(tsim, xsim[1, :])
ylim(0, 1.1)
xlim(0, 1)
yticks([0; 1])
ylabel(L"s")
subplot(3, 1, 2)
plot(tsim, xsim0[2, :], linestyle = "--")
plot(tsim, xsim[2, :])
ylabel(L"p")
xlim(0, 1)
subplot(3, 1, 3)
plot(tsim, usim0[1, :], linestyle = "--")
plot(tsim, usim[1, :])
ylim(-0.2, 2.2)
xlim(0, 1)
yticks([0; 1; 2])
ylabel(L"u")
xlabel(L"t")
gcf()
savefig(foldername * "ScSTO.jl/demo/data/omm_traj.pdf")
=#

################################################################################
# with switching cost
################################################################################

swc_grid = [25.0; 30.0] # 20, 25, 30
ng = length(swc_grid)

objective_swc = Array{Float64}(undef, ng)
swdelta_swc = Array{Float64}(undef, N, ng)
swtau_swc = Array{Float64}(undef, N - 1, ng)
xsim_swc = Array{Float64}(undef, nx, nt, ng)
usim_swc = Array{Float64}(undef, nu, nt, ng)
objf_swc = Array{Any}(undef, ng)

for k = 1:ng
    swc = swc_grid[k]

    # STO problem with switching costs
    prob_swc = scstoproblem(
        x0,
        dynam,
        d_dynam,
        uvec,
        ngrid = ngrid,
        t0 = t0,
        tf = tf,
        Q = Q,
        E = E,
        swc = swc,
        tau0ws = swtau,
    )

    out_swc = solver(prob_swc)
    print(out_swc)

    swdelta_swc[:, k] = out_swc.x
    swtau_swc[:, k], _ = gettau(prob_swc, swdelta_swc[:, k])
    xsim_swc[:, :, k], _, objective_swc[k], _ = simulate(prob_swc, swtau_swc[:, k], tsim)
    usim_swc[:, :, k], _ = simulateinput(prob_swc, swtau_swc[:, k], tsim)
    objf_swc[k] = prob_swc.repo.objf
end

tmpmat = Array{Float64}(undef, nt, 1+3*(ng+1))
tmpmat[:,1] = t
tmpmat[:,2] = xsim[1, :]
tmpmat[:,3] = xsim[2, :]
tmpmat[:,4] = usim[1, :]
for k = 1:ng
    tmpmat[:,3*k+2] = xsim_swc[1, :, k]
    tmpmat[:,3*k+3] = xsim_swc[2, :, k]
    tmpmat[:,3*k+4] = usim_swc[1, :, k]
end
writedlm( foldername * "omm_swc_traj.csv",  tmpmat, ',')

#=
figure()
subplot(3, 1, 1)
plot(tsim, xsim[1, :], linestyle = "--")
for k = 1:ng
    plot(tsim, xsim_swc[1, :, k])
end
ylim(0, 1.1)
xlim(0, 1)
yticks([0; 1])
ylabel(L"s")
subplot(3, 1, 2)
plot(tsim, xsim[2, :], linestyle = "--")
for k = 1:ng
    plot(tsim, xsim_swc[2, :, k])
end
ylabel(L"p")
xlim(0, 1)
subplot(3, 1, 3)
plot(tsim, usim[1, :], linestyle = "--", label = "σ = 0.0")
for k = 1:ng
    plot(tsim, usim_swc[1, :, k], label = "σ = $(swc_grid[k])")
end
#legend(loc = "upper right")
legend(loc = "upper center", bbox_to_anchor=(0.5, -0.2), ncol=5)
ylim(-0.2, 2.2)
xlim(0, 1)
yticks([0; 1; 2])
ylabel(L"u")
xlabel(L"t")
gcf()
savefig(foldername * "ScSTO.jl/demo/data/omm_swc_traj.pdf")
=#

#=
figure()
semilogy(objf_0 .- minimum(objf_0))
for k in 1:ng
    semilogy(objf_swc[k] .- minimum(objf_swc[k]))
end
gcf()
savefig(foldername * "ScSTO.jl/demo/omm_objf.pdf")
=#


text:
using Test

excluded = ["build.jl", "deps.jl", "runtests.jl", "testcommon.jl"]

if Sys.iswindows() && Sys.WORD_SIZE == 32
  push!(excluded, "except.jl")
end

include(joinpath(@__DIR__, "build.jl"))

@testset "CxxWrap tests" begin
  @testset "$f" for f in filter(fname -> fname ∉ excluded, readdir())
    println("Running tests from $f...")
    include(f)
  end
end