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.ExternalComponent.Finish;
  16 import net.sf.openrocket.rocketcomponent.FinSet;
  17 import net.sf.openrocket.rocketcomponent.RocketComponent;
  18 import net.sf.openrocket.rocketcomponent.SymmetricComponent;
  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 AbstractAerodynamicCalculator {
  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
  48         @Override
  49         public BarrowmanCalculator newInstance() {
  50                 return new BarrowmanCalculator();
  51         }
  52
  53
  54         /**
  55          * Calculate the CP according to the extended Barrowman method.
  56          */
  57         @Override
  58         public Coordinate getCP(Configuration configuration, FlightConditions conditions,
  59                         WarningSet warnings) {
  60                 checkCache(configuration);
  61                 AerodynamicForces forces = getNonAxial(configuration, conditions, null, warnings);
  62                 return forces.getCP();
  63         }
  64
  65
  66
  67         @Override
  68         public Map<RocketComponent, AerodynamicForces> getForceAnalysis(Configuration configuration,
  69                         FlightConditions conditions, WarningSet warnings) {
  70                 checkCache(configuration);
  71
  72                 AerodynamicForces f;
  73                 Map<RocketComponent, AerodynamicForces> map =
  74                                 new LinkedHashMap<RocketComponent, AerodynamicForces>();
  75
  76                 // Add all components to the map
  77                 for (RocketComponent c : configuration) {
  78                         f = new AerodynamicForces();
  79                         f.setComponent(c);
  80
  81                         map.put(c, f);
  82                 }
  83
  84
  85                 // Calculate non-axial force data
  86                 AerodynamicForces total = getNonAxial(configuration, conditions, map, warnings);
  87
  88
  89                 // Calculate friction data
  90                 total.setFrictionCD(calculateFrictionDrag(configuration, conditions, map, warnings));
  91                 total.setPressureCD(calculatePressureDrag(configuration, conditions, map, warnings));
  92                 total.setBaseCD(calculateBaseDrag(configuration, conditions, map, warnings));
  93
  94                 total.setComponent(configuration.getRocket());
  95                 map.put(total.getComponent(), total);
  96
  97
  98                 for (RocketComponent c : map.keySet()) {
  99                         f = map.get(c);
 100                         if (Double.isNaN(f.getBaseCD()) && Double.isNaN(f.getPressureCD()) &&
 101                                         Double.isNaN(f.getFrictionCD()))
 102                                 continue;
 103                         if (Double.isNaN(f.getBaseCD()))
 104                                 f.setBaseCD(0);
 105                         if (Double.isNaN(f.getPressureCD()))
 106                                 f.setPressureCD(0);
 107                         if (Double.isNaN(f.getFrictionCD()))
 108                                 f.setFrictionCD(0);
 109                         f.setCD(f.getBaseCD() + f.getPressureCD() + f.getFrictionCD());
 110                         f.setCaxial(calculateAxialDrag(conditions, f.getCD()));
 111                 }
 112
 113                 return map;
 114         }
 115
 116
 117
 118         @Override
 119         public AerodynamicForces getAerodynamicForces(Configuration configuration,
 120                         FlightConditions conditions, WarningSet warnings) {
 121                 checkCache(configuration);
 122
 123                 if (warnings == null)
 124                         warnings = ignoreWarningSet;
 125
 126                 // Calculate non-axial force data
 127                 AerodynamicForces total = getNonAxial(configuration, conditions, null, warnings);
 128
 129                 // Calculate friction data
 130                 total.setFrictionCD(calculateFrictionDrag(configuration, conditions, null, warnings));
 131                 total.setPressureCD(calculatePressureDrag(configuration, conditions, null, warnings));
 132                 total.setBaseCD(calculateBaseDrag(configuration, conditions, null, warnings));
 133
 134                 total.setCD(total.getFrictionCD() + total.getPressureCD() + total.getBaseCD());
 135
 136                 total.setCaxial(calculateAxialDrag(conditions, total.getCD()));
 137
 138                 // Calculate pitch and yaw damping moments
 139                 calculateDampingMoments(configuration, conditions, total);
 140                 total.setCm(total.getCm() - total.getPitchDampingMoment());
 141                 total.setCyaw(total.getCyaw() - total.getYawDampingMoment());
 142
 143                 return total;
 144         }
 145
 146
 147
 148
 149         /*
 150          * Perform the actual CP calculation.
 151          */
 152         private AerodynamicForces getNonAxial(Configuration configuration, FlightConditions conditions,
 153                         Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
 154
 155                 checkCache(configuration);
 156
 157                 AerodynamicForces total = new AerodynamicForces();
 158                 total.zero();
 159
 160                 double radius = 0; // aft radius of previous component
 161                 double componentX = 0; // aft coordinate of previous component
 162                 AerodynamicForces forces = new AerodynamicForces();
 163
 164                 if (warnings == null)
 165                         warnings = ignoreWarningSet;
 166
 167                 if (conditions.getAOA() > 17.5 * Math.PI / 180)
 168                         warnings.add(new Warning.LargeAOA(conditions.getAOA()));
 169
 170
 171                 if (calcMap == null)
 172                         buildCalcMap(configuration);
 173
 174                 for (RocketComponent component : configuration) {
 175
 176                         // Skip non-aerodynamic components
 177                         if (!component.isAerodynamic())
 178                                 continue;
 179
 180                         // Check for discontinuities
 181                         if (component instanceof SymmetricComponent) {
 182                                 SymmetricComponent sym = (SymmetricComponent) component;
 183                                 // TODO:LOW: Ignores other cluster components (not clusterable)
 184                                 double x = component.toAbsolute(Coordinate.NUL)[0].x;
 185
 186                                 // Check for lengthwise discontinuity
 187                                 if (x > componentX + 0.0001) {
 188                                         if (!MathUtil.equals(radius, 0)) {
 189                                                 warnings.add(Warning.DISCONTINUITY);
 190                                                 radius = 0;
 191                                         }
 192                                 }
 193                                 componentX = component.toAbsolute(new Coordinate(component.getLength()))[0].x;
 194
 195                                 // Check for radius discontinuity
 196                                 if (!MathUtil.equals(sym.getForeRadius(), radius)) {
 197                                         warnings.add(Warning.DISCONTINUITY);
 198                                         // TODO: MEDIUM: Apply correction to values to cp and to map
 199                                 }
 200                                 radius = sym.getAftRadius();
 201                         }
 202
 203                         // Call calculation method
 204                         forces.zero();
 205                         calcMap.get(component).calculateNonaxialForces(conditions, forces, warnings);
 206                         forces.setCP(component.toAbsolute(forces.getCP())[0]);
 207                         forces.setCm(forces.getCN() * forces.getCP().x / conditions.getRefLength());
 208                         //                      System.out.println("  CN="+forces.CN+" cp.x="+forces.cp.x+" Cm="+forces.Cm);
 209
 210                         if (map != null) {
 211                                 AerodynamicForces f = map.get(component);
 212
 213                                 f.setCP(forces.getCP());
 214                                 f.setCNa(forces.getCNa());
 215                                 f.setCN(forces.getCN());
 216                                 f.setCm(forces.getCm());
 217                                 f.setCside(forces.getCside());
 218                                 f.setCyaw(forces.getCyaw());
 219                                 f.setCroll(forces.getCroll());
 220                                 f.setCrollDamp(forces.getCrollDamp());
 221                                 f.setCrollForce(forces.getCrollForce());
 222                         }
 223
 224                         total.setCP(total.getCP().average(forces.getCP()));
 225                         total.setCNa(total.getCNa() + forces.getCNa());
 226                         total.setCN(total.getCN() + forces.getCN());
 227                         total.setCm(total.getCm() + forces.getCm());
 228                         total.setCside(total.getCside() + forces.getCside());
 229                         total.setCyaw(total.getCyaw() + forces.getCyaw());
 230                         total.setCroll(total.getCroll() + forces.getCroll());
 231                         total.setCrollDamp(total.getCrollDamp() + forces.getCrollDamp());
 232                         total.setCrollForce(total.getCrollForce() + forces.getCrollForce());
 233                 }
 234
 235                 return total;
 236         }
 237
 238
 239
 240
 241         ////////////////  DRAG CALCULATIONS  ////////////////
 242
 243
 244         private double calculateFrictionDrag(Configuration configuration, FlightConditions conditions,
 245                         Map<RocketComponent, AerodynamicForces> map, WarningSet set) {
 246                 double c1 = 1.0, c2 = 1.0;
 247
 248                 double mach = conditions.getMach();
 249                 double Re;
 250                 double Cf;
 251
 252                 if (calcMap == null)
 253                         buildCalcMap(configuration);
 254
 255                 Re = conditions.getVelocity() * configuration.getLength() /
 256                                 conditions.getAtmosphericConditions().getKinematicViscosity();
 257
 258                 //              System.out.printf("Re=%.3e   ", Re);
 259
 260                 // Calculate the skin friction coefficient (assume non-roughness limited)
 261                 if (configuration.getRocket().isPerfectFinish()) {
 262
 263                         //                      System.out.printf("Perfect finish: Re=%f ",Re);
 264                         // Assume partial laminar layer.  Roughness-limitation is checked later.
 265                         if (Re < 1e4) {
 266                                 // Too low, constant
 267                                 Cf = 1.33e-2;
 268                                 //                              System.out.printf("constant Cf=%f ",Cf);
 269                         } else if (Re < 5.39e5) {
 270                                 // Fully laminar
 271                                 Cf = 1.328 / MathUtil.safeSqrt(Re);
 272                                 //                              System.out.printf("basic Cf=%f ",Cf);
 273                         } else {
 274                                 // Transitional
 275                                 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6) - 1700 / Re;
 276                                 //                              System.out.printf("transitional Cf=%f ",Cf);
 277                         }
 278
 279                         // Compressibility correction
 280
 281                         if (mach < 1.1) {
 282                                 // Below Re=1e6 no correction
 283                                 if (Re > 1e6) {
 284                                         if (Re < 3e6) {
 285                                                 c1 = 1 - 0.1 * pow2(mach) * (Re - 1e6) / 2e6; // transition to turbulent
 286                                         } else {
 287                                                 c1 = 1 - 0.1 * pow2(mach);
 288                                         }
 289                                 }
 290                         }
 291                         if (mach > 0.9) {
 292                                 if (Re > 1e6) {
 293                                         if (Re < 3e6) {
 294                                                 c2 = 1 + (1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25) - 1) * (Re - 1e6) / 2e6;
 295                                         } else {
 296                                                 c2 = 1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25);
 297                                         }
 298                                 }
 299                         }
 300
 301                         //                      System.out.printf("c1=%f c2=%f\n", c1,c2);
 302                         // Applying continuously around Mach 1
 303                         if (mach < 0.9) {
 304                                 Cf *= c1;
 305                         } else if (mach < 1.1) {
 306                                 Cf *= (c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2);
 307                         } else {
 308                                 Cf *= c2;
 309                         }
 310
 311                         //                      System.out.printf("M=%f Cf=%f (smooth)\n",mach,Cf);
 312
 313                 } else {
 314
 315                         // Assume fully turbulent.  Roughness-limitation is checked later.
 316                         if (Re < 1e4) {
 317                                 // Too low, constant
 318                                 Cf = 1.48e-2;
 319                                 //                              System.out.printf("LOW-TURB  ");
 320                         } else {
 321                                 // Turbulent
 322                                 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6);
 323                                 //                              System.out.printf("NORMAL-TURB  ");
 324                         }
 325
 326                         // Compressibility correction
 327
 328                         if (mach < 1.1) {
 329                                 c1 = 1 - 0.1 * pow2(mach);
 330                         }
 331                         if (mach > 0.9) {
 332                                 c2 = 1 / Math.pow(1 + 0.15 * pow2(mach), 0.58);
 333                         }
 334                         // Applying continuously around Mach 1
 335                         if (mach < 0.9) {
 336                                 Cf *= c1;
 337                         } else if (mach < 1.1) {
 338                                 Cf *= c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
 339                         } else {
 340                                 Cf *= c2;
 341                         }
 342
 343                         //                      System.out.printf("M=%f, Cd=%f (turbulent)\n", mach,Cf);
 344
 345                 }
 346
 347                 // Roughness-limited value correction term
 348                 double roughnessCorrection;
 349                 if (mach < 0.9) {
 350                         roughnessCorrection = 1 - 0.1 * pow2(mach);
 351                 } else if (mach > 1.1) {
 352                         roughnessCorrection = 1 / (1 + 0.18 * pow2(mach));
 353                 } else {
 354                         c1 = 1 - 0.1 * pow2(0.9);
 355                         c2 = 1.0 / (1 + 0.18 * pow2(1.1));
 356                         roughnessCorrection = c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
 357                 }
 358
 359                 //              System.out.printf("Cf=%.3f  ", Cf);
 360
 361
 362                 /*
 363                  * Calculate the friction drag coefficient.
 364                  *
 365                  * The body wetted area is summed up and finally corrected with the rocket
 366                  * fineness ratio (calculated in the same iteration).  The fins are corrected
 367                  * for thickness as we go on.
 368                  */
 369
 370                 double finFriction = 0;
 371                 double bodyFriction = 0;
 372                 double maxR = 0, len = 0;
 373
 374                 double[] roughnessLimited = new double[Finish.values().length];
 375                 Arrays.fill(roughnessLimited, Double.NaN);
 376
 377                 for (RocketComponent c : configuration) {
 378
 379                         // Consider only SymmetricComponents and FinSets:
 380                         if (!(c instanceof SymmetricComponent) &&
 381                                         !(c instanceof FinSet))
 382                                 continue;
 383
 384                         // Calculate the roughness-limited friction coefficient
 385                         Finish finish = ((ExternalComponent) c).getFinish();
 386                         if (Double.isNaN(roughnessLimited[finish.ordinal()])) {
 387                                 roughnessLimited[finish.ordinal()] =
 388                                                 0.032 * Math.pow(finish.getRoughnessSize() / configuration.getLength(), 0.2) *
 389                                                                 roughnessCorrection;
 390
 391                                 //                              System.out.printf("roughness["+finish+"]=%.3f  ",
 392                                 //                                              roughnessLimited[finish.ordinal()]);
 393                         }
 394
 395                         /*
 396                          * Actual Cf is maximum of Cf and the roughness-limited value.
 397                          * For perfect finish require additionally that Re > 1e6
 398                          */
 399                         double componentCf;
 400                         if (configuration.getRocket().isPerfectFinish()) {
 401
 402                                 // For perfect finish require Re > 1e6
 403                                 if ((Re > 1.0e6) && (roughnessLimited[finish.ordinal()] > Cf)) {
 404                                         componentCf = roughnessLimited[finish.ordinal()];
 405                                         //                                      System.out.printf("    rl=%f Cf=%f (perfect=%b)\n",
 406                                         //                                      roughnessLimited[finish.ordinal()],
 407                                         //                                      Cf,rocket.isPerfectFinish());
 408
 409                                         //                                      System.out.printf("LIMITED  ");
 410                                 } else {
 411                                         componentCf = Cf;
 412                                         //                                      System.out.printf("NORMAL  ");
 413                                 }
 414
 415                         } else {
 416
 417                                 // For fully turbulent use simple max
 418                                 componentCf = Math.max(Cf, roughnessLimited[finish.ordinal()]);
 419
 420                         }
 421
 422                         //                      System.out.printf("compCf=%.3f  ", componentCf);
 423
 424
 425
 426
 427                         // Calculate the friction drag:
 428                         if (c instanceof SymmetricComponent) {
 429
 430                                 SymmetricComponent s = (SymmetricComponent) c;
 431
 432                                 bodyFriction += componentCf * s.getComponentWetArea();
 433
 434                                 if (map != null) {
 435                                         // Corrected later
 436                                         map.get(c).setFrictionCD(componentCf * s.getComponentWetArea()
 437                                                         / conditions.getRefArea());
 438                                 }
 439
 440                                 double r = Math.max(s.getForeRadius(), s.getAftRadius());
 441                                 if (r > maxR)
 442                                         maxR = r;
 443                                 len += c.getLength();
 444
 445                         } else if (c instanceof FinSet) {
 446
 447                                 FinSet f = (FinSet) c;
 448                                 double mac = ((FinSetCalc) calcMap.get(c)).getMACLength();
 449                                 double cd = componentCf * (1 + 2 * f.getThickness() / mac) *
 450                                                 2 * f.getFinCount() * f.getFinArea();
 451                                 finFriction += cd;
 452
 453                                 if (map != null) {
 454                                         map.get(c).setFrictionCD(cd / conditions.getRefArea());
 455                                 }
 456
 457                         }
 458
 459                 }
 460                 // fB may be POSITIVE_INFINITY, but that's ok for us
 461                 double fB = (len + 0.0001) / maxR;
 462                 double correction = (1 + 1.0 / (2 * fB));
 463
 464                 // Correct body data in map
 465                 if (map != null) {
 466                         for (RocketComponent c : map.keySet()) {
 467                                 if (c instanceof SymmetricComponent) {
 468                                         map.get(c).setFrictionCD(map.get(c).getFrictionCD() * correction);
 469                                 }
 470                         }
 471                 }
 472
 473                 //              System.out.printf("\n");
 474                 return (finFriction + correction * bodyFriction) / conditions.getRefArea();
 475         }
 476
 477
 478
 479         private double calculatePressureDrag(Configuration configuration, FlightConditions conditions,
 480                         Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
 481
 482                 double stagnation, base, total;
 483                 double radius = 0;
 484
 485                 if (calcMap == null)
 486                         buildCalcMap(configuration);
 487
 488                 stagnation = calculateStagnationCD(conditions.getMach());
 489                 base = calculateBaseCD(conditions.getMach());
 490
 491                 total = 0;
 492                 for (RocketComponent c : configuration) {
 493                         if (!c.isAerodynamic())
 494                                 continue;
 495
 496                         // Pressure fore drag
 497                         double cd = calcMap.get(c).calculatePressureDragForce(conditions, stagnation, base,
 498                                         warnings);
 499                         total += cd;
 500
 501                         if (map != null) {
 502                                 map.get(c).setPressureCD(cd);
 503                         }
 504
 505
 506                         // Stagnation drag
 507                         if (c instanceof SymmetricComponent) {
 508                                 SymmetricComponent s = (SymmetricComponent) c;
 509
 510                                 if (radius < s.getForeRadius()) {
 511                                         double area = Math.PI * (pow2(s.getForeRadius()) - pow2(radius));
 512                                         cd = stagnation * area / conditions.getRefArea();
 513                                         total += cd;
 514                                         if (map != null) {
 515                                                 map.get(c).setPressureCD(map.get(c).getPressureCD() + cd);
 516                                         }
 517                                 }
 518
 519                                 radius = s.getAftRadius();
 520                         }
 521                 }
 522
 523                 return total;
 524         }
 525
 526
 527         private double calculateBaseDrag(Configuration configuration, FlightConditions conditions,
 528                         Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
 529
 530                 double base, total;
 531                 double radius = 0;
 532                 RocketComponent prevComponent = null;
 533
 534                 if (calcMap == null)
 535                         buildCalcMap(configuration);
 536
 537                 base = calculateBaseCD(conditions.getMach());
 538                 total = 0;
 539
 540                 for (RocketComponent c : configuration) {
 541                         if (!(c instanceof SymmetricComponent))
 542                                 continue;
 543
 544                         SymmetricComponent s = (SymmetricComponent) c;
 545
 546                         if (radius > s.getForeRadius()) {
 547                                 double area = Math.PI * (pow2(radius) - pow2(s.getForeRadius()));
 548                                 double cd = base * area / conditions.getRefArea();
 549                                 total += cd;
 550                                 if (map != null) {
 551                                         map.get(prevComponent).setBaseCD(cd);
 552                                 }
 553                         }
 554
 555                         radius = s.getAftRadius();
 556                         prevComponent = c;
 557                 }
 558
 559                 if (radius > 0) {
 560                         double area = Math.PI * pow2(radius);
 561                         double cd = base * area / conditions.getRefArea();
 562                         total += cd;
 563                         if (map != null) {
 564                                 map.get(prevComponent).setBaseCD(cd);
 565                         }
 566                 }
 567
 568                 return total;
 569         }
 570
 571
 572
 573         public static double calculateStagnationCD(double m) {
 574                 double pressure;
 575                 if (m <= 1) {
 576                         pressure = 1 + pow2(m) / 4 + pow2(pow2(m)) / 40;
 577                 } else {
 578                         pressure = 1.84 - 0.76 / pow2(m) + 0.166 / pow2(pow2(m)) + 0.035 / pow2(m * m * m);
 579                 }
 580                 return 0.85 * pressure;
 581         }
 582
 583
 584         public static double calculateBaseCD(double m) {
 585                 if (m <= 1) {
 586                         return 0.12 + 0.13 * m * m;
 587                 } else {
 588                         return 0.25 / m;
 589                 }
 590         }
 591
 592
 593
 594         private static final double[] axialDragPoly1, axialDragPoly2;
 595         static {
 596                 PolyInterpolator interpolator;
 597                 interpolator = new PolyInterpolator(
 598                                 new double[] { 0, 17 * Math.PI / 180 },
 599                                 new double[] { 0, 17 * Math.PI / 180 }
 600                                 );
 601                 axialDragPoly1 = interpolator.interpolator(1, 1.3, 0, 0);
 602
 603                 interpolator = new PolyInterpolator(
 604                                 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
 605                                 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
 606                                 new double[] { Math.PI / 2 }
 607                                 );
 608                 axialDragPoly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
 609         }
 610
 611
 612         /**
 613          * Calculate the axial drag from the total drag coefficient.
 614          *
 615          * @param conditions
 616          * @param cd
 617          * @return
 618          */
 619         private double calculateAxialDrag(FlightConditions conditions, double cd) {
 620                 double aoa = MathUtil.clamp(conditions.getAOA(), 0, Math.PI);
 621                 double mul;
 622
 623                 //              double sinaoa = conditions.getSinAOA();
 624                 //              return cd * (1 + Math.min(sinaoa, 0.25));
 625
 626
 627                 if (aoa > Math.PI / 2)
 628                         aoa = Math.PI - aoa;
 629                 if (aoa < 17 * Math.PI / 180)
 630                         mul = PolyInterpolator.eval(aoa, axialDragPoly1);
 631                 else
 632                         mul = PolyInterpolator.eval(aoa, axialDragPoly2);
 633
 634                 if (conditions.getAOA() < Math.PI / 2)
 635                         return mul * cd;
 636                 else
 637                         return -mul * cd;
 638         }
 639
 640
 641         private void calculateDampingMoments(Configuration configuration, FlightConditions conditions,
 642                         AerodynamicForces total) {
 643
 644                 // Calculate pitch and yaw damping moments
 645                 double mul = getDampingMultiplier(configuration, conditions,
 646                                 conditions.getPitchCenter().x);
 647                 double pitch = conditions.getPitchRate();
 648                 double yaw = conditions.getYawRate();
 649                 double vel = conditions.getVelocity();
 650
 651                 //                      double Cm = total.Cm - total.CN * total.cg.x / conditions.getRefLength();
 652                 //                      System.out.printf("Damping pitch/yaw, mul=%.4f pitch rate=%.4f "+
 653                 //                                      "Cm=%.4f / %.4f effect=%.4f aoa=%.4f\n", mul, pitch, total.Cm, Cm,
 654                 //                                      -(mul * MathUtil.sign(pitch) * pow2(pitch/vel)),
 655                 //                                      conditions.getAOA()*180/Math.PI);
 656
 657                 mul *= 3; // TODO: Higher damping yields much more realistic apogee turn
 658
 659                 //                      total.Cm -= mul * pitch / pow2(vel);
 660                 //                      total.Cyaw -= mul * yaw / pow2(vel);
 661                 total.setPitchDampingMoment(mul * MathUtil.sign(pitch) * pow2(pitch / vel));
 662                 total.setYawDampingMoment(mul * MathUtil.sign(yaw) * pow2(yaw / vel));
 663
 664         }
 665
 666         // TODO: MEDIUM: Are the rotation etc. being added correctly?  sin/cos theta?
 667
 668
 669         private double getDampingMultiplier(Configuration configuration, FlightConditions conditions,
 670                         double cgx) {
 671                 if (cacheDiameter < 0) {
 672                         double area = 0;
 673                         cacheLength = 0;
 674                         cacheDiameter = 0;
 675
 676                         for (RocketComponent c : configuration) {
 677                                 if (c instanceof SymmetricComponent) {
 678                                         SymmetricComponent s = (SymmetricComponent) c;
 679                                         area += s.getComponentPlanformArea();
 680                                         cacheLength += s.getLength();
 681                                 }
 682                         }
 683                         if (cacheLength > 0)
 684                                 cacheDiameter = area / cacheLength;
 685                 }
 686
 687                 double mul;
 688
 689                 // Body
 690                 mul = 0.275 * cacheDiameter / (conditions.getRefArea() * conditions.getRefLength());
 691                 mul *= (MathUtil.pow4(cgx) + MathUtil.pow4(cacheLength - cgx));
 692
 693                 // Fins
 694                 // TODO: LOW: This could be optimized a lot...
 695                 for (RocketComponent c : configuration) {
 696                         if (c instanceof FinSet) {
 697                                 FinSet f = (FinSet) c;
 698                                 mul += 0.6 * Math.min(f.getFinCount(), 4) * f.getFinArea() *
 699                                                 MathUtil.pow3(Math.abs(f.toAbsolute(new Coordinate(
 700                                                                                 ((FinSetCalc) calcMap.get(f)).getMidchordPos()))[0].x
 701                                                                                 - cgx)) /
 702                                                                                 (conditions.getRefArea() * conditions.getRefLength());
 703                         }
 704                 }
 705
 706                 return mul;
 707         }
 708
 709
 710
 711         ////////  The calculator map
 712
 713         @Override
 714         protected void voidAerodynamicCache() {
 715                 super.voidAerodynamicCache();
 716
 717                 calcMap = null;
 718                 cacheDiameter = -1;
 719                 cacheLength = -1;
 720         }
 721
 722
 723         private void buildCalcMap(Configuration configuration) {
 724                 Iterator<RocketComponent> iterator;
 725
 726                 calcMap = new HashMap<RocketComponent, RocketComponentCalc>();
 727
 728                 iterator = configuration.getRocket().iterator();
 729                 while (iterator.hasNext()) {
 730                         RocketComponent c = iterator.next();
 731
 732                         if (!c.isAerodynamic())
 733                                 continue;
 734
 735                         calcMap.put(c, (RocketComponentCalc) Reflection.construct(BARROWMAN_PACKAGE,
 736                                         c, BARROWMAN_SUFFIX, c));
 737                 }
 738         }
 739
 740
 741         @Override
 742         public int getModID() {
 743                 // Only cached data is stored, return constant mod ID
 744                 return 0;
 745         }
 746
 747
 748 }