(*     Module #13:  Coupled Harmonic Oscillators *)

(*     Program for Symmetric Initial Conditions      *)

Clear[x1,x2,v1,v2];
ti=0.;tf=32.;delta=0.4;
x1[ti]=1.0; x2[ti]=1.0; v1[ti]=0.0; v2[ti]=0.0;
m1=1.0;m2=1.0;k=1.25; kc=0.25;
Do[ F1[t] = -k*x1[t]+kc*(x1[t]-x2[t]); 
    F2[t] = -k*x2[t]+kc*(x2[t]-x1[t]);
    v1[t+delta] = v1[t] + (F1[t]/m1)*delta; 
    v2[t+delta] = v2[t] + (F2[t]/m2)*delta; 
    x1[t+delta] = x1[t] + v1[t+delta]*delta;    
    x2[t+delta] = x2[t] + v2[t+delta]*delta,
    {t,ti,tf,delta}]

(*  Make Graph of position of particle 1 vs. time  *)

x1data = Table[{t,x1[t]}, {t,ti,tf,delta}];
x1plot=ListPlot[x1data, AxesLabel->{"t","x1"},
           PlotJoined->True]

(*  Make Graph of position of particle 2 vs. time  *)

x2data = Table[{t,x2[t]}, {t,ti,tf,delta}];
x2plot=ListPlot[x2data, AxesLabel->{"t","x2"}]

(*  Show Combined graph of particles 1 and 2       *)

Show[x1plot,x2plot]

(*  Program for Anti-Symmetric Initial Conditions   *)

Clear[x1,x2,v1,v2];
ti=0.;tf=32.;delta=0.4;
x1[ti]=-1.0; x2[ti]=1.0; v1[ti]=0.0; v2[ti]=0.0;
m1=1.0;m2=1.0;k=1.25; kc=0.25;
Do[ F1[t] = -k*x1[t]+kc*(x1[t]-x2[t]); 
    F2[t] = -k*x2[t]+kc*(x2[t]-x1[t]);
    v1[t+delta] = v1[t] + (F1[t]/m1)*delta; 
    v2[t+delta] = v2[t] + (F2[t]/m2)*delta; 
    x1[t+delta] = x1[t] + v1[t+delta]*delta;    
    x2[t+delta] = x2[t] + v2[t+delta]*delta,
    {t,ti,tf,delta}]

x1data = Table[{t,x1[t]}, {t,ti,tf,delta}];
x1plot=ListPlot[x1data, AxesLabel->{"t","x1"},
           PlotJoined->True]

x2data = Table[{t,x2[t]}, {t,ti,tf,delta}];
x2plot=ListPlot[x2data, AxesLabel->{"t","x2"},
                 PlotJoined->True,PlotStyle->Dashing[{0.02}]]

Show[x1plot,x2plot]

(*   Case of More General Initial Initial Conditions   *)

Clear[x1,x2,v1,v2];
ti=0.;tf=32.;delta=0.4;
x1[ti]=1.0; x2[ti]=0.0; v1[ti]=0.0; v2[ti]=0.0;
m1=1.0;m2=1.0;k=1.25; kc=0.25;
Do[ F1[t] = -k*x1[t]+kc*(x1[t]-x2[t]); 
    F2[t] = -k*x2[t]+kc*(x2[t]-x1[t]);
    v1[t+delta] = v1[t] + (F1[t]/m1)*delta; 
    v2[t+delta] = v2[t] + (F2[t]/m2)*delta; 
    x1[t+delta] = x1[t] + v1[t+delta]*delta;    
    x2[t+delta] = x2[t] + v2[t+delta]*delta,
    {t,ti,tf,delta}]

x1data = Table[{t,x1[t]}, {t,ti,tf,delta}];
x1plot=ListPlot[x1data, AxesLabel->{"t","x1"},
           PlotJoined->True]

x2data = Table[{t,x2[t]}, {t,ti,tf,delta}];
x2plot=ListPlot[x2data, AxesLabel->{"t","x2"},
        PlotJoined->True,PlotStyle->Dashing[{0.02}]]

Show[x1plot,x2plot]