git.gag.com Git - debian/openrocket/blob - core/src/net/sf/openrocket/aerodynamics/barrowman/FinSetCalc.java

   1 package net.sf.openrocket.aerodynamics.barrowman;
   2
   3 import static java.lang.Math.pow;
   4 import static net.sf.openrocket.util.MathUtil.pow2;
   5
   6 import java.util.Arrays;
   7
   8 import net.sf.openrocket.aerodynamics.AerodynamicForces;
   9 import net.sf.openrocket.aerodynamics.FlightConditions;
  10 import net.sf.openrocket.aerodynamics.Warning;
  11 import net.sf.openrocket.aerodynamics.WarningSet;
  12 import net.sf.openrocket.rocketcomponent.FinSet;
  13 import net.sf.openrocket.rocketcomponent.RocketComponent;
  14 import net.sf.openrocket.util.BugException;
  15 import net.sf.openrocket.util.Coordinate;
  16 import net.sf.openrocket.util.LinearInterpolator;
  17 import net.sf.openrocket.util.MathUtil;
  18 import net.sf.openrocket.util.PolyInterpolator;
  19
  20
  21 public class FinSetCalc extends RocketComponentCalc {
  22         private static final double STALL_ANGLE = (20 * Math.PI / 180);
  23
  24         /** Number of divisions in the fin chords. */
  25         protected static final int DIVISIONS = 48;
  26
  27
  28
  29         protected double macLength = Double.NaN; // MAC length
  30         protected double macLead = Double.NaN; // MAC leading edge position
  31         protected double macSpan = Double.NaN; // MAC spanwise position
  32         protected double finArea = Double.NaN; // Fin area
  33         protected double ar = Double.NaN; // Fin aspect ratio
  34         protected double span = Double.NaN; // Fin span
  35         protected double cosGamma = Double.NaN; // Cosine of midchord sweep angle
  36         protected double cosGammaLead = Double.NaN; // Cosine of leading edge sweep angle
  37         protected double rollSum = Double.NaN; // Roll damping sum term
  38
  39         protected int interferenceFinCount = -1; // No. of fins in interference
  40
  41         protected double[] chordLead = new double[DIVISIONS];
  42         protected double[] chordTrail = new double[DIVISIONS];
  43         protected double[] chordLength = new double[DIVISIONS];
  44
  45
  46         protected final WarningSet geometryWarnings = new WarningSet();
  47
  48         private double[] poly = new double[6];
  49
  50         private final double thickness;
  51         private final double bodyRadius;
  52         private final int finCount;
  53         private final double baseRotation;
  54         private final double cantAngle;
  55         private final FinSet.CrossSection crossSection;
  56
  57         public FinSetCalc(RocketComponent component) {
  58                 super(component);
  59                 if (!(component instanceof FinSet)) {
  60                         throw new IllegalArgumentException("Illegal component type " + component);
  61                 }
  62
  63                 FinSet fin = (FinSet) component;
  64                 thickness = fin.getThickness();
  65                 bodyRadius = fin.getBodyRadius();
  66                 finCount = fin.getFinCount();
  67                 baseRotation = fin.getBaseRotation();
  68                 cantAngle = fin.getCantAngle();
  69                 span = fin.getSpan();
  70                 finArea = fin.getFinArea();
  71                 crossSection = fin.getCrossSection();
  72
  73                 calculateFinGeometry(fin);
  74                 calculatePoly();
  75                 calculateInterferenceFinCount(fin);
  76         }
  77
  78
  79         /*
  80          * Calculates the non-axial forces produced by the fins (normal and side forces,
  81          * pitch, yaw and roll moments, CP position, CNa).
  82          */
  83         @Override
  84         public void calculateNonaxialForces(FlightConditions conditions,
  85                         AerodynamicForces forces, WarningSet warnings) {
  86
  87
  88                 if (span < 0.001) {
  89                         forces.setCm(0);
  90                         forces.setCN(0);
  91                         forces.setCNa(0);
  92                         forces.setCP(Coordinate.NUL);
  93                         forces.setCroll(0);
  94                         forces.setCrollDamp(0);
  95                         forces.setCrollForce(0);
  96                         forces.setCside(0);
  97                         forces.setCyaw(0);
  98                         return;
  99                 }
 100
 101
 102                 // Add warnings  (radius/2 == diameter/4)
 103                 if (thickness > bodyRadius / 2) {
 104                         warnings.add(Warning.THICK_FIN);
 105                 }
 106                 warnings.addAll(geometryWarnings);
 107
 108
 109
 110                 //////// Calculate CNa.  /////////
 111
 112                 // One fin without interference (both sub- and supersonic):
 113                 double cna1 = calculateFinCNa1(conditions);
 114
 115
 116
 117                 // Multiple fins with fin-fin interference
 118                 double cna;
 119                 double theta = conditions.getTheta();
 120                 double angle = baseRotation;
 121
 122
 123                 // Compute basic CNa without interference effects
 124                 if (finCount == 1 || finCount == 2) {
 125                         // Basic CNa from geometry
 126                         double mul = 0;
 127                         for (int i = 0; i < finCount; i++) {
 128                                 mul += MathUtil.pow2(Math.sin(theta - angle));
 129                                 angle += 2 * Math.PI / finCount;
 130                         }
 131                         cna = cna1 * mul;
 132                 } else {
 133                         // Basic CNa assuming full efficiency
 134                         cna = cna1 * finCount / 2.0;
 135                 }
 136
 137
 138                 // Take into account fin-fin interference effects
 139                 switch (interferenceFinCount) {
 140                 case 1:
 141                 case 2:
 142                 case 3:
 143                 case 4:
 144                         // No interference effect
 145                         break;
 146
 147                 case 5:
 148                         cna *= 0.948;
 149                         break;
 150
 151                 case 6:
 152                         cna *= 0.913;
 153                         break;
 154
 155                 case 7:
 156                         cna *= 0.854;
 157                         break;
 158
 159                 case 8:
 160                         cna *= 0.81;
 161                         break;
 162
 163                 default:
 164                         // Assume 75% efficiency
 165                         cna *= 0.75;
 166                         warnings.add(Warning.PARALLEL_FINS);
 167                         break;
 168                 }
 169
 170                 /*
 171                  * Used in 0.9.5 and earlier.  Takes into account rotation angle for three
 172                  * and four fins, does not take into account interference from other fin sets.
 173                  *
 174                 switch (fins) {
 175                 case 1:
 176                 case 2:
 177                         // from geometry
 178                         double mul = 0;
 179                         for (int i=0; i < fins; i++) {
 180                                 mul += MathUtil.pow2(Math.sin(theta - angle));
 181                                 angle += 2 * Math.PI / fins;
 182                         }
 183                         cna = cna1*mul;
 184                         break;
 185
 186                 case 3:
 187                         // multiplier 1.5, sinusoidal reduction of 15%
 188                         cna = cna1 * 1.5 * (1 - 0.15*pow2(Math.cos(1.5 * (theta-angle))));
 189                         break;
 190
 191                 case 4:
 192                         // multiplier 2.0, sinusoidal reduction of 6%
 193                         cna = cna1 * 2.0 * (1 - 0.06*pow2(Math.sin(2 * (theta-angle))));
 194                         break;
 195
 196                 case 5:
 197                         cna = 2.37 * cna1;
 198                         break;
 199
 200                 case 6:
 201                         cna = 2.74 * cna1;
 202                         break;
 203
 204                 case 7:
 205                         cna = 2.99 * cna1;
 206                         break;
 207
 208                 case 8:
 209                         cna = 3.24 * cna1;
 210                         break;
 211
 212                 default:
 213                         // Assume N/2 * 3/4 efficiency for more fins
 214                         cna = cna1 * fins * 3.0/8.0;
 215                         break;
 216                 }
 217                 */
 218
 219                 // Body-fin interference effect
 220                 double r = bodyRadius;
 221                 double tau = r / (span + r);
 222                 if (Double.isNaN(tau) || Double.isInfinite(tau))
 223                         tau = 0;
 224                 cna *= 1 + tau; // Classical Barrowman
 225                 //              cna *= pow2(1 + tau);   // Barrowman thesis (too optimistic??)
 226
 227
 228
 229                 // TODO: LOW: check for fin tip mach cone interference
 230                 // (Barrowman thesis pdf-page 40)
 231
 232                 // TODO: LOW: fin-fin mach cone effect, MIL-HDBK page 5-25
 233
 234
 235
 236                 // Calculate CP position
 237                 double x = macLead + calculateCPPos(conditions) * macLength;
 238
 239
 240
 241                 // Calculate roll forces, reduce forcing above stall angle
 242
 243                 // Without body-fin interference effect:
 244                 //              forces.CrollForce = fins * (macSpan+r) * cna1 * component.getCantAngle() /
 245                 //                      conditions.getRefLength();
 246                 // With body-fin interference effect:
 247                 forces.setCrollForce(finCount * (macSpan + r) * cna1 * (1 + tau) * cantAngle / conditions.getRefLength());
 248
 249
 250
 251
 252                 if (conditions.getAOA() > STALL_ANGLE) {
 253                         //                      System.out.println("Fin stalling in roll");
 254                         forces.setCrollForce(forces.getCrollForce() * MathUtil.clamp(
 255                                         1 - (conditions.getAOA() - STALL_ANGLE) / (STALL_ANGLE / 2), 0, 1));
 256                 }
 257                 forces.setCrollDamp(calculateDampingMoment(conditions));
 258                 forces.setCroll(forces.getCrollForce() - forces.getCrollDamp());
 259
 260
 261
 262                 //              System.out.printf(component.getName() + ":  roll rate:%.3f  force:%.3f  damp:%.3f  " +
 263                 //                              "total:%.3f\n",
 264                 //                              conditions.getRollRate(), forces.CrollForce, forces.CrollDamp, forces.Croll);
 265
 266                 forces.setCNa(cna);
 267                 forces.setCN(cna * MathUtil.min(conditions.getAOA(), STALL_ANGLE));
 268                 forces.setCP(new Coordinate(x, 0, 0, cna));
 269                 forces.setCm(forces.getCN() * x / conditions.getRefLength());
 270
 271                 /*
 272                  * TODO: HIGH:  Compute actual side force and yaw moment.
 273                  * This is not currently performed because it produces strange results for
 274                  * stable rockets that have two fins in the front part of the fuselage,
 275                  * where the rocket flies at an ever-increasing angle of attack.  This may
 276                  * be due to incorrect computation of pitch/yaw damping moments.
 277                  */
 278                 //              if (fins == 1 || fins == 2) {
 279                 //                      forces.Cside = fins * cna1 * Math.cos(theta-angle) * Math.sin(theta-angle);
 280                 //                      forces.Cyaw = fins * forces.Cside * x / conditions.getRefLength();
 281                 //              } else {
 282                 //                      forces.Cside = 0;
 283                 //                      forces.Cyaw = 0;
 284                 //              }
 285                 forces.setCside(0);
 286                 forces.setCyaw(0);
 287
 288
 289         }
 290
 291
 292         /**
 293          * Returns the MAC length of the fin.  This is required in the friction drag
 294          * computation.
 295          *
 296          * @return  the MAC length of the fin.
 297          */
 298         public double getMACLength() {
 299                 return macLength;
 300         }
 301
 302         public double getMidchordPos() {
 303                 return macLead + 0.5 * macLength;
 304         }
 305
 306
 307
 308         /**
 309          * Pre-calculates the fin geometry values.
 310          */
 311         protected void calculateFinGeometry(FinSet component) {
 312
 313                 span = component.getSpan();
 314                 finArea = component.getFinArea();
 315                 ar = 2 * pow2(span) / finArea;
 316
 317                 Coordinate[] points = component.getFinPoints();
 318
 319                 // Check for jagged edges
 320                 geometryWarnings.clear();
 321                 boolean down = false;
 322                 for (int i = 1; i < points.length; i++) {
 323                         if ((points[i].y > points[i - 1].y + 0.001) && down) {
 324                                 geometryWarnings.add(Warning.JAGGED_EDGED_FIN);
 325                                 break;
 326                         }
 327                         if (points[i].y < points[i - 1].y - 0.001) {
 328                                 down = true;
 329                         }
 330                 }
 331
 332
 333                 // Calculate the chord lead and trail positions and length
 334
 335                 Arrays.fill(chordLead, Double.POSITIVE_INFINITY);
 336                 Arrays.fill(chordTrail, Double.NEGATIVE_INFINITY);
 337                 Arrays.fill(chordLength, 0);
 338
 339                 for (int point = 1; point < points.length; point++) {
 340                         double x1 = points[point - 1].x;
 341                         double y1 = points[point - 1].y;
 342                         double x2 = points[point].x;
 343                         double y2 = points[point].y;
 344
 345                         if (MathUtil.equals(y1, y2))
 346                                 continue;
 347
 348                         int i1 = (int) (y1 * 1.0001 / span * (DIVISIONS - 1));
 349                         int i2 = (int) (y2 * 1.0001 / span * (DIVISIONS - 1));
 350                         i1 = MathUtil.clamp(i1, 0, DIVISIONS - 1);
 351                         i2 = MathUtil.clamp(i2, 0, DIVISIONS - 1);
 352                         if (i1 > i2) {
 353                                 int tmp = i2;
 354                                 i2 = i1;
 355                                 i1 = tmp;
 356                         }
 357
 358                         for (int i = i1; i <= i2; i++) {
 359                                 // Intersection point (x,y)
 360                                 double y = i * span / (DIVISIONS - 1);
 361                                 double x = (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2;
 362                                 if (x < chordLead[i])
 363                                         chordLead[i] = x;
 364                                 if (x > chordTrail[i])
 365                                         chordTrail[i] = x;
 366
 367                                 // TODO: LOW:  If fin point exactly on chord line, might be counted twice:
 368                                 if (y1 < y2) {
 369                                         chordLength[i] -= x;
 370                                 } else {
 371                                         chordLength[i] += x;
 372                                 }
 373                         }
 374                 }
 375
 376                 // Check and correct any inconsistencies
 377                 for (int i = 0; i < DIVISIONS; i++) {
 378                         if (Double.isInfinite(chordLead[i]) || Double.isInfinite(chordTrail[i]) ||
 379                                         Double.isNaN(chordLead[i]) || Double.isNaN(chordTrail[i])) {
 380                                 chordLead[i] = 0;
 381                                 chordTrail[i] = 0;
 382                         }
 383                         if (chordLength[i] < 0 || Double.isNaN(chordLength[i])) {
 384                                 chordLength[i] = 0;
 385                         }
 386                         if (chordLength[i] > chordTrail[i] - chordLead[i]) {
 387                                 chordLength[i] = chordTrail[i] - chordLead[i];
 388                         }
 389                 }
 390
 391
 392                 /* Calculate fin properties:
 393                  *
 394                  * macLength // MAC length
 395                  * macLead   // MAC leading edge position
 396                  * macSpan   // MAC spanwise position
 397                  * ar        // Fin aspect ratio (already set)
 398                  * span      // Fin span (already set)
 399                  */
 400                 macLength = 0;
 401                 macLead = 0;
 402                 macSpan = 0;
 403                 cosGamma = 0;
 404                 cosGammaLead = 0;
 405                 rollSum = 0;
 406                 double area = 0;
 407                 double radius = component.getBodyRadius();
 408
 409                 final double dy = span / (DIVISIONS - 1);
 410                 for (int i = 0; i < DIVISIONS; i++) {
 411                         double length = chordTrail[i] - chordLead[i];
 412                         double y = i * dy;
 413
 414                         macLength += length * length;
 415                         macSpan += y * length;
 416                         macLead += chordLead[i] * length;
 417                         area += length;
 418                         rollSum += chordLength[i] * pow2(radius + y);
 419
 420                         if (i > 0) {
 421                                 double dx = (chordTrail[i] + chordLead[i]) / 2 - (chordTrail[i - 1] + chordLead[i - 1]) / 2;
 422                                 cosGamma += dy / MathUtil.hypot(dx, dy);
 423
 424                                 dx = chordLead[i] - chordLead[i - 1];
 425                                 cosGammaLead += dy / MathUtil.hypot(dx, dy);
 426                         }
 427                 }
 428
 429                 macLength *= dy;
 430                 macSpan *= dy;
 431                 macLead *= dy;
 432                 area *= dy;
 433                 rollSum *= dy;
 434
 435                 macLength /= area;
 436                 macSpan /= area;
 437                 macLead /= area;
 438                 cosGamma /= (DIVISIONS - 1);
 439                 cosGammaLead /= (DIVISIONS - 1);
 440         }
 441
 442
 443         ///////////////  CNa1 calculation  ////////////////
 444
 445         private static final double CNA_SUBSONIC = 0.9;
 446         private static final double CNA_SUPERSONIC = 1.5;
 447         private static final double CNA_SUPERSONIC_B = pow(pow2(CNA_SUPERSONIC) - 1, 1.5);
 448         private static final double GAMMA = 1.4;
 449         private static final LinearInterpolator K1, K2, K3;
 450         private static final PolyInterpolator cnaInterpolator = new PolyInterpolator(
 451                         new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
 452                         new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
 453                         new double[] { CNA_SUBSONIC }
 454                         );
 455         /* Pre-calculate the values for K1, K2 and K3 */
 456         static {
 457                 // Up to Mach 5
 458                 int n = (int) ((5.0 - CNA_SUPERSONIC) * 10);
 459                 double[] x = new double[n];
 460                 double[] k1 = new double[n];
 461                 double[] k2 = new double[n];
 462                 double[] k3 = new double[n];
 463                 for (int i = 0; i < n; i++) {
 464                         double M = CNA_SUPERSONIC + i * 0.1;
 465                         double beta = MathUtil.safeSqrt(M * M - 1);
 466                         x[i] = M;
 467                         k1[i] = 2.0 / beta;
 468                         k2[i] = ((GAMMA + 1) * pow(M, 4) - 4 * pow2(beta)) / (4 * pow(beta, 4));
 469                         k3[i] = ((GAMMA + 1) * pow(M, 8) + (2 * pow2(GAMMA) - 7 * GAMMA - 5) * pow(M, 6) +
 470                                         10 * (GAMMA + 1) * pow(M, 4) + 8) / (6 * pow(beta, 7));
 471                 }
 472                 K1 = new LinearInterpolator(x, k1);
 473                 K2 = new LinearInterpolator(x, k2);
 474                 K3 = new LinearInterpolator(x, k3);
 475
 476                 //              System.out.println("K1[m="+CNA_SUPERSONIC+"] = "+k1[0]);
 477                 //              System.out.println("K2[m="+CNA_SUPERSONIC+"] = "+k2[0]);
 478                 //              System.out.println("K3[m="+CNA_SUPERSONIC+"] = "+k3[0]);
 479         }
 480
 481
 482         protected double calculateFinCNa1(FlightConditions conditions) {
 483                 double mach = conditions.getMach();
 484                 double ref = conditions.getRefArea();
 485                 double alpha = MathUtil.min(conditions.getAOA(),
 486                                 Math.PI - conditions.getAOA(), STALL_ANGLE);
 487
 488                 // Subsonic case
 489                 if (mach <= CNA_SUBSONIC) {
 490                         return 2 * Math.PI * pow2(span) / (1 + MathUtil.safeSqrt(1 + (1 - pow2(mach)) *
 491                                         pow2(pow2(span) / (finArea * cosGamma)))) / ref;
 492                 }
 493
 494                 // Supersonic case
 495                 if (mach >= CNA_SUPERSONIC) {
 496                         return finArea * (K1.getValue(mach) + K2.getValue(mach) * alpha +
 497                                         K3.getValue(mach) * pow2(alpha)) / ref;
 498                 }
 499
 500                 // Transonic case, interpolate
 501                 double subV, superV;
 502                 double subD, superD;
 503
 504                 double sq = MathUtil.safeSqrt(1 + (1 - pow2(CNA_SUBSONIC)) * pow2(span * span / (finArea * cosGamma)));
 505                 subV = 2 * Math.PI * pow2(span) / ref / (1 + sq);
 506                 subD = 2 * mach * Math.PI * pow(span, 6) / (pow2(finArea * cosGamma) * ref *
 507                                 sq * pow2(1 + sq));
 508
 509                 superV = finArea * (K1.getValue(CNA_SUPERSONIC) + K2.getValue(CNA_SUPERSONIC) * alpha +
 510                                 K3.getValue(CNA_SUPERSONIC) * pow2(alpha)) / ref;
 511                 superD = -finArea / ref * 2 * CNA_SUPERSONIC / CNA_SUPERSONIC_B;
 512
 513                 //              System.out.println("subV="+subV+" superV="+superV+" subD="+subD+" superD="+superD);
 514
 515                 return cnaInterpolator.interpolate(mach, subV, superV, subD, superD, 0);
 516         }
 517
 518
 519
 520
 521         private double calculateDampingMoment(FlightConditions conditions) {
 522                 double rollRate = conditions.getRollRate();
 523
 524                 if (Math.abs(rollRate) < 0.1)
 525                         return 0;
 526
 527                 double mach = conditions.getMach();
 528                 double absRate = Math.abs(rollRate);
 529
 530
 531                 /*
 532                  * At low speeds and relatively large roll rates (i.e. near apogee) the
 533                  * fin tips rotate well above stall angle.  In this case sum the chords
 534                  * separately.
 535                  */
 536                 if (absRate * (bodyRadius + span) / conditions.getVelocity() > 15 * Math.PI / 180) {
 537                         double sum = 0;
 538                         for (int i = 0; i < DIVISIONS; i++) {
 539                                 double dist = bodyRadius + span * i / DIVISIONS;
 540                                 double aoa = Math.min(absRate * dist / conditions.getVelocity(), 15 * Math.PI / 180);
 541                                 sum += chordLength[i] * dist * aoa;
 542                         }
 543                         sum = sum * (span / DIVISIONS) * 2 * Math.PI / conditions.getBeta() /
 544                                         (conditions.getRefArea() * conditions.getRefLength());
 545
 546                         //                      System.out.println("SPECIAL: " +
 547                         //                                      (MathUtil.sign(rollRate) *component.getFinCount() * sum));
 548                         return MathUtil.sign(rollRate) * finCount * sum;
 549                 }
 550
 551
 552
 553                 if (mach <= CNA_SUBSONIC) {
 554                         //                      System.out.println("BASIC:   "+
 555                         //                                      (component.getFinCount() * 2*Math.PI * rollRate * rollSum /
 556                         //                      (conditions.getRefArea() * conditions.getRefLength() *
 557                         //                                      conditions.getVelocity() * conditions.getBeta())));
 558
 559                         return finCount * 2 * Math.PI * rollRate * rollSum /
 560                                         (conditions.getRefArea() * conditions.getRefLength() *
 561                                                         conditions.getVelocity() * conditions.getBeta());
 562                 }
 563                 if (mach >= CNA_SUPERSONIC) {
 564
 565                         double vel = conditions.getVelocity();
 566                         double k1 = K1.getValue(mach);
 567                         double k2 = K2.getValue(mach);
 568                         double k3 = K3.getValue(mach);
 569
 570                         double sum = 0;
 571
 572                         for (int i = 0; i < DIVISIONS; i++) {
 573                                 double y = i * span / (DIVISIONS - 1);
 574                                 double angle = rollRate * (bodyRadius + y) / vel;
 575
 576                                 sum += (k1 * angle + k2 * angle * angle + k3 * angle * angle * angle)
 577                                                 * chordLength[i] * (bodyRadius + y);
 578                         }
 579
 580                         return finCount * sum * span / (DIVISIONS - 1) /
 581                                         (conditions.getRefArea() * conditions.getRefLength());
 582                 }
 583
 584                 // Transonic, do linear interpolation
 585
 586                 FlightConditions cond = conditions.clone();
 587                 cond.setMach(CNA_SUBSONIC - 0.01);
 588                 double subsonic = calculateDampingMoment(cond);
 589                 cond.setMach(CNA_SUPERSONIC + 0.01);
 590                 double supersonic = calculateDampingMoment(cond);
 591
 592                 return subsonic * (CNA_SUPERSONIC - mach) / (CNA_SUPERSONIC - CNA_SUBSONIC) +
 593                                 supersonic * (mach - CNA_SUBSONIC) / (CNA_SUPERSONIC - CNA_SUBSONIC);
 594         }
 595
 596
 597
 598
 599         /**
 600          * Return the relative position of the CP along the mean aerodynamic chord.
 601          * Below mach 0.5 it is at the quarter chord, above mach 2 calculated using an
 602          * empirical formula, between these two using an interpolation polynomial.
 603          *
 604          * @param cond   Mach speed used
 605          * @return               CP position along the MAC
 606          */
 607         private double calculateCPPos(FlightConditions cond) {
 608                 double m = cond.getMach();
 609                 if (m <= 0.5) {
 610                         // At subsonic speeds CP at quarter chord
 611                         return 0.25;
 612                 }
 613                 if (m >= 2) {
 614                         // At supersonic speeds use empirical formula
 615                         double beta = cond.getBeta();
 616                         return (ar * beta - 0.67) / (2 * ar * beta - 1);
 617                 }
 618
 619                 // In between use interpolation polynomial
 620                 double x = 1.0;
 621                 double val = 0;
 622
 623                 for (int i = 0; i < poly.length; i++) {
 624                         val += poly[i] * x;
 625                         x *= m;
 626                 }
 627
 628                 return val;
 629         }
 630
 631         /**
 632          * Calculate CP position interpolation polynomial coefficients from the
 633          * fin geometry.  This is a fifth order polynomial that satisfies
 634          *
 635          * p(0.5)=0.25
 636          * p'(0.5)=0
 637          * p(2) = f(2)
 638          * p'(2) = f'(2)
 639          * p''(2) = 0
 640          * p'''(2) = 0
 641          *
 642          * where f(M) = (ar*sqrt(M^2-1) - 0.67) / (2*ar*sqrt(M^2-1) - 1).
 643          *
 644          * The values were calculated analytically in Mathematica.  The coefficients
 645          * are used as poly[0] + poly[1]*x + poly[2]*x^2 + ...
 646          */
 647         private void calculatePoly() {
 648                 double denom = pow2(1 - 3.4641 * ar); // common denominator
 649
 650                 poly[5] = (-1.58025 * (-0.728769 + ar) * (-0.192105 + ar)) / denom;
 651                 poly[4] = (12.8395 * (-0.725688 + ar) * (-0.19292 + ar)) / denom;
 652                 poly[3] = (-39.5062 * (-0.72074 + ar) * (-0.194245 + ar)) / denom;
 653                 poly[2] = (55.3086 * (-0.711482 + ar) * (-0.196772 + ar)) / denom;
 654                 poly[1] = (-31.6049 * (-0.705375 + ar) * (-0.198476 + ar)) / denom;
 655                 poly[0] = (9.16049 * (-0.588838 + ar) * (-0.20624 + ar)) / denom;
 656         }
 657
 658
 659         //      @SuppressWarnings("null")
 660         //      public static void main(String arg[]) {
 661         //              Rocket rocket = TestRocket.makeRocket();
 662         //              FinSet finset = null;
 663         //
 664         //              Iterator<RocketComponent> iter = rocket.deepIterator();
 665         //              while (iter.hasNext()) {
 666         //                      RocketComponent c = iter.next();
 667         //                      if (c instanceof FinSet) {
 668         //                              finset = (FinSet)c;
 669         //                              break;
 670         //                      }
 671         //              }
 672         //
 673         //              ((TrapezoidFinSet)finset).setHeight(0.10);
 674         //              ((TrapezoidFinSet)finset).setRootChord(0.10);
 675         //              ((TrapezoidFinSet)finset).setTipChord(0.10);
 676         //              ((TrapezoidFinSet)finset).setSweep(0.0);
 677         //
 678         //
 679         //              FinSetCalc calc = new FinSetCalc(finset);
 680         //
 681         //              calc.calculateFinGeometry();
 682         //              FlightConditions cond = new FlightConditions(new Configuration(rocket));
 683         //              for (double m=0; m < 3; m+=0.05) {
 684         //                      cond.setMach(m);
 685         //                      cond.setAOA(0.0*Math.PI/180);
 686         //                      double cna = calc.calculateFinCNa1(cond);
 687         //                      System.out.printf("%5.2f "+cna+"\n", m);
 688         //              }
 689         //
 690         //      }
 691
 692
 693         @Override
 694         public double calculatePressureDragForce(FlightConditions conditions,
 695                         double stagnationCD, double baseCD, WarningSet warnings) {
 696
 697                 double mach = conditions.getMach();
 698                 double drag = 0;
 699
 700                 // Pressure fore-drag
 701                 if (crossSection == FinSet.CrossSection.AIRFOIL ||
 702                                 crossSection == FinSet.CrossSection.ROUNDED) {
 703
 704                         // Round leading edge
 705                         if (mach < 0.9) {
 706                                 drag = Math.pow(1 - pow2(mach), -0.417) - 1;
 707                         } else if (mach < 1) {
 708                                 drag = 1 - 1.785 * (mach - 0.9);
 709                         } else {
 710                                 drag = 1.214 - 0.502 / pow2(mach) + 0.1095 / pow2(pow2(mach));
 711                         }
 712
 713                 } else if (crossSection == FinSet.CrossSection.SQUARE) {
 714                         drag = stagnationCD;
 715                 } else {
 716                         throw new UnsupportedOperationException("Unsupported fin profile: " + crossSection);
 717                 }
 718
 719                 // Slanted leading edge
 720                 drag *= pow2(cosGammaLead);
 721
 722                 // Trailing edge drag
 723                 if (crossSection == FinSet.CrossSection.SQUARE) {
 724                         drag += baseCD;
 725                 } else if (crossSection == FinSet.CrossSection.ROUNDED) {
 726                         drag += baseCD / 2;
 727                 }
 728                 // Airfoil assumed to have zero base drag
 729
 730
 731                 // Scale to correct reference area
 732                 drag *= finCount * span * thickness / conditions.getRefArea();
 733
 734                 return drag;
 735         }
 736
 737
 738         private void calculateInterferenceFinCount(FinSet component) {
 739                 RocketComponent parent = component.getParent();
 740                 if (parent == null) {
 741                         throw new IllegalStateException("fin set without parent component");
 742                 }
 743
 744                 double lead = component.toRelative(Coordinate.NUL, parent)[0].x;
 745                 double trail = component.toRelative(new Coordinate(component.getLength()),
 746                                         parent)[0].x;
 747
 748                 /*
 749                  * The counting fails if the fin root chord is very small, in that case assume
 750                  * no other fin interference than this fin set.
 751                  */
 752                 if (trail - lead < 0.007) {
 753                         interferenceFinCount = finCount;
 754                 } else {
 755                         interferenceFinCount = 0;
 756                         for (RocketComponent c : parent.getChildren()) {
 757                                 if (c instanceof FinSet) {
 758                                         double finLead = c.toRelative(Coordinate.NUL, parent)[0].x;
 759                                         double finTrail = c.toRelative(new Coordinate(c.getLength()), parent)[0].x;
 760
 761                                         // Compute overlap of the fins
 762
 763                                         if ((finLead < trail - 0.005) && (finTrail > lead + 0.005)) {
 764                                                 interferenceFinCount += ((FinSet) c).getFinCount();
 765                                         }
 766                                 }
 767                         }
 768                 }
 769                 if (interferenceFinCount < component.getFinCount()) {
 770                         throw new BugException("Counted " + interferenceFinCount + " parallel fins, " +
 771                                                 "when component itself has " + component.getFinCount() +
 772                                                 ", fin points=" + Arrays.toString(component.getFinPoints()));
 773                 }
 774         }
 775
 776 }