text:
using TrussMorph
tm = TrussMorph
import JSON
using Dates
using Printf

function parse_truss_json(file_path::String; parse_morph=false)
    data = Dict()
    open(file_path, "r") do f
        data_txt = read(f, String)
        data = JSON.parse(data_txt)
    end
    dim = data["dimension"]
    n_nodes = data["node_num"]
    n_elements = data["element_num"]

    # get material properties
    # pressure: kN/cm^2 -> kN/m^2
    # density: kN/m^3
    mp = tm.MaterialProperties(data["material_properties"]["material_name"],
                            data["material_properties"]["youngs_modulus"] * 1e4,
                            data["material_properties"]["shear_modulus"] * 1e4,
                            data["material_properties"]["poisson_ratio"],
                            data["material_properties"]["density"])

    X = Matrix{Float64}(undef, n_nodes,2)
    T = Matrix{Int}(undef, n_elements,2)
    fix_node_ids = []

    # get node coord
    for i=1:n_nodes
        X[i,:] = hcat(data["node_list"][i]["point"]["X"],
                      data["node_list"][i]["point"]["Y"])
                      # data["node_list"][i]["point"]["Z"])
        if 1 == data["node_list"][i]["is_grounded"]
            push!(fix_node_ids, i)
        end
    end

    # get fixities
    fix_dof = [1,2,6]
    if 2 != dim
        fix_dof = 1:1:6
    end

    S = Matrix{Int}(undef, length(fix_node_ids),length(fix_dof)+1)
    for i=1:length(fix_node_ids)
        S[i,1] = fix_node_ids[i]
        S[i,2:end] = data["node_list"][fix_node_ids[i]]["fixities"][fix_dof]'
    end

    # get element node ids
    for i=1:n_elements
        T[i,:] = (data["element_list"][i]["end_node_ids"] .+ 1)'
    end

    morph_data = Dict()
    if parse_morph
        morph_data = data["morph_data"]
        # @show morph_data
    end

    return Truss(X, T, S, mp), morph_data
end

function parse_load_json(file_path::String, node_dof::Int)
    data = Dict()
    open(file_path, "r") do f
        data_txt = read(f, String)
        data = JSON.parse(data_txt)
    end

    dim = data["dimension"]
    n_load_nodes = length(data["point_load_list"])
    @assert(dim == 2)

    Load = zeros(n_load_nodes, 1+node_dof)
    for i=1:n_load_nodes
        Load[i,1] = data["point_load_list"][i]["applied_node_id"] + 1
        if 2 == dim
            Load[i,2] = data["point_load_list"][i]["Fx"]
            Load[i,3] = data["point_load_list"][i]["Fy"]
            if 3 == node_dof
                Load[i,4] = data["point_load_list"][i]["Mz"]
            end
        end
    end

    @assert(n_load_nodes > 0)
    return Load
    # TODO: include_self_weight
end

function save_morph_path_json(morph_path::Array{Matrix{Float64}}, file_dir::String, st_file_name::String, end_file_name::String, t0::TrussMorph.Truss, sE::Vector{Float64}, wE::Vector{Float64}, totE::Vector{Float64}, sumE::Float64, parm_smooth::Float64, parm_weight::Float64)
    pure_st_file_name = SubString(st_file_name, 1:length(st_file_name)-length(".json"))
    pure_end_file_name = SubString(end_file_name, 1:length(end_file_name)-length(".json"))
    result_file_name = pure_st_file_name * "-" *  pure_end_file_name
    f_file_dir = joinpath(file_dir, result_file_name)
    if !ispath(f_file_dir)
        mkpath(f_file_dir)
    end

    # same material
    mp_data = Dict()
    mp_data["material_name"] = t0.mp.name
    mp_data["youngs_modulus"] = t0.mp.E * 1e-4
    mp_data["youngs_modulus_unit"] = "kN/cm2"
    mp_data["shear_modulus"] = t0.mp.G * 1e-4
    mp_data["shear_modulus_unit"] = "kN/cm2"
    mp_data["tensile_yeild_stress"] = "N/A"
    mp_data["tensile_yeild_stress_unit"] = "kN/cm2"
    mp_data["density"] = t0.mp.ρ
    mp_data["density_unit"] = "kN/m3"
    mp_data["poisson_ratio"] = t0.mp.μ

    mp_data["radius"] = "N/A"
    mp_data["radius_unit"] = "centimeter"

    # same topology
    topo_data = Dict[]
    for e=1:size(t0.T,1)
        e_data = Dict()
        e_data["end_node_ids"] = t0.T[e,:]
        e_data["element_id"] = e-1
        e_data["layer_id"] = 0
        push!(topo_data, e_data)
    end

    # morph data
    morph_data = Dict()
    morph_data["smooth_parameter"] = parm_smooth
    morph_data["weight_parameter"] = parm_weight
    morph_data["time_in_path"] = string(0,"/",0)
    morph_data["weight_energy"] = 0.0
    morph_data["smoothness_energy"] = 0.0
    morph_data["tot_energy"] = 0.0

    plen = length(morph_path)
    for i=1:plen
        data = Dict()
        data["model_name"] = result_file_name * "_mp" * string(i,"-",plen)
        data["model_type"]= "2D_frame"
        data["unit"] = "meter"
        data["generate_time"] = string(Dates.now())
        data["dimension"] = size(morph_path[i],2)

        i_morph_data = morph_data #deepcopy
        i_morph_data["time_in_path"] = string(i,"/",plen)
        i_morph_data["smoothness_energy"] = sE[i]
        i_morph_data["weight_energy"] = wE[i]
        i_morph_data["tot_energy"] = totE[i]
        i_morph_data["sum_energy"] = sumE
        data["morph_data"] = i_morph_data

        data["node_num"] = size(morph_path[i],1)
        data["element_num"] = size(t0.T,1)

        data["material_properties"] = mp_data
        data["node_list"] = Dict[]
        for j=1:size(morph_path[i],1)
            pt_data = Dict()
            pt_data["point"] = Dict()
            pt_data["point"]["X"] = morph_path[i][j,1]
            pt_data["point"]["Y"] = morph_path[i][j,2]
            pt_data["node_id"] = j-1

            pt_fix = findall(x->x==j, t0.S[:,1])
            pt_data["is_grounded"] = !isempty(pt_fix)
            if !isempty(pt_fix)
                @assert(length(pt_fix) == 1)
                pt_data["fixities"] = ones(Int, 6)
                if 2 == size(morph_path[i],1)
                    pt_data["fixities"][1] = t0.S[pt_fix,1]
                    pt_data["fixities"][2] = t0.S[pt_fix,2]
                    pt_data["fixities"][6] = t0.S[pt_fix,3]
                end
            else
                pt_data["fixities"] = []
            end
            push!(data["node_list"], pt_data)
        end
        data["element_list"] = topo_data

        # write
        stringdata = JSON.json(data)
        result_json_name = result_file_name * "_" * string(@sprintf("%03d",i),"-",plen) * ".json"
        open(joinpath(f_file_dir,result_json_name), "w") do f
                write(f, stringdata)
        end
    end
end


text:
using TypeClasses

# MonoidAlternative
# =================

function create_tuple_monoid_definitions(n)
  # unary and type arguments for `neutral`
  unary_typevars = Symbol.(:A, 1:n)
  unary_args = [:(a::Type{Tuple{$(unary_typevars...)}})]
  unary_return(unary) = :(tuple($(unary.(unary_typevars)...)))

  # binary and value arguments, for both `combine` and `orelse`
  binary_typevars_1 = Symbol.(:A, 1:n)
  binary_typevars_2 = Symbol.(:B, 1:n)
  binary_args = [:(a::Tuple{$(binary_typevars_1...)}), :(b::Tuple{$(binary_typevars_2...)})]
  binary_return(binary) = :(tuple($(binary.([:(getfield(a, $i)) for i in 1:n], [:(getfield(b, $i)) for i in 1:n])...)))
  quote
    function TypeClasses.neutral($(unary_args...)) where {$(unary_typevars...)}
      $(unary_return(A -> :(neutral($A))))
    end
    function TypeClasses.combine($(binary_args...)) where {$(binary_typevars_1...), $(binary_typevars_2...)}
      $(binary_return((a, b) -> :($a ⊕ $b)))
    end

    # we don't implement orelse, as it is commonly meant on container level, but there is no obvious failure semantics here
  end
end

macro create_tuple_monoid_definitions(n::Int)
  esc(Expr(:block, [create_tuple_monoid_definitions(i) for i ∈ 1:n]...))
end

@create_tuple_monoid_definitions 20  # 20 is an arbitrary number

ext: jl


text:
using FLOWFarm
using CCBlade

const ff=FLOWFarm

"""
    multiple_components_op(U, V, W, Omega, r, precone, yaw, tilt, azimuth, rho, mu=1.81206e-05, asound=1.0)

Return the operating points along the blade considering all the inflow velocity components.

# Arguments
- `U::Array{Float}`: u velocity component of the inflow at each r
- `V::Array{Float}`: v velocity component of the inflow at each r
- `W::Array{Float}`: w velocity component of the inflow at each r
- `Omega::Float`: rotor rotational speed
- `r::Array{Float}`: radial locations of interest
- `precone::Float`: rotor precone angle
- `yaw::Float`: rotor yaw angle
- `tilt::Float`: rotor tilt angle
- `azimuth::Float`: blade azimuth angle
- `rho::Float`: air density

# Keyword Arguments
- `mu::Float`: air viscocity (can usually use the default)
- `asound::Float`: speed of sound (can usually use the default)
"""
function multiple_components_op(U, V, W, Omega, r, precone, yaw, tilt, azimuth, rho, mu=1.81206e-05, asound=1.0)

    sy = sin(yaw)
    cy = cos(yaw)
    st = sin(tilt)
    ct = cos(tilt)
    sa = sin(azimuth)
    ca = cos(azimuth)
    sc = sin(precone)
    cc = cos(precone)

    magnitude = sqrt.(U.^2 + V.^2 + W.^2)
    x0 = U./magnitude
    y0 = V./magnitude
    z0 = W./magnitude

    #first
    x1 = x0.*cy + y0.*sy
    y1 = -x0.*sy + y0.*cy
    z1 = z0

    #second
    x2 = -z1.*st + x1.*ct
    y2 = y1
    z2 = z1.*ct + x1.*st

    #third
    x3 = x2
    y3 = y2.*ca + z2.*sa
    z3 = -y2.*sa + z2.*ca

    #fourth
    x4 = z3.*sc + x3.*cc
    y4 = y3
    z4 = z3.*cc - x3.*sc

    # transform wind to blade c.s.
    Vwind_x = magnitude.*x4
    Vwind_y = magnitude.*y4
    Vwind_z = magnitude.*z4

    # coordinate in azimuthal coordinate system
    z_az = r.*cos(precone)

    # wind from rotation to blade c.s.
    Vrot_x = -Omega*sc
    Vrot_y = Omega.*z_az

    # total velocity
    Vx = Vwind_x .+ Vrot_x
    Vy = Vwind_y + Vrot_y

    # operating point
    return OperatingPoint(Vx, Vy, rho, mu, asound)
end

"""
    distributed_velocity_op(V, Omega, r, precone, yaw, tilt, azimuth, rho, mu=1.81206e-05, asound=1.0)

Return the operating points along the blade considering varied inflow along the blade.

# Arguments
- `V::Array{Float}`: velocity inflow at each r
- `Omega::Float`: rotor rotational speed
- `r::Array{Float}`: radial locations of interest
- `precone::Float`: rotor precone angle
- `yaw::Float`: rotor yaw angle
- `tilt::Float`: rotor tilt angle
- `azimuth::Float`: blade azimuth angle
- `rho::Float`: air density

# Keyword Arguments
- `mu::Float`: air viscocity (can usually use the default)
- `asound::Float`: speed of sound (can usually use the default)
"""
function distributed_velocity_op(V, Omega, r, precone, yaw, tilt, azimuth, rho, mu=1.81206e-05, asound=1.0)

    sy = sin(yaw)
    cy = cos(yaw)
    st = sin(tilt)
    ct = cos(tilt)
    sa = sin(azimuth)
    ca = cos(azimuth)
    sc = sin(precone)
    cc = cos(precone)

    # coordinate in azimuthal coordinate system
    x_az = -r.*sin(precone)
    z_az = r.*cos(precone)

    # get section heights in wind-aligned coordinate system
    heightFromHub = (z_az*ca)*ct - x_az*st

    # transform wind to blade c.s.
    Vwind_x = V * ((cy*st*ca + sy*sa)*sc + cy*ct*cc)
    Vwind_y = V * (cy*st*sa - sy*ca)

    # wind from rotation to blade c.s.
    Vrot_x = -Omega*sc
    Vrot_y = Omega*z_az

    # total velocity
    Vx = Vwind_x + Vrot_x
    Vy = Vwind_y + Vrot_y

    # operating point
    return OperatingPoint(Vx, Vy, rho, mu, asound)
end

"""
    find_xyz_simple(x_hub,y_hub,z_hub,r,yaw,azimuth)

Find the xyz locations of points along a blade given it's location and azimuth angle.
Currently doesn't consider precone or tilt.

# Arguments
- `x_hub::Float`: x location of hub
- `y_hub::Float`: y location of hub
- `z_hub::Float`: z location of hub (hub height if no topology)
- `r::Array{Float}`: radial locations of interest
- `precone::Float`: rotor precone angle
- `yaw::Float`: rotor yaw angle
- `azimuth::Float`: blade azimuth angle
"""
function find_xyz_simple(x_hub,y_hub,z_hub,r,yaw,azimuth)

        sy = sin(yaw)
        cy = cos(yaw)
        sa = sin(azimuth)
        ca = cos(azimuth)

        x_locs = x_hub .- r.*(sy*sa)
        y_locs = y_hub .- r.*(cy*sa)
        z_locs = z_hub .+ r.*(ca)

        return x_locs, y_locs, z_locs
end

"""
    calc_moment_stress(mx,my,dx,dy,Rcyl=1.771,tcyl=0.06)

Calculates stresses from bending moments on a hollow cylinder

# Arguments
- `mx::Float`: x moment
- `my::Float`: y moment
- `dx::Float`: x distance to the location of interest
- `dy::Float`: y distance to the location of interest

# Keyword Arguments
- `Rcyl::Float`: radius of the cylinder
- `tcyl::Float`: thickenss of the cylinder
"""
function calc_moment_stress(mx,my,dx,dy,Rcyl=1.771,tcyl=0.06)

        I = 0.25*pi*(Rcyl^4-(Rcyl-tcyl)^4)
        stress = mx*dy/I + my*dx/I

        return stress
end

function get_speeds(turbine_x,turbine_y,turbine_z,turb_index,hub_height,r,turbine_yaw,azimuth,turbine_ct,turbine_ai,rotor_diameter,turbine_local_ti,
                    sorted_turbine_index,wtvelocities,wind_resource,model_set;wind_farm_state_id=1,downwind_turbine_id=0)
        x_locs, y_locs, z_locs = find_xyz_simple(turbine_x[turb_index],turbine_y[turb_index],hub_height[turb_index],r,turbine_yaw[turb_index],azimuth)
        npts = length(x_locs)
        point_velocities = zeros(npts)
            for i in 1:npts
                loc = [x_locs[i], y_locs[i], z_locs[i]]
                # point_velocities[i] = FLOWFarm.point_velocity(loc, windfarm, windfarmstate, windresource, wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, 0)

                # point_velocities[i] = FLOWFarm.point_velocity(loc, model_set, problem_description)
                point_velocities[i] = FLOWFarm.point_velocity(loc, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
                                    wind_resource, model_set;
                                    wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=downwind_turbine_id)
            end
        return point_velocities
end

"""

    rainflow(array_ext,uc_mult=0.5)
Rainflow counting of a signal's turning points

# Arguments
        array_ext (numpy.ndarray): array of turning points

# Keyword Arguments
        uc_mult (float): partial-load scaling [opt, default=0.5]

# Returns
        array_out (numpy.ndarray): (3 x n_cycle) array of rainflow values:
                                    1) load range
                                    2) range mean
                                    3) cycle count
"""
function rainflow(array_ext,uc_mult=0.5)

    tot_num = length(array_ext)             # total size of input array
    array_out = zeros(3, tot_num-1)         # initialize output array

    pr = 1                                  # index of input array
    po = 1                                  # index of output array
    j = 0                                   # index of temporary array "a"
    a  = zeros(tot_num)                     # temporary array for algorithm

    # loop through each turning point stored in input array
    for i = 1:tot_num
        j += 1                  # increment "a" counter
        a[j] = array_ext[pr]    # put turning point into temporary array
        pr += 1                 # increment input array pointer

        while j >= 3 && abs( a[j-1] - a[j-2]) <= abs(a[j] - a[j-1])

            lrange = abs( a[j-1] - a[j-2] )

            # partial range
            if j == 3
                mean      = ( a[1] + a[2] ) / 2.0
                a[1]=a[2]
                a[2]=a[3]
                j=2
                if lrange > 0
                    array_out[1,po] = lrange
                    array_out[2,po] = mean
                    array_out[3,po] = uc_mult
                    po += 1
                end

            # full range
            else
                mean      = ( a[j-1] + a[j-2] ) / 2.0
                a[j-2]=a[j]
                j=j-2
                if (lrange > 0)
                    array_out[1,po] = lrange
                    array_out[2,po] = mean
                    array_out[3,po] = 1.00
                    po += 1
                end
            end
        end
    end

    # partial range
    for i = 1:j-1
        lrange    = abs( a[i] - a[i+1] )
        mean      = ( a[i] + a[i+1] ) / 2.0
        if lrange > 0
            array_out[1,po] = lrange
            array_out[2,po] = mean
            array_out[3,po] = uc_mult
            po += 1
        end
    end

    # get rid of unused entries
    out = array_out[:,1:po-1]

    return out
end

"""
    get_peaks(array)

get the turning point values of a signal

# Arguments
- `array::Array{Float}`: the signal to find the turning points, or peaks
"""
function get_peaks(array)
    A = array[:]
    # get rid of any zero slope in the beginning
    while A[2] == A[1]
        A = A[2:length(A)]
    end
    peaks = [A[1]]
    if A[2] > A[1]
        slope = "p"
    elseif A[2] < A[1]
        slope = "m"
    end
    for i = 1:length(A)-2
        ind = i+1
        if slope == "p"
            if A[ind+1] < A[ind]
                peaks = append!(peaks,A[ind])
                slope = "m"
            end
        elseif slope == "m"
            if A[ind+1] > A[ind]
                peaks = append!(peaks,A[ind])
                slope = "p"
            end
        end
    end
    peaks = append!(peaks,A[length(A)])
    return peaks
end

"""
    get_peaks_indices(array)

return the indices of the signal peaks

# Arguments
- `array::Array{Float}`: the signal to find the turning points, or peaks
"""
function get_peaks_indices(array)
    A = array[:]
    # get rid of any zero slope in the beginning
    while A[2] == A[1]
        A = A[2:length(A)]
    end
    peaks = [1]
    if A[2] > A[1]
        slope = "p"
    elseif A[2] < A[1]
        slope = "m"
    end
    for i = 1:length(A)-2
        ind = i+1
        if slope == "p"
            if A[ind+1] < A[ind]
                peaks = append!(peaks,ind)
                slope = "m"
            end
        elseif slope == "m"
            if A[ind+1] > A[ind]
                peaks = append!(peaks,ind)
                slope = "p"
            end
        end
    end
    peaks = append!(peaks,length(A))
    return peaks
end

function calc_TI(constant,ai,TI_free,initial,sep,downstream)

    ti_calculation = constant * (1.0/3.0)^ai * TI_free^initial * sep^downstream

    TI = sqrt(ti_calculation^2 + TI_free^2)

    return TI
end

"""
    get_moments(out,Rhub,Rtip,r,az,precone,tilt)

Trapezoidal integration to find the blade root bending moment using the loads distribution

# Arguments
- `out::CCBlade dict: output from running CCBlade solve
- `Rhub::Float`: radius of the rotor hub
- `Rtip::Float`: radius of the blade tip
- `r::Array{Float}`: radial locations of interest
- `az::Float`: blade azimuth angle
- `precone::Float`: rotor precone angle
- `tilt::Float`: rotor tilt angle
"""
function get_moments(out,Rhub,Rtip,r,az,precone,tilt)
    nr = length(r)+2

    loads_flap = zeros(nr)
    loads_edge = zeros(nr)
    r_arr = zeros(nr)

    loads_flap[2:end-1] = out.Np/1000.0
    loads_edge[2:end-1] = out.Tp/1000.0
    r_arr[2:end-1] = r
    #
    #approximate loads at r = Rhub
    dL_flap = loads_flap[3]-loads_flap[2]
    dL_edge = loads_edge[3]-loads_edge[2]
    dr = r_arr[3]-r_arr[2]

    m_flap = dL_flap/dr
    m_edge = dL_edge/dr

    Lhub_flap = loads_flap[2] + m_flap*(Rhub-r_arr[2])
    Lhub_edge = loads_edge[2] + m_edge*(Rhub-r_arr[2])

    loads_flap[1] = Lhub_flap
    loads_edge[1] = Lhub_edge
    r_arr[1] = Rhub

    #approximate loads at r = Rtip
    dL_flap = loads_flap[end-1]-loads_flap[end-2]
    dL_edge = loads_edge[end-1]-loads_edge[end-2]
    dr = r_arr[end-1]-r_arr[end-2]

    m_flap = dL_flap/dr
    m_edge = dL_edge/dr

    Lhub_flap = loads_flap[end] + m_flap*(Rtip-r_arr[end-1])
    Lhub_edge = loads_edge[end] + m_edge*(Rtip-r_arr[end-1])

    loads_flap[end] = Lhub_flap
    loads_edge[end] = Lhub_edge
    r_arr[end] = Rtip

    M_flap = trapz(r_arr,loads_flap.*(r_arr.-Rhub))
    M_edge = trapz(r_arr,loads_edge.*(r_arr.-Rhub))

    # add gravity loads
    blade_mass=17536.617
    blade_cm=20.650
    grav=9.81
    M_edge += sin(az)*cos(precone)*cos(tilt)*blade_mass*grav*blade_cm/1000.0
    return M_flap, M_edge
end

function get_single_damage(model_set,problem_description,turbine_ID,state_ID,nCycles,az_arr,
    turb_samples,points_x,points_y,omega_func,pitch_func,turbulence_func,r,rotor,sections,Rhub,Rtip,precone,tilt,rho;
    Nlocs=20,fos=1.15)

        naz = length(az_arr)

        turbine_x = problem_description.wind_farm_states[state_ID].turbine_x
        turbine_y = problem_description.wind_farm_states[state_ID].turbine_y
        turbine_z = problem_description.wind_farm_states[state_ID].turbine_z
        nturbines = length(turbine_x)

        naz = length(az_arr)
        ws = problem_description.wind_resource.wind_speeds[state_ID]
        # println("ws: ", ws)
        measurementheight = problem_description.wind_resource.measurement_heights
        air_density = problem_description.wind_resource.air_density
        wind_shear_model = problem_description.wind_resource.wind_shear_model
        ambient_tis = problem_description.wind_resource.ambient_tis
        windfarm = problem_description.wind_farm
        windfarmstate = problem_description.wind_farm_states[state_ID]
        turb_type = problem_description.wind_farm.turbine_definition_ids[turbine_ID]
        hub_height = problem_description.wind_farm.turbine_definitions[turb_type].hub_height
        yaw = problem_description.wind_farm_states[state_ID].turbine_yaw[turbine_ID]

        # flap = []
        # edge = []
        # oms = []
        flap = zeros(nCycles*naz)
        edge = zeros(nCycles*naz)
        oms = zeros(nCycles*naz)

        turbine_inflow_velcities = zeros(nturbines) .+ ws
        temp_resource = ff.DiscretizedWindResource([3*pi/2], [ws], [1.0], measurementheight, air_density,ambient_tis, wind_shear_model)
        temp_pd = ff.WindFarmProblemDescription(windfarm, temp_resource, [windfarmstate])
        ff.turbine_velocities_one_direction!(points_x, points_y, model_set, temp_pd)

        wfs_temp = temp_pd.wind_farm_states[1]
        for i = 1:nCycles*naz
            az = az_arr[(i+1)%naz+1]

            """need to figure this out"""
            x_locs, y_locs, z_locs = ff.find_xyz_simple(turbine_x[turbine_ID],turbine_y[turbine_ID],turbine_z[turbine_ID].+hub_height,r,yaw,az)

            TI_arr = zeros(length(r))
            for k = 1:length(r)
                loc = [x_locs[k],y_locs[k],z_locs[k]]
                # TI_arr[k] = turbulence_func(loc,windfarm,windfarmstate,ambient_tis[state_ID])
                TI_arr[k] = turbulence_func(loc,windfarm,wfs_temp,ambient_tis[state_ID])
            end

            TI_inst = 0.0
            for k = 1:length(r)
                if k == length(r)
                    TI_inst += (Rtip-r[end])*TI_arr[end]/Rtip
                else
                    TI_inst += (r[k+1]-r[k])*TI_arr[k]/Rtip
                end
            end

            turbine_inflow_velcities = temp_pd.wind_farm_states[state_ID].turbine_inflow_velcities
            turbine_inflow_velcities = turbine_inflow_velcities + turb_samples[i]*TI_inst.*turbine_inflow_velcities
            Omega_rpm = omega_func(turbine_inflow_velcities[turbine_ID]) #in rpm
            Omega = Omega_rpm*0.10471975512 #convert to rad/s
            pitch_deg = pitch_func(turbine_inflow_velcities[turbine_ID]) #in degrees
            pitch = pitch_deg*pi/180.0 #convert to rad
            U = get_speeds(turbine_x,turbine_y,turbine_ID,hub_height,r,yaw,az,model_set,temp_pd)
            U = U + turb_samples[i]*TI_inst.*U
            op = distributed_velocity_op.(U, Omega, r, precone, yaw, tilt, az, rho)
            out = CCBlade.solve.(Ref(rotor), sections, op)
            flap[i],edge[i] = ff.get_moments(out,Rhub,Rtip,r,az,precone,tilt)
            oms[i] = Omega

        end

        xlocs = zeros(Nlocs)
        ylocs = zeros(Nlocs)
        angles = range(-pi/2.0,stop=pi/2.0,length=Nlocs)
        root_rad=3.542/2.0
        for i in 1:Nlocs
            xlocs[i] = cos(angles[i])*root_rad
            ylocs[i] = sin(angles[i])*root_rad
        end

        avg_omega = sum(oms)/length(oms)
        total_time = nCycles/(avg_omega/(2.0*pi))

        d_arr = zeros(Nlocs)
        damage = 0.0
        ind = 0
        for i  = 1:Nlocs
            sigma = calc_moment_stress.(edge,flap,xlocs[i],ylocs[i])
            peaks = get_peaks(sigma)
            out = rainflow(peaks)

            alternate = out[1,:]/2.
            mean = out[2,:]
            count = out[3,:]

            su = 70000.0
            m = 10.0
            years = 25.0
            freq = 1.0

            mar = alternate./(1.0.-mean./su)

            npts = length(mar)
            #damage calculations
            d = 0.0
            for i = 1:npts
                Nfail = ((su)/(mar[i]*fos))^m
                mult = years*365.25*24.0*3600.0*freq/total_time
                d += count[i]*mult/Nfail
            end
            d_arr[i] = d
            if d > damage
                damage = d
                ind = i
            end
        end
        # println("angle: ", rad2deg.(angles[ind]))
        return damage
end

function get_single_state_damage(model_set,problem_description,state_ID,nCycles,az_arr,
    turb_samples,points_x,points_y,omega_func,pitch_func,turbulence_func,r,rotor,sections,Rhub,Rtip,precone,tilt,rho;
    Nlocs=20,fos=1.15)

    turbine_x = problem_description.wind_farm_states[state_ID].turbine_x
    nturbines = length(turbine_x)
    damage = zeros(nturbines)
    for k = 1:nturbines
        damage[k] = ff.get_single_damage(model_set,problem_description,k,state_ID,nCycles,az_arr,
            turb_samples,points_x,points_y,omega_func,pitch_func,turbulence_func,r,rotor,sections,Rhub,Rtip,precone,tilt,rho,
            Nlocs=Nlocs,fos=fos)
    end

    return damage
end

function get_total_farm_damage(model_set,problem_description,nCycles,az_arr,
    turb_samples,points_x,points_y,omega_func,pitch_func,turbulence_func,r,rotor,sections,Rhub,Rtip,precone,tilt,rho;
    Nlocs=20,fos=1.15)

    wind_farm = problem_description.wind_farm
    wind_resource = problem_description.wind_resource
    turbine_x = problem_description.wind_farm_states[1].turbine_x
    nturbines = length(turbine_x)

    frequencies = problem_description.wind_resource.wind_probabilities
    ndirections = length(frequencies)

    damage = zeros(nturbines)
    for k = 1:ndirections

        problem_description.wind_farm_states[k].turbine_x[:],problem_description.wind_farm_states[k].turbine_y[:] =
                ff.rotate_to_wind_direction(wind_farm.turbine_x, wind_farm.turbine_y, wind_resource.wind_directions[k])

        problem_description.wind_farm_states[k].sorted_turbine_index[:] = sortperm(problem_description.wind_farm_states[k].turbine_x)

        state_damage = ff.get_single_state_damage(model_set,problem_description,k,nCycles,az_arr,
            turb_samples,points_x,points_y,omega_func,pitch_func,turbulence_func,r,rotor,sections,Rhub,Rtip,precone,tilt,rho,
            Nlocs=Nlocs,fos=fos)
        damage = damage + state_damage.*frequencies[k]
    end

    return damage
end

struct NpTp{T}
    Np::T
    Tp::T
end

function get_single_damage_surr(turbine_x,turbine_y,turbine_z,rotor_diameter,hub_height,turbine_yaw,turbine_ai,sorted_turbine_index,ct_model,turbine_ID,state_ID,nCycles,az_arr,
    turb_samples,flap_func,edge_func,omega_func,pitch_func,turbulence_func,r,rotor,sections,Rhub,Rtip,precone,tilt,rho,wind_resource,model_set;
    Nlocs=20,fos=1.15,rotor_sample_points_y=[0.69,0.69,-0.69,-0.69],rotor_sample_points_z=[0.69,-0.69,0.69,-0.69],div_sigma=2.5,div_ti=1.2)

        naz = length(az_arr)

        nturbines = length(turbine_x)

        naz = length(az_arr)
        ws = wind_resource.wind_speeds[state_ID]
        # println("ws: ", ws)
        measurementheight = wind_resource.measurement_heights
        air_density = wind_resource.air_density
        wind_shear_model = wind_resource.wind_shear_model
        ambient_tis = wind_resource.ambient_tis

        flap = zeros(nCycles*naz)
        edge = zeros(nCycles*naz)
        oms = zeros(nCycles*naz)

        temp_resource = ff.DiscretizedWindResource([3*pi/2], [ws], [1.0], measurementheight, air_density, ambient_tis, wind_shear_model)

        turbine_velocities, turbine_ct, turbine_local_ti = ff.turbine_velocities_one_direction(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
                            turbine_ai, sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                            model_set; wind_farm_state_id=state_ID)

        for i = 1:nCycles*naz
            az = az_arr[(i+1)%naz+1]

            """need to figure this out"""
            x_locs, y_locs, z_locs = ff.find_xyz_simple(turbine_x[turbine_ID],turbine_y[turbine_ID],turbine_z[turbine_ID].+hub_height[turbine_ID],r,turbine_yaw[turbine_ID],az)

            TI_arr = zeros(length(r))
            for k = 1:length(r)
                loc = [x_locs[k],y_locs[k],z_locs[k]]
                # TI_arr[k] = turbulence_func(loc,windfarm,windfarmstate,ambient_tis[state_ID])
                # TI_arr[k] = turbulence_func(loc,windfarm,wfs_temp,ambient_tis[state_ID])
                TI_arr[k] = turbulence_func(loc,turbine_x, turbine_y, rotor_diameter, hub_height, turbine_ct, sorted_turbine_index, ambient_tis[state_ID], div_sigma=div_sigma, div_ti=div_ti)
            end

            TI_inst = 0.0
            for k = 1:length(r)
                if k == length(r)
                    TI_inst += (Rtip-r[end])*TI_arr[end]/Rtip
                else
                    TI_inst += (r[k+1]-r[k])*TI_arr[k]/Rtip
                end
            end

            time1 = time()
            turbine_velocity_with_turb = turbine_velocities[turbine_ID] + turb_samples[i]*TI_inst*turbine_velocities[turbine_ID]
            Omega_rpm = omega_func(turbine_velocity_with_turb) #in rpm
            Omega = Omega_rpm*0.10471975512 #convert to rad/s
            pitch_deg = pitch_func(turbine_velocity_with_turb) #in degrees
            pitch = pitch_deg*pi/180.0 #convert to rad
            U = get_speeds(turbine_x,turbine_y,turbine_z,turbine_ID,hub_height,r,turbine_yaw,az,turbine_ct,turbine_ai,rotor_diameter,turbine_local_ti,
                                sorted_turbine_index,turbine_velocities,wind_resource,model_set;wind_farm_state_id=state_ID,downwind_turbine_id=turbine_ID)
            U = U + turb_samples[i]*TI_inst.*U

            # op = distributed_velocity_op.(U, Omega, r, precone, turbine_yaw[turbine_ID], tilt, az, rho)
            # out_CC = CCBlade.solve.(Ref(rotor), sections, op)
            # println("CC: ", out_CC.Np[1], " ", out_CC.Tp[1])

            np = flap_func(U,turbine_velocity_with_turb)
            tp = edge_func(U,turbine_velocity_with_turb)
            out = NpTp(np,tp)
            # println("spline: ", out.Np[1], " ", out.Tp[1])
            # println("___________________________")
            flap[i],edge[i] = ff.get_moments(out,Rhub,Rtip,r,az,precone,tilt)
            oms[i] = Omega
            println("loop time: ", time()-time1)
        end

        xlocs = zeros(Nlocs)
        ylocs = zeros(Nlocs)
        angles = range(-pi/2.0,stop=pi/2.0,length=Nlocs)
        root_rad=3.542/2.0
        for i in 1:Nlocs
            xlocs[i] = cos(angles[i])*root_rad
            ylocs[i] = sin(angles[i])*root_rad
        end

        avg_omega = sum(oms)/length(oms)
        total_time = nCycles/(avg_omega/(2.0*pi))

        d_arr = zeros(Nlocs)
        damage = 0.0
        ind = 0
        for i  = 1:Nlocs
            sigma = calc_moment_stress.(edge,flap,xlocs[i],ylocs[i])
            peaks = get_peaks(sigma)
            out = rainflow(peaks)

            alternate = out[1,:]/2.
            mean = out[2,:]
            count = out[3,:]

            su = 70000.0
            m = 10.0
            years = 25.0
            freq = 1.0

            mar = alternate./(1.0.-mean./su)

            npts = length(mar)
            #damage calculations
            d = 0.0
            for i = 1:npts
                Nfail = ((su)/(mar[i]*fos))^m
                mult = years*365.25*24.0*3600.0*freq/total_time
                d += count[i]*mult/Nfail
            end
            d_arr[i] = d
            if d > damage
                damage = d
                ind = i
            end
        end
        # println("angle: ", rad2deg.(angles[ind]))
        return damage
end


text:
## Deprecated in CUDA 11.1

struct nvmlDeviceAttributesV1_st
    multiprocessorCount::UInt32
    sharedCopyEngineCount::UInt32
    sharedDecoderCount::UInt32
    sharedEncoderCount::UInt32
    sharedJpegCount::UInt32
    sharedOfaCount::UInt32
end

const nvmlDeviceAttributesV1_t = nvmlDeviceAttributesV1_st

@checked function nvmlDeviceGetAttributes(device, attributes)
    initialize_context()
    ccall((:nvmlDeviceGetAttributes, libnvml()), nvmlReturn_t,
                   (nvmlDevice_t, Ptr{nvmlDeviceAttributesV1_t}),
                   device, attributes)
end

struct nvmlProcessInfoV1_st
    pid::UInt32
    usedGpuMemory::Culonglong
end

const nvmlProcessInfoV1_t = nvmlProcessInfoV1_st

@checked function nvmlDeviceGetComputeRunningProcesses(device, infoCount, infos)
    initialize_context()
    ccall((:nvmlDeviceGetComputeRunningProcesses, libnvml()), nvmlReturn_t,
                   (nvmlDevice_t, Ptr{UInt32}, Ptr{nvmlProcessInfoV1_t}),
                   device, infoCount, infos)
end

@checked function nvmlDeviceGetGraphicsRunningProcesses(device, infoCount, infos)
    initialize_context()
    ccall((:nvmlDeviceGetGraphicsRunningProcesses, libnvml()), nvmlReturn_t,
                   (nvmlDevice_t, Ptr{UInt32}, Ptr{nvmlProcessInfoV1_t}),
                   device, infoCount, infos)
end

## Superseded in CUDA 11.2

struct nvmlComputeInstanceInfoV1_st
    device::nvmlDevice_t
    gpuInstance::nvmlGpuInstance_t
    id::UInt32
    profileId::UInt32
end

const nvmlComputeInstanceInfoV1_t = nvmlComputeInstanceInfoV1_st

@checked function nvmlComputeInstanceGetInfo(computeInstance, info)
    initialize_context()
    ccall((:nvmlComputeInstanceGetInfo, libnvml()), nvmlReturn_t,
                   (nvmlComputeInstance_t, Ptr{nvmlComputeInstanceInfoV1_t}),
                   computeInstance, info)
end

##


text:
using HPAT
using MPI

@acc hpat function calcKmeans(k::Int64, file_name)
    points = DataSource(Matrix{Float64},HDF5,"/points",file_name)
    clusters = Kmeans(points, k,10)
    return clusters
end

K = 3

if (length(ARGS) > 0)
	K = parse(Int, ARGS[1])
end

c2 = calcKmeans(K, HPAT.getDefaultDataPath()*"/kmeans_small.hdf5")
rank = MPI.Comm_rank(MPI.COMM_WORLD)
pes = MPI.Comm_size(MPI.COMM_WORLD)

MPI.Barrier(MPI.COMM_WORLD)
t1 = time_ns()
c2 = calcKmeans(K, HPAT.getDefaultDataPath()*"/kmeans_large.hdf5")
t2 = time_ns()

if rank==0
	println("nodes: ", pes, "\ntime: ", (t2-t1)/1.0e9, "\nCentroids: ", c2)
end

MPI.Barrier(MPI.COMM_WORLD)