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

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