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