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