001/*
002*   NACA4ModUncambered -- An arbitrary uncambered NACA 4 digit modified 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.awt.geom.Point2D;
025import java.text.DecimalFormat;
026import java.text.NumberFormat;
027
028
029/**
030*  <p> This class represents an arbitrary uncambered NACA
031*      modified 4 digit airfoil section such as a NACA 0012-34
032*      airfoil.  The 1st two digits indicate a symmetric airfoil,
033*      the 2nd two, the thickness-chord ratio.  The 1st digit
034*      after the hyphen is the leading edge radius index and
035*      the last digit after the hyphen is the position of
036*      maximum thickness in tenths of chord.
037*  </p>
038*
039*  <p> Ported from FORTRAN "NACA4.FOR" to Java by:
040*                Joseph A. Huwaldt, October 19, 2000     </p>
041*
042*  <p> Original FORTRAN "NACA4" code had the following note:  </p>
043*
044*  <pre>
045*         AUTHORS - Charles L.Ladson and Cuyler W. Brooks, NASA Langley
046*                   Liam Hardy, NASA Ames
047*                   Ralph Carmichael, Public Domain Aeronautical Software
048*         Last FORTRAN version:  8Aug95  1.7   RLC
049*
050*         NOTES - This program has been known by the names ANALIN, FOURDIGIT and NACA4.
051*         REFERENCES-  NASA Technical Memorandum TM X-3284 (November, 1975),
052*                      "Development of a Computer Program to Obtain Ordinates for
053*                      NACA 4-Digit, 4-Digit Modified, 5-Digit and 16-Digit Airfoils",
054*                      by Charles L. Ladson and Cuyler W. Brooks, Jr.,
055*                      NASA Langley Research Center.
056*
057*                      NASA Technical Memorandum TM 4741 (December 1996),
058*                      "Computer Program to Obtain Ordinates for NACA Airfoils",
059*                      by Charles L. Ladson, Cuyler W. Brooks, Jr., and Acquilla S. Hill,
060*                      NASA Langley Research Center and
061*                      Darrell W. Sproles, Computer Sciences Corporation, Hampton, VA.
062*
063*                      "Theory of Wing Sections", by Ira Abbott and Albert Von Doenhoff.
064*  </pre>
065*
066*  <p>  Modified by:  Joseph A. Huwaldt  </p>
067*
068*  @author  Joseph A. Huwaldt   Date:  October 19, 2000
069*  @version February 22, 2025
070*/
071public class NACA4ModUncambered extends NACA4Uncambered {
072
073        /**
074        *  Table relating t/c ratio and LE index to LE radius.
075        */
076        private static final double[][] tableRle = {
077                { 0.000077, 0.000306, 0.000689, 0.001224, 0.001913, 0.002755, 0.003749, 0.004897 },
078                { 0.000110, 0.000441, 0.000992, 0.001763, 0.002755, 0.003967, 0.005399, 0.007052 },
079                { 0.000150, 0.000600, 0.001350, 0.002400, 0.003749, 0.005399, 0.007349, 0.009599 },
080                { 0.000196, 0.000784, 0.001763, 0.003134, 0.004897, 0.007052, 0.009599, 0.012537 },
081                { 0.000248, 0.000992, 0.002231, 0.003967, 0.006198, 0.008925, 0.012148, 0.015867 },
082                { 0.000306, 0.001224, 0.002755, 0.004897, 0.007652, 0.011019, 0.014998, 0.019589 },
083                { 0.000370, 0.001481, 0.003333, 0.005925, 0.009259, 0.013333, 0.018147, 0.023703 },
084                { 0.000441, 0.001763, 0.003967, 0.007052, 0.011019, 0.015867, 0.021597, 0.028208 },
085                { 0.000517, 0.002069, 0.004655, 0.008276, 0.012932, 0.018622, 0.025346, 0.033105 },
086                { 0.000600, 0.002400, 0.005399, 0.009599, 0.014998, 0.021597, 0.029395, 0.038394 },
087                { 0.000689, 0.002755, 0.006198, 0.011019, 0.017217, 0.024792, 0.033745, 0.044075 },
088                { 0.000784, 0.003134, 0.007052, 0.012537, 0.019589, 0.028208, 0.038394, 0.050147 },
089                { 0.000885, 0.003438, 0.007961, 0.014153, 0.022114, 0.031844, 0.043343, 0.056612 },
090                { 0.000992, 0.003967, 0.008925, 0.015867, 0.024792, 0.035701, 0.048592, 0.063468 },
091                { 0.001105, 0.004420, 0.009944, 0.017679, 0.027623, 0.039778, 0.054142, 0.070716 },
092                { 0.001224, 0.004897, 0.011019, 0.019589, 0.030608, 0.044075, 0.059991, 0.078355 },
093                { 0.001350, 0.005399, 0.012148, 0.021597, 0.033745, 0.048592, 0.066140, 0.086387 } };
094        
095        private static final double[] toc = {0.05, 0.06, 0.07, 0.08, 0.09,
096                                0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19,
097                                0.20, 0.21 };
098        
099        /**
100        *  The leading edge radius in percent chord length.
101        */
102        protected double RLE;
103        
104        /**
105        *  The leading edge radius index.
106        */
107        protected int leIndex = 0;
108        
109        /**
110        *  The coefficients of the profile equation.
111        */
112        private double a0, a1, a2, a3;
113        private double d0, d1, d2, d3;                  //      d1 = trailing edge slope.
114        
115        
116        /**
117        *  Create an uncambered NACA modified 4 digit airfoil with the
118        *  specified parameters.  For example:  a NACA 0012-34 airfoil
119        *  translates to:  thickness = 0.12, Rle = 0.003967, xThickness = 0.40.
120        *
121        *  @param  thickness  The thickness to chord ratio (e.g.: 0.09 == 9% t/c).
122        *  @param  Rle        The radius of the leading edge of the airfoil in
123        *                     percent chord.
124        *  @param  xThickness The position of maximum thickness in percent chord
125        *                     (e.g.:  0.40 == 40% t/c).
126        *  @param  length     The chord length.
127        */
128        public NACA4ModUncambered(double thickness, double Rle, double xThickness, double length) {
129                super(thickness, length);
130                
131                xMaxT = xThickness;
132                RLE = Rle;
133                
134                calcCoefficients();
135                
136        }
137        
138        /**
139        *  Create an uncambered NACA modified 4 digit airfoil with the
140        *  specified parameters.  For example:  a NACA 0012-34 airfoil
141        *  translates to:  thickness = 0.12, LEindex = 3, xThickness = 0.40.
142        *
143        *  @param  thickness  The thickness to chord ratio (e.g.: 0.09 == 9% t/c).
144        *  @param  LEindex    The leading edge radius index for this airfoil.  Should
145        *                     be a number from 1 to 8.
146        *  @param  xThickness The position of maximum thickness in percent chord
147        *                     (e.g.:  0.40 == 40% t/c).
148        *  @param  length     The chord length.
149        */
150        public NACA4ModUncambered(double thickness, int LEindex, double xThickness, double length) {
151                super(thickness, length);
152                
153                xMaxT = xThickness;
154                RLE = lookupRle(thickness, LEindex);
155                leIndex = LEindex;
156                
157                calcCoefficients();
158                
159        }
160        
161
162        /**
163        *  Returns a string representation of this airfoil
164        *  (the NACA designation of this airfoil).
165        */
166    @Override
167        public String toString() {
168                StringBuilder buffer = new StringBuilder("NACA 00");
169                
170                if (TOC < 0.1)
171                        buffer.append("0");
172                buffer.append((int)(TOC*100));
173                
174                buffer.append("-");
175                
176                if (leIndex == 0)
177                        buffer.append("?");
178                else
179                        buffer.append(leIndex);
180                
181                buffer.append((int)(xMaxT*10));
182                
183                if (chord != 1.0) {
184                        buffer.append(" c=");
185                        buffer.append(chord);
186                }
187                return buffer.toString();
188        }
189        
190        /**
191        *  Method to calculate the ordinates of the uncambered airfoil
192        *  forward of the maximum thickness point.  This implementation
193        *  computes a set of ordinates that are appropriate for a modified
194        *  NACA 4 digit airfoil.
195        *
196        *  @param  x  The x/c location currently being calculated.
197        *  @param  o  The ordinate data structure. The following are set: y, yp, ypp.
198        */
199    @Override
200        protected void calcOrdinateForward(double x, Ordinate o) {
201                //      Uncambered ordinate.
202                double sqrtx = Math.sqrt(x);
203                double x2 = x*x;
204                o.y = a0*sqrtx + a1*x + a2*x2 + a3*x2*x;
205                o.yp = 0.5*a0/sqrtx + a1 + 2*a2*x + 3*a3*x2;
206                o.ypp = -0.5*0.5*a0/Math.sqrt(x2*x) + 2*a2 + 3*2*a3*x;
207        }
208        
209        /**
210        *  Method to calculate the ordinates of the uncambered airfoil
211        *  aft of the maximum thickness point.  This implementation
212        *  computes a set of ordinates that are appropriate for a modified
213        *  NACA 4 digit airfoil.
214        *
215        *  @param  x  The x/c location currently being calculated.
216        *  @param  o  The ordinate data structure. The following are set: y, yp, ypp.
217        */
218    @Override
219        protected void calcOrdinateAft(double x, Ordinate o) {
220                double oneMx = 1 - x;
221                double oneMx2 = oneMx*oneMx;
222                o.y = d0 + d1*oneMx + d2*oneMx2 + d3*oneMx2*oneMx;
223                o.yp = -d1 - 2*d2*oneMx - 3*d3*oneMx2;
224                o.ypp = 2*d2 + 6*d3*oneMx;
225        }
226        
227        //-----------------------------------------------------------------------------------
228
229        /**
230        *  Calculate the ordinate equation coefficients required
231        *  for a modified, uncambered, NACA 4 digit airfoil.
232        */
233        private void calcCoefficients() {
234                a0 = Math.sqrt(2*RLE)*0.2/TOC;
235                d0 = 0.002;
236                double xM2 = xMaxT*xMaxT;
237                d1 = 0.1*(2.24 - 5.42*xMaxT + 12.3*xM2)/(1 - 0.878*xMaxT);
238                double oneMxM = 1 - xMaxT;
239                double oneMxM2 = oneMxM*oneMxM;
240                d3 = (3*d1 - 0.588/oneMxM)/(3*oneMxM2);
241                d2 = -1.5*oneMxM*d3 - 0.5*d1/oneMxM;
242                a3 = 0.1/xM2/xMaxT + (2*d1*oneMxM - 0.588)/(2*xMaxT*oneMxM2) -
243                                                3*a0/(8*Math.pow(xMaxT, 2.5));
244                a2 = -0.10/xM2 + 0.5*a0/Math.pow(xMaxT, 1.5) - 2*xMaxT*a3;
245                a1 = -0.5*a0/Math.sqrt(xMaxT) - 2*xMaxT*a2 - 3*xM2*a3;
246        }
247
248        /**
249        *  Lookup the appropriate RLE value from the LE index and t/c.
250        */
251        private static double lookupRle(double thickness, int LEindex) {
252                --LEindex;
253                
254                //      Do a range check.
255                if (LEindex < 0)
256                        LEindex = 0;
257                else if (LEindex > 7)
258                        LEindex = 7;
259                
260                //      Find the indexes of the toc values that bound the thickness.
261                int length = toc.length;
262                int i=0;
263                for (; i < length; ++i)
264                        if (toc[i] > thickness)
265                                break;
266                int indexL = 0;
267                if (i == length)
268                        indexL = i - 2;
269                else if (i != 0)
270                        indexL = i - 1;
271                int indexH = indexL + 1;
272                
273                //      Do the interpolation into the table.
274                double value = lineInterp(toc[indexH], tableRle[indexH][LEindex],
275                                                                toc[indexL], tableRle[indexL][LEindex], thickness);
276                
277                return value;
278        }
279        
280        /**
281        *  Simple linear 1D interpolation between two points.
282        *
283        *  @param   x1,y1  Coordinates of the 1st point (the high point).
284        *  @param   x2,y2  Coordinates of the 2nd point (the low point).
285        *  @param   x      The X coordinate of the point for which we want to
286        *                  interpolate to determine a Y coordinate.  Will
287        *                  extrapolate if X is outside of the bounds of the 
288        *                  point arguments.
289        *  @return  The interpolated Y value corresponding to the input X
290        *           value is returned.
291        */
292        private static double lineInterp( double x1, double y1, double x2, double y2, double x ) {
293                return ((y2 - y1)/(x2 - x1)*(x - x1) + y1);
294        }
295
296
297        /**
298        *  Simple method to test this class.
299        *
300        * @param args the command-line arguments
301        */
302        public static void main(String[] args) {
303        
304                DecimalFormat nf = (DecimalFormat)NumberFormat.getInstance();
305                nf.setMaximumFractionDigits(5);
306                nf.setMinimumFractionDigits(5);
307                
308                System.out.println("Start NACA4ModUncambered...");
309                
310                System.out.println("Creating a NACA 0010-64 airfoil...");
311                Airfoil af = new NACA4ModUncambered(0.10, 6, 0.4, 1);
312                
313                System.out.println("Airfoil = " + af.toString());
314                
315                //      Output the upper surface of the airfoil.
316                List<Point2D> upper = af.getUpper();
317                List<Double> ypArr = af.getUpperYp();
318                
319                System.out.println("        X    \t    Y    \t    dy/dx");
320                int length = upper.size();
321                for (int i=0; i < length; ++i) {
322                        Point2D o = upper.get(i);
323                        System.out.println("    " + nf.format(o.getX()) + "\t" + nf.format(o.getY()) +
324                                                                        "\t" + nf.format(ypArr.get(i)));
325                }
326                
327                System.out.println("# ordinates = " + length);
328                System.out.println("Done!");
329        }
330        
331}
332
333