001/* 002* NACA6Series -- An arbitrary NACA 6 or 6*A series airfoil. 003* 004* Copyright (C) 2010-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; 026 027import static java.lang.Math.*; 028 029 030/** 031* <p> This class represents an arbitrary NACA 6 or 6*A series 032* airfoil section such as a NACA 63-020 airfoil. 033* </p> 034* 035* <p> Ported from FORTRAN "NACA6.FOR" to Java by: 036* Joseph A. Huwaldt, June 3, 2010 </p> 037* 038* <p> Original FORTRAN "NACA4" code had the following note: </p> 039* 040* <pre> 041* AUTHORS - Charles L.Ladson and Cuyler W. Brooks, NASA Langley 042* Liam Hardy, NASA Ames 043* Ralph Carmichael, Public Domain Aeronautical Software 044* Last FORTRAN version: 23Nov96 2.0 RLC 045* 046* NOTES - This program has also been known as LADSON and SIXSERIES and 047* as SIXSERIE on systems with a 8-letter name limit. 048* REFERENCES- NASA Technical Memorandum TM X-3069 (September, 1974), 049* by Charles L. Ladson and Cuyler W. Brooks, Jr., NASA Langley Research Center. 050* 051* "Theory of Wing Sections", by Ira Abbott and Albert Von Doenhoff. 052* </pre> 053* 054* <p> Modified by: Joseph A. Huwaldt </p> 055* 056* @author Joseph A. Huwaldt Date: June 3, 2010 057* @version February 22, 2025 058*/ 059public abstract class NACA6Series implements Airfoil { 060 061 // Constants 062 private static final int NX = 200; 063 private static final int MXCOMB = 10; 064 private static final double TINY = 1.0e-10; 065 private static final double EPS = 0.1e-10; 066 private static final double DX = 0.01; 067 068 // Ordinate equation coefficients. 069 070 /** 071 * Thickness-to-Chord Ratio 072 */ 073 private double toc = 0.20; 074 075 /** 076 * Chord length. Ordinates are scaled to this chord length. 077 */ 078 private double chord = 1; 079 080 private double cli = 0; // Design lift coefficient. 081 082 private double aa = 1; // Mean line chordwise loading (use 0.8 for 6A-series). 083 084 085 // Buffers for the output data. 086 private transient List<Point2D> upper, lower; 087 private transient List<Point2D> camberLine; 088 private transient List<Double> yUp, yLp; 089 090 //----------------------------------------------------------------------------------- 091 092 /** 093 * Create a NACA 6 series airfoil with the specified parameters. 094 * 095 * @param CLi Design lift coefficient (e.g.: 63-206 has CLi = 0.2). 096 * @param thickness The thickness to chord ratio (e.g.: 0.20 == 20% t/c). 097 * @param length The chord length. 098 */ 099 public NACA6Series(double CLi, double thickness, double length) { 100 cli = CLi; 101 chord = length; 102 toc = thickness; 103 } 104 105 //----------------------------------------------------------------------------------- 106 107 /** 108 * Returns a list of points containing the abscissas (X coordinate) and 109 * ordinates (Y coordinate) of the points defining the upper surface of the airfoil. 110 */ 111 @Override 112 public List<Point2D> getUpper() { 113 if (upper == null) 114 calcOrdinatesSlopes(); 115 116 return upper; 117 } 118 119 /** 120 * Returns a list of points containing the abscissas (X coordinate) and 121 * ordinates (Y coordinate) of the points defining the lower surface of the airfoil. 122 */ 123 @Override 124 public List<Point2D> getLower() { 125 if (lower == null) 126 calcOrdinatesSlopes(); 127 128 return lower; 129 } 130 131 /** 132 * Returns a list of points containing the camber line of the airfoil. 133 */ 134 @Override 135 public List<Point2D> getCamber() { 136 if (camberLine == null) 137 calcOrdinatesSlopes(); 138 139 return camberLine; 140 } 141 142 /** 143 * Returns a list containing the slope (dy/dx) of the upper 144 * surface of the airfoil at each ordinate. 145 */ 146 @Override 147 public List<Double> getUpperYp() { 148 if (yUp == null) 149 calcOrdinatesSlopes(); 150 151 return yUp; 152 } 153 154 /** 155 * Returns a list containing the slope (dy/dx) of the lower 156 * surface of the airfoil at each ordinate. 157 */ 158 @Override 159 public List<Double> getLowerYp() { 160 if (yLp == null) 161 calcOrdinatesSlopes(); 162 163 return yLp; 164 } 165 166 /** 167 * Returns a string representation of this airfoil 168 * (the NACA designation of this airfoil). 169 */ 170 @Override 171 public String toString() { 172 StringBuilder buffer = new StringBuilder("NACA "); 173 buffer.append(getProfile()); 174 buffer.append("-"); 175 buffer.append((int)(cli*10)); 176 if (toc < 0.1F) 177 buffer.append("0"); 178 buffer.append((int)(toc*100)); 179 180 if (chord != 1.0F) { 181 buffer.append(" c="); 182 buffer.append(chord); 183 } 184 return buffer.toString(); 185 } 186 187 /** 188 * Returns a String that represents the profile type of this airfoil. 189 * 190 * @return A String that represents the profile type of this airfoil. 191 */ 192 public abstract String getProfile(); 193 194 //----------------------------------------------------------------------------------- 195 196 //----------------------------------------------------------------------------------- 197 198 /** 199 * A method that is called when a new airfoil is instantiated 200 * that calculates the ordinates and slopes of the airfoil. 201 * This method is general to all NACA 4 digit, 5 digit, 16 series 202 * and mod-4 digit airfoils. 203 */ 204 private void calcOrdinatesSlopes() { 205 206 // Allocate memory for buffer arrays. 207 upper = new ArrayList(); 208 lower = new ArrayList(); 209 yUp = new ArrayList(); 210 yLp = new ArrayList(); 211 camberLine = new ArrayList(); 212 213 upper.add(new Point2D.Double(0,0)); 214 lower.add(new Point2D.Double(0,0)); 215 camberLine.add(new Point2D.Double(0,0)); 216 217 boolean if6xa = false; 218 double frac = 1; 219 double clis = cli; 220 int l = 0; 221 222// l = l + 1; 223 double y; 224 double xc = 0; 225 double yc = 0; 226 double xuc; 227 double yuc; 228 double xlc; 229 double ylc; 230 double xl; 231 232 double[] xau = new double[NX]; 233 double[] yau = new double[NX]; 234 double[] xal = new double[NX]; 235 double[] yal = new double[NX]; 236 double[] phi = new double[NX+1]; 237 double[] eps = new double[NX+1]; 238 double[] psi = new double[NX+1]; 239 double[] xt = new double[NX+1]; 240 double[] yt = new double[NX+1]; 241 double[] ytp = new double[NX+1]; 242 double[] bspline = new double[NX+1]; 243 double[] cspline = new double[NX+1]; 244 double[] dspline = new double[NX+1]; 245 double[] sout = new double[4]; 246 247 248 double u = 0.005; 249 double v= -(aa-u)/abs(aa-u); 250 double omxl = (1.-u)*log(1.-u); 251 double amxl = (aa-u)*log(abs(aa-u)); 252 double omxl1 = -log(1.-u) - 1.; 253 double amxl1 = -log(abs(aa-u)) + v; 254 double omxl2 = 1./(1.-u); 255 double amxl2 = -v/abs(aa-u); 256 double ali = 0; 257 258 double h=0, g, q=0, z, z1=0, z2; 259 if ( aa >= EPS && abs(1.-aa) >= EPS ) { 260 g = -(aa*aa*(.5*log(aa)-0.25)+0.25)/(1.-aa); 261 q = 1; 262 double omaa2 = 1.-aa; 263 omaa2 *= omaa2; 264 h = (0.5*omaa2*log(1.-aa)-0.25*omaa2)/(1.-aa) + g; 265 double aamu = aa-u; 266 double omu = 1-u; 267 z = .5*aamu*amxl - .5*omu*omxl - .25*aamu*aamu + .25*omu*omu; 268 z1 = .5*(aamu*amxl1-amxl-omu*omxl1+omxl+aamu-omu); 269 z2 = .5*aamu*amxl2 - amxl1 - .5*omu*omxl2 + omxl1; 270 } 271 272 if ( aa < EPS ) { 273 h = -0.5; 274 q = 1; 275 z1 = u*log(u) - .5*u - .5*(1.-u)*omxl1 + .5*omxl - .5; 276 277 } else if (abs(aa - 1) < EPS) { 278 h = 0; 279 q = 0; 280 z1 = -omxl1; 281 } 282 283 double tanth0 = cli*(z1/(1.-q*aa)-1.-log(u)-h)/PI/(aa+1.)/2.; 284 285 // Slope of profile at origin, Upper and Lower: 286 double yp; 287 double yup=0, ylp=0; 288 if (tanth0 != 0.) { 289 yup = -1./tanth0; 290 ylp = -1./tanth0; 291 } else { 292 yup = ylp = 0; 293 } 294 295 // First station aft of origin on uncambered profile: 296 double x = 0.00025; 297 int i = 0; 298 299 // Start loop for X increment: 300 while ( x <= 1.0000000001) { 301 302 // Skip thickness computation after first pass: 303 if ( i == 0) { 304 phep(phi,eps); 305 phps(phi,psi); 306 307 double rat = 1; 308 int it = 0; 309 double acrat = 1; 310 311 while (true) { 312 // Loop start for thickness iteration: 313 ++it; 314 acrat = acrat*rat; 315 double ymax=0, xym=0; 316 317 for (int j=0; j < NX+1; ++j) { 318 xt[j] = -2.0*cosh(psi[j]*acrat)*cos(phi[j]-eps[j]*acrat); 319 yt[j] = 2.0*sinh(psi[j]*acrat)*sin(phi[j]-eps[j]*acrat); 320 if ( yt[j] > ymax ) { 321 xym = xt[j]; 322 ymax = yt[j]; 323 } 324 } 325 326 // Estimate first,second and third derivatives by spline fit. 327 spline(NX+1, xt, yt, bspline, cspline, dspline); 328 for (int j=0; j < NX+1; ++j) { 329 seval2(NX+1, xt[j], xt, yt, bspline, cspline, dspline, sout); 330 ytp[j] = sout[1]; 331 } 332 333 // Estimate the location of maximum thickness. Look for the 334 // interval where dt/dx changes sign. XTP,YM are x,y of max thickness. 335 double xtp = 1; 336 for (int j=2; j < NX; ++j) { 337 if ( ytp[j] < 0.0 && ytp[j-1] >= 0.0) { 338 xtp = xt[j-1] + ytp[j-1]*(xt[j]-xt[j-1])/(ytp[j-1]-ytp[j]); 339 } 340 } 341 342 double ym = seval(NX+1, xtp, xt, yt, bspline, cspline, dspline); 343 double xo = xt[0]; 344 xl = xt[NX-1]; 345 double tr = 2.*ym/(xl-xo); 346 rat = toc/tr; 347 double sf = rat; 348 349 if (toc <= EPS || abs(rat-1) <= 0.0001 || it > 10) { 350 351 if ( i == 0 ) { 352 for (int j = 0; j < NX+1; ++j) { 353 xt[j] = (xt[j]-xo)/(xl-xo); 354 // Scale linearly to exact thickness: 355 yt[j] = sf*yt[j]/(xl-xo); 356 ytp[j] = sf*ytp[j]; 357 } 358 } 359 360 // Now the XT,YT,YTP array are loaded and never change 361 xtp = (xtp-xo)/(xl-xo); 362 ymax = ymax*sf/(xl-xo); 363 ym = ym*sf/(xl-xo); 364 xym = (xym-xo)/(xl-xo); 365 xl = 0; 366 367 if (toc > EPS) { 368 // Fit tilted ellipse at eleventh profile point: 369 double cn = 2.*ytp[10] - yt[10]/xt[10] + 0.1; 370 double an = xt[10]*(ytp[10]*xt[10]-yt[10])/(xt[10]*(2.*ytp[10]-cn)-yt[10]); 371 double ytmcnxt = yt[10] - cn*xt[10]; 372 double xtman = xt[10] - an; 373 double an2 = an*an; 374 double bn = sqrt(ytmcnxt*ytmcnxt/(1.-xtman*xtman/an2)); 375 double bn2 = bn*bn; 376 for (int j = 0; j < 10; ++j) { 377 xtman = xt[j] - an; 378 yt[j] = bn*sqrt(1.-xtman*xtman/an2) + cn*xt[j]; 379 if ( xt[j] > EPS ) { 380 ytmcnxt = yt[j]-cn*xt[j]; 381 ytp[j] = bn2*(an-xt[j])/an2/ytmcnxt + cn; 382 } 383 } 384 385 // Update spline. 386 spline(NX+1, xt, yt, bspline, cspline, dspline); 387 } 388 389 x = 0.00025; 390 ali = abs(cli); 391 xl = 0; 392 break; 393 } 394 } // End while(true) 395 } // end if (i == 0) 396 397 // Interpolate for thickness and derivatives at desired values of X: 398 spline(NX+1,xt,yt,bspline,cspline,dspline); 399 seval2(NX+1,x,xt,yt,bspline,cspline,dspline,sout); 400 y = sout[0]; 401 yp = sout[1]; 402 403 // Compute camberline: 404 cli = clis; 405 l = 0; 406 407 xc = x*chord; 408 yc = y*chord; 409 double xll = x*log(x); 410 q = 1; 411 if ( abs(1.-aa) < EPS && abs(1.-x) < EPS) { 412 g = h = q = z = 0; 413 z1 = z2 = -10.e10; 414 415 } else if (aa < EPS && (1-x) < EPS ) { 416 g = -0.25; 417 h = -.5; 418 q = 1; 419 z = -.25; 420 z1 = 0; 421 z2 = -10.e10; 422 423 } else if ( abs(aa - x) < EPS) { 424// System.out.println("Here 1"); 425 double omx = 1.-x; 426 z = -.5*omx*omx*log(omx) + 0.25*omx*omx; 427 z1 = -.5*omx*(-log(omx)-1.) + .5*omx*log(omx) - .5*(omx); 428 z2 = -log(omx) - 0.5; 429 g = -(aa*aa*(.5*log(aa)-0.25)+0.25)/(1.-aa); 430 double omaa = 1.-aa; 431 h = (0.5*omaa*omaa*log(omaa)-0.25*omaa*omaa)/omaa + g; 432 433 } else if ( abs(1.-x) < EPS ) { 434// System.out.println("Here 2"); 435 g = -(aa*aa*(.5*log(aa)-0.25)+0.25)/(1.-aa); 436 double omaa = 1.-aa; 437 h = (0.5*omaa*omaa*log(omaa)-0.25*omaa*omaa)/omaa + g; 438 double aam1 = aa-1.; 439 z = .5*aam1*aam1*log(abs(aam1)) - 0.25*aam1*aam1; 440 z1 = -aam1*log(abs(aam1)); 441 z2 = -10.e10; 442 443 } else if ( abs(aa-1.) < EPS ) { 444// System.out.println("Here 3"); 445 g = h = q = 0; 446 double omx = 1.-x; 447 z = -omx*log(omx); 448 z1 = log(omx) + 1.; 449 z2 = -1./omx; 450 451 } else { 452// System.out.println("Here 4"); 453 double aamx = aa-x; 454 double omx = 1.-x; 455 omxl = omx*log(omx); 456 amxl = aamx*log(abs(aamx)); 457 omxl1 = -log(omx) - 1.; 458 amxl1 = -log(abs(aamx)) - 1.; 459 omxl2 = 1./omx; 460 amxl2 = 1./aamx; 461 z = .5*aamx*amxl - .5*omx*omxl - .25*aamx*aamx + .25*omx*omx; 462 z1 = .5*(aamx*amxl1-amxl-omx*omxl1+omxl+aamx-omx); 463 z2 = .5*aamx*amxl2 - amxl1 - .5*omx*omxl2 + omxl1; 464 if ( aa <= EPS ) { 465 g = -0.25; 466 h = -.5; 467 468 } else { 469 double omaa = 1.-aa; 470 g = -(aa*aa*(.5*log(aa)-0.25)+0.25)/omaa; 471 h = (0.5*omaa*omaa*log(omaa)-0.25*omaa*omaa)/omaa + g; 472 } 473 } 474 475 double ycmb = cli*(z/(1.-q*aa)-xll+g-h*x)/PI/(aa+1.)/2.; 476 477 double xsv = x; 478 if ( x < 0.005 ) x = 0.005; 479 double tanth = cli*(z1/(1.-q*aa)-1.-log(x)-h)/PI/(aa+1.)/2.0; 480 x = xsv; 481 if (if6xa) tanth = -2; 482 double ycp2; 483 if (x <= 0.005) 484 ycp2 = 0; 485 else if ( abs(1.-x) <= EPS ) 486 ycp2 = 1./EPS; 487 else { 488 double pia = PI*(aa+1)*2; 489 ycp2 = cli*(z2/(1.-q*aa)-1./x)/pia; 490 } 491 492 // Modified camberline option: 493 if6xa = modCamberLine(x, if6xa, ycmb, tanth, ycp2, cli, sout); 494 ycmb = sout[0]; 495 tanth = sout[1]; 496 ycp2 = sout[2]; 497 498 double f = sqrt(1.+tanth*tanth); 499 double thp = ycp2/(f*f); 500 double sinth = tanth/f; 501 double costh = 1./f; 502 503 // Camberline and derivatives computed: 504 505 // Combine thickness distributuion and camberline: 506 double xu = x - y*sinth; 507 double yu = ycmb + y*costh; 508 xl = x + y*sinth; 509 double yl = ycmb - y*costh; 510 modTrailingEdge(x, xu, yu, xl, yl, sout); 511 yu = sout[0]; 512 yl = sout[1]; 513 514 // Multiply by chord: 515 xuc = xu*chord; 516 yuc = yu*chord; 517 xlc = xl*chord; 518 ylc = yl*chord; 519 520 if (ali > EPS) { 521 // Find local slope of cambered profile: 522 yup = (tanth*f+yp-tanth*y*thp)/(f-yp*tanth-y*thp); 523 ylp = (tanth*f-yp+tanth*y*thp)/(f+yp*tanth+y*thp); 524 } 525 526 ++i; 527 xau[i] = xuc; 528 yau[i] = yuc; 529 xal[i] = xlc; 530 yal[i] = ylc; 531 532 533 // Store scaled profile in appropriate arrays 534 upper.add(new Point2D.Double(xuc,yuc)); 535 if (ali >= TINY) 536 lower.add(new Point2D.Double(xlc,ylc)); 537 else 538 lower.add(new Point2D.Double(xuc,-yuc)); 539 camberLine.add(new Point2D.Double(x*chord,ycmb*chord)); 540 541 542 // Find X increment: 543 if ( x <= 0.012251) frac = 0.025; 544 else if ( x <= 0.09751 ) frac = 0.25; 545 546 // Increment X and return to start of X loop: 547 x += frac*DX; 548 frac = 1; 549 } // end while( x <= 1.0 ) 550 551 // Calculate derivatives of final coordinates; tabulate and save results. 552 ++i; 553 if ( ali >= EPS ) { 554 // output for cambered airfoil 555 spline(i,xau,yau,bspline,cspline,dspline); 556 for (int j=0; j < i; ++j) { 557 seval2(i,xau[j],xau,yau,bspline,cspline,dspline,sout); 558 yUp.add(sout[1]); 559 } 560 spline(i,xal,yal,bspline,cspline,dspline); 561 for (int j=0; j < i; ++j) { 562 seval2(i,xal[j],xal,yal,bspline,cspline,dspline,sout); 563 yLp.add(sout[1]); 564 } 565 566 } else { 567 // uncambered 568 spline(i,xau,yau,bspline,cspline,dspline); 569 for (int j=0; j < i; ++j) { 570 seval2(i,xau[j],xau,yau,bspline,cspline,dspline,sout); 571 yUp.add(sout[1]); 572 yLp.add(-sout[1]); 573 } 574 } 575 576 } 577 578 579 /** 580 * Method used by 6*A series airfoils to modify the camber line. 581 * The default implementation simply returns the inputs as outputs 582 * and returns false. 583 * 584 * @param x The current x/c being calculated. 585 * @param if6xa Flag indicating if 6*A specific calculations have been enabled. 586 * @param ycmb The y position of the camber line. 587 * @param tanth The slope of the camber line. 588 * @param ycp2 ? 589 * @param cli The design lift coefficient. 590 * @param outputs An existing 3-element array filled in with 591 * outputs: outputs[ycmb, tanth, ycp2]. 592 * @return true if this is a 6*A series airfoil, false if it is not. 593 */ 594 protected boolean modCamberLine(double x, boolean if6xa, double ycmb, double tanth, double ycp2, double cli, double[] outputs) { 595 outputs[0] = ycmb; 596 outputs[1] = tanth; 597 outputs[2] = ycp2; 598 return if6xa; 599 } 600 601 602 /** 603 * Method used by 6*A series airfoils to modify the airfoil trailing edge. 604 * The default implementation simply returns the input yu and yl values. 605 * 606 * @param x The x/c location being calculated. 607 * @param xu The x/c location of the upper surface ordinate. 608 * @param yu The upper surface ordinate. 609 * @param xl The x/c location of hte lower surface ordinate. 610 * @param yl The lower surface ordinate. 611 * @param outputs An existing 2-element array filled in with 612 * outputs: outputs[yu, yl]. 613 */ 614 protected void modTrailingEdge(double x, double xu, double yu, double xl, double yl, double[] outputs) { 615 outputs[0] = yu; 616 outputs[1] = yl; 617 } 618 619 620 /** 621 * Fill in phi, eps vectors for 63 series airfoil. 622 * 623 * @param phi An existing array with 201 elements to be filled in 624 * by this method. 625 * @param eps An existing array with 201 elements to be filled in 626 * by this method. 627 */ 628 protected abstract void phep(double[] phi, double[] eps); 629 630 631 /** 632 * Fill in the psi vector for a 63 series airfoil. 633 * 634 * @param phi An array filled in by phep(). 635 * @param psi An existing array with 201 elements to be filled in 636 * by this method. 637 */ 638 protected abstract void phps(double[] phi, double[] psi); 639 640 641 /** 642 * Compute the coefficients of a cubic spline 643 * NOTES - From "Computer Methods for Mathematical Computations", by 644 * Forsythe, Malcolm, and Moser. Prentice-Hall, 1977. 645 * 646 * @param n # of data points (n>1) 647 * @param x abscissas of knots: x[0..n-1] 648 * @param y ordinates of knots: x[0..n-1] 649 * @param b Existing array filled with linear coeff.: b[0..n-1] 650 * @param c Existing array filled with quadratic coeff.: c[0..n-1] 651 * @param d Existing array filled with cubic coeff.: d[0..n-1] 652 */ 653 protected final void spline(int n, double[] x, double[] y, double[] b, double[] c, double[] d) { 654 655 if (n < 3) { 656 // special treatment for n < 3 657 b[0] = 0; 658 if (n == 2) b[0] = (y[1] - y[0])/(x[1] - x[0]); 659 c[0] = 0; 660 d[0] = 0; 661 b[1] = b[0]; 662 c[1] = 0; 663 d[1] = 0; 664 return; 665 } 666 667 // Set up tridiagonal system 668 // b=diagonal, d=offdiagonal, c=right-hand side 669 d[0] = x[1] - x[0]; 670 c[1] = (y[1]-y[0])/d[0]; 671 for (int i=1; i < n-1; ++i) { 672 d[i] = x[i+1] - x[i]; 673 b[i] = 2.0*(d[i-1]+d[i]); 674 c[i+1] = (y[i+1]-y[i])/d[i]; 675 c[i] = c[i+1] - c[i]; 676 } 677 678 // End conditions. third derivatives at x[1] and x[n] obtained 679 // from divided differences. 680 b[0] = -d[0]; 681 b[n-1] = -d[n-2]; 682 c[0] = 0.0; 683 c[n-1] = 0.0; 684 if (n > 3) { 685 c[0] = c[2]/(x[3]-x[1]) - c[1]/(x[2]-x[0]); 686 c[n-1] = c[n-2]/(x[n-1]-x[n-3]) - c[n-3]/(x[n-2]-x[n-4]); 687 c[0] = c[0]*d[0]*d[0]/(x[3]-x[0]); 688 c[n-1] = -c[n-1]*d[n-2]*d[n-2]/(x[n-1]-x[n-4]); 689 } 690 691 for (int i=1; i < n; ++i) { 692 // forward elimination 693 double t = d[i-1]/b[i-1]; 694 b[i] = b[i] - t*d[i-1]; 695 c[i] = c[i] - t*c[i-1]; 696 } 697 698 c[n-1] = c[n-1]/b[n-1]; // back substitution ( makes C the sigma of text) 699 for (int i=n-2; i >= 0; --i) { 700 c[i] = (c[i]-d[i]*c[i+1])/b[i]; 701 } 702 703 // Compute polynomial coefficients 704 b[n-1] = (y[n-1]-y[n-2])/d[n-2] + d[n-2]*(c[n-2]+c[n-1]+c[n-1]); 705 for (int i=0; i < n-1; ++i) { 706 b[i] = (y[i+1]-y[i])/d[i] - d[i]*(c[i+1]+c[i]+c[i]); 707 d[i] = (c[i+1]-c[i])/d[i]; 708 c[i] = 3.0*c[i]; 709 } 710 c[n-1] = 3.0*c[n-1]; 711 d[n-1] = d[n-2]; 712 713 } 714 715 716 /** 717 * Evaluate the cubic spline function. 718 * seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3 719 * where: x(i)=<u<x(i+1) 720 * NOTES - From "Computer Methods for Mathematical Computations", by 721 * Forsythe, Malcolm, and Moser. Prentice-Hall, 1977. 722 * 723 * @param n # of data points (n>1) 724 * @param u abscissa at which the spline is to be evaluated. 725 * @param x abscissas of knots: x[0..n-1] 726 * @param y ordinates of knots: x[0..n-1] 727 * @param b linear coeff.: b[0..n-1] 728 * @param c quadratic coeff.: c[0..n-1] 729 * @param d cubic coeff.: d[0..n-1] 730 * @return The spline function value evaluated at u. 731 */ 732 protected final double seval(int n, double u, double[] x, double[] y, double[] b, double[] c, double[] d) { 733 734 if (n <= 1 ) return y[0]; 735 736 int i = 0; 737 int j = n-1; 738 while (true) { 739 int k = (i+j)/2; 740 if (u < x[k]) 741 j = k; 742 else 743 i = k; 744 if ( j <= i+1) { 745 double dx = u - x[i]; 746 return (y[i] + dx*(b[i] + dx*(c[i] + dx*d[i]))); 747 } 748 } 749 } 750 751 /** 752 * Evaluate the cubic spline function and it's 1st, 2nd and 3rd derivatives. 753 * seval=y(i)+b(i)*(u-x(i))+c(i)*(u-x(i))**2+d(i)*(u-x(i))**3 754 * where: x(i)=<u<x(i+1) 755 * NOTES- if u < x(1), i=0 is used;if u > x(n), i=n-1 is used 756 * 757 * NOTES - From "Computer Methods for Mathematical Computations", by 758 * Forsythe, Malcolm, and Moser. Prentice-Hall, 1977. 759 * 760 * @param n # of data points (n>1) 761 * @param u abscissa at which the spline is to be evaluated. 762 * @param x abscissas of knots: x[0..n-1] 763 * @param y ordinates of knots: x[0..n-1] 764 * @param b linear coeff.: b[0..n-1] 765 * @param c quadratic coeff.: c[0..n-1] 766 * @param d cubic coeff.: d[0..n-1] 767 * @param outputs An existing 2-element array that will be filed with the 768 * the spline function value and it's 1st derivative: 769 * outputs[seval..deriv]. 770 */ 771 private void seval2(int n, double u, double[] x, double[] y, double[] b, double[] c, double[] d, 772 double[] outputs) { 773 774 if (n <= 1) { 775 outputs[0] = y[0]; 776 outputs[1] = 0; 777 return; 778 } 779 780 int i = 0; 781 int j = n-1; 782 while (true) { 783 int k = (i+j)/2; 784 if (u < x[k]) 785 j = k; 786 else 787 i = k; 788 if (j <= i+1) { 789 double dx = u - x[i]; 790 outputs[0] = y[i] + dx*(b[i] + dx*(c[i] + dx*d[i])); 791 outputs[1] = b[i] + dx*(c[i]+c[i]+dx*(d[i]+d[i]+d[i])); 792 return; 793 } 794 } 795 } 796 797} 798 799