JuliaCode_complete

function IFR(surf::TwoDSurf,x,z,t::Float64)
    # Given body frame and kinematics, find global positions
    alpha = surf.kindef.alpha(t)
    X0 = -surf.kindef.u(t)*t
    Z0 = surf.kindef.h(t)
    pvt = surf.pvt

    X = (x .- pvt).*cos(alpha) + z.*sin(alpha) .+ X0 .+ pvt
    Z = -(x .- pvt).*sin(alpha) + z.*cos(alpha) .+ Z0

    return X,Z
end

function BFR(surf::TwoDSurf,X,Z,t::Float64)
    # Given inertial frame and kinematics, find body positions
    alpha = surf.kindef.alpha(t)
    X0 = -surf.kindef.u(t)*t
    Z0 = surf.kindef.h(t)
    pvt = surf.pvt

    x = (X .- X0 .- pvt).*cos(alpha) - (Z .- Z0).*sin(alpha) .+ pvt
    z = (X .- X0 .- pvt).*sin(alpha) + (Z .- Z0).*cos(alpha)

    return x,z
end


function refresh_vortex(surf::TwoDSurf,vor_loc)
    # Updates vortex locations to vor_loc
    for i=1:length(surf.bv)
        surf.bv[i].x = vor_loc[i,1]
        surf.bv[i].z = vor_loc[i,2]
    end
    return surf
end

function place_tev2(surf::TwoDSurf,field::TwoDFlowField,dt)
    vorx = map(q -> q.x,surf.bv)
    vorz = map(q -> q.z,surf.bv)
    if size(field.tev,1) == 0
        xloc = surf.bv[end].x + 0.5*surf.kinem.u*dt
        zloc = surf.bv[end].z
    else
        xloc = surf.bv[end].x + (1. /3.)*(field.tev[end].x - surf.bv[end].x)
        zloc = surf.bv[end].z + (1. /3.)*(field.tev[end].z - surf.bv[end].z)
    end
    push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))
    return field
end

function vor_BFR(surf::TwoDSurf,t,IFR_field::TwoDFlowField,curfield::TwoDFlowField)
    # finds and updates the BFR positions of the vortexes in IFR_field to
    #   curfield
    ntev = size(IFR_field.tev,1)
    nlev = size(IFR_field.lev,1)
    vorX = [map(q -> q.x,IFR_field.tev); map(q -> q.x,IFR_field.lev)]
    vorZ = [map(q -> q.z,IFR_field.tev); map(q -> q.z,IFR_field.lev)]

    # Body frame conversion
    vorx, vorz = BFR(surf,vorX,vorZ,t)

    # Update bfr TEV
    for i = 1:ntev
        curfield.tev[i].x = vorx[i]
        curfield.tev[i].z = vorz[i]
    end
    # Update bfr LEV
    for i = 1:nlev
        curfield.lev[i].x = vorx[i+ntev]
        curfield.lev[i].z = vorz[i+ntev]
    end
    return curfield
end

function influence_coeff(surf::TwoDSurf,curfield::TwoDFlowField,coll_loc,n,dt)
    # With the surface and flowfield, determine the influence matrix "a"
    a = zeros(surf.ndiv,surf.ndiv)
    a[end,:] .= 1. # for wake portion
    vc = surf.bv[1].vc

    if size(curfield.tev,1) == 0 # calculate TEV placement location
        xloc = surf.bv[end].x + 0.5*surf.kinem.u*dt*cos(surf.kinem.alpha)
        zloc = surf.bv[end].z + 0.5*surf.kinem.u*dt*sin(surf.kinem.alpha)
    else
        xloc = surf.bv[end].x + (1. /3.)*(curfield.tev[end].x - surf.bv[end].x)
        zloc = surf.bv[end].z + (1. /3.)*(curfield.tev[end].z - surf.bv[end].z)
    end

    for i = 1:surf.ndiv-1, j = 1:surf.ndiv
        # i is test location, j is vortex source
        t_x = coll_loc[i,1]
        t_z = coll_loc[i,2]

        if j < surf.ndiv # Bound vortices (a_ij)
            src = surf.bv[j]
            xdist = src.x .- t_x
            zdist = src.z .- t_z
        else # Wake vorticy (a_iw)
            xdist = xloc .- t_x
            zdist = zloc .- t_z
        end

        distsq = xdist.*xdist + zdist.*zdist
        u = -zdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))
        w = xdist./(2*pi*sqrt.(vc*vc*vc*vc .+ distsq.*distsq))

        a[i,j] = u*n[i,1] + w*n[i,2]
    end
    return a
end

function LVE(surf::TwoDSurf,curfield::TwoDFlowField,nsteps::Int64 = 500,dtstar::Float64 = 0.015,frameCnt = 20,view = "square")
    ## Initialize variables
    mat = zeros(8 ,0) # loads output matrix
    prevCirc = 0
    vor_loc = zeros(surf.ndiv-1,2) # (x1,z1 ; x2,z2 ...) Inertial frame
    coll_loc = zeros(surf.ndiv-1,2) # (x1,z1 ; x2,z2 ...)
    global RHS = zeros(surf.ndiv)
    IFR_vor = 0 .* vor_loc # Global frame
    IFR_field = TwoDFlowField()
    IFR_surf = TwoDSurf(surf.coord_file,surf.pvt,surf.kindef,surf.lespcrit; ndiv = surf.ndiv, camberType = "linear")
    dp = zeros(surf.ndiv-1,1) # change in panel pressure
    aniStep = round((nsteps-1) / frameCnt) # Frames to capture animation on
    aniStep < 1 ? aniStep = 1 : aniStep
    frames = 0
    prevNow = 0
    tevCirc = []

    #   length of each panel
    ds = sqrt.(( surf.x[1:end-1]-surf.x[2:end]).^2 + (surf.cam[1:end-1]-surf.cam[2:end]).^2)
    #   Surf cam_slope does not give the correct panel locations
    cam_slope = asind.((surf.cam[2:end]-surf.cam[1:end-1])./ds) # [deg]

    # Normal vectors for each panel
    n = hcat(-sind.(cam_slope),cosd.(cam_slope))
    # Tangential vectors for each panel
    tau = hcat(cosd.(cam_slope),sind.(cam_slope))

    ## Vortex Locations
    # Located at 25% of each panel
    #   Vortex loactions (at 25% of panel length)
    vor_loc[:,1] = surf.x[1:end-1] + .25 .* ds .* cosd.(cam_slope)
    vor_loc[:,2] = surf.cam[1:end-1] + .25 .* ds .* sind.(cam_slope)
    #   Collocation points (at 75% panel chordwise)
    coll_loc[:,1] = surf.x[1:end-1] + .75 .* ds .* cosd.(cam_slope)
    coll_loc[:,2] = surf.cam[1:end-1] + .75 .* ds .* sind.(cam_slope)

    ## Refresh vortex values
    refresh_vortex(surf,vor_loc)

    ## Time Stepping Loop
    t = -dtstar

    for i = 0:nsteps
        t += dtstar

        # Update Kinematics
        update_kinem(surf, t)

        ## Influence Coeffcients
        a = influence_coeff(surf,curfield,coll_loc,n,dtstar)

        ## RHS Vector
        RHS[end] = prevCirc # previous total circulation
        # u_w, w_w
        if length(curfield.tev) > 0 # Trailing edge vortexs exist
            surf.uind[1:end-1], surf.wind[1:end-1] = ind_vel([curfield.tev; curfield.lev],coll_loc[:,1],coll_loc[:,2])
            #surf.uind[end-1] = surf.uind[end-2]
            #surf.wind[end-1] = surf.wind[end-2]
        end
        for j = 1:surf.ndiv-1
            alpha = surf.kinem.alpha
            # relVel = [U(t) ; W(t)]
            global relVel = [cos(alpha) -sin(alpha) ; sin(alpha) cos(alpha)] * [surf.kinem.u ; -surf.kinem.hdot] + [-surf.kinem.alphadot*coll_loc[j,2] ; surf.kinem.alphadot * (coll_loc[j,1] - surf.pvt)]

            # RHS = -[U_t + u_w , W_t + w_w] dot n
            global RHS[j] = -((relVel[1] + surf.uind[j])*n[j,1] + (relVel[2] + surf.wind[j])*n[j,2])

        end

        ## Vortex Strenghths a*circ = RHS
        global circ = a\RHS

        circChange = sum(circ[1:end-1]) - prevCirc #bound vorticy circ change
        prevCirc = sum(circ[1:end-1])

        for j = 1:surf.ndiv-1 # Set bv circs
            surf.bv[j].s = circ[j]
        end
        # Inertial reference frame
        IFR_vor[:,1],IFR_vor[:,2] = IFR(surf,vor_loc[:,1],vor_loc[:,2],t)
        IFR_surf = refresh_vortex(IFR_surf,IFR_vor)

        ## TEV setup
        place_tev2(IFR_surf,IFR_field,dtstar)
        if length(IFR_field.tev) > 1
            IFR_field.tev[end].s = circ[end]
        end

        # Position mesh
        if i % aniStep == 0.

            temp = Int(i / aniStep) + 1
            now = Dates.format(Dates.now(), "HH:MM")
            if prevNow == now
                println("$temp / $frameCnt")
            else
                println("$temp / $frameCnt,  $now")
            end
            prevNow = now
            # Velocities at mesh points
            vorts = IFR_field.tev[:]
            global mesh = meshgrid(surf,vorts,.25,t,100,view)
            vorts = [vorts ; IFR_surf.bv]
            for j = 1:size(vorts,1)
                u,w = ind_vel(vorts[j],mesh.x,mesh.z)
                global mesh.uMat = mesh.uMat .+ u
                global mesh.wMat = mesh.wMat .+ w
                global mesh.velMag = sqrt.(mesh.uMat.^2 + mesh.wMat.^2)
                global mesh.t = t
                global mesh.circ = circ
            end
            if frames == 0
                frames = [mesh;]
            else
                push!(frames,mesh)
            end
        end
        if i > 0
            ## Pressures (Change in pressure for each panel)
            for j = 1:length(dp)
                dp[j] = surf.rho * ( ( (relVel[1] + surf.uind[j])*tau[j,1] + (relVel[2] + surf.wind[j])*tau[j,2] ) * circ[j]/ds[j] + circChange )
            end
            ## Loads
            # Cn normal
            ndem = .5 * surf.rho * surf.kinem.u^2 * surf.c # non-dimensionalization constant
            cn = sum( dp[j]*ds[j] for j = 1:length(dp) ) / ndem
            # LESP
            x = 1 - 2*ds[1] / surf.c
            surf.a0[1] = 1.13 * circ[1] / ( surf.kinem.u * surf.c * ( acos(x) + sin(acos(x) )))
            # Cs suction
            cs = 2*pi*surf.a0[1]^2
            # Cl lift
            alpha = surf.kinem.alpha
            cl = cn*cos(alpha) + cs*sin(alpha)
            # Cd drag
            cd = cn*sin(alpha) - cs*cos(alpha)
            # Cm moment
            ndem = ndem * surf.c # change ndem to moment ndem
            cm = -sum( dp[j]*ds[j]*(vor_loc[j,1]-surf.pvt) for j = 1:length(dp) ) / ndem

            mat = hcat(mat,[t, surf.kinem.alpha*180/pi, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm]) # Pressures and loads

        end
        ## Wake Rollup
        IFR_field = wakeroll(IFR_surf,IFR_field,dtstar)

        # Convert IFR_field values to curfield (IFR to BFR)
        push!(curfield.tev,TwoDVort(0,0,IFR_field.tev[end].s,0.02*surf.c,0.,0.))
        vor_BFR(surf,t,IFR_field,curfield)

        tevCirc = [tevCirc ; prevCirc]

    end
    mat = mat
    test = tevCirc
    return frames,IFR_field,mat,test
end

#### MAIN File
include("src/UnsteadyFlowSolvers.jl")
using Revise
## Define Geometry and Kinematics
# Kinematics
case = 2 # kinematics case to run
if case == 1
    amp = 35 #degrees
    k = .005
    tstart = 60 # non-dimensional time
elseif case == 2
    amp = 45
    k = .05
    tstart = 20
elseif case == 3
    amp = 90
    k = .4
    tstart = 10
elseif case == 4
    amp = 25
    k = .11
    a = 11
    tstart = 1
end
#amp = amp * pi / 180
if case != 4
a = (pi^2 * k*180 )/(2*amp*pi *(1-0.1))
end

alpha = 10
alphadef = UnsteadyFlowSolvers.ConstDef(alpha *pi/180)#EldUptstartDef(amp,k,a,1)#tstart)
hdef = UnsteadyFlowSolvers.ConstDef(0)
udef = UnsteadyFlowSolvers.ConstDef(1)
full_kinem = UnsteadyFlowSolvers.KinemDef(alphadef, hdef, udef)
# Geometry
pvt = 0.25 ;
geometry = "FlatPlate"#"bin/sd7003.dat"
lespcrit = [10000;]
surf = UnsteadyFlowSolvers.TwoDSurf(geometry,pvt,full_kinem,lespcrit; ndiv = 70, camberType = "linear")
curfield = UnsteadyFlowSolvers.TwoDFlowField()
# Iteration Parameters
dtstar = UnsteadyFlowSolvers.find_tstep(alphadef)

#=Solve for time needed to acheive maximum AoA with a proportional loop
global t = tstart
global alpha = alphadef(t) *180/pi
P = 10
while alpha <= amp - .05 || alpha >= amp + .05 # Timestep is outside .05 degrees
    global alpha = alphadef(t) *180/pi
    error = amp - alpha
    global t += P*error
    #println(alpha )#* 180 / pi)
    if t > 10000
        break
    end
end
t_tot = ceil(t)
=#
t_tot = 9
nsteps = Int(round(t_tot/dtstar))
nsteps = 40 #DEBUG

println("Running LVE")
frames, ifr_field, mat, test = UnsteadyFlowSolvers.LVE(surf,curfield,nsteps,dtstar,40,"longwake")
#circ,a,RHS = UnsteadyFlowSolvers.LVE(surf,curfield,nsteps,dtstar)

#=
println("Running LDVM")
startflag = 0
writeflag = 0
writeInterval = t_tot/18.
surf2 = UnsteadyFlowSolvers.TwoDSurf(geometry,pvt,full_kinem,lespcrit; ndiv = 4)
curfield2 = UnsteadyFlowSolvers.TwoDFlowField()
delvort = UnsteadyFlowSolvers.delNone()

mat2, surf2, curfield2 = UnsteadyFlowSolvers.ldvm(surf2, curfield2, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)
#mat2 = UnsteadyFlowSolvers.ldvm(surf2, curfield2, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)
=#
# mat = [t , alpha , h , u , LESP , cl , cd , cm] for each timestep

# Velocity plot
using Plots
pyplot()
single = "false" ;
if single == "true"
    mesh = frames[end]

    contour(mesh.x[1,:],mesh.z[:,1], mesh.velMag,fill=true,levels = 200,c = :lightrainbow_r)
    scatter!(mesh.vorX, mesh.vorZ,markersize = 1.5, markerstrokestyle = :dot, markerstrokecolor = :white)
    plot!(mesh.camX,mesh.camZ,color = :black, linewidth = 3, aspect_ratio = :equal,legend = :none,axis = false)
elseif single == "false" # Make full list of images in folder
    dir = "case_$case steps_$nsteps"
    mkdir(dir)
    for i = 1:length(frames)
        mesh = frames[i]

        contour(mesh.x[1,:],mesh.z[:,1], mesh.velMag,fill=true,levels = 200,c = :lightrainbow_r)
        scatter!(mesh.vorX, mesh.vorZ,markersize = 1.5, markerstrokestyle = :dot, markerstrokecolor = :white)
        plot!(mesh.camX,mesh.camZ,color = :black, linewidth = 3, aspect_ratio = :equal,legend = :none,axis = false)

        savefig("$dir/$i")
    end
end
#
#= Plot circ vs alpha
circ = map(q -> q.circ[1],frames)
alpha = map(q -> q.alpha,frames)
plot(alpha,circ)
circ1 = frames[1].circ
x = (1:length(circ1))/length(circ1)
plot(x,circ1)
=#

#= Force plots
#UnsteadyFlowSolvers.makeForcePlots2D()

plot(mat[:,1],mat[:,6],color = :black) # cl vs AoA
plot!(mat2[:,1],mat2[:,6],color = :red)
=#
alpha = Int(alpha)
#plot(mat2[:,1],mat2[:,2], xlabel = "t*", ylabel = "AoA")
#savefig("RameshData_Motion_$alpha")

UnsteadyFlowSolvers.subPlotForces(mat,1,"t*","MyData_$alpha")

plot(mat[:,1],test[1:end-1], xlabel = "t*", ylabel = "TEV Strength")