#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
// #define ONE_PERIOD_ONLY // uncomment this to do only one period, then print stats and exit.
void get_forces(double *x, double *y, double *Fx, double *Fy, int N, double k, double r0) // pointers here are array parameters
{
int i;
double r;
Fx[0]=Fy[0]=Fx[N]=Fy[N]=0; // have to set the end forces properly to avoid possible uninitialized memory shenanigans
double U=0,E=0,T=0;
for (i=1;i<N;i++)
{
// left force
r = hypot(x[i]-x[i-1],y[i]-y[i-1]);
Fx[i] = -(x[i]-x[i-1]) * k * (r-r0)/r;
Fy[i] = -(y[i]-y[i-1]) * k * (r-r0)/r;
// right force
r = hypot(x[i]-x[i+1],y[i]-y[i+1]);
Fx[i] += -(x[i]-x[i+1]) * k * (r-r0)/r;
Fy[i] += -(y[i]-y[i+1]) * k * (r-r0)/r;
}
}
void evolve_euler(double *x, double *y, double *vx, double *vy, int N, double k, double m, double r0, double dt)
{
int i;
double Fx[N+1],Fy[N+1]; // this could be made faster by mallocing this once, but we leave it this way for students
// to avoid having to deal with malloc(). In any case memory allocation is faster than a bunch
// of square root calls in hypot().
get_forces(x,y,Fx,Fy,N,k,r0);
for (i=1;i<N;i++)
{
x[i] += vx[i]*dt;
y[i] += vy[i]*dt;
vx[i] += Fx[i]/m*dt;
vy[i] += Fy[i]/m*dt;
}
}
// this function is around to go from Euler-Cromer to leapfrog, if we want second-order precision
void evolve_velocity_half(double *x, double *y, double *vx, double *vy, int N, double k, double m, double r0, double dt)
{
int i;
double Fx[N+1],Fy[N+1]; // this could be made faster by mallocing this once, but we leave it this way for students
get_forces(x,y,Fx,Fy,N,k,r0);
for (i=1;i<N;i++)
{
vx[i] += Fx[i]/m*dt/2;
vy[i] += Fy[i]/m*dt/2;
}
}
// Students might not be familiar with pass-by-reference as a trick for returning multiple values yet.
// Ideally they should be coding this anyway, and there are a number of workarounds, in particular
// just not using a function for this.
void get_energy(double *x, double *y, double *vx, double *vy, int N, double k, double m, double r0, double *E, double *T, double *U)
{
*E=*T=*U=0;
int i;
double r;
for (i=0;i<N;i++)
{
*T+=0.5*m*(vx[i]*vx[i] + vy[i]*vy[i]);
r = hypot(x[i]-x[i+1],y[i]-y[i+1]);
*U+=0.5*k*(r-r0)*(r-r0);
}
*E=*T+*U;
}
// does what it says on the tin
void evolve_euler_cromer(double *x, double *y, double *vx, double *vy, int N, double k, double m, double r0, double dt)
{
int i;
double Fx[N+1],Fy[N+1]; // this could be made faster by mallocing this once, but we leave it this way for students
for (i=1;i<N;i++)
{
x[i] += vx[i]*dt;
y[i] += vy[i]*dt;
}
get_forces(x,y,Fx,Fy,N,k,r0);
for (i=1;i<N;i++)
{
vx[i] += Fx[i]/m*dt;
vy[i] += Fy[i]/m*dt;
}
}
// does what it says on the tin
void evolve_rk2(double *x, double *y, double *vx, double *vy, int N, double k, double m, double r0, double dt)
{
int i;
double Fx[N+1],Fy[N+1]; // this could be made faster by mallocing this once, but we leave it this way for students
double xh[N+1],yh[N+1],vxh[N+1],vyh[N+1];
vxh[0]=vyh[0]=vxh[N]=vyh[N]=0;
get_forces(x,y,Fx,Fy,N,k,r0);
for (i=0;i<=N;i++)
{
xh[i] = x[i] + vx[i]*dt/2;
yh[i] = y[i] + vy[i]*dt/2;
vxh[i] = vx[i] + Fx[i]/m*dt/2;
vyh[i] = vy[i] + Fy[i]/m*dt/2;
}
get_forces(xh,yh,Fx,Fy,N,k,r0);
for (i=0;i<=N;i++) // need two for loops -- can't interleave halfstep/fullstep updates (students have trouble with this sometimes!)
{
x[i] = x[i] + vx[i]*dt;
y[i] = y[i] + vy[i]*dt;
vx[i] = vx[i] + Fx[i]/m*dt;
vy[i] = vy[i] + Fy[i]/m*dt;
}
}
// function to encapsulate determining whether we need to shovel another frame to the animator. delay is the delay in msec.
// I usually aim to hide library calls, like clock() and CLOCKS_PER_SEC, from students in the beginning, since it's not really
// relevant to their development of computational skills, which are what I really care about.
int istime(int delay)
{
static int nextdraw=0;
if (clock() > nextdraw)
{
nextdraw = clock() + delay * CLOCKS_PER_SEC/1000.0;
return 1;
}
return 0;
}
int main(int argc, char **argv)
{
int i,N=80; // number of links, not number of nodes!! Careful for the lurking fencepost errors
int modenumber=3; // put in some defaults just in case
double t, dt=2e-6, amplitude=0.1;
double stiffness=10, density=1, length=1; // unstretched properties of original string
double k, m, r0; // properties of single string
double tension=1,Ls;
if (argc < 9) // if they've not given me the parameters I need, don't just segfault -- tell the user what to do, then let them try again
{
printf("!Usage: <this> <N> <modenumber> <dt> <amplitude> <stiffness> <density> <length> <tension>\n");
exit(0);
}
N=atoi(argv[1]);
modenumber=atoi(argv[2]);
dt=atof(argv[3]);
amplitude=atof(argv[4]);
stiffness=atof(argv[5]);
density=atof(argv[6]);
length=atof(argv[7]);
tension=atof(argv[8]);
double x[N+1], y[N+1], vx[N+1], vy[N+1], E, T, U;
// compute microscopic properties from macroscopic ones
r0=length/N;
m=density*length/N;
k=stiffness*N/length;
// figure out stretched length
Ls=length + tension * length / stiffness;
// make predictions based on what our freshman mechanics class taught us
double density_stretched = density * length / Ls;
double wavespeed = sqrt(tension/density_stretched);
double period_predict = 2 * Ls / wavespeed / modenumber;
double vym_last=0;
int monitor_node = N/modenumber/2; // this is the node that we'll be watching to see when a period has elapsed.
int nperiods=0;
for (i=0;i<=N;i++) // remember, we have N+1 of these
{
x[i] = Ls*i/N;
y[i] = amplitude*sin(modenumber * M_PI * x[i] / Ls);
vx[i]=0;
vy[i]=0;
}
// now, loop over time forever...
printf("!N\tmode\tgamma\t\tL\t\tamp.\t\tdens.\tTension\tT\tT_pred\tDelta\t\tFreq (Hz)\n"); // print header
for (t=0;1;t+=dt)
{
vym_last=vy[monitor_node];
evolve_euler_cromer(x,y,vx,vy,N,k,m,r0,dt);
// "if we were going up, but now we're going down, then a period is complete"
// this is crude and will fail if it has enough "wobble", but it's sufficient for this project.
if (vym_last > 0 && vy[monitor_node] < 0)
{
nperiods++;
printf("!%d\t%d\t%.2e\t%.2e\t%.4e\t%.1e\t%.1e\t%.4f\t%.4f\t%.4e\t%.4e\n",N,
modenumber,stiffness,length,amplitude,density,tension,t,period_predict*nperiods,1-t/period_predict/nperiods,nperiods/t);
#ifdef ONE_PERIOD_ONLY
printf("Q\n"); // kill anim before we die ourselves
exit(1);
#endif
}
if (istime(15)) // wait 15ms between frames sent to anim; this gives us about 60fps.
{
printf("C 1 0 0\nc %f %f 0.02\nC 1 1 1\n",x[monitor_node],y[monitor_node]); // draw a big red blob around the node we're watching
for (i=0;i<=N;i++)
{
printf("C %e 0.5 %e\n",0.5+y[i]/amplitude,0.5-y[i]/amplitude); // use red/blue shading; this will make vibrations visible even if amp<<1
printf("c %f %f %f\n",x[i],y[i],length/N/3); // draw circles, with radius scaled to separation
if (i<N) printf("l %f %f %f %f\n",x[i],y[i],x[i+1],y[i+1]); // the if call ensures we don't drive off the array
}
printf("T -0.5 -0.7\ntime = %f\n",t);
get_energy(x,y,vx,vy,N,k,m,r0,&E,&T,&U);
printf("T -0.5 -0.63\nenergy = %e + %e = %e\n",T,U,E);
printf("F\n"); // flush frame
}
}
}