git.gag.com Git - debian/openrocket/blob - core/src/net/sf/openrocket/rocketcomponent/SymmetricComponent.java

   1 package net.sf.openrocket.rocketcomponent;
   2
   3 import static net.sf.openrocket.util.MathUtil.pow2;
   4
   5 import java.util.ArrayList;
   6 import java.util.Collection;
   7 import java.util.List;
   8
   9 import net.sf.openrocket.preset.ComponentPreset;
  10 import net.sf.openrocket.util.Coordinate;
  11 import net.sf.openrocket.util.MathUtil;
  12
  13
  14 /**
  15  * Class for an axially symmetric rocket component generated by rotating
  16  * a function y=f(x) >= 0 around the x-axis (eg. tube, cone, etc.)
  17  *
  18  * @author Sampo Niskanen <sampo.niskanen@iki.fi>
  19  */
  20
  21 public abstract class SymmetricComponent extends BodyComponent implements RadialParent {
  22         public static final double DEFAULT_RADIUS = 0.025;
  23         public static final double DEFAULT_THICKNESS = 0.002;
  24
  25         private static final int DIVISIONS = 100; // No. of divisions when integrating
  26
  27         protected boolean filled = false;
  28         protected double thickness = DEFAULT_THICKNESS;
  29
  30
  31         // Cached data, default values signify not calculated
  32         private double wetArea = -1;
  33         private double planArea = -1;
  34         private double planCenter = -1;
  35         private double volume = -1;
  36         private double fullVolume = -1;
  37         private double longitudinalInertia = -1;
  38         private double rotationalInertia = -1;
  39         private Coordinate cg = null;
  40
  41
  42
  43         public SymmetricComponent() {
  44                 super();
  45         }
  46
  47
  48         /**
  49          * Return the component radius at position x.
  50          * @param x Position on x-axis.
  51          * @return  Radius of the component at the given position, or 0 if outside
  52          *          the component.
  53          */
  54         public abstract double getRadius(double x);
  55
  56         @Override
  57         public abstract double getInnerRadius(double x);
  58
  59         public abstract double getForeRadius();
  60
  61         public abstract boolean isForeRadiusAutomatic();
  62
  63         public abstract double getAftRadius();
  64
  65         public abstract boolean isAftRadiusAutomatic();
  66
  67
  68         // Implement the Radial interface:
  69         @Override
  70         public final double getOuterRadius(double x) {
  71                 return getRadius(x);
  72         }
  73
  74
  75         @Override
  76         public final double getRadius(double x, double theta) {
  77                 return getRadius(x);
  78         }
  79
  80         @Override
  81         public final double getInnerRadius(double x, double theta) {
  82                 return getInnerRadius(x);
  83         }
  84
  85
  86
  87         /**
  88          * Return the component wall thickness.
  89          */
  90         public double getThickness() {
  91                 if (filled)
  92                         return Math.max(getForeRadius(), getAftRadius());
  93                 return Math.min(thickness, Math.max(getForeRadius(), getAftRadius()));
  94         }
  95
  96
  97         /**
  98          * Set the component wall thickness.  Values greater than the maximum radius are not
  99          * allowed, and will result in setting the thickness to the maximum radius.
 100          */
 101         public void setThickness(double thickness) {
 102                 if ((this.thickness == thickness) && !filled)
 103                         return;
 104                 this.thickness = MathUtil.clamp(thickness, 0, Math.max(getForeRadius(), getAftRadius()));
 105                 filled = false;
 106                 fireComponentChangeEvent(ComponentChangeEvent.MASS_CHANGE);
 107                 clearPreset();
 108         }
 109
 110
 111         /**
 112          * Returns whether the component is set as filled.  If it is set filled, then the
 113          * wall thickness will have no effect.
 114          */
 115         public boolean isFilled() {
 116                 return filled;
 117         }
 118
 119
 120         /**
 121          * Sets whether the component is set as filled.  If the component is filled, then
 122          * the wall thickness will have no effect.
 123          */
 124         public void setFilled(boolean filled) {
 125                 if (this.filled == filled)
 126                         return;
 127                 this.filled = filled;
 128                 fireComponentChangeEvent(ComponentChangeEvent.MASS_CHANGE);
 129                 clearPreset();
 130         }
 131
 132
 133         /**
 134          * Adds component bounds at a number of points between 0...length.
 135          */
 136         @Override
 137         public Collection<Coordinate> getComponentBounds() {
 138                 List<Coordinate> list = new ArrayList<Coordinate>(20);
 139                 for (int n = 0; n <= 5; n++) {
 140                         double x = n * length / 5;
 141                         double r = getRadius(x);
 142                         addBound(list, x, r);
 143                 }
 144                 return list;
 145         }
 146
 147
 148
 149         @Override
 150         protected void loadFromPreset(ComponentPreset preset) {
 151                 if ( preset.has(ComponentPreset.THICKNESS) ) {
 152                         this.thickness = preset.get(ComponentPreset.THICKNESS);
 153                         this.filled = false;
 154                 }
 155                 if ( preset.has(ComponentPreset.FILLED)) {
 156                         this.filled = true;
 157                 }
 158
 159                 super.loadFromPreset(preset);
 160         }
 161
 162
 163
 164
 165         /**
 166          * Calculate volume of the component by integrating over the length of the component.
 167          * The method caches the result, so subsequent calls are instant.  Subclasses may
 168          * override this method for simple shapes and use this method as necessary.
 169          *
 170          * @return  The volume of the component.
 171          */
 172         @Override
 173         public double getComponentVolume() {
 174                 if (volume < 0)
 175                         integrate();
 176                 return volume;
 177         }
 178
 179
 180         /**
 181          * Calculate full (filled) volume of the component by integrating over the length
 182          * of the component.  The method caches the result, so subsequent calls are instant.
 183          * Subclasses may override this method for simple shapes and use this method as
 184          * necessary.
 185          *
 186          * @return  The filled volume of the component.
 187          */
 188         public double getFullVolume() {
 189                 if (fullVolume < 0)
 190                         integrate();
 191                 return fullVolume;
 192         }
 193
 194
 195         /**
 196          * Calculate the wetted area of the component by integrating over the length
 197          * of the component.  The method caches the result, so subsequent calls are instant.
 198          * Subclasses may override this method for simple shapes and use this method as
 199          * necessary.
 200          *
 201          * @return  The wetted area of the component.
 202          */
 203         public double getComponentWetArea() {
 204                 if (wetArea < 0)
 205                         integrate();
 206                 return wetArea;
 207         }
 208
 209
 210         /**
 211          * Calculate the planform area of the component by integrating over the length of
 212          * the component.  The method caches the result, so subsequent calls are instant.
 213          * Subclasses may override this method for simple shapes and use this method as
 214          * necessary.
 215          *
 216          * @return  The planform area of the component.
 217          */
 218         public double getComponentPlanformArea() {
 219                 if (planArea < 0)
 220                         integrate();
 221                 return planArea;
 222         }
 223
 224
 225         /**
 226          * Calculate the planform center X-coordinate of the component by integrating over
 227          * the length of the component.  The planform center is defined as
 228          * <pre>   integrate(x*2*r(x)) / planform area  </pre>
 229          * The method caches the result, so subsequent calls are instant.  Subclasses may
 230          * override this method for simple shapes and use this method as necessary.
 231          *
 232          * @return  The planform center of the component.
 233          */
 234         public double getComponentPlanformCenter() {
 235                 if (planCenter < 0)
 236                         integrate();
 237                 return planCenter;
 238         }
 239
 240
 241         /**
 242          * Calculate CG of the component by integrating over the length of the component.
 243          * The method caches the result, so subsequent calls are instant.  Subclasses may
 244          * override this method for simple shapes and use this method as necessary.
 245          *
 246          * @return  The CG+mass of the component.
 247          */
 248         @Override
 249         public Coordinate getComponentCG() {
 250                 if (cg == null)
 251                         integrate();
 252                 return cg;
 253         }
 254
 255
 256         @Override
 257         public double getLongitudinalUnitInertia() {
 258                 if (longitudinalInertia < 0) {
 259                         if (getComponentVolume() > 0.0000001) // == 0.1cm^3
 260                                 integrateInertiaVolume();
 261                         else
 262                                 integrateInertiaSurface();
 263                 }
 264                 return longitudinalInertia;
 265         }
 266
 267
 268         @Override
 269         public double getRotationalUnitInertia() {
 270                 if (rotationalInertia < 0) {
 271                         if (getComponentVolume() > 0.0000001)
 272                                 integrateInertiaVolume();
 273                         else
 274                                 integrateInertiaSurface();
 275                 }
 276                 return rotationalInertia;
 277         }
 278
 279
 280
 281         /**
 282          * Performs integration over the length of the component and updates the cached variables.
 283          */
 284         private void integrate() {
 285                 double x, r1, r2;
 286                 double cgx;
 287
 288                 // Check length > 0
 289                 if (length <= 0) {
 290                         wetArea = 0;
 291                         planArea = 0;
 292                         planCenter = 0;
 293                         volume = 0;
 294                         cg = Coordinate.NUL;
 295                         return;
 296                 }
 297
 298
 299                 // Integrate for volume, CG, wetted area and planform area
 300
 301                 final double step = length / DIVISIONS;
 302                 final double pi3 = Math.PI / 3.0;
 303                 r1 = getRadius(0);
 304                 x = 0;
 305                 wetArea = 0;
 306                 planArea = 0;
 307                 planCenter = 0;
 308                 fullVolume = 0;
 309                 volume = 0;
 310                 cgx = 0;
 311
 312                 for (int n = 1; n <= DIVISIONS; n++) {
 313                         /*
 314                          * r1 and r2 are the two radii
 315                          * x is the position of r1
 316                          * hyp is the length of the hypotenuse from r1 to r2
 317                          * height if the y-axis height of the component if not filled
 318                          */
 319                         /*
 320                          * l is the step size for the current loop.  Could also be called delta-x.
 321                          *
 322                          * to account for accumulated errors in the x position during the loop
 323                          * during the last iteration (n== DIVISIONS) we recompute l to be
 324                          * whatever is left.
 325                          */
 326                         double l = (n==DIVISIONS) ? length -x : step;
 327
 328                         // Further to prevent round off error from the previous statement,
 329                         // we clamp r2 to length at the last iteration.
 330                         r2 = getRadius((n==DIVISIONS) ? length : x + l);
 331
 332                         final double hyp = MathUtil.hypot(r2 - r1, l);
 333
 334                         // Volume differential elements
 335                         final double dV;
 336                         final double dFullV;
 337
 338                         dFullV = pi3 * l * (r1 * r1 + r1 * r2 + r2 * r2);
 339
 340                         if ( filled ) {
 341                                 dV = dFullV;
 342                         } else {
 343                                 // hollow
 344                                 // Thickness is normal to the surface of the component
 345                                 // here we use simple trig to project the Thickness
 346                                 // on to the y dimension (radius).
 347                                 double height = thickness * hyp / l;
 348                                 if (r1 < height || r2 < height) {
 349                                         // Filled portion of piece
 350                                         dV = dFullV;
 351                                 } else {
 352                                         // Hollow portion of piece
 353                                         dV = MathUtil.max(Math.PI* l * height * (r1 + r2 - height), 0);
 354                                 }
 355                         }
 356
 357                         // Add to the volume-related components
 358                         volume += dV;
 359                         fullVolume += dFullV;
 360                         cgx += (x + l / 2) * dV;
 361
 362                         // Wetted area ( * PI at the end)
 363                         wetArea += hyp * (r1 + r2);
 364
 365                         // Planform area & center
 366                         final double p = l * (r1 + r2);
 367                         planArea += p;
 368                         planCenter += (x + l / 2) * p;
 369
 370                         // Update for next iteration
 371                         r1 = r2;
 372                         x += l;
 373                 }
 374
 375                 wetArea *= Math.PI;
 376
 377                 if (planArea > 0)
 378                         planCenter /= planArea;
 379
 380                 if (volume < 0.0000000001) { // 0.1 mm^3
 381                         volume = 0;
 382                         cg = new Coordinate(length / 2, 0, 0, 0);
 383                 } else {
 384                         // the mass of this shape is the material density * volume.
 385                         // it cannot come from super.getComponentMass() since that
 386                         // includes the shoulders
 387                         cg = new Coordinate(cgx / volume, 0, 0, getMaterial().getDensity() * volume );
 388                 }
 389         }
 390
 391
 392         /**
 393          * Integrate the longitudinal and rotational inertia based on component volume.
 394          * This method may be used only if the total volume is zero.
 395          */
 396         private void integrateInertiaVolume() {
 397                 double x, r1, r2;
 398
 399                 final double l = length / DIVISIONS;
 400                 final double pil = Math.PI * l; // PI * l
 401                 final double pil3 = Math.PI * l / 3; // PI * l/3
 402
 403                 r1 = getRadius(0);
 404                 x = 0;
 405                 longitudinalInertia = 0;
 406                 rotationalInertia = 0;
 407
 408                 double vol = 0;
 409
 410                 for (int n = 1; n <= DIVISIONS; n++) {
 411                         /*
 412                          * r1 and r2 are the two radii, outer is their average
 413                          * x is the position of r1
 414                          * hyp is the length of the hypotenuse from r1 to r2
 415                          * height if the y-axis height of the component if not filled
 416                          */
 417                         r2 = getRadius(x + l);
 418                         final double outer = (r1 + r2) / 2;
 419
 420
 421                         // Volume differential elements
 422                         final double inner;
 423                         final double dV;
 424
 425                         if (filled || r1 < thickness || r2 < thickness) {
 426                                 inner = 0;
 427                                 dV = pil3 * (r1 * r1 + r1 * r2 + r2 * r2);
 428                         } else {
 429                                 final double hyp = MathUtil.hypot(r2 - r1, l);
 430                                 final double height = thickness * hyp / l;
 431                                 dV = pil * height * (r1 + r2 - height);
 432                                 inner = Math.max(outer - height, 0);
 433                         }
 434
 435                         rotationalInertia += dV * (pow2(outer) + pow2(inner)) / 2;
 436                         longitudinalInertia += dV * ((3 * (pow2(outer) + pow2(inner)) + pow2(l)) / 12
 437                                         + pow2(x + l / 2));
 438
 439                         vol += dV;
 440
 441                         // Update for next iteration
 442                         r1 = r2;
 443                         x += l;
 444                 }
 445
 446                 if (MathUtil.equals(vol, 0)) {
 447                         integrateInertiaSurface();
 448                         return;
 449                 }
 450
 451                 rotationalInertia /= vol;
 452                 longitudinalInertia /= vol;
 453
 454                 // Shift longitudinal inertia to CG
 455                 longitudinalInertia = Math.max(longitudinalInertia - pow2(getComponentCG().x), 0);
 456         }
 457
 458
 459
 460         /**
 461          * Integrate the longitudinal and rotational inertia based on component surface area.
 462          * This method may be used only if the total volume is zero.
 463          */
 464         private void integrateInertiaSurface() {
 465                 double x, r1, r2;
 466
 467                 final double l = length / DIVISIONS;
 468
 469                 r1 = getRadius(0);
 470                 //System.out.println(r1);
 471                 x = 0;
 472
 473                 longitudinalInertia = 0;
 474                 rotationalInertia = 0;
 475
 476                 double surface = 0;
 477
 478                 for (int n = 1; n <= DIVISIONS; n++) {
 479                         /*
 480                          * r1 and r2 are the two radii, outer is their average
 481                          * x is the position of r1
 482                          * hyp is the length of the hypotenuse from r1 to r2
 483                          * height if the y-axis height of the component if not filled
 484                          */
 485                         r2 = getRadius(x + l);
 486                         final double hyp = MathUtil.hypot(r2 - r1, l);
 487                         final double outer = (r1 + r2) / 2;
 488
 489                         final double dS = hyp * (r1 + r2) * Math.PI;
 490
 491                         rotationalInertia += dS * pow2(outer);
 492                         longitudinalInertia += dS * ((6 * pow2(outer) + pow2(l)) / 12 + pow2(x + l / 2));
 493
 494                         surface += dS;
 495
 496                         // Update for next iteration
 497                         r1 = r2;
 498                         x += l;
 499                 }
 500
 501                 if (MathUtil.equals(surface, 0)) {
 502                         longitudinalInertia = 0;
 503                         rotationalInertia = 0;
 504                         return;
 505                 }
 506
 507                 longitudinalInertia /= surface;
 508                 rotationalInertia /= surface;
 509
 510                 // Shift longitudinal inertia to CG
 511                 longitudinalInertia = Math.max(longitudinalInertia - pow2(getComponentCG().x), 0);
 512         }
 513
 514
 515
 516
 517         /**
 518          * Invalidates the cached volume and CG information.
 519          */
 520         @Override
 521         protected void componentChanged(ComponentChangeEvent e) {
 522                 super.componentChanged(e);
 523                 if (!e.isOtherChange()) {
 524                         wetArea = -1;
 525                         planArea = -1;
 526                         planCenter = -1;
 527                         volume = -1;
 528                         fullVolume = -1;
 529                         longitudinalInertia = -1;
 530                         rotationalInertia = -1;
 531                         cg = null;
 532                 }
 533         }
 534
 535
 536
 537         ///////////   Auto radius helper methods
 538
 539
 540         /**
 541          * Returns the automatic radius for this component towards the
 542          * front of the rocket.  The automatics will not search towards the
 543          * rear of the rocket for a suitable radius.  A positive return value
 544          * indicates a preferred radius, a negative value indicates that a
 545          * match was not found.
 546          */
 547         protected abstract double getFrontAutoRadius();
 548
 549         /**
 550          * Returns the automatic radius for this component towards the
 551          * end of the rocket.  The automatics will not search towards the
 552          * front of the rocket for a suitable radius.  A positive return value
 553          * indicates a preferred radius, a negative value indicates that a
 554          * match was not found.
 555          */
 556         protected abstract double getRearAutoRadius();
 557
 558
 559
 560         /**
 561          * Return the previous symmetric component, or null if none exists.
 562          * NOTE: This method currently assumes that there are no external
 563          * "pods".
 564          *
 565          * @return      the previous SymmetricComponent, or null.
 566          */
 567         protected final SymmetricComponent getPreviousSymmetricComponent() {
 568                 RocketComponent c;
 569                 for (c = this.getPreviousComponent(); c != null; c = c.getPreviousComponent()) {
 570                         if (c instanceof SymmetricComponent) {
 571                                 return (SymmetricComponent) c;
 572                         }
 573                         if (!(c instanceof Stage) &&
 574                                         (c.relativePosition == RocketComponent.Position.AFTER))
 575                                 return null; // Bad component type as "parent"
 576                 }
 577                 return null;
 578         }
 579
 580         /**
 581          * Return the next symmetric component, or null if none exists.
 582          * NOTE: This method currently assumes that there are no external
 583          * "pods".
 584          *
 585          * @return      the next SymmetricComponent, or null.
 586          */
 587         protected final SymmetricComponent getNextSymmetricComponent() {
 588                 RocketComponent c;
 589                 for (c = this.getNextComponent(); c != null; c = c.getNextComponent()) {
 590                         if (c instanceof SymmetricComponent) {
 591                                 return (SymmetricComponent) c;
 592                         }
 593                         if (!(c instanceof Stage) &&
 594                                         (c.relativePosition == RocketComponent.Position.AFTER))
 595                                 return null; // Bad component type as "parent"
 596                 }
 597                 return null;
 598         }
 599
 600 }