attachment:main.cpp of NumericalSimulation

   1 // Calculation of the movement of the launch loop during launch
   2 // By Kristian Andresen
   3 // To do:
   4 // - Conservation of energy
   5 // - Planetary rotation
   6 // - Pressure as function of height
   7 // - Wind loading
   8 // - Spaceship launch
   9 
  10 #include <iostream>
  11 #include <fstream>
  12 #include <stdlib.h>
  13 #include <math.h>
  14 #include <sstream> 
  15 #include <vector>
  16 using namespace std;
  17 
  18 // Utilities
  19 namespace transient {
  20 #define swap(T,a,b) { T buff = b; b = a; a = buff; }
  21 double temp[3]; // Temporary holder of vectors
  22 #define plusMinus(i) for(int i = -1; i <= 1; i = i+2)
  23 #define for2(i) for(int i = 0; i < 2; i++)
  24 #define for3(i) for(int i = 0; i < 3; i++)
  25 #define forN(N,i) for(int i = 0; i < N; i++)
  26 inline double sqr(double x) { return x*x; }
  27 inline double dot(double a[], double b[]) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
  28 inline double sqr(double a[]) { return dot(a, a); }
  29 inline double* sub(double res[], double a[], double b[]) { for3(i) res[i] = a[i] - b[i]; return res; }
  30 inline double* sub(double a[], double b[]) { for3(i) temp[i] = a[i] - b[i]; return temp; }
  31 inline void normalize(double a[]) { double length = sqrt(sqr(a)); for3(i) a[i] /= length; } 
  32 inline double dist(double a[], double b[]) { return sqrt(sqr(sub(a, b))); }
  33 inline int wrap(int i) { return i >= 3? i-3 : i; }
  34 inline double abs(double a) { return a < 0? -a : a; } 
  35 inline double* cross(double c[], double a[], double b[]) { for3(i) c[i] = a[wrap(1+i)]*b[wrap(2+i)] - a[wrap(2+i)]*b[wrap(1+i)]; return c; }
  36 int getColor(int r, int g, int b) { return r * 65536 + g * 256 + b; }
  37 const int rotorColor[4] = {getColor(255,0,0), getColor(0,255,0), getColor(255,255,0), getColor(0,255,255)};
  38 }
  39 
  40 // Constants and globals
  41 namespace transient {
  42 // Physical constants
  43 const double earthRadius = 6378100; // [m]
  44 const double earthMassTimesG = 3.98734836e14; // [N*m^2/kg]
  45 const double earthMass = 5.9736e24; // [kg]
  46 const double pi = 3.1415926;
  47 
  48 // Track design parameters
  49 const double distance = 100e3; // [m] Distance between ground stations
  50 const double rotorDensity = 1; // [kg/m] Rotor density at apex
  51 const double apex[3] = {0.0, 0.0, earthRadius + 40e3}; // [x,y,z] 
  52 const double setupApexSpeed = 2000; // [m/s] Speed of rotor at apex during initialization
  53 const double simulationApexSpeed = 2000; // [m/s] Speed of rotor at apex during simulation
  54 const double soundInTrack = 4000; // [m/s] Speed of sound in track
  55 const double specificStrength = 1e5; // [N/(kg/m)] Maximum tensile force in track depends on density
  56 const double trackArea = 0.001; // [m^2] This influences track resistance to bending
  57 const double shearModulus = 60e9; // [Pa] Shear modulus of steel
  58 const double youngModulus = 200e9; // [Pa] Youngs modulus of steel
  59 const int tracks = 4; // Number of tracks
  60 const bool eastbound[4] = {true, false, true, false}; // Direction of each track
  61 const bool helicalTrack = true; // If true, track is initialized as a spiral
  62 const int guyWires = 5; // Number of guy wires attached to each track segment, uneven number for balance
  63 const bool guyWire = (tracks > 1) && helicalTrack; // If true, employ guy wires
  64 const double helixRadius = 1000; // [m]
  65 const double helixes = 7; // Number of twists along the track, need not be integer
  66 const double centrifugalPeriod = distance * 2 / simulationApexSpeed / helixes; // [s]
  67 const double centrifugalAcceleration = 4 * sqr(pi) * helixRadius / sqr(centrifugalPeriod); // [m/s^2]
  68 
  69 // Simulation parameters
  70 const int segments = 200; // Defines the resolution of the calculation, even number
  71 const int timesteps = 20000; // Length of simulation
  72 const int sampling = 200; // Data is sampled after this number of timesteps
  73 const double setupTick = distance / segments / setupApexSpeed; // [s] Duration of one tick during initialization
  74 // If simulationTick is too large, simulation will be numerically unstable, but not for reasons which have to do with physics
  75 const double simulationTick = setupTick / 100; // [s] A method is needed to calculate a suitable simulationTick
  76 
  77 // Globals
  78 double tick; // [s] Duration of one time step
  79 double apexSpeed; // [m/s] Speed of rotor at apex
  80 double trackDensity; // [kg/m]
  81 }
  82 
  83 namespace transient {
  84 
  85 struct Spring; // Forward declaration
  86 
  87 // Main data structure
  88 struct Segment
  89 {
  90   // Segment attributes
  91   double position[3]; // [m]
  92   double velocity[3]; // [m/s]
  93   double rotorMass; // [kg]
  94   double trackMass; // [kg] 
  95   double massRatio; // Ratio between gravitational and inertial mass, used to apply gravity during setup
  96   double direction[3]; // Normalized vector to next track segment
  97 
  98   // Mass and momentum transfer variables
  99   double designSpeed; // [m/s]
 100   double rotorSpeed; // [m/s] Actual speed
 101   double stableSpeed; // [m/s]
 102   double incomingMass; // [kg]
 103   double incomingVelocity[3]; // [m/s]
 104   double incomingPosition[3]; // [m]
 105 
 106   // Spring variables
 107   Spring* shearStress[2][2]; // Springs to model bending stresses [0 is upstream, 1 is downstram][0 is down, 1 is sideways]
 108   Spring* track; // Connection to downstream track segment
 109   Spring* wire[tracks][guyWires]; // Guy wires
 110 
 111   Segment() { 
 112     for2(i) for2(j) shearStress[i][j] = 0;
 113     track = 0;
 114     forN(tracks,d) forN(guyWires,w) wire[d][w] = 0;
 115   };
 116   double mass() { return rotorMass + trackMass; };
 117   double length() { return dist(position, (this+1)->position); };
 118   double speed() { return sqrt(sqr(velocity)); };
 119   double height() { return sqrt(sqr(position)) - earthRadius; }; // Height above ground
 120   void findDirection() { normalize(sub(direction, (this+1)->position, position)); };
 121   void applyGravity();
 122   void handleMomentum() { for3(a) position[a] += velocity[a] * tick; };
 123   void handleArrival();
 124   void activeStabilization();
 125   void handleDeparture();
 126   void connect(Spring** c, Segment* s, double* vector, double length, bool setup);
 127   void connectTrack(Segment* s, bool setup);
 128   void connectWire(Segment* s, int i, int j, int w, bool setup);
 129   void shearForces(bool setup);
 130   double* guyForce();
 131 };
 132 
 133 // Data structure for guy wires, track connections and shear stress, modelled as Kelvin-Voigt materials
 134 struct Spring // Each spring has design length, present length, and old length from previous timestep
 135 {
 136   bool bendable; // True for guy wires, false for track
 137   double designLength; // [m]
 138   double length; // [m] The length can be negative when modelling shear forces with springs
 139   double oldLength; // [m] Used to find damping
 140   double designTension; // [N]
 141   double guyTension; // [N] Only used by guy wires
 142   Segment* end[2]; // Remember that pointers can't handle being copied
 143   double dampingForce; // [N]
 144   double springConstant; // [kg/s^2]
 145   double vector[3]; // [m] Projection of vector from end[0] to end[1] onto spring direction
 146   void contract();
 147   double dampingConstant() { return 2 * sqrt(springConstant * (end[0]->mass() + end[1]->mass())) / 100; }; // [kg/s] Guess
 148   double springForce() { return - springConstant * (length - designLength) - designTension; };
 149   Spring() { guyTension = 0; };
 150 };
 151 
 152 // Data
 153 Segment track[tracks][segments+1];
 154 vector<Spring*> spring;
 155 
 156 // Member functions
 157 void Segment::applyGravity() {
 158   double down[3];
 159   for3(a) down[a] = -position[a];
 160   double g = earthMassTimesG * massRatio / sqr(down);
 161   normalize(down);
 162   for3(a) velocity[a] += down[a] * g * tick;
 163   for3(a) position[a] += down[a] * g * sqr(tick) / 2;
 164 }
 165 void Segment::handleArrival() {
 166   for3(a) velocity[a] = (incomingMass * incomingVelocity[a] + mass() * velocity[a]) / (mass() + incomingMass); 
 167   for3(a) position[a] = (incomingMass * incomingPosition[a] + mass() * position[a]) / (mass() + incomingMass);
 168   rotorMass += incomingMass;
 169 }
 170 void Segment::activeStabilization() { // This doesn't perform much better than simply setting stableSpeed = designSpeed
 171   stableSpeed = mass() * dot(velocity, direction) / rotorMass / sqr(direction); // Minimizing kinetic energy of track
 172   double minMaxSpeed[2]; // Index 0 is not necessarily the minimum
 173   for2(i) { // Alternatively, choose rotor speed such that potential plus kinetic energy of rotor is constant
 174     double g = earthMassTimesG / sqr((this-1+i)->position);
 175     minMaxSpeed[i] = sqrt(sqr((this-1+2*i)->rotorSpeed)
 176       - 2 * (this-1+i)->mass() / (this-1+i)->rotorMass * g * (this->height() - (this-1+2*i)->height()));
 177   }
 178   if ((stableSpeed > minMaxSpeed[0] && stableSpeed < minMaxSpeed[1]) || (stableSpeed < minMaxSpeed[0] && stableSpeed > minMaxSpeed[1])) {
 179 
 180   } else stableSpeed = designSpeed;
 181 }
 182 void Segment::handleDeparture() { 
 183   rotorSpeed = stableSpeed;
 184   // First we find the mass to be moved based on the speed and density
 185   (this+1)->incomingMass = rotorSpeed * rotorMass / length() * tick;
 186   rotorMass -= (this+1)->incomingMass;
 187   // Next, find the outgoing velocity
 188   normalize(sub((this+1)->incomingVelocity, (this+1)->position, position));
 189   for3(a) (this+1)->incomingVelocity[a] *= rotorSpeed;
 190   // Next, conserve momentum
 191   for3(a) velocity[a] = ((mass() + (this+1)->incomingMass) * velocity[a] - (this+1)->incomingMass * (this+1)->incomingVelocity[a]) / mass();
 192   // Lastly, adjust center of mass
 193   for3(a) (this+1)->incomingPosition[a] = (position[a] + (this+1)->position[a]) / 2;
 194   for3(a) position[a] = ((mass() + (this+1)->incomingMass) * position[a] - (this+1)->incomingMass * (this+1)->incomingPosition[a]) / mass();
 195 }
 196 void Segment::connect(Spring** c, Segment* s, double* vector, double length, bool setup) {
 197   if (setup) {
 198     spring.push_back((*c) = new Spring());
 199     (*c)->designLength = length;
 200     (*c)->oldLength = length;
 201     (*c)->end[0] = this;
 202     (*c)->end[1] = s;
 203     (*c)->bendable = false;
 204     (*c)->designTension = 0;
 205   }
 206   for3(a) (*c)->vector[a] = vector[a];
 207   (*c)->length = length;
 208 }
 209 void Segment::connectTrack(Segment* s, bool setup) {
 210   Spring** c = &track;
 211   connect(c, s, sub(s->position, position), dist(s->position, position), setup);
 212   if (setup) {
 213     (*c)->springConstant = sqr(soundInTrack) * trackDensity / (*c)->length;
 214     (*c)->designTension = 1e4; 
 215   } else (*c)->contract();
 216 }
 217 void Segment::connectWire(Segment* s, int i, int j, int w, bool setup) {
 218   Spring** c = &wire[j][w];
 219   connect(c, s, sub(s->position, position), dist(s->position, position), setup);
 220   if (setup) {
 221     (*c)->springConstant = sqr(soundInTrack) * trackDensity / (*c)->length / 10; // Artificially high, active tensioning
 222     (*c)->bendable = true;
 223     (*c)->designTension = 1; // Temporary value
 224     s->wire[i][w] = (*c); // Give opposing segment access to spring data
 225   } else (*c)->contract();
 226 }
 227 void Segment::shearForces(bool setup) {
 228   // First, find guy wire and track direction - these change with time and cannot be cached
 229   double axis[3][3]; // [0 is the track direction, 1 is down, 2 is sideways][x,y,z]
 230   for3(i) for3(j) axis[i][j] = 0;
 231   sub(axis[0], (this+1)->position, (this-1)->position);
 232   if (guyWire) { // Use guy wires to define the local "down" direction - should perhaps also use neighbour segments
 233     forN(tracks,d) forN(guyWires,w) if (wire[d][w]) for3(a) axis[1][a] += (wire[d][w]->end[0] == this? 1 : -1) * wire[d][w]->vector[a];
 234     cross(axis[2], axis[0], axis[1]);
 235   } else { // Use north to define sideways
 236     axis[2][1] = 1;
 237     cross(axis[1], axis[0], axis[2]);
 238     cross(axis[2], axis[0], axis[1]);
 239   }
 240   for3(a) normalize(axis[a]);
 241   // Next, do projections, reducing to one dimensional strains
 242   double projection[3][3]; // Indices are segment and axis 
 243   for3(s) for3(a) projection[s][a] = dot((this-1+s)->position, axis[a]) / sqr(axis[a]);
 244   // Next, find the deflection
 245   double deflection[2]; // [0 is guy wire direction, 1 is sideways]
 246   for2(i) deflection[i] = projection[1][i+1] - (projection[0][i+1] + projection[2][i+1])/2;
 247   // Move the data to the spring
 248   for2(i) for2(j) {
 249     Spring** c = &shearStress[i][j];
 250     connect(c, this+2*i-1, axis[1+j], -deflection[j], setup);
 251     if (setup) (*c)->springConstant = shearModulus * trackArea / length();
 252     else (*c)->contract();
 253   }
 254 }
 255 double* Segment::guyForce() { // Find the combined force of all guy wires
 256   for3(a) temp[a] = 0;
 257   forN(tracks,j) forN(guyWires,w) if (wire[j][w]) {
 258     normalize(wire[j][w]->vector);
 259     for3(a) temp[a] += (wire[j][w]->end[0] == this? -1 : 1) * wire[j][w]->springForce() * wire[j][w]->vector[a];
 260   }
 261   return temp;
 262 }
 263 
 264 void Spring::contract() {
 265   double springSpeed = (length - oldLength) / tick;
 266   oldLength = length;
 267   dampingForce = - dampingConstant() * springSpeed;
 268   normalize(vector); // Get direction of acceleration
 269   double force = springForce() + dampingForce;
 270   if (bendable && force > 0) force = 0; // A bendable spring can only provide a contracting (negative) force, or it buckles
 271   for2(i) for3(a) end[i]->velocity[a] += (2*i-1) * force * tick * vector[a] / end[i]->mass(); // Newtons law
 272 }
 273 
 274 // Initialization
 275 void traceTrack(Segment* t) { 
 276   // Define initial parabola - disperse track segments so the slugs pass by them at constant time intervals
 277   forN(segments+1,s) t[s].rotorMass = rotorDensity * apexSpeed * tick; // Is this balanced right?
 278   Segment* middle = &t[segments/2];
 279   for3(a) middle->position[a] = apex[a];
 280   middle->velocity[0] = apexSpeed;
 281   middle->trackMass = trackDensity * apexSpeed * tick;
 282   middle->massRatio = (middle->trackMass + middle->rotorMass) / middle->rotorMass;
 283   for (int d = -1; d <= 1; d += 2) { // 2 different directions
 284     for3(a) middle->velocity[a] *= -1; // Switch direction
 285     for (int i = 1; i <= segments/2; i++) { 
 286       Segment* s = t + segments/2 + d*i;
 287       for3(a) s->position[a] = (s-d)->position[a] + (s-d)->velocity[a] * tick;
 288       for3(a) s->velocity[a] = (s-d)->velocity[a];
 289       s->trackMass = trackDensity * s->speed() * tick;
 290       s->massRatio = s->mass() / s->rotorMass;
 291       s->applyGravity();
 292     }
 293   }
 294   forN(segments+1,s) t[s].position[1] = (double(rand())/RAND_MAX-0.5)/100; // Centimeter precision in startup
 295 }
 296 void findTrackDensity(Segment* t) {
 297   // Adjust the weight of the track to keep the rotor down
 298   trackDensity = 10 * rotorDensity; // Guess, to be improved upon
 299   for (double adjustment = trackDensity/2; adjustment > 0.01; adjustment /= 2) {
 300     traceTrack(t);
 301     trackDensity += adjustment * (t[0].height() > 0? 1 : -1);
 302   }
 303   cout << "Track density " << trackDensity << endl;
 304   cout << "Rotor speed at ground " << t[0].speed() << endl;
 305 }
 306 void applyHelix(Segment* t, int d)
 307 {
 308   double north[3] = {0,1,0};
 309   double up[segments+1][3];
 310   forN(segments+1,s) {
 311     normalize(cross(up[s], t[s].direction, north));
 312     // Equal time between segments leads to equal periods of centrifugal rotation
 313     for3(a) t[s].position[a] += north[a] * helixRadius * sin(helixes*pi*s/segments + d*2*pi/tracks);
 314     for3(a) t[s].position[a] += up[s][a] * helixRadius * cos(helixes*pi*s/segments + d*2*pi/tracks);
 315   }
 316 }
 317 void handleSprings(bool setup) {
 318   if (guyWire) forN(segments+1,i) forN(tracks,d) forN(tracks,j) forN(guyWires,w) {
 319     int opposite = i + w - (guyWires-1)/2;
 320     if (eastbound[d] != eastbound[j]) opposite = segments - opposite; 
 321     Segment* s = &track[d][i];
 322     if (opposite >= 0 && opposite <= segments && d != j) { // Stay within track, don't connect to self
 323       if ((setup && !s->wire[j][w] && !track[j][opposite].wire[d][w]) // Check that this isn't a repeat
 324         || (!setup && s->wire[j][w] && s->wire[j][w]->end[0] == s)) { // Check that wire exists
 325         s->connectWire(&track[j][opposite], d, j, w, setup);
 326       }
 327     }
 328   }
 329   // Fix guy wire tension so it balances centrifugal force
 330   if (setup && guyWire) forN(segments+1,i) forN(tracks,d) {
 331     Segment* s = &track[d][i];
 332     double* sum = s->guyForce();
 333     double factor = centrifugalAcceleration * s->rotorMass / s->length() * s->designSpeed * tick / 2 / sqrt(sqr(sum));
 334     forN(tracks,j) forN(guyWires,w) if (s->wire[j][w]) s->wire[j][w]->guyTension += factor;
 335   }
 336   if (setup && guyWire) forN(segments+1,s) forN(tracks,i) forN(tracks,j) forN(guyWires,w)
 337     if (track[i][s].wire[j][w]) track[i][s].wire[j][w]->designTension = track[i][s].wire[j][w]->guyTension;
 338   // Track connection
 339   forN(tracks,d) forN(segments,s) track[d][s].connectTrack(&track[d][s+1], setup);
 340   // Shear forces
 341   forN(tracks,d) forN(segments-1,s) track[d][s+1].shearForces(setup);
 342 }
 343 void setupLoop() {
 344   findTrackDensity(track[0]);
 345   forN(tracks,d) forN(segments+1,s) track[d][s] = track[0][s]; // Copy initial parabola
 346   forN(tracks,d) {
 347     forN(segments,s) track[d][s].findDirection();
 348     for3(a) track[d][segments].direction[a] = track[d][segments-1].direction[a]; // This is an approximation
 349     if (helicalTrack) applyHelix(track[d], d);
 350     // The westbound track must move opposite the eastbound track, so we swap
 351     if (!eastbound[d]) forN(segments/2,s) swap(Segment,track[d][s],track[d][segments-s]); 
 352     forN(segments,s) track[d][s].findDirection(); // Direction has changed now
 353     forN(segments+1,s) track[d][s].designSpeed = track[d][s].speed() * (simulationApexSpeed/setupApexSpeed);
 354     forN(segments+1,s) track[d][s].rotorSpeed = track[d][s].designSpeed; 
 355     forN(segments+1,s) for3(a) track[d][s].velocity[a] = track[d][s].designSpeed * track[d][s].direction[a] / track[d][s].massRatio;
 356     forN(segments+1,s) track[d][s].massRatio = 1; 
 357   }
 358   handleSprings(true); // Setup of springs
 359 }
 360 
 361 // Output
 362 void saveToFile(int t) {
 363   ofstream file;
 364   stringstream filename;
 365   filename << "data/t" << 1000000 + t << ".txt"; // Inelegant method for getting alphabetical order
 366   file.open(filename.str());
 367   forN(tracks,d) forN(segments-1,i)  {
 368     Segment* s = &track[d][i+1];
 369     int color = (i + int((timesteps-t) * (simulationApexSpeed/setupApexSpeed) / (setupTick/simulationTick))) % (segments/2);
 370     if (color > segments/4) color = rotorColor[d]; else color = getColor(0,0,0);
 371     if (i == 0) color = getColor(255,255,255); // Invisible jump from eastbound to westbound
 372     file << s->position[0] / 1000 << " ";
 373     file << s->position[1] << " ";
 374     file << s->height() / 1000 << " ";
 375     file << color << " ";
 376     file << log(abs(s->track->springForce())+10e-20) << " ";
 377     file << log(sqrt(sqr(s->guyForce()))+10e-20) << " ";
 378     file << log(s->rotorMass/track[d][i].length()*sqr(s->designSpeed)+10e-20) << "\n";
 379   }
 380   file.close(); 
 381 }
 382 void evaluate() {
 383   // Rotor mass
 384   double maxRotor[4];
 385   for2(i) maxRotor[i] = track[0][0].rotorMass; for2(i) maxRotor[2+i] = 0;
 386   forN(tracks,t) forN(segments,s) for2(i) if (track[t][s].rotorMass * (2*i-1) > (2*i-1) * maxRotor[i]) {
 387     maxRotor[i] = track[t][s].rotorMass;
 388     maxRotor[2+i] = s;
 389   }
 390   cout << "Rotor mass extremum "; for2(i) cout << 100*maxRotor[i]/track[0][0].rotorMass << "% ";
 391   cout << "Index "; for2(i) cout << maxRotor[2+i] << " "; cout << endl;
 392   // Rotor speed
 393   for2(i) maxRotor[i] = 1; for2(i) maxRotor[2+i] = 0;
 394   forN(tracks,t) forN(segments-1,s) for2(i) if (track[t][s+1].rotorSpeed/track[t][s+1].designSpeed * (2*i-1) > (2*i-1) * maxRotor[i]) {
 395     maxRotor[i] = track[t][s+1].rotorSpeed/track[t][s+1].designSpeed;
 396     maxRotor[2+i] = s+1;
 397   }
 398   cout << "Rotor speed extremum "; for2(i) cout << 100*maxRotor[i] << "% ";
 399   cout << "Index "; for2(i) cout << maxRotor[2+i] << " "; cout << endl;
 400 }
 401 
 402 // Simulation
 403 void simulateTrack(Segment* t) {
 404   for (int i = 1; i < segments; i++) t[i].findDirection();
 405   // Handle the injection of new slugs
 406   t[1].incomingMass = t[0].rotorMass / t[0].length() * t[0].designSpeed * tick;
 407   for3(a) t[1].incomingVelocity[a] = t[0].designSpeed * t[0].direction[a];
 408   for3(a) t[1].incomingPosition[a] = (t[0].position[a] + t[1].position[a]) / 2;
 409   // Handle the main track
 410   for (int i = 1; i < segments; i++) t[i].applyGravity();
 411   for (int i = 1; i < segments; i++) t[i].activeStabilization();
 412   for (int i = 1; i < segments; i++) t[i].handleDeparture();
 413   for (int i = 1; i < segments; i++) t[i].handleArrival();
 414   for (int i = 1; i < segments; i++) t[i].handleMomentum();
 415 }
 416 void simulateLoop() {
 417   cout << "Remaining timesteps ";
 418   forN(timesteps+1,t) {
 419     if (t % (timesteps/10) == 0) { cout << 100 - 100*t/timesteps << "%" << endl; if (t > 0) evaluate(); }
 420     else if (t % (timesteps/100) == 0) cout << ".";
 421     forN(tracks,d) simulateTrack(track[d]);
 422     handleSprings(false);
 423     if (t % sampling == 0) saveToFile(t);
 424   }
 425 }
 426 
 427 } // End of namespace
 428 
 429 using namespace transient;
 430 
 431 int main () {
 432   // Initialize
 433   tick = setupTick; 
 434   apexSpeed = setupApexSpeed;
 435   setupLoop();
 436   // Simulate
 437   tick = simulationTick;
 438   apexSpeed = simulationApexSpeed;
 439   simulateLoop();
 440   return 0;
 441 }