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