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 = MathUtil.max(pil*height*(r1+r2-height), 0);
 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.0000000001) {  // 0.1 mm^3
 338                         volume = 0;
 339                         cg = new Coordinate(length/2, 0, 0, 0);
 340                 } else {
 341                         // getComponentMass is safe now
 342                         // Use super.getComponentMass() to ensure only the transition shape mass
 343                         // is used, not the shoulders
 344                         cg = new Coordinate(cgx/volume,0,0,super.getComponentMass());
 345                 }
 346         }
 347
 348
 349         /**
 350          * Integrate the longitudal and rotational inertia based on component volume.
 351          * This method may be used only if the total volume is zero.
 352          */
 353         private void integrateInertiaVolume() {
 354                 double x, r1, r2;
 355
 356                 final double l = length/DIVISIONS;
 357                 final double pil  = Math.PI*l;    // PI * l
 358                 final double pil3 = Math.PI*l/3;  // PI * l/3
 359
 360                 r1 = getRadius(0);
 361                 x = 0;
 362                 longitudalInertia = 0;
 363                 rotationalInertia = 0;
 364
 365                 double volume = 0;
 366
 367                 for (int n=1; n<=DIVISIONS; n++) {
 368                         /*
 369                          * r1 and r2 are the two radii, outer is their average
 370                          * x is the position of r1
 371                          * hyp is the length of the hypotenuse from r1 to r2
 372                          * height if the y-axis height of the component if not filled
 373                          */
 374                         r2 = getRadius(x+l);
 375                         final double outer = (r1 + r2)/2;
 376
 377
 378                         // Volume differential elements
 379                         final double inner;
 380                         final double dV;
 381
 382                         if (filled || r1<thickness || r2<thickness) {
 383                                 inner = 0;
 384                                 dV = pil3*(r1*r1 + r1*r2 + r2*r2);
 385                         } else {
 386                                 final double hyp = MathUtil.hypot(r2-r1, l);
 387                                 final double height = thickness*hyp/l;
 388                                 dV = pil*height*(r1+r2-height);
 389                                 inner = Math.max(outer-height, 0);
 390                         }
 391
 392                         rotationalInertia += dV * (pow2(outer) + pow2(inner))/2;
 393                         longitudalInertia += dV * ((3 * (pow2(outer) + pow2(inner)) + pow2(l))/12
 394                                         + pow2(x+l/2));
 395
 396                         volume += dV;
 397
 398                         // Update for next iteration
 399                         r1 = r2;
 400                         x += l;
 401                 }
 402
 403                 if (MathUtil.equals(volume,0)) {
 404                         integrateInertiaSurface();
 405                         return;
 406                 }
 407
 408                 rotationalInertia /= volume;
 409                 longitudalInertia /= volume;
 410
 411                 // Shift longitudal inertia to CG
 412                 longitudalInertia = Math.max(longitudalInertia - pow2(getComponentCG().x), 0);
 413         }
 414
 415
 416
 417         /**
 418          * Integrate the longitudal and rotational inertia based on component surface area.
 419          * This method may be used only if the total volume is zero.
 420          */
 421         private void integrateInertiaSurface() {
 422                 double x, r1, r2;
 423
 424                 final double l = length/DIVISIONS;
 425
 426                 r1 = getRadius(0);
 427                 x = 0;
 428                 longitudalInertia = 0;
 429                 rotationalInertia = 0;
 430
 431                 double surface = 0;
 432
 433                 for (int n=1; n<=DIVISIONS; n++) {
 434                         /*
 435                          * r1 and r2 are the two radii, outer is their average
 436                          * x is the position of r1
 437                          * hyp is the length of the hypotenuse from r1 to r2
 438                          * height if the y-axis height of the component if not filled
 439                          */
 440                         r2 = getRadius(x+l);
 441                         final double hyp = MathUtil.hypot(r2-r1, l);
 442                         final double outer = (r1 + r2)/2;
 443
 444                         final double dS = hyp * (r1+r2) * Math.PI;
 445
 446                         rotationalInertia += dS * pow2(outer);
 447                         longitudalInertia += dS * ((6 * pow2(outer) + pow2(l))/12 + pow2(x+l/2));
 448
 449                         surface += dS;
 450
 451                         // Update for next iteration
 452                         r1 = r2;
 453                         x += l;
 454                 }
 455
 456                 if (MathUtil.equals(surface, 0)) {
 457                         longitudalInertia = 0;
 458                         rotationalInertia = 0;
 459                         return;
 460                 }
 461
 462                 longitudalInertia /= surface;
 463                 rotationalInertia /= surface;
 464
 465                 // Shift longitudal inertia to CG
 466                 longitudalInertia = Math.max(longitudalInertia - pow2(getComponentCG().x), 0);
 467         }
 468
 469
 470
 471
 472         /**
 473          * Invalidates the cached volume and CG information.
 474          */
 475         @Override
 476         protected void componentChanged(ComponentChangeEvent e) {
 477                 super.componentChanged(e);
 478                 if (!e.isOtherChange()) {
 479                         wetArea = -1;
 480                         planArea = -1;
 481                         planCenter = -1;
 482                         volume = -1;
 483                         fullVolume = -1;
 484                         longitudalInertia = -1;
 485                         rotationalInertia = -1;
 486                         cg = null;
 487                 }
 488         }
 489
 490
 491
 492         ///////////   Auto radius helper methods
 493
 494
 495         /**
 496          * Returns the automatic radius for this component towards the
 497          * front of the rocket.  The automatics will not search towards the
 498          * rear of the rocket for a suitable radius.  A positive return value
 499          * indicates a preferred radius, a negative value indicates that a
 500          * match was not found.
 501          */
 502         protected abstract double getFrontAutoRadius();
 503
 504         /**
 505          * Returns the automatic radius for this component towards the
 506          * end of the rocket.  The automatics will not search towards the
 507          * front of the rocket for a suitable radius.  A positive return value
 508          * indicates a preferred radius, a negative value indicates that a
 509          * match was not found.
 510          */
 511         protected abstract double getRearAutoRadius();
 512
 513
 514
 515         /**
 516          * Return the previous symmetric component, or null if none exists.
 517          * NOTE: This method currently assumes that there are no external
 518          * "pods".
 519          *
 520          * @return      the previous SymmetricComponent, or null.
 521          */
 522         protected final SymmetricComponent getPreviousSymmetricComponent() {
 523                 RocketComponent c;
 524                 for (c = this.getPreviousComponent(); c != null; c = c.getPreviousComponent()) {
 525                         if (c instanceof SymmetricComponent) {
 526                                 return (SymmetricComponent)c;
 527                         }
 528                         if (!(c instanceof Stage) &&
 529                                  (c.relativePosition == RocketComponent.Position.AFTER))
 530                                 return null;   // Bad component type as "parent"
 531                 }
 532                 return null;
 533         }
 534
 535         /**
 536          * Return the next symmetric component, or null if none exists.
 537          * NOTE: This method currently assumes that there are no external
 538          * "pods".
 539          *
 540          * @return      the next SymmetricComponent, or null.
 541          */
 542         protected final SymmetricComponent getNextSymmetricComponent() {
 543                 RocketComponent c;
 544                 for (c = this.getNextComponent(); c != null; c = c.getNextComponent()) {
 545                         if (c instanceof SymmetricComponent) {
 546                                 return (SymmetricComponent)c;
 547                         }
 548                         if (!(c instanceof Stage) &&
 549                                  (c.relativePosition == RocketComponent.Position.AFTER))
 550                                 return null;   // Bad component type as "parent"
 551                 }
 552                 return null;
 553         }
 554
 555 }