1 (* -*- mode: tuareg; -*- *) 2 3 type naca4_num = string 4 type naca4_params = { m : float; p : float; t : float } 5 type coord = Coordinate of { x : float; y : float } 6 type naca4_coords = { c : coord array; u : coord array; l : coord array } 7 type naca4_arfoil = { naca_num : naca4_num; coords : naca4_coords } 8 9 let get_params naca_num = 10 if String.length naca_num <> 4 then 11 raise (Invalid_argument "NACA number must be a 4-digit string"); 12 (* Maximum camber *) 13 let m = String.sub naca_num 0 1 |> float_of_string |> ( *. ) 0.01 in 14 (* Location of maximum camber *) 15 let p = String.sub naca_num 1 1 |> float_of_string |> ( *. ) 0.10 in 16 (* Maximum thickness *) 17 let t = String.sub naca_num 2 2 |> float_of_string |> ( *. ) 0.01 in 18 { m; p; t } 19 20 let y_t t x = 21 5. *. t 22 *. ((0.2969 *. sqrt x) 23 -. (0.1260 *. x) 24 -. (0.3516 *. (x ** 2.)) 25 +. (0.2843 *. (x ** 3.)) 26 -. (0.1015 *. (x ** 4.))) 27 28 let array_of_amt amt = 29 Array.init (amt + 1) (fun x -> float_of_int x /. float_of_int amt) 30 31 let densified_array_of_amt amt = 32 let pct = 0.25 in 33 let dense_amt = int_of_float (pct *. float_of_int amt) in 34 let factor = 5. in 35 let dense_array = 36 Array.init 37 (int_of_float (pct *. amt *. factor)) 38 (fun x -> float_of_int x /. (10. *. factor)) 39 in 40 let sparse_array = 41 let pct = 1. -. pct in 42 Array.init 43 (int_of_float (pct *. (amt))) 44 (fun x -> (float_of_int x +. (pct *. factor)) /. factor) 45 in 46 Array.append dense_array sparse_array 47 48 (* Symmetrical airfoil *) 49 module Symmetrical = struct 50 let center_coord x = Coordinate { x; y = 0. } 51 52 let upper_coord t x = 53 let y = y_t x t in 54 Coordinate { x; y } 55 56 let lower_coord t x = 57 let y = -.y_t x t in 58 Coordinate { x; y } 59 60 let get_coords t xs = 61 { 62 c = Array.map center_coord xs; 63 u = Array.map (upper_coord t) xs; 64 l = Array.map (lower_coord t) xs; 65 } 66 end 67 68 (* Cambered airfoil *) 69 module Cambered = struct 70 let y_c m p x = 71 if 0. <= x && x <= p then m /. (p ** 2.) *. ((2. *. p *. x) -. (x ** 2.)) 72 else if p <= x && x <= 1. then 73 m /. ((1. -. p) ** 2.) *. (1. -. (2. *. p) +. (2. *. p *. x) -. (x ** 2.)) 74 else failwith "x value out of bounds" 75 76 let theta m p x = 77 let dy_c_over_dx = 78 if 0. <= x && x <= p then 2. *. m /. (p ** 2.) *. (p -. x) 79 else if p <= x && x <= 1. then 2. *. m /. ((1. -. p) ** 2.) *. (p -. x) 80 else failwith "x value out of bounds" 81 in 82 atan dy_c_over_dx 83 84 let camber_coord m p x = 85 let y = y_c m p x in 86 Coordinate { x; y } 87 88 let upper_coord m p t x = 89 let y_t = y_t t x in 90 let y_c = y_c m p x in 91 let theta = theta m p x in 92 let x = x -. (y_t *. sin theta) in 93 let y = y_c +. (y_t *. cos theta) in 94 Coordinate { x; y } 95 96 let lower_coord m p t x = 97 let y_t = y_t t x in 98 let y_c = y_c m p x in 99 let theta = theta m p x in 100 let x = x +. (y_t *. sin theta) in 101 let y = y_c -. (y_t *. cos theta) in 102 Coordinate { x; y } 103 104 let get_coords m p t xs = 105 { 106 c = Array.map (camber_coord m p) xs; 107 u = Array.map (upper_coord m p t) xs; 108 l = Array.map (lower_coord m p t) xs; 109 } 110 end 111 112 let get_coords naca_num ?(amt = 100) () = 113 let { m; p; t } = get_params naca_num in 114 let x_coords = densified_array_of_amt amt in 115 let airfoil_is_symmetrical = m = 0. && p = 0. in 116 if airfoil_is_symmetrical then Symmetrical.get_coords t x_coords 117 else Cambered.get_coords m p t x_coords 118 119 let extract_coords coords = 120 let get_x (Coordinate c) = c.x in 121 let get_y (Coordinate c) = c.y in 122 let xs = Array.map get_x coords in 123 let ys = Array.map get_y coords in 124 (xs, ys) 125 126 let plot_curve coords = 127 let get_x (Coordinate c) = c.x in 128 let get_y (Coordinate c) = c.y in 129 let xs = Array.map get_x coords in 130 let ys = Array.map get_y coords in 131 (xs, ys) 132 133 let plot_airfoil naca_num ?(amt = 100) () = 134 let naca_coords = get_coords naca_num ~amt () in 135 let xs, ys = extract_coords naca_coords.c in 136 Plplot.plcol0 1; 137 Plplot.plpoin xs ys 1; 138 let xs, ys = extract_coords naca_coords.u in 139 Plplot.plcol0 6; 140 Plplot.plline xs ys; 141 let xs, ys = extract_coords naca_coords.l in 142 Plplot.plcol0 6; 143 Plplot.plline xs ys 144