(* -*- mode: tuareg; -*- *)

type naca4_num = string
type naca4_params = { m : float; p : float; t : float }
type coord = Coordinate of { x : float; y : float }
type naca4_coords = { c : coord array; u : coord array; l : coord array }
type naca4_arfoil = { naca_num : naca4_num; coords : naca4_coords }

let get_params naca_num =
  if String.length naca_num <> 4 then
    raise (Invalid_argument "NACA number must be a 4-digit string");
  (* Maximum camber *)
  let m = String.sub naca_num 0 1 |> float_of_string |> ( *. ) 0.01 in
  (* Location of maximum camber *)
  let p = String.sub naca_num 1 1 |> float_of_string |> ( *. ) 0.10 in
  (* Maximum thickness *)
  let t = String.sub naca_num 2 2 |> float_of_string |> ( *. ) 0.01 in
  { m; p; t }

let y_t t x =
  5. *. t
  *. ((0.2969 *. sqrt x)
     -. (0.1260 *. x)
     -. (0.3516 *. (x ** 2.))
     +. (0.2843 *. (x ** 3.))
     -. (0.1015 *. (x ** 4.)))

let array_of_amt amt =
  Array.init (amt + 1) (fun x -> float_of_int x /. float_of_int amt)

let densified_array_of_amt amt =
  let pct = 0.25 in
  let dense_amt = int_of_float (pct *. float_of_int amt) in
  let factor = 5. in
  let dense_array =
    Array.init
      (int_of_float (pct *. amt *. factor))
      (fun x -> float_of_int x /. (10. *. factor))
  in
  let sparse_array =
    let pct = 1. -. pct in
    Array.init
      (int_of_float (pct *. (amt)))
      (fun x -> (float_of_int x +. (pct *. factor)) /. factor)
  in
  Array.append dense_array sparse_array

(* Symmetrical airfoil *)
module Symmetrical = struct
  let center_coord x = Coordinate { x; y = 0. }

  let upper_coord t x =
    let y = y_t x t in
    Coordinate { x; y }

  let lower_coord t x =
    let y = -.y_t x t in
    Coordinate { x; y }

  let get_coords t xs =
    {
      c = Array.map center_coord xs;
      u = Array.map (upper_coord t) xs;
      l = Array.map (lower_coord t) xs;
    }
end

(* Cambered airfoil *)
module Cambered = struct
  let y_c m p x =
    if 0. <= x && x <= p then m /. (p ** 2.) *. ((2. *. p *. x) -. (x ** 2.))
    else if p <= x && x <= 1. then
      m /. ((1. -. p) ** 2.) *. (1. -. (2. *. p) +. (2. *. p *. x) -. (x ** 2.))
    else failwith "x value out of bounds"

  let theta m p x =
    let dy_c_over_dx =
      if 0. <= x && x <= p then 2. *. m /. (p ** 2.) *. (p -. x)
      else if p <= x && x <= 1. then 2. *. m /. ((1. -. p) ** 2.) *. (p -. x)
      else failwith "x value out of bounds"
    in
    atan dy_c_over_dx

  let camber_coord m p x =
    let y = y_c m p x in
    Coordinate { x; y }

  let upper_coord m p t x =
    let y_t = y_t t x in
    let y_c = y_c m p x in
    let theta = theta m p x in
    let x = x -. (y_t *. sin theta) in
    let y = y_c +. (y_t *. cos theta) in
    Coordinate { x; y }

  let lower_coord m p t x =
    let y_t = y_t t x in
    let y_c = y_c m p x in
    let theta = theta m p x in
    let x = x +. (y_t *. sin theta) in
    let y = y_c -. (y_t *. cos theta) in
    Coordinate { x; y }

  let get_coords m p t xs =
    {
      c = Array.map (camber_coord m p) xs;
      u = Array.map (upper_coord m p t) xs;
      l = Array.map (lower_coord m p t) xs;
    }
end

let get_coords naca_num ?(amt = 100) () =
  let { m; p; t } = get_params naca_num in
  let x_coords = densified_array_of_amt amt in
  let airfoil_is_symmetrical = m = 0. && p = 0. in
  if airfoil_is_symmetrical then Symmetrical.get_coords t x_coords
  else Cambered.get_coords m p t x_coords

let extract_coords coords =
  let get_x (Coordinate c) = c.x in
  let get_y (Coordinate c) = c.y in
  let xs = Array.map get_x coords in
  let ys = Array.map get_y coords in
  (xs, ys)

let plot_curve coords =
  let get_x (Coordinate c) = c.x in
  let get_y (Coordinate c) = c.y in
  let xs = Array.map get_x coords in
  let ys = Array.map get_y coords in
  (xs, ys)

let plot_airfoil naca_num ?(amt = 100) () =
  let naca_coords = get_coords naca_num ~amt () in
  let xs, ys = extract_coords naca_coords.c in
  Plplot.plcol0 1;
  Plplot.plpoin xs ys 1;
  let xs, ys = extract_coords naca_coords.u in
  Plplot.plcol0 6;
  Plplot.plline xs ys;
  let xs, ys = extract_coords naca_coords.l in
  Plplot.plcol0 6;
  Plplot.plline xs ys