public
Created

Outputs SVG image of swinging Atwood's machine trajectory given the mass ratio, initial distance from the centre, initial angle, initial radial velocity, and initial angular velocity. The numerical scheme used is the fourth order Runge Kutta. This generated the trajectory images on https://en.wikipedia.org/wiki/Swinging_Atwood's_Machine

  • Download Gist
sam.cpp
C++
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
/*
Swinging Atwood's Machine simulator by Daniel L. Lu
Outputs an SVG image showing the orbit for initial conditions specified via standard input.
The trajectory is normalized such that the furthest distance from the origin is always 500 pixels.
 
Copyright (c) 2013, Daniel L. Lu
All rights reserved.
 
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
 
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
 
#include <iostream>
#include <fstream>
#include <cmath>
 
using namespace std;
 
const long double h = 0.001L, mug = 45.0L, PI = 3.1415926535897932384626433832795028L;
 
void derivatives(long double *y, long double *dy, long double mu) {
long double g = mug/mu;
dy[0] = y[2]/(1+mu); //r
dy[1] = y[3]/powl(y[0],2); //theta
dy[2] = powl(y[3],2)/powl(y[0],3) - g*(mu-cosl(y[1])); //pr
dy[3] = -g*y[0]*sinl(y[1]);
}
 
void rk4(long double *x, long double mu) {
/* fourth order Runge Kutta integrator */
long double dx[4], yt[4], dym[4], dyt[4];
derivatives(x, dx, mu); // k1
for(int i=0; i<4; i++) {
yt[i] = x[i] + h/2 * dx[i];
}
derivatives(yt, dyt, mu); // k2
for(int i=0; i<4; i++) {
yt[i] = x[i] + h/2 * dyt[i];
}
derivatives(yt, dym, mu); // k3
for(int i=0; i<4; i++) {
yt[i] = x[i] + h * dym[i];
dym[i] += dyt[i];
}
derivatives(yt, dyt, mu); // k4
for(int i=0; i<4; i++) {
x[i] += h/6*(dx[i]+dyt[i]+2*dym[i]);
}
}
 
int main() {
long double mu, r0, t0, dr0, dt0;
while(cin >> mu >> r0 >> t0 >> dr0 >> dt0) {
/* open a file for writing */
char filename[200];
sprintf(filename, "sam_mu_%Lf_r0_%Lf_t0_%Lf_dr0_%Lf_dt0_%Lf.svg", mu, r0, t0, dr0, dt0);
fstream fout(filename, fstream::out);
long double *x = new long double[4];
/* set the initial conditions */
x[0] = r0; // initial radius
x[1] = t0; // initial angle
x[2] = dr0; // initial radial velocity
x[3] = dt0; // initial angular velocity
/* yay svg */
fout << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" width=\"1200px\" height=\"1200px\">" << endl;
/* draw the x and y axes */
fout << "<line x1=\"0\" y1=\"600\" x2=\"1200\" y2=\"600\" style=\"fill:none;stroke:#eeeeff;stroke-width:1\"/>" << endl;
fout << "<line x1=\"600\" y1=\"0\" x2=\"600\" y2=\"1200\" style=\"fill:none;stroke:#eeeeff;stroke-width:1\"/>" << endl;
/* draw the concentric circles (gridlines) */
for(int i=10; i<600; i+= 10) {
fout << "<circle cx=\"600\" cy=\"600\" r=\"" << i << "\" style=\"fill:none;stroke:#eeeeff;stroke-width:1\"/>" << endl;
}
/* draw a polyline that is the entire trajectory */
fout << "<polyline points=\"";
/* simulate for 50000 points at most... arbitrarily chosen so that the svg doesn't take too long to load */
long double r[50000], t[50000], rr=0;
int imax = 50000; // not to be confused with the cinema technology
for(int i=0; i<imax; i++) {
/* break if it becomes too close to singularity (although sometimes it gets imprecise even when r > 3)*/
if(i>10 && x[0]<3) {
imax = i;
break;
}
/* keep track of the largest possible radius to scale it such that it doesn't go out of bounds */
if(x[0]>rr) rr = x[0];
/* copy simulation results to an array */
r[i] = x[0];
t[i] = x[1];
/* actually in the simulation we do 100 steps of simulation for each point drawn to the svg */
for(int j=0; j<100; j++) rk4(x,mu);
}
/* output the simulation results to the svg */
for(int i=0; i<imax; i++) {
fout << 600+500*r[i]/rr*sinl(t[i]) << ',' << 600+500*r[i]/rr*cosl(t[i]); if(i < imax-1) fout << ' ';
}
fout << "\" style=\"fill:none;stroke:black;stroke-width:1\" />\n</svg>" << endl;
fout.close();
cout << filename << endl;
}
return 0;
}

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.