Code-Eli  0.3.6
four_digit.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2014 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef four_digit_hpp
14 #define four_digit_hpp
15 
16 #include <string> // std::string
17 #include <sstream> // std::ostringstream, std::istringstream
18 #include <iomanip> // std::setw
19 #include <algorithm> // std::transform
20 
21 #include "eli/code_eli.hpp"
22 
23 namespace eli
24 {
25  namespace geom
26  {
27  namespace curve
28  {
29  namespace pseudo
30  {
31  template<typename data__>
32  class four_digit
33  {
34  public:
35  typedef data__ data_type;
36  typedef Eigen::Matrix<data_type, 1, 2> point_type;
37  typedef Eigen::Matrix<data_type, 5, 1> coefficient_type;
38 
39  public:
40  four_digit() : thickness(00), camber(0), camber_loc(0), sharp_te(false)
41  {
42  recalc_params();
44  }
45 
48  sharp_te(fs.sharp_te)
49  {
50  recalc_params();
52  }
53 
54  coefficient_type get_thickness_coefficients() const {return a;}
55 
56  data_type get_u_min() const {return static_cast<data_type>(-1);}
57  data_type get_u_max() const {return static_cast<data_type>(1);}
58 
60  {
61  sharp_te=fl;
63  }
64  bool sharp_trailing_edge() const {return sharp_te;}
65 
66  bool set_thickness(const data_type &t)
67  {
68  if (t<100)
69  {
70  thickness=t;
71  recalc_params();
72  return true;
73  }
74  return false;
75  }
76  data_type get_thickness() const {return thickness;}
77 
78  bool set_camber(const data_type &cam, const data_type &cam_loc)
79  {
80  if ( (cam<50) && (cam_loc<10) && (cam_loc>0))
81  {
82  camber=cam;
83  camber_loc=cam_loc;
84  recalc_params();
85  return true;
86  }
87 
88  return false;
89  }
90  data_type get_maximum_camber() const {return camber;}
91  data_type get_maximum_camber_location() const {return camber_loc;}
92 
93  bool set_name(const std::string &name)
94  {
95  std::istringstream istr(name);
96  std::string buffer;
97 
98  // get first character and eat preceding white space
99  istr >> std::setw(4) >> buffer;
100  if ((buffer == "NACA") || (buffer == "naca") || (buffer=="Naca"))
101  {
102  unsigned short mm, pp;
103  char buf[3] = "\n\n";
104 
105  istr >> std::setw(4) >> buffer;
106  buf[0]=buffer[0]; mm=std::atoi(buf);
107  buf[0]=buffer[1]; pp=std::atoi(buf);
108 
109  if (set_camber(mm, pp))
110  {
111  buf[0]=buffer[2];
112  buf[1]=buffer[3];
113  if (set_thickness(std::atoi(buf)))
114  return true;
115  }
116  }
117 
118  return false;
119  }
120 
121  std::string get_name() const
122  {
123  std::ostringstream str;
124 
125  str.fill('0');
126  str << "NACA " << std::setw(1) << camber << std::setw(1) << camber_loc << std::setw(2) << thickness;
127  return str.str();
128  }
129 
130  point_type f(const data_type &xi) const
131  {
132  point_type x, xp, xpp;
133 
134  evaluate(x, xp, xpp, xi);
135 
136  return x;
137  }
138 
139  point_type fp(const data_type &xi) const
140  {
141  point_type x, xp, xpp;
142 
143  evaluate(x, xp, xpp, xi);
144 
145  return xp;
146  }
147 
148  point_type fpp(const data_type &xi) const
149  {
150  point_type x, xp, xpp;
151 
152  evaluate(x, xp, xpp, xi);
153 
154  return xpp;
155  }
156 
157  void evaluate(point_type &x, point_type &xp, point_type &xpp, const data_type &xi) const
158  {
159  // check to make sure given valid parametric value
160  assert((xi>=-1) && (xi<=1));
161 
162  data_type xc, xcp, yc, ycp, ycpp, ycppp, yt, ytp, ytpp;
163  const data_type one(1), two(2), three(3);
164  bool lower;
165 
166  // calculate the lower surface
167  if (xi<0)
168  {
169  xcp=-one;
170  xc=-xi;
171  lower=true;
172  }
173  // calculate the upper surface
174  else
175  {
176  if (xi==0)
177  {
178  xcp=0;
179  }
180  else
181  {
182  xcp=one;
183  }
184 
185  xc=xi;
186  lower=false;
187  }
188 
189  // calculate the supporting quantities needed
190  calc_camber(yc, ycp, ycpp, ycppp, xc, lower);
191  calc_thickness(yt, ytp, ytpp, xc, lower);
192 
193  data_type tmp1, cos_theta, cos2_theta, cos3_theta, sin_theta;
194 
195  tmp1=std::sqrt(one+ycp*ycp);
196  cos_theta=one/tmp1;
197  sin_theta=ycp/tmp1;
198  cos2_theta=cos_theta*cos_theta;
199  cos3_theta=cos2_theta*cos_theta;
200 
201  // calculate the info
202  x(0)=xc-yt*sin_theta;
203  x(1)=yc+yt*cos_theta;
204  xp(0)=xcp-ytp*sin_theta-yt*cos3_theta*ycpp;
205  xp(1)=ycp+ytp*cos_theta-yt*cos2_theta*sin_theta*ycpp;
206  xpp(0)=-ytpp*sin_theta-two*ytp*cos3_theta*ycpp
207  +yt*cos2_theta*(three*cos2_theta*sin_theta*ycpp*ycpp-ycppp);
208  xpp(1)=ycpp+ytpp*cos_theta-two*ytp*sin_theta*cos2_theta*ycpp
209  -yt*cos_theta*((two-cos2_theta)*cos2_theta*ycpp*ycpp-sin_theta*ycppp);
210  }
211 
212  point_type tangent(const data_type &xi) const
213  {
214  point_type tgt(fp(xi));
215 
216  tgt.normalize();
217  return tgt;
218  }
219 
220  protected:
222  {
223  data_type ten(10), one_hundred(100);
224 
225  p=camber/one_hundred;
226  m=camber_loc/ten;
227  t=thickness/one_hundred;
228  }
229 
231  {
232  typedef Eigen::Matrix<data_type, 5, 5> coefficient_matrix_type;
233 
234  coefficient_matrix_type coef_mat;
235  coefficient_type rhs;
236  coefficient_type orig_a;
237 
238  // set the specified coefficients
239  orig_a << static_cast<data_type>(0.2969),
240  static_cast<data_type>(-0.1260),
241  static_cast<data_type>(-0.3516),
242  static_cast<data_type>(0.2843),
243  static_cast<data_type>(-0.1015);
244 
245  // if blunt trailing edge, then use specified coefficients
246  if (!sharp_trailing_edge())
247  {
248  a=orig_a;
249  return;
250  }
251 
252  // if want sharp trailing edge then find "actual" constraints that
253  // are placed on thickness distribution
254 
255  // calculate the constraint coefficients
256  calc_four_digit_args(coef_mat, 0, static_cast<data_type>(1)); // (1) trailing edge location
257  calc_four_digit_der_args(coef_mat, 1, static_cast<data_type>(1)); // (2) trailing edge slope
258  calc_four_digit_args(coef_mat, 2, static_cast<data_type>(0.1)); // (3) leading edge shape
259  calc_four_digit_args(coef_mat, 3, static_cast<data_type>(0.3)); // (4) thickness at x/c=0.3 (should be max thickness location, but isn't)
260  calc_four_digit_der_args(coef_mat, 4, static_cast<data_type>(0.3)); // (5) slope at x/c=0.3 (should be zero slope at max thickness, but isn't)
261 
262  // calculate the corresponding constraints for the blunt trailing edge
263  rhs=coef_mat*orig_a;
264 
265  // correct the trailing edge thickness constraint to zero while leaving the rest un changed
266  rhs(0)=static_cast<data_type>(0);
267  a=coef_mat.lu().solve(rhs);
268  }
269 
270  template<typename Derived1>
271  static void calc_four_digit_args(Eigen::MatrixBase<Derived1> &A, const typename Derived1::Index &i, const data_type &xi)
272  {
273  data_type xi2(xi*xi), xi3(xi2*xi), xi4(xi3*xi);
274 
275  A(i,0)=std::sqrt(xi);
276  A(i,1)=xi;
277  A(i,2)=xi2;
278  A(i,3)=xi3;
279  A(i,4)=xi4;
280  }
281 
282  template<typename Derived1>
283  static void calc_four_digit_der_args(Eigen::MatrixBase<Derived1> &A, const typename Derived1::Index &i, const data_type &xi)
284  {
285  data_type xi2(xi*xi), xi3(xi2*xi);
286 
287  A(i,0)=static_cast<data_type>(0.5)/std::sqrt(xi);
288  A(i,1)=static_cast<data_type>(1);
289  A(i,2)=static_cast<data_type>(2)*xi;
290  A(i,3)=static_cast<data_type>(3)*xi2;
291  A(i,4)=static_cast<data_type>(4)*xi3;
292  }
293  void calc_camber(data_type &y, data_type &yp, data_type &ypp, data_type &yppp, const data_type &xi, bool lower) const
294  {
295  // check to make sure given valid parametric value
296  assert((xi>=0) && (xi<=1));
297 
298  data_type zero(0), one(1), two(2);
299 
300  // short circuit if no camber
301  if (camber==0)
302  {
303  y=zero;
304  yp=zero;
305  ypp=zero;
306  yppp=zero;
307 
308  return;
309  }
310 
311  if (xi<=m)
312  {
313  data_type pm2=p/(m*m);
314 
315  y=pm2*(xi*(two*m-xi));
316  yp=two*pm2*(m-xi);
317  ypp=-two*pm2;
318  yppp=zero;
319  }
320  else
321  {
322  data_type p1m2=p/(one+m*(m-two));
323 
324  y=p1m2*(one-two*m+xi*(two*m-xi));
325  yp=two*p1m2*(m-xi);
326  ypp=-two*p1m2;
327  yppp=zero;
328 
329  if (lower)
330  {
331  yp*=-one;
332  yppp*=-one;
333  }
334  }
335  }
336 
337  void calc_thickness(data_type &y, data_type &yp, data_type &ypp, const data_type &xi, bool lower) const
338  {
339  // check to make sure given valid parametric value
340  assert((xi>=0) && (xi<=1));
341 
342  const data_type zero(0), one(1), two(2), three(3), four(4), six(6), twelve(12), half(one/two), quarter(one/four);
343  const data_type xi2(xi*xi), xi3(xi*xi2), xi4(xi2*xi2), sqrtxi(std::sqrt(xi));
344  const data_type trat(t/static_cast<data_type>(0.20));
345 
346  // short circuit for no thickness
347  if (thickness==0)
348  {
349  y=zero;
350  yp=zero;
351  ypp=zero;
352  return;
353  }
354 
355  if (xi==0)
356  {
357  y=zero;
358  yp=one/std::numeric_limits<data_type>::epsilon();
359  ypp=yp;
360  return;
361  }
362  else if ((xi==1) && sharp_trailing_edge())
363  {
364  y=zero;
365  yp=trat*(a.sum()-half*a(0));
366  ypp=trat*(-quarter*a(0)+two*a(2)+six*a(3)+twelve*a(4));
367  if (lower)
368  {
369  ypp*=-one;
370  }
371  return;
372  }
373 
374  y=trat*(a(0)*sqrtxi+a(1)*xi+a(2)*xi2+a(3)*xi3+a(4)*xi4);
375  yp=trat*(half*a(0)/sqrtxi+a(1)+two*a(2)*xi+three*a(3)*xi2+four*a(4)*xi3);
376  ypp=trat*(-quarter*a(0)/sqrtxi/xi+two*a(2)+six*a(3)*xi+twelve*a(4)*xi2);
377 
378  if (lower)
379  {
380  y*=-one;
381  ypp*=-one;
382  }
383  }
384 
385  private:
386  data_type thickness; // thickness index (integer [00,99] with practical limit of [00, 30].
387  // Index is interpreted as 100 times the maximum thickness.
388  data_type camber; // max camber index (integer [0,9]). Index will be
389  // interpreted as 100 times the maximum camber.
390  data_type camber_loc; // chord-wise location of maximum camber index (integer [0,9]).
391  // Index is interpreted as 10 times the maximum camber location.
392 
393  data_type m; // location of maximum camber
394  data_type p; // maximum camber
395  data_type t; // maximum thickness
396  bool sharp_te; // flag to indicate if the trailing edge should be sharp
397  coefficient_type a; // coefficients for thickness distribution
398  };
399  }
400  }
401  }
402 }
403 #endif
point_type f(const data_type &xi) const
Definition: four_digit.hpp:130
data_type get_u_max() const
Definition: four_digit.hpp:57
data_type t
Definition: four_digit.hpp:395
Definition: math.hpp:20
data_type m
Definition: four_digit.hpp:393
static void calc_four_digit_der_args(Eigen::MatrixBase< Derived1 > &A, const typename Derived1::Index &i, const data_type &xi)
Definition: four_digit.hpp:283
point_type tangent(const data_type &xi) const
Definition: four_digit.hpp:212
four_digit(const four_digit< data_type > &fs)
Definition: four_digit.hpp:46
data_type thickness
Definition: four_digit.hpp:386
void calc_thickness(data_type &y, data_type &yp, data_type &ypp, const data_type &xi, bool lower) const
Definition: four_digit.hpp:337
four_digit()
Definition: four_digit.hpp:40
void recalc_coefficients()
Definition: four_digit.hpp:230
bool sharp_te
Definition: four_digit.hpp:396
coefficient_type a
Definition: four_digit.hpp:397
bool set_camber(const data_type &cam, const data_type &cam_loc)
Definition: four_digit.hpp:78
data_type camber_loc
Definition: four_digit.hpp:390
data_type get_thickness() const
Definition: four_digit.hpp:76
Eigen::Matrix< data_type, 5, 1 > coefficient_type
Definition: four_digit.hpp:37
std::string get_name() const
Definition: four_digit.hpp:121
coefficient_type get_thickness_coefficients() const
Definition: four_digit.hpp:54
point_type fpp(const data_type &xi) const
Definition: four_digit.hpp:148
data_type get_maximum_camber() const
Definition: four_digit.hpp:90
data_type camber
Definition: four_digit.hpp:388
bool set_name(const std::string &name)
Definition: four_digit.hpp:93
data_type get_u_min() const
Definition: four_digit.hpp:56
Eigen::Matrix< data_type, 1, 2 > point_type
Definition: four_digit.hpp:36
static void calc_four_digit_args(Eigen::MatrixBase< Derived1 > &A, const typename Derived1::Index &i, const data_type &xi)
Definition: four_digit.hpp:271
void calc_camber(data_type &y, data_type &yp, data_type &ypp, data_type &yppp, const data_type &xi, bool lower) const
Definition: four_digit.hpp:293
void recalc_params()
Definition: four_digit.hpp:221
void evaluate(point_type &x, point_type &xp, point_type &xpp, const data_type &xi) const
Definition: four_digit.hpp:157
Definition: four_digit.hpp:32
void set_sharp_trailing_edge(bool fl)
Definition: four_digit.hpp:59
data__ data_type
Definition: four_digit.hpp:35
data_type get_maximum_camber_location() const
Definition: four_digit.hpp:91
point_type fp(const data_type &xi) const
Definition: four_digit.hpp:139
bool set_thickness(const data_type &t)
Definition: four_digit.hpp:66
bool sharp_trailing_edge() const
Definition: four_digit.hpp:64
data_type p
Definition: four_digit.hpp:394