001/*
002*   NACA4Uncambered -- An arbitrary uncambered NACA 4 digit airfoil.
003*
004*   Copyright (C) 2000-2013, by Joseph A. Huwaldt
005*   All rights reserved.
006*   
007*   This library is free software; you can redistribute it and/or
008*   modify it under the terms of the GNU Lesser General Public
009*   License as published by the Free Software Foundation; either
010*   version 2.1 of the License, or (at your option) any later version.
011*   
012*   This library is distributed in the hope that it will be useful,
013*   but WITHOUT ANY WARRANTY; without even the implied warranty of
014*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
015*   Lesser General Public License for more details.
016*
017*   You should have received a copy of the GNU Lesser General Public License
018*   along with this program; if not, write to the Free Software
019*   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
020*   Or visit:  http://www.gnu.org/licenses/lgpl.html
021**/
022package jahuwaldt.aero.airfoils;
023
024import java.util.List;
025import java.util.ArrayList;
026import java.awt.geom.Point2D;
027import java.text.DecimalFormat;
028import java.text.NumberFormat;
029
030
031/**
032*  <p> This class represents an arbitrary uncambered NACA 4
033*      digit airfoil section such as a NACA 0012 airfoil.
034*      The 1st two digits indicate a symmetric airfoil,
035*      the 2nd two, the thickness to chord ratio.  All NACA 4
036*      and 5 digit (and 16 series and mod-4 digit) airfoils
037*      inherit from this class the generic code for generating
038*      airfoil ordinates.
039*  </p>
040*
041*  <p> Ported from FORTRAN "NACA4.FOR" to Java by:
042*                Joseph A. Huwaldt, October 8, 2000     </p>
043*
044*  <p> Original FORTRAN "NACA4" code had the following note:  </p>
045*
046*  <pre>
047*         AUTHORS - Charles L.Ladson and Cuyler W. Brooks, NASA Langley
048*                   Liam Hardy, NASA Ames
049*                   Ralph Carmichael, Public Domain Aeronautical Software
050*         Last FORTRAN version:  8Aug95  1.7   RLC
051*
052*         NOTES - This program has been known by the names ANALIN, FOURDIGIT and NACA4.
053*         REFERENCES-  NASA Technical Memorandum TM X-3284 (November, 1975),
054*                      "Development of a Computer Program to Obtain Ordinates for
055*                      NACA 4-Digit, 4-Digit Modified, 5-Digit and 16-Digit Airfoils",
056*                      by Charles L. Ladson and Cuyler W. Brooks, Jr.,
057*                      NASA Langley Research Center.
058*
059*                      NASA Technical Memorandum TM 4741 (December 1996),
060*                      "Computer Program to Obtain Ordinates for NACA Airfoils",
061*                      by Charles L. Ladson, Cuyler W. Brooks, Jr., and Acquilla S. Hill,
062*                      NASA Langley Research Center and
063*                      Darrell W. Sproles, Computer Sciences Corporation, Hampton, VA.
064*
065*                      "Theory of Wing Sections", by Ira Abbott and Albert Von Doenhoff.
066*  </pre>
067*
068*  <p>  Modified by:  Joseph A. Huwaldt  </p>
069*
070*  @author  Joseph A. Huwaldt   Date:  October 8, 2000
071*  @version March 7, 2013
072**/
073public class NACA4Uncambered implements Airfoil {
074
075        // Constants
076        protected static final double EPS = 1.E-10;
077        protected static final double BIG = 1./EPS;
078        private static final double kDX = 0.01;
079        
080        // Ordinate equation coefficients.
081        private static final double a0=0.2969, a1=-0.1260, a2=-0.3516;
082        private static final double a3=0.2843, a4=-0.1015;
083
084        /**
085        *  Chord location of maximum thickness.
086        **/
087        protected double xMaxT = 0.50F;
088        
089        
090        // Inputs required by an uncambered NACA 4 digit airfoil.
091        /**
092        *  Thickness-to-Chord Ratio
093        **/
094        protected double TOC = 0.20F;
095        
096        /**
097        *  Chord length.  Ordinates are scaled to this chord length.
098        **/
099        protected double chord = 1;
100        
101        
102        // Buffers for the output data.
103        protected transient List<Point2D> upper, lower;
104        protected transient List<Point2D> camberLine;
105        protected transient List<Double> yUp, yLp;
106        
107        //-----------------------------------------------------------------------------------
108        
109        /**
110        *  Creates an uncambered NACA 4 digit airfoil with a
111        *  thickness to chord ratio of 20% and a chord length of 1.0.
112        **/
113        public NACA4Uncambered() { }
114        
115        /**
116        *  Create an uncambered NACA 4 digit airfoil with the
117        *  specified parameters.
118        *
119        *  @param  thickness  The thickness to chord ratio (e.g.: 0.20 ==> 20% t/c).
120        *  @param  length     The chord length.
121        **/
122        public NACA4Uncambered(double thickness, double length) {
123                chord = length;
124                TOC = thickness;
125        }
126        
127        //-----------------------------------------------------------------------------------
128        
129        /**
130        *  Returns a list of points containing the abscissas (X coordinate) and
131        *  ordinates (Y coordinate) of the points defining the upper surface of the airfoil.
132        **/
133    @Override
134        public List<Point2D> getUpper() {
135                if (upper == null)
136                        calcOrdinatesSlopes();
137
138                return upper;
139        }
140        
141        /**
142        *  Returns a list of points containing the abscissas (X coordinate) and
143        *  ordinates (Y coordinate) of the points defining the lower surface of the airfoil.
144        **/
145    @Override
146        public List<Point2D> getLower() {
147                if (lower == null)
148                        calcOrdinatesSlopes();
149
150                return lower;
151        }
152        
153        /**
154        *  Returns a list of points containing the camber line of the airfoil.
155        **/
156    @Override
157        public List<Point2D> getCamber() {
158                if (camberLine == null)
159                        calcOrdinatesSlopes();
160
161                return camberLine;
162        }
163        
164        /**
165        *  Returns a list containing the slope (dy/dx) of the upper
166        *  surface of the airfoil at each ordinate.
167        **/
168    @Override
169        public List<Double> getUpperYp() {
170                if (yUp == null)
171                        calcOrdinatesSlopes();
172
173                return yUp;
174        }
175        
176        /**
177        *  Returns a list containing the slope (dy/dx) of the lower
178        *  surface of the airfoil at each ordinate.
179        **/
180    @Override
181        public List<Double> getLowerYp() {
182                if (yLp == null)
183                        calcOrdinatesSlopes();
184
185                return yLp;
186        }
187        
188        /**
189        *  Returns a string representation of this airfoil
190        *  (the NACA designation of this airfoil).
191        **/
192    @Override
193        public String toString() {
194                StringBuilder buffer = new StringBuilder("NACA 00");
195                
196                if (TOC < 0.1F)
197                        buffer.append("0");
198                buffer.append((int)(TOC*100));
199                
200                if (chord != 1.0F) {
201                        buffer.append(" c=");
202                        buffer.append(chord);
203                }
204                return buffer.toString();
205        }
206        
207
208        //-----------------------------------------------------------------------------------
209
210        /**
211        *  Method to determine the local slope of the airfoil at
212        *  the leading edge.  This implementation sets the leading
213        *  edge slope to a very large number (essentially infinity)
214        *  as is appropriate for uncambered NACA 4 airfoils.  Sub-
215        *  classes should override this method to provide an appropriate
216        *  leading edge slope.
217        *
218        *  @param  o  The ordinate data structure.  The following are set:  tanth, yp, ypp.
219        **/
220        protected void calcLESlope(Ordinate o) {
221                o.tanth = EPS;
222                o.yp = BIG;
223                o.ypp = BIG;
224        }
225        
226        /**
227        *  Method to calculate the ordinate of the uncambered airfoil
228        *  forward of the maximum thickness point.  This implementation
229        *  is appropriate for uncambered NACA 4 airfoils.  Sub-classes
230        *  should override this method to compute the ordinates forward
231        *  of the max thickness point.
232        *
233        *  @param  x  The x/c location currently being calculated.
234        *  @param  o  The ordinate data structure. The following are set: y, yp, ypp.
235        **/
236        protected void calcOrdinateForward(double x, Ordinate o) {
237                //      Uncambered ordinate.
238                double sqrtx = Math.sqrt(x);
239                double x2 = x*x;
240                o.y = a0*sqrtx + a1*x + a2*x2 + a3*x2*x + a4*x2*x2;
241                o.yp = 0.5*a0/sqrtx + a1 + 2*a2*x + 3*a3*x2 + 4*a4*x2*x;
242                o.ypp = -0.5*0.5*a0/Math.sqrt(x2*x) + 2*a2 + 3*2*a3*x + 4*3*a4*x2;
243        }
244        
245        /**
246        *  Method to calculate the ordinate of the uncambered airfoil
247        *  aft of the maximum thickness point.  This implementation
248        *  simply calls calcOrdinateForward().  Sub-classes should
249        *  override this method to compute the ordinates aft of the
250        *  max thickness point.
251        *
252        *  @param  x  The x/c location currently being calculated.
253        *  @param  o  The ordinate data structure. The following are set: y, yp, ypp.
254        **/
255        protected void calcOrdinateAft(double x, Ordinate o) {
256                calcOrdinateForward(x,o);
257        }
258        
259        /**
260        *  Method to calculate the camber distribution for the airfoil.
261        *  This implementation simply sets camber ralated variables
262        *  to zero (no camber).  Sub-classes should override this method
263        *  to calculate camber.
264        *
265        *  @param  x  The x/c location currently being calculated.
266        *  @param  o  The ordinate data structure. The following are set: yCMB, tanth, thp.
267        *  @return The following function value is returned: sqrt(1 + tanth^2).
268        **/
269        protected double calcCamber(double x, Ordinate o) {
270                //      Uncambered airfoil.
271                o.yCMB = 0;
272                o.tanth = 0;
273                o.thp = 0;
274                
275                double func = 1;
276                return func;
277        }
278
279        //-----------------------------------------------------------------------------------
280        
281        /**
282        *  A method that is called when a new airfoil is instantiated
283        *  that calculates the ordinates and slopes of the airfoil.
284        *  This method is general to all NACA 4 digit, 5 digit, 16 series
285        *  and mod-4 digit airfoils.
286        **/
287        private void calcOrdinatesSlopes() {
288        
289                //      Allocate memory for buffer arrays.
290                upper = new ArrayList();
291                lower = new ArrayList();
292                yUp = new ArrayList();
293                yLp = new ArrayList();
294                camberLine = new ArrayList();
295                
296                //      Create an ordinate data structure for passing around data.
297                Ordinate ord = new Ordinate();
298
299                // Calculate the slope of the leading edge point.
300                calcLESlope(ord);
301                double tanth = ord.tanth;
302                double y;
303                double yp;
304                double ypp = ord.ypp;
305                
306                // Start filling in the data arrays.
307                int i = 0;
308                upper.add(new Point2D.Double(0,0));
309                lower.add(new Point2D.Double(0,0));
310                yUp.add(-1./tanth);
311                yLp.add(-1./tanth);
312                
313                // Generate airfoil upstream of max thickness point.
314                double x = 0.00025;
315                while (x < xMaxT && Math.abs(x-xMaxT) >= EPS) {
316                
317                        // Calculate the ordinate of an uncambered airfoil.
318                        calcOrdinateForward(x,ord);
319                        y = ord.y*TOC/0.2;
320                        yp = ord.yp*TOC/0.2;
321                        ypp = ord.ypp*TOC/0.2;
322                        double xC = x*chord;
323                        double yC = y*chord;
324                        
325                        // Calculate the camber profile offset.
326                        double func = calcCamber(x, ord);
327                        tanth = ord.tanth;
328                        double yCMB = ord.yCMB;
329                        double thp = ord.thp;
330                        double sinth = tanth/func;
331                        double costh = 1./func;
332                        
333                        // Add camber to uncambered ordinate and output.
334                        ++i;
335                        double xUi = (x - y*sinth)*chord;
336                        double yUi = (yCMB + y*costh)*chord;
337                        double xLi = (x + y*sinth)*chord;
338                        double yLi = (yCMB - y*costh)*chord;
339                        upper.add(new Point2D.Double(xUi, yUi));
340                        lower.add(new Point2D.Double(xLi, yLi));
341                        camberLine.add(new Point2D.Double(x*chord, yCMB*chord));
342                        
343                        // Calculate the slope of the airfoil upper and lower surface.
344                        double yUpr;
345                        double yLpr;
346                        if (Math.abs(tanth) >= EPS) {
347                                yUpr = (tanth*func + yp - tanth*y*thp)/(func - yp*tanth - y*thp);
348                                yLpr = (tanth*func - yp + tanth*y*thp)/(func + yp*tanth + y*thp);
349                        } else 
350                                yUpr = yLpr = yp;
351                        
352                        yUp.add(yUpr);
353                        yLp.add(yLpr);
354                        
355                        // Prepare to calculate the next point on the airfoil.
356                        double frac = 1;
357                        if (x <= 0.00225)
358                                frac = 0.025;
359                        else if (x <= 0.0124)
360                                frac = 0.025;
361                        else if (x <= 0.0975)
362                                frac = 0.25;
363                        
364                        x += frac*kDX;
365                }
366                
367                //      Ensure the x = xMaxT is a computed point.
368                x = xMaxT;
369                
370                //      Generate airfoil downstream of max thickness point.
371                while (x <= 1.0) {
372                        calcOrdinateAft(x,ord);
373                        y = ord.y*TOC/0.2;
374                        yp = ord.yp*TOC/0.2;
375                        ypp = ord.ypp*TOC/0.2;
376                        double xC = x*chord;
377                        double yC = y*chord;
378                
379                        double func = calcCamber(x, ord);
380                        tanth = ord.tanth;
381                        double yCMB = ord.yCMB;
382                        double thp = ord.thp;
383                        double sinth = tanth/func;
384                        double costh = 1./func;
385                        
386                        ++i;
387                        double xUi = (x - y*sinth)*chord;
388                        double yUi = (yCMB + y*costh)*chord;
389                        double xLi = (x + y*sinth)*chord;
390                        double yLi = (yCMB - y*costh)*chord;
391                        upper.add(new Point2D.Double(xUi, yUi));
392                        lower.add(new Point2D.Double(xLi, yLi));
393                        camberLine.add(new Point2D.Double(x*chord, yCMB*chord));
394                        
395                        double yUpr;
396                        double yLpr;
397                        if (Math.abs(tanth) >= EPS) {
398                                yUpr = tanth*(func + yp/tanth - y*thp)/(func - yp*tanth - y*thp);
399                                yLpr = tanth*(func - yp/tanth + y*thp)/(func + yp*tanth + y*thp);
400                        } else
401                                yUpr = yLpr = yp;
402                        
403                        yUp.add(yUpr);
404                        yLp.add(yLpr);
405                        
406                        x += kDX;
407                        
408                        //      Ensure that x=1 point is calculated.
409                        if (x > 1 && x < 1 + kDX)       x = 1;
410                }
411        }
412        
413        /**
414        *  Class that serves as a simple container for airfoil ordinate data.
415        **/
416        protected class Ordinate {
417                /**
418                *  Ordinate (Y coordinate) of the uncambered airfoil.
419                **/
420                public double y;
421                
422                /**
423                *  The local slope of the uncambered airfoil:  yp = dy/dx.
424                **/
425                public double yp;
426                
427                /**
428                *  The rate of change of the local slope of the
429                *  uncambered airfoil: ypp = dy^2/dx^2.
430                **/
431                public double ypp;
432                
433                /**
434                *  The offset in the ordinate from the uncambered airfoil
435                *  to the cambered airfoil.  This is the camber profile.
436                **/
437                public double yCMB;
438                
439                /**
440                *  The tangent of the local surface angle due to camber.
441                **/
442                public double tanth;
443                
444                /**
445                *  The local surface angle due to camber.
446                **/
447                public double thp;
448        }
449
450        /**
451        *  Simple method to test this class.
452        **/
453        public static void main(String[] args) {
454        
455                DecimalFormat nf = (DecimalFormat)NumberFormat.getInstance();
456                nf.setMaximumFractionDigits(5);
457                nf.setMinimumFractionDigits(5);
458                
459                System.out.println("Start NACA4Uncambered...");
460                
461                System.out.println("Creating a NACA 0012 airfoil...");
462                Airfoil af = new NACA4Uncambered(0.12, 1);
463                
464                System.out.println("Airfoil = " + af.toString());
465                
466                //      Output the upper surface of the airfoil.
467                List<Point2D> upper = af.getUpper();
468                List<Double> ypArr = af.getUpperYp();
469                
470                System.out.println("        X    \t    Y    \t    dy/dx");
471                int length = upper.size();
472                for (int i=0; i < length; ++i) {
473                        Point2D o = upper.get(i);
474                        System.out.println("    " + nf.format(o.getX()) + "\t" + nf.format(o.getY()) +
475                                                                        "\t" + nf.format(ypArr.get(i)));
476                }
477                
478                System.out.println("# ordinates = " + length);
479                System.out.println("Done!");
480        }
481}
482
483