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

   1 package net.sf.openrocket.aerodynamics;
   2
   3 import static net.sf.openrocket.util.MathUtil.pow2;
   4
   5 import java.util.Arrays;
   6 import java.util.HashMap;
   7 import java.util.Iterator;
   8 import java.util.LinkedHashMap;
   9 import java.util.Map;
  10
  11 import net.sf.openrocket.aerodynamics.barrowman.FinSetCalc;
  12 import net.sf.openrocket.aerodynamics.barrowman.RocketComponentCalc;
  13 import net.sf.openrocket.rocketcomponent.Configuration;
  14 import net.sf.openrocket.rocketcomponent.ExternalComponent;
  15 import net.sf.openrocket.rocketcomponent.FinSet;
  16 import net.sf.openrocket.rocketcomponent.RocketComponent;
  17 import net.sf.openrocket.rocketcomponent.SymmetricComponent;
  18 import net.sf.openrocket.rocketcomponent.ExternalComponent.Finish;
  19 import net.sf.openrocket.util.Coordinate;
  20 import net.sf.openrocket.util.MathUtil;
  21 import net.sf.openrocket.util.PolyInterpolator;
  22 import net.sf.openrocket.util.Reflection;
  23
  24 /**
  25  * An aerodynamic calculator that uses the extended Barrowman method to
  26  * calculate the CP of a rocket.
  27  *
  28  * @author Sampo Niskanen <sampo.niskanen@iki.fi>
  29  */
  30 public class BarrowmanCalculator extends AerodynamicCalculator {
  31
  32         private static final String BARROWMAN_PACKAGE = "net.sf.openrocket.aerodynamics.barrowman";
  33         private static final String BARROWMAN_SUFFIX = "Calc";
  34
  35
  36         private Map<RocketComponent, RocketComponentCalc> calcMap = null;
  37
  38         private double cacheDiameter = -1;
  39         private double cacheLength = -1;
  40
  41
  42
  43         public BarrowmanCalculator() {
  44
  45         }
  46
  47         public BarrowmanCalculator(Configuration config) {
  48                 super(config);
  49         }
  50
  51
  52         @Override
  53         public BarrowmanCalculator newInstance() {
  54                 return new BarrowmanCalculator();
  55         }
  56
  57
  58         /**
  59          * Calculate the CP according to the extended Barrowman method.
  60          */
  61         @Override
  62         public Coordinate getCP(FlightConditions conditions, WarningSet warnings) {
  63                 AerodynamicForces forces = getNonAxial(conditions, null, warnings);
  64                 return forces.cp;
  65         }
  66
  67
  68
  69         @Override
  70         public Map<RocketComponent, AerodynamicForces> getForceAnalysis(FlightConditions conditions,
  71                         WarningSet warnings) {
  72
  73                 AerodynamicForces f;
  74                 Map<RocketComponent, AerodynamicForces> map =
  75                         new LinkedHashMap<RocketComponent, AerodynamicForces>();
  76
  77                 // Add all components to the map
  78                 for (RocketComponent c: configuration) {
  79                         f = new AerodynamicForces();
  80                         f.component = c;
  81
  82                         // Calculate CG
  83                         f.cg = Coordinate.NUL;
  84                         for (Coordinate coord: c.toAbsolute(c.getCG())) {
  85                                 f.cg = f.cg.average(coord);
  86                         }
  87
  88                         map.put(c, f);
  89                 }
  90
  91
  92                 // Calculate non-axial force data
  93                 AerodynamicForces total = getNonAxial(conditions, map, warnings);
  94
  95
  96                 // Calculate friction data
  97                 total.frictionCD = calculateFrictionDrag(conditions, map, warnings);
  98                 total.pressureCD = calculatePressureDrag(conditions, map, warnings);
  99                 total.baseCD = calculateBaseDrag(conditions, map, warnings);
 100                 total.cg = getCG(0);
 101
 102                 total.component = rocket;
 103                 map.put(rocket, total);
 104
 105
 106                 for (RocketComponent c: map.keySet()) {
 107                         f = map.get(c);
 108                         if (Double.isNaN(f.baseCD) && Double.isNaN(f.pressureCD) &&
 109                                         Double.isNaN(f.frictionCD))
 110                                 continue;
 111                         if (Double.isNaN(f.baseCD))
 112                                 f.baseCD = 0;
 113                         if (Double.isNaN(f.pressureCD))
 114                                 f.pressureCD = 0;
 115                         if (Double.isNaN(f.frictionCD))
 116                                 f.frictionCD = 0;
 117                         f.CD = f.baseCD + f.pressureCD + f.frictionCD;
 118                         f.Caxial = calculateAxialDrag(conditions, f.CD);
 119                 }
 120
 121                 return map;
 122         }
 123
 124
 125
 126         @Override
 127         public AerodynamicForces getAerodynamicForces(double time, FlightConditions conditions,
 128                         WarningSet warnings) {
 129
 130                 if (warnings == null)
 131                         warnings = ignoreWarningSet;
 132
 133                 // Calculate non-axial force data
 134                 AerodynamicForces total = getNonAxial(conditions, null, warnings);
 135
 136                 // Calculate friction data
 137                 total.frictionCD = calculateFrictionDrag(conditions, null, warnings);
 138                 total.pressureCD = calculatePressureDrag(conditions, null, warnings);
 139                 total.baseCD = calculateBaseDrag(conditions, null, warnings);
 140
 141                 total.CD = total.frictionCD + total.pressureCD + total.baseCD;
 142
 143                 total.Caxial = calculateAxialDrag(conditions, total.CD);
 144
 145
 146                 // Calculate CG and moments of inertia
 147                 total.cg = this.getCG(time);
 148                 total.longitudalInertia = this.getLongitudalInertia(time);
 149                 total.rotationalInertia = this.getRotationalInertia(time);
 150
 151
 152                 // Calculate pitch and yaw damping moments
 153                 calculateDampingMoments(conditions, total);
 154                 total.Cm -= total.pitchDampingMoment;
 155                 total.Cyaw -= total.yawDampingMoment;
 156
 157
 158 //              System.out.println("Conditions are "+conditions + "
 159 //              pitch rate="+conditions.getPitchRate());
 160 //              System.out.println("Total Cm="+total.Cm+" damping effect="+
 161 //                              (12 * Math.signum(conditions.getPitchRate()) *
 162 //                              MathUtil.pow2(conditions.getPitchRate()) /
 163 //                              MathUtil.pow2(conditions.getVelocity())));
 164
 165 //              double ef = Math.abs(12 *
 166 //                              MathUtil.pow2(conditions.getPitchRate()) /
 167 //                              MathUtil.pow2(conditions.getVelocity()));
 168 //
 169 ////            System.out.println("maxEffect="+maxEffect);
 170 //              total.Cm -= 12 * Math.signum(conditions.getPitchRate()) *
 171 //                              MathUtil.pow2(conditions.getPitchRate()) /
 172 //                              MathUtil.pow2(conditions.getVelocity());
 173 //
 174 //              total.Cyaw -= 0.06 * Math.signum(conditions.getYawRate()) *
 175 //                                              MathUtil.pow2(conditions.getYawRate()) /
 176 //                                              MathUtil.pow2(conditions.getVelocity());
 177
 178                 return total;
 179         }
 180
 181
 182
 183         @Override
 184         public AerodynamicForces getAxialForces(double time,
 185                         FlightConditions conditions, WarningSet warnings) {
 186
 187                 if (warnings == null)
 188                         warnings = ignoreWarningSet;
 189
 190                 AerodynamicForces total = new AerodynamicForces();
 191                 total.zero();
 192
 193                 // Calculate friction data
 194                 total.frictionCD = calculateFrictionDrag(conditions, null, warnings);
 195                 total.pressureCD = calculatePressureDrag(conditions, null, warnings);
 196                 total.baseCD = calculateBaseDrag(conditions, null, warnings);
 197
 198                 total.CD = total.frictionCD + total.pressureCD + total.baseCD;
 199
 200                 total.Caxial = calculateAxialDrag(conditions, total.CD);
 201
 202                 // Calculate CG and moments of inertia
 203                 total.cg = this.getCG(time);
 204                 total.longitudalInertia = this.getLongitudalInertia(time);
 205                 total.rotationalInertia = this.getRotationalInertia(time);
 206
 207                 return total;
 208         }
 209
 210
 211
 212
 213
 214         /*
 215          * Perform the actual CP calculation.
 216          */
 217         private AerodynamicForces getNonAxial(FlightConditions conditions,
 218                         Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
 219
 220                 AerodynamicForces total = new AerodynamicForces();
 221                 total.zero();
 222
 223                 double radius = 0;      // aft radius of previous component
 224                 double componentX = 0;  // aft coordinate of previous component
 225                 AerodynamicForces forces = new AerodynamicForces();
 226
 227                 if (warnings == null)
 228                         warnings = ignoreWarningSet;
 229
 230                 if (conditions.getAOA() > 17.5*Math.PI/180)
 231                         warnings.add(new Warning.LargeAOA(conditions.getAOA()));
 232
 233                 checkCache();
 234
 235                 if (calcMap == null)
 236                         buildCalcMap();
 237
 238                 for (RocketComponent component: configuration) {
 239
 240                         // Skip non-aerodynamic components
 241                         if (!component.isAerodynamic())
 242                                 continue;
 243
 244                         // Check for discontinuities
 245                         if (component instanceof SymmetricComponent) {
 246                                 SymmetricComponent sym = (SymmetricComponent) component;
 247                                 // TODO:LOW: Ignores other cluster components (not clusterable)
 248                                 double x = component.toAbsolute(Coordinate.NUL)[0].x;
 249
 250                                 // Check for lengthwise discontinuity
 251                                 if (x > componentX + 0.0001){
 252                                         if (!MathUtil.equals(radius, 0)) {
 253                                                 warnings.add(Warning.DISCONTINUITY);
 254                                                 radius = 0;
 255                                         }
 256                                 }
 257                                 componentX = component.toAbsolute(new Coordinate(component.getLength()))[0].x;
 258
 259                                 // Check for radius discontinuity
 260                                 if (!MathUtil.equals(sym.getForeRadius(), radius)) {
 261                                         warnings.add(Warning.DISCONTINUITY);
 262                                         // TODO: MEDIUM: Apply correction to values to cp and to map
 263                                 }
 264                                 radius = sym.getAftRadius();
 265                         }
 266
 267                         // Call calculation method
 268                         forces.zero();
 269                         calcMap.get(component).calculateNonaxialForces(conditions, forces, warnings);
 270                         forces.cp = component.toAbsolute(forces.cp)[0];
 271                         forces.Cm = forces.CN * forces.cp.x / conditions.getRefLength();
 272 //                      System.out.println("  CN="+forces.CN+" cp.x="+forces.cp.x+" Cm="+forces.Cm);
 273
 274                         if (map != null) {
 275                                 AerodynamicForces f = map.get(component);
 276
 277                                 f.cp = forces.cp;
 278                                 f.CNa = forces.CNa;
 279                                 f.CN = forces.CN;
 280                                 f.Cm = forces.Cm;
 281                                 f.Cside = forces.Cside;
 282                                 f.Cyaw = forces.Cyaw;
 283                                 f.Croll = forces.Croll;
 284                                 f.CrollDamp = forces.CrollDamp;
 285                                 f.CrollForce = forces.CrollForce;
 286                         }
 287
 288                         total.cp = total.cp.average(forces.cp);
 289                         total.CNa += forces.CNa;
 290                         total.CN += forces.CN;
 291                         total.Cm += forces.Cm;
 292                         total.Cside += forces.Cside;
 293                         total.Cyaw += forces.Cyaw;
 294                         total.Croll += forces.Croll;
 295                         total.CrollDamp += forces.CrollDamp;
 296                         total.CrollForce += forces.CrollForce;
 297                 }
 298
 299                 return total;
 300         }
 301
 302
 303
 304
 305         ////////////////  DRAG CALCULATIONS  ////////////////
 306
 307
 308         private double calculateFrictionDrag(FlightConditions conditions,
 309                         Map<RocketComponent, AerodynamicForces> map, WarningSet set) {
 310                 double c1=1.0, c2=1.0;
 311
 312                 double mach = conditions.getMach();
 313                 double Re;
 314                 double Cf;
 315
 316                 if (calcMap == null)
 317                         buildCalcMap();
 318
 319                 Re = conditions.getVelocity() * configuration.getLength() /
 320                 conditions.getAtmosphericConditions().getKinematicViscosity();
 321
 322 //              System.out.printf("Re=%.3e   ", Re);
 323
 324                 // Calculate the skin friction coefficient (assume non-roughness limited)
 325                 if (configuration.getRocket().isPerfectFinish()) {
 326
 327 //                      System.out.printf("Perfect finish: Re=%f ",Re);
 328                         // Assume partial laminar layer.  Roughness-limitation is checked later.
 329                         if (Re < 1e4) {
 330                                 // Too low, constant
 331                                 Cf = 1.33e-2;
 332 //                              System.out.printf("constant Cf=%f ",Cf);
 333                         } else if (Re < 5.39e5) {
 334                                 // Fully laminar
 335                                 Cf = 1.328 / Math.sqrt(Re);
 336 //                              System.out.printf("basic Cf=%f ",Cf);
 337                         } else {
 338                                 // Transitional
 339                                 Cf = 1.0/pow2(1.50 * Math.log(Re) - 5.6) - 1700/Re;
 340 //                              System.out.printf("transitional Cf=%f ",Cf);
 341                         }
 342
 343                         // Compressibility correction
 344
 345                         if (mach < 1.1) {
 346                                 // Below Re=1e6 no correction
 347                                 if (Re > 1e6) {
 348                                         if (Re < 3e6) {
 349                                                 c1 = 1 - 0.1*pow2(mach)*(Re-1e6)/2e6;  // transition to turbulent
 350                                         } else {
 351                                                 c1 = 1 - 0.1*pow2(mach);
 352                                         }
 353                                 }
 354                         }
 355                         if (mach > 0.9) {
 356                                 if (Re > 1e6) {
 357                                         if (Re < 3e6) {
 358                                                 c2 = 1 + (1.0 / Math.pow(1+0.045*pow2(mach), 0.25) -1) * (Re-1e6)/2e6;
 359                                         } else {
 360                                                 c2 = 1.0 / Math.pow(1+0.045*pow2(mach), 0.25);
 361                                         }
 362                                 }
 363                         }
 364
 365 //                      System.out.printf("c1=%f c2=%f\n", c1,c2);
 366                         // Applying continuously around Mach 1
 367                         if (mach < 0.9) {
 368                                 Cf *= c1;
 369                         } else if (mach < 1.1) {
 370                                 Cf *= (c2 * (mach-0.9)/0.2 + c1 * (1.1-mach)/0.2);
 371                         } else {
 372                                 Cf *= c2;
 373                         }
 374
 375 //                      System.out.printf("M=%f Cf=%f (smooth)\n",mach,Cf);
 376
 377                 } else {
 378
 379                         // Assume fully turbulent.  Roughness-limitation is checked later.
 380                         if (Re < 1e4) {
 381                                 // Too low, constant
 382                                 Cf = 1.48e-2;
 383 //                              System.out.printf("LOW-TURB  ");
 384                         } else {
 385                                 // Turbulent
 386                                 Cf = 1.0/pow2(1.50 * Math.log(Re) - 5.6);
 387 //                              System.out.printf("NORMAL-TURB  ");
 388                         }
 389
 390                         // Compressibility correction
 391
 392                         if (mach < 1.1) {
 393                                 c1 = 1 - 0.1*pow2(mach);
 394                         }
 395                         if (mach > 0.9) {
 396                                 c2 = 1/Math.pow(1 + 0.15*pow2(mach), 0.58);
 397                         }
 398                         // Applying continuously around Mach 1
 399                         if (mach < 0.9) {
 400                                 Cf *= c1;
 401                         } else if (mach < 1.1) {
 402                                 Cf *= c2 * (mach-0.9)/0.2 + c1 * (1.1-mach)/0.2;
 403                         } else {
 404                                 Cf *= c2;
 405                         }
 406
 407 //                      System.out.printf("M=%f, Cd=%f (turbulent)\n", mach,Cf);
 408
 409                 }
 410
 411                 // Roughness-limited value correction term
 412                 double roughnessCorrection;
 413                 if (mach < 0.9) {
 414                         roughnessCorrection = 1 - 0.1*pow2(mach);
 415                 } else if (mach > 1.1) {
 416                         roughnessCorrection = 1/(1 + 0.18*pow2(mach));
 417                 } else {
 418                         c1 = 1 - 0.1*pow2(0.9);
 419                         c2 = 1.0/(1+0.18 * pow2(1.1));
 420                         roughnessCorrection = c2 * (mach-0.9)/0.2 + c1 * (1.1-mach)/0.2;
 421                 }
 422
 423 //              System.out.printf("Cf=%.3f  ", Cf);
 424
 425
 426                 /*
 427                  * Calculate the friction drag coefficient.
 428                  *
 429                  * The body wetted area is summed up and finally corrected with the rocket
 430                  * fineness ratio (calculated in the same iteration).  The fins are corrected
 431                  * for thickness as we go on.
 432                  */
 433
 434                 double finFriction = 0;
 435                 double bodyFriction = 0;
 436                 double maxR=0, len=0;
 437
 438                 double[] roughnessLimited = new double[Finish.values().length];
 439                 Arrays.fill(roughnessLimited, Double.NaN);
 440
 441                 for (RocketComponent c: configuration) {
 442
 443                         // Consider only SymmetricComponents and FinSets:
 444                         if (!(c instanceof SymmetricComponent) &&
 445                                         !(c instanceof FinSet))
 446                                 continue;
 447
 448                         // Calculate the roughness-limited friction coefficient
 449                         Finish finish = ((ExternalComponent)c).getFinish();
 450                         if (Double.isNaN(roughnessLimited[finish.ordinal()])) {
 451                                 roughnessLimited[finish.ordinal()] =
 452                                         0.032 * Math.pow(finish.getRoughnessSize()/configuration.getLength(), 0.2) *
 453                                         roughnessCorrection;
 454
 455 //                              System.out.printf("roughness["+finish+"]=%.3f  ",
 456 //                                              roughnessLimited[finish.ordinal()]);
 457                         }
 458
 459                         /*
 460                          * Actual Cf is maximum of Cf and the roughness-limited value.
 461                          * For perfect finish require additionally that Re > 1e6
 462                          */
 463                         double componentCf;
 464                         if (configuration.getRocket().isPerfectFinish()) {
 465
 466                                 // For perfect finish require Re > 1e6
 467                                 if ((Re > 1.0e6) && (roughnessLimited[finish.ordinal()] > Cf)) {
 468                                         componentCf = roughnessLimited[finish.ordinal()];
 469 //                                      System.out.printf("    rl=%f Cf=%f (perfect=%b)\n",
 470 //                                      roughnessLimited[finish.ordinal()],
 471 //                                      Cf,rocket.isPerfectFinish());
 472
 473 //                                      System.out.printf("LIMITED  ");
 474                                 } else {
 475                                         componentCf = Cf;
 476 //                                      System.out.printf("NORMAL  ");
 477                                 }
 478
 479                         } else {
 480
 481                                 // For fully turbulent use simple max
 482                                 componentCf = Math.max(Cf, roughnessLimited[finish.ordinal()]);
 483
 484                         }
 485
 486 //                      System.out.printf("compCf=%.3f  ", componentCf);
 487
 488
 489
 490
 491                         // Calculate the friction drag:
 492                         if (c instanceof SymmetricComponent) {
 493
 494                                 SymmetricComponent s = (SymmetricComponent)c;
 495
 496                                 bodyFriction += componentCf * s.getComponentWetArea();
 497
 498                                 if (map != null) {
 499                                         // Corrected later
 500                                         map.get(c).frictionCD = componentCf * s.getComponentWetArea()
 501                                                 / conditions.getRefArea();
 502                                 }
 503
 504                                 double r = Math.max(s.getForeRadius(), s.getAftRadius());
 505                                 if (r > maxR)
 506                                         maxR = r;
 507                                 len += c.getLength();
 508
 509                         } else if (c instanceof FinSet) {
 510
 511                                 FinSet f = (FinSet)c;
 512                                 double mac = ((FinSetCalc)calcMap.get(c)).getMACLength();
 513                                 double cd = componentCf * (1 + 2*f.getThickness()/mac) *
 514                                         2*f.getFinCount() * f.getFinArea();
 515                                 finFriction += cd;
 516
 517                                 if (map != null) {
 518                                         map.get(c).frictionCD = cd / conditions.getRefArea();
 519                                 }
 520
 521                         }
 522
 523                 }
 524                 // fB may be POSITIVE_INFINITY, but that's ok for us
 525                 double fB = (len+0.0001) / maxR;
 526                 double correction = (1 + 1.0/(2*fB));
 527
 528                 // Correct body data in map
 529                 if (map != null) {
 530                         for (RocketComponent c: map.keySet()) {
 531                                 if (c instanceof SymmetricComponent) {
 532                                         map.get(c).frictionCD *= correction;
 533                                 }
 534                         }
 535                 }
 536
 537 //              System.out.printf("\n");
 538                 return (finFriction + correction*bodyFriction) / conditions.getRefArea();
 539         }
 540
 541
 542
 543         private double calculatePressureDrag(FlightConditions conditions,
 544                         Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
 545
 546                 double stagnation, base, total;
 547                 double radius = 0;
 548
 549                 if (calcMap == null)
 550                         buildCalcMap();
 551
 552                 stagnation = calculateStagnationCD(conditions.getMach());
 553                 base = calculateBaseCD(conditions.getMach());
 554
 555                 total = 0;
 556                 for (RocketComponent c: configuration) {
 557                         if (!c.isAerodynamic())
 558                                 continue;
 559
 560                         // Pressure fore drag
 561                         double cd = calcMap.get(c).calculatePressureDragForce(conditions, stagnation, base,
 562                                         warnings);
 563                         total += cd;
 564
 565                         if (map != null) {
 566                                 map.get(c).pressureCD = cd;
 567                         }
 568
 569
 570                         // Stagnation drag
 571                         if (c instanceof SymmetricComponent) {
 572                                 SymmetricComponent s = (SymmetricComponent)c;
 573
 574                                 if (radius < s.getForeRadius()) {
 575                                         double area = Math.PI*(pow2(s.getForeRadius()) - pow2(radius));
 576                                         cd = stagnation * area / conditions.getRefArea();
 577                                         total += cd;
 578                                         if (map != null) {
 579                                                 map.get(c).pressureCD += cd;
 580                                         }
 581                                 }
 582
 583                                 radius = s.getAftRadius();
 584                         }
 585                 }
 586
 587                 return total;
 588         }
 589
 590
 591         private double calculateBaseDrag(FlightConditions conditions,
 592                         Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
 593
 594                 double base, total;
 595                 double radius = 0;
 596                 RocketComponent prevComponent = null;
 597
 598                 if (calcMap == null)
 599                         buildCalcMap();
 600
 601                 base = calculateBaseCD(conditions.getMach());
 602                 total = 0;
 603
 604                 for (RocketComponent c: configuration) {
 605                         if (!(c instanceof SymmetricComponent))
 606                                 continue;
 607
 608                         SymmetricComponent s = (SymmetricComponent)c;
 609
 610                         if (radius > s.getForeRadius()) {
 611                                 double area = Math.PI*(pow2(radius) - pow2(s.getForeRadius()));
 612                                 double cd = base * area / conditions.getRefArea();
 613                                 total += cd;
 614                                 if (map != null) {
 615                                         map.get(prevComponent).baseCD = cd;
 616                                 }
 617                         }
 618
 619                         radius = s.getAftRadius();
 620                         prevComponent = c;
 621                 }
 622
 623                 if (radius > 0) {
 624                         double area = Math.PI*pow2(radius);
 625                         double cd = base * area / conditions.getRefArea();
 626                         total += cd;
 627                         if (map != null) {
 628                                 map.get(prevComponent).baseCD = cd;
 629                         }
 630                 }
 631
 632                 return total;
 633         }
 634
 635
 636
 637         public static double calculateStagnationCD(double m) {
 638                 double pressure;
 639                 if (m <=1) {
 640                         pressure = 1 + pow2(m)/4 + pow2(pow2(m))/40;
 641                 } else {
 642                         pressure = 1.84 - 0.76/pow2(m) + 0.166/pow2(pow2(m)) + 0.035/pow2(m*m*m);
 643                 }
 644                 return 0.85 * pressure;
 645         }
 646
 647
 648         public static double calculateBaseCD(double m) {
 649                 if (m <= 1) {
 650                         return 0.12 + 0.13 * m*m;
 651                 } else {
 652                         return 0.25 / m;
 653                 }
 654         }
 655
 656
 657
 658         private static final double[] axialDragPoly1, axialDragPoly2;
 659         static {
 660                 PolyInterpolator interpolator;
 661                 interpolator = new PolyInterpolator(
 662                                 new double[] { 0, 17*Math.PI/180 },
 663                                 new double[] { 0, 17*Math.PI/180 }
 664                 );
 665                 axialDragPoly1 = interpolator.interpolator(1, 1.3, 0, 0);
 666
 667                 interpolator = new PolyInterpolator(
 668                                 new double[] { 17*Math.PI/180, Math.PI/2 },
 669                                 new double[] { 17*Math.PI/180, Math.PI/2 },
 670                                 new double[] { Math.PI/2 }
 671                 );
 672                 axialDragPoly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
 673         }
 674
 675
 676         /**
 677          * Calculate the axial drag from the total drag coefficient.
 678          *
 679          * @param conditions
 680          * @param cd
 681          * @return
 682          */
 683         private double calculateAxialDrag(FlightConditions conditions, double cd) {
 684                 double aoa = MathUtil.clamp(conditions.getAOA(), 0, Math.PI);
 685                 double mul;
 686
 687 //              double sinaoa = conditions.getSinAOA();
 688 //              return cd * (1 + Math.min(sinaoa, 0.25));
 689
 690
 691                 if (aoa > Math.PI/2)
 692                         aoa = Math.PI - aoa;
 693                 if (aoa < 17*Math.PI/180)
 694                         mul = PolyInterpolator.eval(aoa, axialDragPoly1);
 695                 else
 696                         mul = PolyInterpolator.eval(aoa, axialDragPoly2);
 697
 698                 if (conditions.getAOA() < Math.PI/2)
 699                         return mul * cd;
 700                 else
 701                         return -mul * cd;
 702         }
 703
 704
 705         private void calculateDampingMoments(FlightConditions conditions,
 706                         AerodynamicForces total) {
 707
 708                 // Calculate pitch and yaw damping moments
 709                 if (conditions.getPitchRate() > 0.1 || conditions.getYawRate() > 0.1 || true) {
 710                         double mul = getDampingMultiplier(conditions, total.cg.x);
 711                         double pitch = conditions.getPitchRate();
 712                         double yaw = conditions.getYawRate();
 713                         double vel = conditions.getVelocity();
 714
 715 //                      double Cm = total.Cm - total.CN * total.cg.x / conditions.getRefLength();
 716 //                      System.out.printf("Damping pitch/yaw, mul=%.4f pitch rate=%.4f "+
 717 //                                      "Cm=%.4f / %.4f effect=%.4f aoa=%.4f\n", mul, pitch, total.Cm, Cm,
 718 //                                      -(mul * MathUtil.sign(pitch) * pow2(pitch/vel)),
 719 //                                      conditions.getAOA()*180/Math.PI);
 720
 721                         mul *= 3;   // TODO: Higher damping yields much more realistic apogee turn
 722
 723 //                      total.Cm -= mul * pitch / pow2(vel);
 724 //                      total.Cyaw -= mul * yaw / pow2(vel);
 725                         total.pitchDampingMoment = mul * MathUtil.sign(pitch) * pow2(pitch/vel);
 726                         total.yawDampingMoment = mul * MathUtil.sign(yaw) * pow2(yaw/vel);
 727                 } else {
 728                         total.pitchDampingMoment = 0;
 729                         total.yawDampingMoment = 0;
 730                 }
 731
 732         }
 733
 734         // TODO: MEDIUM: Are the rotation etc. being added correctly?  sin/cos theta?
 735
 736
 737         private double getDampingMultiplier(FlightConditions conditions, double cgx) {
 738                 if (cacheDiameter < 0) {
 739                         double area = 0;
 740                         cacheLength = 0;
 741                         cacheDiameter = 0;
 742
 743                         for (RocketComponent c: configuration) {
 744                                 if (c instanceof SymmetricComponent) {
 745                                         SymmetricComponent s = (SymmetricComponent)c;
 746                                         area += s.getComponentPlanformArea();
 747                                         cacheLength += s.getLength();
 748                                 }
 749                         }
 750                         if (cacheLength > 0)
 751                                 cacheDiameter = area / cacheLength;
 752                 }
 753
 754                 double mul;
 755
 756                 // Body
 757                 mul = 0.275 * cacheDiameter / (conditions.getRefArea() * conditions.getRefLength());
 758                 mul *= (MathUtil.pow4(cgx) + MathUtil.pow4(cacheLength - cgx));
 759
 760                 // Fins
 761                 // TODO: LOW: This could be optimized a lot...
 762                 for (RocketComponent c: configuration) {
 763                         if (c instanceof FinSet) {
 764                                 FinSet f = (FinSet)c;
 765                                 mul += 0.6 * Math.min(f.getFinCount(), 4) * f.getFinArea() *
 766                                                 MathUtil.pow3(Math.abs(f.toAbsolute(new Coordinate(
 767                                                                                 ((FinSetCalc)calcMap.get(f)).getMidchordPos()))[0].x
 768                                                                                 - cgx)) /
 769                                                                                 (conditions.getRefArea() * conditions.getRefLength());
 770                         }
 771                 }
 772
 773                 return mul;
 774         }
 775
 776
 777
 778         ////////  The calculator map
 779
 780         @Override
 781         protected void voidAerodynamicCache() {
 782                 super.voidAerodynamicCache();
 783
 784                 calcMap = null;
 785                 cacheDiameter = -1;
 786                 cacheLength = -1;
 787         }
 788
 789
 790         private void buildCalcMap() {
 791                 Iterator<RocketComponent> iterator;
 792
 793                 calcMap = new HashMap<RocketComponent, RocketComponentCalc>();
 794
 795                 iterator = rocket.deepIterator();
 796                 while (iterator.hasNext()) {
 797                         RocketComponent c = iterator.next();
 798
 799                         if (!c.isAerodynamic())
 800                                 continue;
 801
 802                         calcMap.put(c, (RocketComponentCalc) Reflection.construct(BARROWMAN_PACKAGE,
 803                                         c, BARROWMAN_SUFFIX, c));
 804                 }
 805         }
 806
 807
 808
 809
 810 //      public static void main(String[] arg) {
 811 //
 812 //              PolyInterpolator interpolator;
 813 //
 814 //              interpolator = new PolyInterpolator(
 815 //                              new double[] { 0, 17*Math.PI/180 },
 816 //                              new double[] { 0, 17*Math.PI/180 }
 817 //              );
 818 //              double[] poly1 = interpolator.interpolator(1, 1.3, 0, 0);
 819 //
 820 //              interpolator = new PolyInterpolator(
 821 //                              new double[] { 17*Math.PI/180, Math.PI/2 },
 822 //                              new double[] { 17*Math.PI/180, Math.PI/2 },
 823 //                              new double[] { Math.PI/2 }
 824 //              );
 825 //              double[] poly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
 826 //
 827 //
 828 //              for (double a=0; a<=180.1; a++) {
 829 //                      double r = a*Math.PI/180;
 830 //                      if (r > Math.PI/2)
 831 //                              r = Math.PI - r;
 832 //
 833 //                      double value;
 834 //                      if (r < 18*Math.PI/180)
 835 //                              value = PolyInterpolator.eval(r, poly1);
 836 //                      else
 837 //                              value = PolyInterpolator.eval(r, poly2);
 838 //
 839 //                      System.out.println(""+a+" "+value);
 840 //              }
 841 //
 842 //              System.exit(0);
 843 //
 844 //
 845 //              Rocket normal = TestRocket.makeRocket();
 846 //              Rocket perfect = TestRocket.makeRocket();
 847 //              normal.setPerfectFinish(false);
 848 //              perfect.setPerfectFinish(true);
 849 //
 850 //              Configuration confNormal = new Configuration(normal);
 851 //              Configuration confPerfect = new Configuration(perfect);
 852 //
 853 //              for (RocketComponent c: confNormal) {
 854 //                      if (c instanceof ExternalComponent) {
 855 //                              ((ExternalComponent)c).setFinish(Finish.NORMAL);
 856 //                      }
 857 //              }
 858 //              for (RocketComponent c: confPerfect) {
 859 //                      if (c instanceof ExternalComponent) {
 860 //                              ((ExternalComponent)c).setFinish(Finish.NORMAL);
 861 //                      }
 862 //              }
 863 //
 864 //
 865 //              confNormal.setToStage(0);
 866 //              confPerfect.setToStage(0);
 867 //
 868 //
 869 //
 870 //              BarrowmanCalculator calcNormal = new BarrowmanCalculator(confNormal);
 871 //              BarrowmanCalculator calcPerfect = new BarrowmanCalculator(confPerfect);
 872 //
 873 //              FlightConditions conditions = new FlightConditions(confNormal);
 874 //
 875 //              for (double mach=0; mach < 3; mach += 0.1) {
 876 //                      conditions.setMach(mach);
 877 //
 878 //                      Map<RocketComponent, AerodynamicForces> data =
 879 //                              calcNormal.getForceAnalysis(conditions, null);
 880 //                      AerodynamicForces forcesNormal = data.get(normal);
 881 //
 882 //                      data = calcPerfect.getForceAnalysis(conditions, null);
 883 //                      AerodynamicForces forcesPerfect = data.get(perfect);
 884 //
 885 //                      System.out.printf("%f %f %f %f %f %f %f\n",mach,
 886 //                                      forcesNormal.pressureCD, forcesPerfect.pressureCD,
 887 //                                      forcesNormal.frictionCD, forcesPerfect.frictionCD,
 888 //                                      forcesNormal.CD, forcesPerfect.CD);
 889 //              }
 890 //
 891 //
 892 //
 893 //      }
 894
 895 }