## Attachment 'main.cpp'

```   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
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; // 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*b + a*b + a*b; }
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 = {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 = {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 = {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; // [m]
92   double velocity; // [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; // 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; // [m/s]
104   double incomingPosition; // [m]
105
106   // Spring variables
107   Spring* shearStress; // 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; // Remember that pointers can't handle being copied
143   double dampingForce; // [N]
144   double springConstant; // [kg/s^2]
145   double vector; // [m] Projection of vector from end to end onto spring direction
146   void contract();
147   double dampingConstant() { return 2 * sqrt(springConstant * (end->mass() + end->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;
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; // 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 && stableSpeed < minMaxSpeed) || (stableSpeed < minMaxSpeed && stableSpeed > minMaxSpeed)) {
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 = this;
202     (*c)->end = 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; // [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, (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[a] += (wire[d][w]->end == this? 1 : -1) * wire[d][w]->vector[a];
234     cross(axis, axis, axis);
235   } else { // Use north to define sideways
236     axis = 1;
237     cross(axis, axis, axis);
238     cross(axis, axis, axis);
239   }
240   for3(a) normalize(axis[a]);
241   // Next, do projections, reducing to one dimensional strains
242   double projection; // 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; // [0 is guy wire direction, 1 is sideways]
246   for2(i) deflection[i] = projection[i+1] - (projection[i+1] + projection[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 == 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 = 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 = (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
300     traceTrack(t);
301     trackDensity += adjustment * (t.height() > 0? 1 : -1);
302   }
303   cout << "Track density " << trackDensity << endl;
304   cout << "Rotor speed at ground " << t.speed() << endl;
305 }
306 void applyHelix(Segment* t, int d)
307 {
308   double north = {0,1,0};
309   double up[segments+1];
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 == 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);
345   forN(tracks,d) forN(segments+1,s) track[d][s] = track[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 / 1000 << " ";
373     file << s->position << " ";
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;
385   for2(i) maxRotor[i] = track.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.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.incomingMass = t.rotorMass / t.length() * t.designSpeed * tick;
407   for3(a) t.incomingVelocity[a] = t.designSpeed * t.direction[a];
408   for3(a) t.incomingPosition[a] = (t.position[a] + t.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 }
```

## Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
• [get | view] (2011-05-24 15:10:38, 1.9 KB) [[attachment:LaunchLoopScript]]
• [get | view] (2011-05-24 14:43:53, 228.0 KB) [[attachment:NumericalSimulation1]]
• [get | view] (2011-05-24 15:00:36, 228.0 KB) [[attachment:NumericalSimulation1.gif]]
• [get | view] (2011-05-24 14:50:37, 557.1 KB) [[attachment:NumericalSimulation2.gif]]
• [get | view] (2011-05-24 15:06:56, 20.3 KB) [[attachment:main.cpp]]
All files | Selected Files: delete move to page