1 files changed, 143 insertions, 0 deletions
diff --git a/lib/naca/naca4.ml b/lib/naca/naca4.ml
new file mode 100644
index 0000000..d65e581
--- /dev/null
+++ b/lib/naca/naca4.ml
@@ -0,0 +1,143 @@
+(* -*- 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