git.gag.com Git - debian/openrocket/blob - src/net/sf/openrocket/simulation/RK4Simulator.java

   1 package net.sf.openrocket.simulation;
   2
   3 import java.util.Collection;
   4
   5 import net.sf.openrocket.aerodynamics.AerodynamicCalculator;
   6 import net.sf.openrocket.aerodynamics.AerodynamicForces;
   7 import net.sf.openrocket.aerodynamics.AtmosphericConditions;
   8 import net.sf.openrocket.aerodynamics.FlightConditions;
   9 import net.sf.openrocket.aerodynamics.GravityModel;
  10 import net.sf.openrocket.aerodynamics.Warning;
  11 import net.sf.openrocket.aerodynamics.WarningSet;
  12 import net.sf.openrocket.aerodynamics.WindSimulator;
  13 import net.sf.openrocket.rocketcomponent.Configuration;
  14 import net.sf.openrocket.rocketcomponent.LaunchLug;
  15 import net.sf.openrocket.rocketcomponent.RocketComponent;
  16 import net.sf.openrocket.simulation.exception.SimulationException;
  17 import net.sf.openrocket.util.Coordinate;
  18 import net.sf.openrocket.util.MathUtil;
  19 import net.sf.openrocket.util.Quaternion;
  20 import net.sf.openrocket.util.Rotation2D;
  21
  22
  23 public class RK4Simulator extends FlightSimulator {
  24
  25         /**
  26          * A recommended reasonably accurate time step.
  27          */
  28         public static final double RECOMMENDED_TIME_STEP = 0.05;
  29
  30         /**
  31          * A recommended maximum angle step value.
  32          */
  33         public static final double RECOMMENDED_ANGLE_STEP = 3*Math.PI/180;
  34
  35         /**
  36          * Maximum roll step allowed.  This is selected as an uneven division of the full
  37          * circle so that the simulation will sample the most wind directions
  38          */
  39         private static final double MAX_ROLL_STEP_ANGLE = 28.32 * Math.PI/180;
  40 //      private static final double MAX_ROLL_STEP_ANGLE = 8.32 * Math.PI/180;
  41
  42         private static final double MAX_ROLL_RATE_CHANGE = 2 * Math.PI/180;
  43         private static final double MAX_PITCH_CHANGE = 2 * Math.PI/180;
  44
  45
  46         private static final boolean DEBUG = false;
  47
  48
  49         /* Single instance so it doesn't have to be created each semi-step. */
  50         private final FlightConditions flightConditions = new FlightConditions(null);
  51
  52
  53         private Coordinate linearAcceleration;
  54         private Coordinate angularAcceleration;
  55
  56         // set by calculateFlightConditions and calculateAcceleration:
  57         private double timestep;
  58         private double oldTimestep;
  59         private AerodynamicForces forces;
  60         private double windSpeed;
  61         private double thrustForce, dragForce;
  62         private double lateralPitchRate = 0;
  63
  64         private double rollAcceleration = 0;
  65         private double lateralPitchAcceleration = 0;
  66
  67         private double maxVelocityZ = 0;
  68         private double startWarningTime = -1;
  69
  70         private Rotation2D thetaRotation;
  71
  72
  73         public RK4Simulator() {
  74                 super();
  75         }
  76
  77
  78         public RK4Simulator(AerodynamicCalculator calculator) {
  79                 super(calculator);
  80         }
  81
  82
  83
  84
  85
  86         @Override
  87         protected RK4SimulationStatus initializeSimulation(Configuration configuration,
  88                         SimulationConditions simulation) {
  89
  90                 RK4SimulationStatus status = new RK4SimulationStatus();
  91
  92                 status.startConditions = simulation;
  93
  94                 status.configuration = configuration;
  95                 // TODO: LOW: Branch names
  96                 status.flightData = new FlightDataBranch("Main", FlightDataBranch.TYPE_TIME);
  97                 status.launchRod = true;
  98                 status.time = 0.0;
  99                 status.simulationStartTime = System.nanoTime();
 100
 101                 status.launchRodDirection = new Coordinate(
 102                                 Math.sin(simulation.getLaunchRodAngle()) *
 103                                 Math.cos(simulation.getLaunchRodDirection()),
 104                                 Math.sin(simulation.getLaunchRodAngle()) *
 105                                 Math.sin(simulation.getLaunchRodDirection()),
 106                                 Math.cos(simulation.getLaunchRodAngle())
 107                 );
 108                 status.launchRodLength = simulation.getLaunchRodLength();
 109
 110                 // Take into account launch lug positions
 111                 double lugPosition = Double.NaN;
 112                 for (RocketComponent c: configuration) {
 113                         if (c instanceof LaunchLug) {
 114                                 double pos = c.toAbsolute(new Coordinate(c.getLength()))[0].x;
 115                                 if (Double.isNaN(lugPosition) || pos > lugPosition) {
 116                                         lugPosition = pos;
 117                                 }
 118                         }
 119                 }
 120                 if (!Double.isNaN(lugPosition)) {
 121                         double maxX = 0;
 122                         for (Coordinate c: configuration.getBounds()) {
 123                                 if (c.x > maxX)
 124                                         maxX = c.x;
 125                         }
 126                         if (maxX >= lugPosition) {
 127                                 status.launchRodLength = Math.max(0,
 128                                                 status.launchRodLength - (maxX - lugPosition));
 129                         }
 130                 }
 131
 132
 133                 Quaternion o = new Quaternion();
 134                 o.multiplyLeft(Quaternion.rotation(
 135                                 new Coordinate(0, simulation.getLaunchRodAngle(), 0)));
 136                 o.multiplyLeft(Quaternion.rotation(
 137                                 new Coordinate(0, 0, simulation.getLaunchRodDirection())));
 138                 status.orientation = o;
 139                 status.position = Coordinate.NUL;
 140                 status.velocity = Coordinate.NUL;
 141                 status.rotation = Coordinate.NUL;
 142
 143                 /*
 144                  * Force a very small deviation to the wind speed to avoid insanely
 145                  * perfect conditions (rocket dropping at exactly 180 deg AOA).
 146                  */
 147                 status.windSimulator = new WindSimulator();
 148                 status.windSimulator.setAverage(simulation.getWindSpeedAverage());
 149                 status.windSimulator.setStandardDeviation(
 150                                 Math.max(simulation.getWindSpeedDeviation(), 0.005));
 151 //              status.windSimulator.reset();
 152
 153                 status.gravityModel = new GravityModel(simulation.getLaunchLatitude());
 154
 155                 rollAcceleration = 0;
 156                 lateralPitchAcceleration = 0;
 157                 oldTimestep = -1;
 158                 maxVelocityZ = 0;
 159                 startWarningTime = -1;
 160
 161                 return status;
 162         }
 163
 164
 165
 166         @Override
 167         protected Collection<FlightEvent> step(SimulationConditions simulation,
 168                         SimulationStatus simulationStatus) throws SimulationException {
 169
 170                 RK4SimulationStatus status = (RK4SimulationStatus)simulationStatus;
 171
 172                 ////////  Perform RK4 integration:  ////////
 173
 174                 Coordinate k1a, k1v, k1ra, k1rv; // Acceleration, velocity, rot.acc, rot.vel
 175                 Coordinate k2a, k2v, k2ra, k2rv;
 176                 Coordinate k3a, k3v, k3ra, k3rv;
 177                 Coordinate k4a, k4v, k4ra, k4rv;
 178                 RK4SimulationStatus status2;
 179
 180
 181                 // Calculate time step and store data after first call to calculateFlightConditions
 182                 calculateFlightConditions(status);
 183
 184
 185                 /*
 186                  * Select the time step to use.  It is the minimum of the following:
 187                  *  1. the user-specified time step
 188                  *  2. the maximum pitch step angle limit
 189                  *  3. the maximum roll step angle limit
 190                  *  4. the maximum roll rate change limit (using previous step acceleration)
 191                  *  5. the maximum pitch change limit (using previous step acceleration)
 192                  *
 193                  * The last two are required since near the steady-state roll rate the roll rate
 194                  * may oscillate significantly even between the sub-steps of the RK4 integration.
 195                  *
 196                  * Additionally a low-pass filter is applied to the time step selectively
 197                  * if the new time step is longer than the previous time step.
 198                  */
 199                 double dt1 = simulation.getTimeStep();
 200                 double dt2 = simulation.getMaximumStepAngle() / lateralPitchRate;
 201                 double dt3 = Math.abs(MAX_ROLL_STEP_ANGLE / flightConditions.getRollRate());
 202                 double dt4 = Math.abs(MAX_ROLL_RATE_CHANGE / rollAcceleration);
 203                 double dt5 = Math.abs(MAX_PITCH_CHANGE / lateralPitchAcceleration);
 204                 timestep = MathUtil.min(dt1,dt2,dt3);
 205                 timestep = MathUtil.min(timestep,dt4,dt5);
 206
 207                 if (oldTimestep > 0 && oldTimestep < timestep) {
 208                         timestep = 0.3*timestep + 0.7*oldTimestep;
 209                 }
 210
 211                 if (timestep < 0.001)
 212                         timestep = 0.001;
 213
 214                 oldTimestep = timestep;
 215                 if (DEBUG)
 216                         System.out.printf("Time step: %.3f  dt1=%.3f dt2=%.3f dt3=%.3f dt4=%.3f dt5=%.3f\n",
 217                         timestep,dt1,dt2,dt3,dt4,dt5);
 218
 219
 220
 221                 //// First position, k1 = f(t, y)
 222
 223                 //calculateFlightConditions already called
 224                 calculateAcceleration(status);
 225                 k1a = linearAcceleration;
 226                 k1ra = angularAcceleration;
 227                 k1v = status.velocity;
 228                 k1rv = status.rotation;
 229
 230
 231                 // Store the flight information
 232                 storeData(status);
 233
 234
 235
 236                 //// Second position, k2 = f(t + h/2, y + k1*h/2)
 237
 238                 status2 = status.clone();
 239                 status2.time = status.time + timestep/2;
 240                 status2.position = status.position.add(k1v.multiply(timestep/2));
 241                 status2.velocity = status.velocity.add(k1a.multiply(timestep/2));
 242                 status2.orientation.multiplyLeft(Quaternion.rotation(k1rv.multiply(timestep/2)));
 243                 status2.rotation = status.rotation.add(k1ra.multiply(timestep/2));
 244
 245                 calculateFlightConditions(status2);
 246                 calculateAcceleration(status2);
 247                 k2a = linearAcceleration;
 248                 k2ra = angularAcceleration;
 249                 k2v = status2.velocity;
 250                 k2rv = status2.rotation;
 251
 252
 253                 //// Third position, k3 = f(t + h/2, y + k2*h/2)
 254
 255                 status2.orientation = status.orientation.clone();  // All others are set explicitly
 256                 status2.position = status.position.add(k2v.multiply(timestep/2));
 257                 status2.velocity = status.velocity.add(k2a.multiply(timestep/2));
 258                 status2.orientation.multiplyLeft(Quaternion.rotation(k2rv.multiply(timestep/2)));
 259                 status2.rotation = status.rotation.add(k2ra.multiply(timestep/2));
 260
 261                 calculateFlightConditions(status2);
 262                 calculateAcceleration(status2);
 263                 k3a = linearAcceleration;
 264                 k3ra = angularAcceleration;
 265                 k3v = status2.velocity;
 266                 k3rv = status2.rotation;
 267
 268
 269
 270                 //// Fourth position, k4 = f(t + h, y + k3*h)
 271
 272                 status2.orientation = status.orientation.clone();  // All others are set explicitly
 273                 status2.time = status.time + timestep;
 274                 status2.position = status.position.add(k3v.multiply(timestep));
 275                 status2.velocity = status.velocity.add(k3a.multiply(timestep));
 276                 status2.orientation.multiplyLeft(Quaternion.rotation(k3rv.multiply(timestep)));
 277                 status2.rotation = status.rotation.add(k3ra.multiply(timestep));
 278
 279                 calculateFlightConditions(status2);
 280                 calculateAcceleration(status2);
 281                 k4a = linearAcceleration;
 282                 k4ra = angularAcceleration;
 283                 k4v = status2.velocity;
 284                 k4rv = status2.rotation;
 285
 286
 287
 288                 //// Sum all together,  y(n+1) = y(n) + h*(k1 + 2*k2 + 2*k3 + k4)/6
 289
 290                 Coordinate deltaV, deltaP, deltaR, deltaO;
 291                 deltaV = k2a.add(k3a).multiply(2).add(k1a).add(k4a).multiply(timestep/6);
 292                 deltaP = k2v.add(k3v).multiply(2).add(k1v).add(k4v).multiply(timestep/6);
 293                 deltaR = k2ra.add(k3ra).multiply(2).add(k1ra).add(k4ra).multiply(timestep/6);
 294                 deltaO = k2rv.add(k3rv).multiply(2).add(k1rv).add(k4rv).multiply(timestep/6);
 295
 296                 if (DEBUG)
 297                         System.out.println("Rot.Acc: "+deltaR+"  k1:"+k1ra+"  k2:"+k2ra+"  k3:"+k3ra+
 298                                         "  k4:"+k4ra);
 299
 300                 status.velocity = status.velocity.add(deltaV);
 301                 status.position = status.position.add(deltaP);
 302                 status.rotation = status.rotation.add(deltaR);
 303                 status.orientation.multiplyLeft(Quaternion.rotation(deltaO));
 304
 305
 306                 status.orientation.normalizeIfNecessary();
 307
 308                 status.time = status.time + timestep;
 309
 310
 311                 return null;
 312         }
 313
 314
 315
 316         /**
 317          * Calculate the linear and angular acceleration at the given status.  The results
 318          * are stored in the fields {@link #linearAcceleration} and {@link #angularAcceleration}.
 319          *
 320          * @param status   the status of the rocket.
 321          * @throws SimulationException
 322          */
 323         private void calculateAcceleration(RK4SimulationStatus status) throws SimulationException {
 324
 325                 /**
 326                  * Check whether to store warnings or not.  Warnings are ignored when on the
 327                  * launch rod or 0.25 seconds after departure, and when the velocity has dropped
 328                  * below 20% of the max. velocity.
 329                  */
 330                 WarningSet warnings = status.warnings;
 331                 maxVelocityZ = MathUtil.max(maxVelocityZ, status.velocity.z);
 332                 if (status.launchRod) {
 333                         warnings = null;
 334                 } else {
 335                         if (status.velocity.z < 0.2 * maxVelocityZ)
 336                                 warnings = null;
 337                         if (startWarningTime < 0)
 338                                 startWarningTime = status.time + 0.25;
 339                 }
 340                 if (status.time < startWarningTime)
 341                         warnings = null;
 342
 343
 344                 // Calculate aerodynamic forces  (only axial if still on launch rod)
 345                 calculator.setConfiguration(status.configuration);
 346
 347                 if (status.launchRod) {
 348                         forces = calculator.getAxialForces(status.time, flightConditions, warnings);
 349                 } else {
 350                         forces = calculator.getAerodynamicForces(status.time, flightConditions, warnings);
 351                 }
 352
 353
 354                 // Allow listeners to modify the forces
 355                 int mod = flightConditions.getModCount();
 356                 SimulationListener[] list = listeners.toArray(new SimulationListener[0]);
 357                 for (SimulationListener l: list) {
 358                         l.forceCalculation(status, flightConditions, forces);
 359                 }
 360                 if (flightConditions.getModCount() != mod) {
 361                         status.warnings.add(Warning.LISTENERS_AFFECTED);
 362                 }
 363
 364
 365                 assert(!Double.isNaN(forces.CD));
 366                 assert(!Double.isNaN(forces.CN));
 367                 assert(!Double.isNaN(forces.Caxial));
 368                 assert(!Double.isNaN(forces.Cm));
 369                 assert(!Double.isNaN(forces.Cyaw));
 370                 assert(!Double.isNaN(forces.Cside));
 371                 assert(!Double.isNaN(forces.Croll));
 372
 373
 374                 ////////  Calculate forces and accelerations  ////////
 375
 376                 double dynP = (0.5 * flightConditions.getAtmosphericConditions().getDensity() *
 377                                 MathUtil.pow2(flightConditions.getVelocity()));
 378                 double refArea = flightConditions.getRefArea();
 379                 double refLength = flightConditions.getRefLength();
 380
 381
 382                 // Linear forces
 383                 thrustForce = calculateThrust(status, timestep);
 384                 dragForce = forces.Caxial * dynP * refArea;
 385                 double fN = forces.CN * dynP * refArea;
 386                 double fSide = forces.Cside * dynP * refArea;
 387
 388 //              double sin = Math.sin(flightConditions.getTheta());
 389 //              double cos = Math.cos(flightConditions.getTheta());
 390
 391 //              double forceX = - fN * cos - fSide * sin;
 392 //              double forceY = - fN * sin - fSide * cos;
 393                 double forceZ = thrustForce - dragForce;
 394
 395
 396 //              linearAcceleration = new Coordinate(forceX / forces.cg.weight,
 397 //                              forceY / forces.cg.weight, forceZ / forces.cg.weight);
 398                 linearAcceleration = new Coordinate(-fN / forces.cg.weight, -fSide / forces.cg.weight,
 399                                 forceZ / forces.cg.weight);
 400
 401                 linearAcceleration = thetaRotation.rotateZ(linearAcceleration);
 402                 linearAcceleration = status.orientation.rotate(linearAcceleration);
 403
 404                 linearAcceleration = linearAcceleration.sub(0, 0, status.gravityModel.getGravity());
 405
 406
 407                 // If still on launch rod, project acceleration onto launch rod direction and
 408                 // set angular acceleration to zero.
 409                 if (status.launchRod) {
 410                         linearAcceleration = status.launchRodDirection.multiply(
 411                                         linearAcceleration.dot(status.launchRodDirection));
 412                         angularAcceleration = Coordinate.NUL;
 413                         rollAcceleration = 0;
 414                         lateralPitchAcceleration = 0;
 415                         return;
 416                 }
 417
 418
 419                 // Convert momenta
 420                 double Cm = forces.Cm - forces.CN * forces.cg.x / refLength;
 421                 double Cyaw = forces.Cyaw - forces.Cside * forces.cg.x / refLength;
 422
 423 //              double momX = (-Cm * sin - Cyaw * cos) * dynP * refArea * refLength;
 424 //              double momY = ( Cm * cos - Cyaw * sin) * dynP * refArea * refLength;
 425                 double momX = -Cyaw * dynP * refArea * refLength;
 426                 double momY = Cm * dynP * refArea * refLength;
 427
 428                 double momZ = forces.Croll * dynP * refArea * refLength;
 429                 if (DEBUG)
 430                         System.out.printf("Croll:  %.3f  dynP=%.3f  momZ=%.3f\n",forces.Croll,dynP,momZ);
 431
 432                 assert(!Double.isNaN(momX));
 433                 assert(!Double.isNaN(momY));
 434                 assert(!Double.isNaN(momZ));
 435                 assert(!Double.isNaN(forces.longitudalInertia));
 436                 assert(!Double.isNaN(forces.rotationalInertia));
 437
 438                 angularAcceleration = new Coordinate(momX / forces.longitudalInertia,
 439                                 momY / forces.longitudalInertia, momZ / forces.rotationalInertia);
 440
 441                 rollAcceleration = angularAcceleration.z;
 442                 // TODO: LOW: This should be hypot, but does it matter?
 443                 lateralPitchAcceleration = MathUtil.max(Math.abs(angularAcceleration.x),
 444                                 Math.abs(angularAcceleration.y));
 445
 446                 if (DEBUG)
 447                         System.out.println("rot.inertia = "+forces.rotationalInertia);
 448
 449                 angularAcceleration = thetaRotation.rotateZ(angularAcceleration);
 450
 451                 angularAcceleration = status.orientation.rotate(angularAcceleration);
 452         }
 453
 454
 455         /**
 456          * Calculate the flight conditions for the current rocket status.  The conditions
 457          * are stored in the field {@link #flightConditions}.  Additional information that
 458          * is calculated and will be stored in the flight data is also computed into the
 459          * suitable fields.
 460          * @throws SimulationException
 461          */
 462         private void calculateFlightConditions(RK4SimulationStatus status) throws SimulationException {
 463
 464                 flightConditions.setReference(status.configuration);
 465
 466
 467                 //// Atmospheric conditions
 468                 AtmosphericConditions atmosphere = status.startConditions.getAtmosphericModel().
 469                         getConditions(status.position.z + status.startConditions.getLaunchAltitude());
 470                 flightConditions.setAtmosphericConditions(atmosphere);
 471
 472
 473                 //// Local wind speed and direction
 474                 windSpeed = status.windSimulator.getWindSpeed(status.time);
 475                 Coordinate airSpeed = status.velocity.add(windSpeed, 0, 0);
 476                 airSpeed = status.orientation.invRotate(airSpeed);
 477
 478
 479         // Lateral direction:
 480         double len = MathUtil.hypot(airSpeed.x, airSpeed.y);
 481         if (len > 0.0001) {
 482             thetaRotation = new Rotation2D(airSpeed.y/len, airSpeed.x/len);
 483             flightConditions.setTheta(Math.atan2(airSpeed.y, airSpeed.x));
 484         } else {
 485             thetaRotation = Rotation2D.ID;
 486             flightConditions.setTheta(0);
 487         }
 488
 489                 double velocity = airSpeed.length();
 490         flightConditions.setVelocity(velocity);
 491         if (velocity > 0.01) {
 492             // aoa must be calculated from the monotonous cosine
 493             // sine can be calculated by a simple division
 494             flightConditions.setAOA(Math.acos(airSpeed.z / velocity), len / velocity);
 495         } else {
 496             flightConditions.setAOA(0);
 497         }
 498
 499
 500                 // Roll, pitch and yaw rate
 501                 Coordinate rot = status.orientation.invRotate(status.rotation);
 502                 rot = thetaRotation.invRotateZ(rot);
 503
 504                 flightConditions.setRollRate(rot.z);
 505                 if (len < 0.001) {
 506                         flightConditions.setPitchRate(0);
 507                         flightConditions.setYawRate(0);
 508                         lateralPitchRate = 0;
 509                 } else {
 510                         flightConditions.setPitchRate(rot.y);
 511                         flightConditions.setYawRate(rot.x);
 512                         // TODO: LOW: set this as power of two?
 513                         lateralPitchRate = MathUtil.hypot(rot.x, rot.y);
 514                 }
 515
 516
 517                 // Allow listeners to modify the conditions
 518                 int mod = flightConditions.getModCount();
 519                 SimulationListener[] list = listeners.toArray(new SimulationListener[0]);
 520                 for (SimulationListener l: list) {
 521                         l.flightConditions(status, flightConditions);
 522                 }
 523                 if (mod != flightConditions.getModCount()) {
 524                         // Re-calculate cached values
 525                         thetaRotation = new Rotation2D(flightConditions.getTheta());
 526                         lateralPitchRate = MathUtil.hypot(flightConditions.getPitchRate(),
 527                                         flightConditions.getYawRate());
 528                         status.warnings.add(Warning.LISTENERS_AFFECTED);
 529                 }
 530
 531         }
 532
 533
 534
 535         private void storeData(RK4SimulationStatus status) {
 536                 FlightDataBranch data = status.flightData;
 537                 boolean extra = status.startConditions.getCalculateExtras();
 538
 539                 data.addPoint();
 540                 data.setValue(FlightDataBranch.TYPE_TIME, status.time);
 541                 data.setValue(FlightDataBranch.TYPE_ALTITUDE, status.position.z);
 542                 data.setValue(FlightDataBranch.TYPE_POSITION_X, status.position.x);
 543                 data.setValue(FlightDataBranch.TYPE_POSITION_Y, status.position.y);
 544
 545                 if (extra) {
 546                         data.setValue(FlightDataBranch.TYPE_POSITION_XY,
 547                                         MathUtil.hypot(status.position.x, status.position.y));
 548                         data.setValue(FlightDataBranch.TYPE_POSITION_DIRECTION,
 549                                         Math.atan2(status.position.y, status.position.x));
 550
 551                         data.setValue(FlightDataBranch.TYPE_VELOCITY_XY,
 552                                         MathUtil.hypot(status.velocity.x, status.velocity.y));
 553                         data.setValue(FlightDataBranch.TYPE_ACCELERATION_XY,
 554                                         MathUtil.hypot(linearAcceleration.x, linearAcceleration.y));
 555
 556                         data.setValue(FlightDataBranch.TYPE_ACCELERATION_TOTAL,linearAcceleration.length());
 557
 558                         double Re = flightConditions.getVelocity() *
 559                                         calculator.getConfiguration().getLength() /
 560                                         flightConditions.getAtmosphericConditions().getKinematicViscosity();
 561                         data.setValue(FlightDataBranch.TYPE_REYNOLDS_NUMBER, Re);
 562                 }
 563
 564                 data.setValue(FlightDataBranch.TYPE_VELOCITY_Z, status.velocity.z);
 565                 data.setValue(FlightDataBranch.TYPE_ACCELERATION_Z, linearAcceleration.z);
 566
 567                 data.setValue(FlightDataBranch.TYPE_VELOCITY_TOTAL, flightConditions.getVelocity());
 568                 data.setValue(FlightDataBranch.TYPE_MACH_NUMBER, flightConditions.getMach());
 569
 570                 if (!status.launchRod) {
 571                         data.setValue(FlightDataBranch.TYPE_CP_LOCATION, forces.cp.x);
 572                         data.setValue(FlightDataBranch.TYPE_CG_LOCATION, forces.cg.x);
 573                         data.setValue(FlightDataBranch.TYPE_STABILITY,
 574                                         (forces.cp.x - forces.cg.x) / flightConditions.getRefLength());
 575                 }
 576                 data.setValue(FlightDataBranch.TYPE_MASS, forces.cg.weight);
 577
 578                 data.setValue(FlightDataBranch.TYPE_THRUST_FORCE, thrustForce);
 579                 data.setValue(FlightDataBranch.TYPE_DRAG_FORCE, dragForce);
 580
 581                 if (!status.launchRod) {
 582                         data.setValue(FlightDataBranch.TYPE_PITCH_MOMENT_COEFF,
 583                                         forces.Cm - forces.CN * forces.cg.x / flightConditions.getRefLength());
 584                         data.setValue(FlightDataBranch.TYPE_YAW_MOMENT_COEFF,
 585                                         forces.Cyaw - forces.Cside * forces.cg.x / flightConditions.getRefLength());
 586                         data.setValue(FlightDataBranch.TYPE_NORMAL_FORCE_COEFF, forces.CN);
 587                         data.setValue(FlightDataBranch.TYPE_SIDE_FORCE_COEFF, forces.Cside);
 588                         data.setValue(FlightDataBranch.TYPE_ROLL_MOMENT_COEFF, forces.Croll);
 589                         data.setValue(FlightDataBranch.TYPE_ROLL_FORCING_COEFF, forces.CrollForce);
 590                         data.setValue(FlightDataBranch.TYPE_ROLL_DAMPING_COEFF, forces.CrollDamp);
 591                         data.setValue(FlightDataBranch.TYPE_PITCH_DAMPING_MOMENT_COEFF,
 592                                         forces.pitchDampingMoment);
 593                 }
 594
 595                 data.setValue(FlightDataBranch.TYPE_DRAG_COEFF, forces.CD);
 596                 data.setValue(FlightDataBranch.TYPE_AXIAL_DRAG_COEFF, forces.Caxial);
 597                 data.setValue(FlightDataBranch.TYPE_FRICTION_DRAG_COEFF, forces.frictionCD);
 598                 data.setValue(FlightDataBranch.TYPE_PRESSURE_DRAG_COEFF, forces.pressureCD);
 599                 data.setValue(FlightDataBranch.TYPE_BASE_DRAG_COEFF, forces.baseCD);
 600
 601                 data.setValue(FlightDataBranch.TYPE_REFERENCE_LENGTH, flightConditions.getRefLength());
 602                 data.setValue(FlightDataBranch.TYPE_REFERENCE_AREA, flightConditions.getRefArea());
 603
 604
 605                 data.setValue(FlightDataBranch.TYPE_PITCH_RATE, flightConditions.getPitchRate());
 606                 data.setValue(FlightDataBranch.TYPE_YAW_RATE, flightConditions.getYawRate());
 607
 608
 609
 610                 if (extra) {
 611                         Coordinate c = status.orientation.rotateZ();
 612                         double theta = Math.atan2(c.z, MathUtil.hypot(c.x, c.y));
 613                         double phi = Math.atan2(c.y, c.x);
 614                         if (phi < -(Math.PI-0.0001))
 615                                 phi = Math.PI;
 616                         data.setValue(FlightDataBranch.TYPE_ORIENTATION_THETA, theta);
 617                         data.setValue(FlightDataBranch.TYPE_ORIENTATION_PHI, phi);
 618                 }
 619
 620                 data.setValue(FlightDataBranch.TYPE_AOA, flightConditions.getAOA());
 621                 data.setValue(FlightDataBranch.TYPE_ROLL_RATE, flightConditions.getRollRate());
 622
 623                 data.setValue(FlightDataBranch.TYPE_WIND_VELOCITY, windSpeed);
 624                 data.setValue(FlightDataBranch.TYPE_AIR_TEMPERATURE,
 625                                 flightConditions.getAtmosphericConditions().temperature);
 626                 data.setValue(FlightDataBranch.TYPE_AIR_PRESSURE,
 627                                 flightConditions.getAtmosphericConditions().pressure);
 628                 data.setValue(FlightDataBranch.TYPE_SPEED_OF_SOUND,
 629                                 flightConditions.getAtmosphericConditions().getMachSpeed());
 630
 631
 632                 data.setValue(FlightDataBranch.TYPE_TIME_STEP, timestep);
 633                 data.setValue(FlightDataBranch.TYPE_COMPUTATION_TIME,
 634                                 (System.nanoTime() - status.simulationStartTime)/1000000000.0);
 635
 636
 637 //              data.setValue(FlightDataBranch.TYPE_, 0);
 638
 639         }
 640
 641
 642 }