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