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