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