# ==============================================================================
# Setup
# ==============================================================================

import FLOWUnsteady as uns
import FLOWUnsteady: vlm, vpm, gt, Im

fidelity = :low
delete_previous_sim = true
nsteps_restart = -1
n_revs = 2 # number of revolutions to simulate
restart_file = nothing
run_name = "sim_12"

nsteps_save = 5
if fidelity == :low
    nsteps_save *= 1
elseif fidelity == :mid
    nsteps_save *= 5
else
    nsteps_save *= 10
end

stride_fluiddomain = max(25, nsteps_save)


maneuver_plot = true
paraview = true # Whether to visualize with Paraview

dir_base = @__DIR__
save_path = joinpath("out", run_name)

if delete_previous_sim
    if isdir(save_path)
        rm(save_path; force=true, recursive=true)
    end
end
mkpath(save_path)

# ==============================================================================
# Operating Point
# ==============================================================================

Vinf(X, t) = [-18.5, 0.0001, 0.0001]

ρ = 1.2
μ = 0.000018
c = 340


Vcruise = 0.0 # (m/s) cruise speed (reference)
RPMh_w = 2000 # (RPM) hover RPM

ttot = n_revs / (RPMh_w / 60)  # (s) total time to perform maneuver

# ==============================================================================
# Vehicle
# ==============================================================================

# Drone parameters
pos_LE = [0.0, 0.0, 0.0]
pos_G_LE = [-0.093, 0.0, -0.055] # center of gravity w.r.t. LE
tiphub_expansion = 1 / 10
m = 22.636 # (kg) mass


# Wing
AR = 5.0
b = 4.95
c_tip = 0.19200
c_hub = 0.32
AR = b / c_tip
n_span = 40
pos = [0, 0.8251649999999999, b / 2] / (b / 2)
clen = [c_hub / c_tip, c_hub / c_tip, c_tip / c_tip]
twist = -[-2.540, -2.540, -2.650]
sweep = [0.0, 0.0]
dihed = [0.0, 0.0]
wing = uns.vlm.complexWing(b, AR, n_span, pos, clen, twist, sweep, dihed; symmetric=true)
wing.Oaxis = gt.rotation_matrix(180, 0, 0)
wing.O = -pos_G_LE


system = uns.vlm.WingSystem()

vlm_system = uns.vlm.WingSystem()
wake_system = uns.vlm.WingSystem()
uns.vlm.addwing(vlm_system, "Wing", wing)
uns.vlm.addwing(wake_system, "SolveVLM", vlm_system)

grids = []

vehicle = uns.VLMVehicle(
    system;
    vlm_system=vlm_system,
    wake_system=wake_system,
    grids=grids
)


# ==============================================================================
# Maneuver
# ==============================================================================

angle_vehicle_x(t) = 20
angle_vehicle_y(t) = 10
angle_vehicle_z(t) = 30
angle_vehicle(t) = [angle_vehicle_x(t), angle_vehicle_y(t), angle_vehicle_z(t)]

Vvehicle(t) = [0, 0, 0]


RPM = ()
angle_tilts = ()


maneuver = uns.KinematicManeuver(angle_tilts, RPM, Vvehicle, angle_vehicle)


# ==============================================================================
# Simulation
# ==============================================================================

use_variable_pitch = true # Whether to use variable pitch in cruise


tstart = 0.00 * ttot # (s) start simulation at this point in time
tquit = 1.00 * ttot # (s) end simulation at this point in time

start_kinmaneuver = true # If true, it starts the maneuver with the
# velocity and angles of tstart.
# If false, starts with velocity=0
# and angles as initiated by the geometry

# -------------------------------------
# Solver
# -------------------------------------

# Aerodynamic solver
VehicleType = uns.UVLMVehicle # Unsteady solver
# VehicleType = uns.QVLMVehicle # Quasi-steady solver

# Time parameters
if fidelity == :low
    nsteps_per_rev = 36 # 10 °/step
elseif fidelity == :mid
    nsteps_per_rev = 72 # 5 °/step
else
    nsteps_per_rev = 360 # 1 °/step
end
resolution_deg = 360 / nsteps_per_rev

nsteps = Int(round(RPMh_w / 60 * 360 * ttot / resolution_deg))
# nsteps = 250 # Time steps for entire maneuver
dt = ttot / nsteps # (s) time step

# VPM particle shedding
p_per_step = 4 # Sheds per time step
shed_starting = fidelity == :low || fidelity == :mid ? false : true # Whether to shed starting vortex
shed_unsteady = true # Whether to shed vorticity from unsteady loading
unsteady_shedcrit = 0.001 # Shed unsteady loading whenever circulation
# fluctuates by more than this ratio

# Regularization of embedded vorticity
r_tip = 0.2
sigma_rotor_surf = fidelity == :low || fidelity == :mid ? r_tip / 10 : r_tip / 80
sigma_vlm_surf = b / 400 # sigma wing 

lambda_vpm = 2.125 # VPM core overlap
# VPM smoothing radius
sigma_vpm_overwrite = lambda_vpm * 2 * pi * r_tip / (nsteps_per_rev * p_per_step)
sigmafactor_vpmonvlm = fidelity == :low || fidelity == :mid ? 1 : 5.5 # Shrink particles by this factor when
# calculating VPM-on-VLM/Rotor induced velocities

# Rotor solver
vlm_rlx = 0.2 # VLM relaxation <-- this also applied to rotors
hubtiploss_correction = vlm.hubtiploss_correction_prandtl # Hub and tip correction


# ==============================================================================
# Wing Solver Settings
# ==============================================================================

# Wing solver: actuator surface model (ASM)
vlm_vortexsheet = fidelity != :low ? true : false # Whether to spread the wing surface vorticity as a vortex sheet (activates ASM)
vlm_vortexsheet_overlap = 2.125 # Overlap of the particles that make the vortex sheet
vlm_vortexsheet_distribution = uns.g_pressure # Distribution of the vortex sheet
chordwing_mean = 0.3
chord = chordwing_mean
thickness_norm = 0.13698
vlm_vortexsheet_sigma_tbv = thickness_norm * chord / 100 # Size of particles in trailing bound vortices
# vlm_vortexsheet_sigma_tbv = sigma_vpm_overwrite
# How many particles to preallocate for the vortex sheet
vlm_vortexsheet_maxstaticparticle = vlm_vortexsheet == false ? nothing : 6000000

# Wing solver: force calculation
KJforce_type = "regular" # KJ force evaluated at middle of bound vortices_vortexsheet also true
include_trailingboundvortex = false # Include trailing bound vortices in force calculations

include_unsteadyforce = true # Include unsteady force
add_unsteadyforce = false # Whether to add the unsteady force to Ftot or to simply output it

include_parasiticdrag = true # Include parasitic-drag force
add_skinfriction = true # If false, the parasitic drag is purely parasitic, meaning no skin friction
calc_cd_from_cl = false # Whether to calculate cd from cl or effective AOA


# VPM solver
vpm_integration = fidelity == :mid || fidelity == :high ? vpm.rungekutta3 : vpm.euler


# Viscosous particles
vpm_viscous = vpm.Inviscid() # VPM viscous diffusion scheme
# vpm_viscous = vpm.CoreSpreading(-1, -1, vpm.zeta_fmm; beta=100.0, itmax=20, tol=1e-1)


# VPM LES subfilter-scale model
if fidelity == :high
    # vpm_SFS = vpm.DynamicSFS(vpm.Estr_fmm, vpm.pseudo3level_positive;
    #     alpha=0.999, maxC=1.0,
    #     clippings=[vpm.clipping_backscatter],
    #     controls=[vpm.control_directional, vpm.control_magnitude]
    # )
    vpm_SFS = vpm.SFS_Cd_twolevel_nobackscatter()
else
    vpm_SFS = vpm.SFS_none
end


if VehicleType == uns.QVLMVehicle
    uns.vlm.VLMSolver._mute_warning(true)
end


# -------------------------------------
# Acoustics 
# -------------------------------------
save_wopwop_input = false
save_init_plots = false



if maneuver_plot
    uns.plot_maneuver(maneuver; ti=0, tf=tquit, save_path=save_path, size_factor=2 / 3)
end


# ==============================================================================
# Monitors
# ==============================================================================

# -------------------- WING MONITORS ---------------------------------------
# Reference parameters for calculating coefficients
# NOTE: make b, ar, and qinf equals to 1.0 to obtain dimensional force
b_ref, ar_ref = 1.0, 1.0
qinf = 1.0
Jref = 1.0

# Force axis labels
CL_lbl = "Lift (N)"
CD_lbl = "Drag (N)"

# Directions of force components
L_dir = [0, 0, 1]
D_dir = [-1, 0, 0]

calc_aerodynamicforce_fun_wing2 = uns.generate_calc_aerodynamicforce(;
    add_parasiticdrag=include_parasiticdrag,
    add_skinfriction=add_skinfriction,
    airfoilpolar="n0012.csv",
)

monitor_name = "Wing"
wingL = vlm.get_wing(vehicle.vlm_system, "Wing")
monitor_wing = uns.generate_monitor_wing(
    wingL, Vinf, b_ref, ar_ref,
    ρ, qinf, nsteps;
    calc_aerodynamicforce_fun=calc_aerodynamicforce_fun_wing2,
    save_path=save_path,
    run_name=monitor_name,
    figname=monitor_name,
    CL_lbl=CL_lbl,
    CD_lbl=CD_lbl,
    L_dir=L_dir,
    D_dir=D_dir,
    # monitor_optargs_wing...
)

runtime_function = monitor_wing


# -------------------- OTHER MONITORS --------------------------------------

# State-variable monitor
statevariable_monitor = uns.generate_monitor_statevariables(; save_path=save_path)
push!(monitors, statevariable_monitor)

# Global enstrophy monitor (numerical stability)
monitor_enstrophy = uns.generate_monitor_enstrophy(; save_path=save_path)
push!(monitors, monitor_enstrophy)

# Monitor of SFS model coefficient Cd
if vpm_SFS != vpm.SFS_none
    monitor_Cd = uns.generate_monitor_Cd(; save_path=save_path)
    push!(monitors, monitor_Cd)
end

# -------------------- CONCATENATE MONITORS --------------------------------
monitors = uns.concatenate(monitors)


# Reference parameters
Vref = Vcruise # Reference velocity to scale maneuver by
RPMref = RPMh_w # Reference RPM to scale maneuver by

# Initial conditions
t0 = tstart / ttot * start_kinmaneuver # Time at start of simulation
Vinit = Vref * maneuver.Vvehicle(t0) # Initial vehicle velocity
Winit = pi / 180 * (maneuver.anglevehicle(t0   1e-12) - maneuver.anglevehicle(t0)) / (ttot * 1e-12) # Initial angular velocity

# Define simulation
simulation = uns.Simulation(
    vehicle,
    maneuver,
    Vref,
    RPMref,
    ttot;
    Vinit=Vinit,
    Winit=Winit,
    t=tstart
);

# Maximum number of particles (for pre-allocating memory)
max_particles = ceil(Int, (nsteps   2) * (2 * vlm.get_m(vehicle.wake_system) * (p_per_step   1)   p_per_step))
max_particles = tquit != Inf ? ceil(Int, max_particles * (tquit - tstart) / ttot) : max_particles
max_particles = min(10000000, max_particles)
max_particles = VehicleType == uns.QVLMVehicle ? 10000 : max_particles

omit_shedding = []


# ==============================================================================
# RUN SIMULATION
# ==============================================================================

@time uns.run_simulation(
    simulation, nsteps;
    # ----- SIMULATION OPTIONS -------------
    Vinf=Vinf,
    rho=ρ, mu=μ, tquit=tquit, sound_spd=c,
    # ----- SOLVERS OPTIONS ----------------
    p_per_step=p_per_step,
    max_particles=max_particles,
    max_static_particles=vlm_vortexsheet_maxstaticparticle,
    vpm_integration=vpm_integration,
    vpm_viscous=vpm_viscous,
    vpm_SFS=vpm_SFS,
    sigma_vlm_surf=sigma_vlm_surf,
    sigma_rotor_surf=sigma_rotor_surf,
    sigma_vpm_overwrite=sigma_vpm_overwrite,
    sigmafactor_vpmonvlm=sigmafactor_vpmonvlm,
    vlm_vortexsheet=vlm_vortexsheet,
    vlm_vortexsheet_overlap=vlm_vortexsheet_overlap,
    vlm_vortexsheet_distribution=vlm_vortexsheet_distribution,
    vlm_vortexsheet_sigma_tbv=vlm_vortexsheet_sigma_tbv,
    vlm_rlx=vlm_rlx,
    hubtiploss_correction=hubtiploss_correction,
    shed_starting=shed_starting,
    shed_unsteady=shed_unsteady,
    unsteady_shedcrit=unsteady_shedcrit,
    omit_shedding=omit_shedding,
    extra_runtime_function=runtime_function,
    # ----- RESTART OPTIONS -----------------
    restart_vpmfile=restart_file,
    # ----- OUTPUT OPTIONS ------------------
    save_path=save_path,
    run_name=run_name,
    save_wopwopin=save_wopwop_input, # <--- Generates input files for PSU-WOPWOP noise analysis
    # save_code=splitdir(@__FILE__)[1]
    nsteps_save=nsteps_save,
    nsteps_restart=nsteps_restart,
    prompt=false,
    create_savepath=false,
);