summaryrefslogtreecommitdiff
path: root/lib/naca/naca4.ml
blob: d65e58147f172929569c0921aab756ce442e6d8f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
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
Copyright 2019--2025 Marius PETER