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