(* Set the values of the various parameters *)
c = 0.2;
w = 0.694;
F = 0.52;
(* Increase F in steps of 0.1 to see the different solutions in the chapter notes. *)
pi = N[Pi];
cycles = 50;
steps = 30;


(* Solve the differential equations with the given set of initial conditions. *)
sol = NDSolve[
{x'[t]==v[t],
v'[t] == -c*v[t]- Sin[x[t]] + F*Cos[w*t],
x[0] == 0.8,
v[0] == 0.8},
{x,v},{t,0,(cycles*(2Pi)/w)},MaxSteps -> 20000];


(* Create a function "reduce" that translate all angles back to the region between -Pi and +Pi *)
reduce[x_] := Mod[x,2pi] /;Mod[x, 2pi] <= pi;
reduce[x_] := (Mod[x,2pi] - 2pi) /;Mod[x, 2pi] > pi;


(* Create a table with the angles mapped to the -Pi to Pi region. *)
points = Flatten[
Table[
{x[t], v[t]}/.sol,
{t, 0, (cycles*(2Pi)/w),(1/steps)(2Pi/w) // N}
], 1
];
xposition = Table[{i, 1},{i, Length[points]}];
newpoints = MapAt[reduce, points,xposition];


(* Plot the solution of the differential equations: v[t] vs t. *)
ParametricPlot[
Evaluate[{t,v[t]} /. sol],{t, 0, (cycles*(2Pi)/w)}, 
PlotRange -> {{((cycles-20)*(2Pi)/w),(cycles*(2Pi)/w)},{-3,3}},
AxesLabel -> {"t/t0","v (rad/s)"},
AspectRatio->1/2
]


(* Plot the phase diagram. *)
ParametricPlot[
Evaluate[{x[t],v[t]} /. sol],{t, ((cycles-20)*(2Pi)/w), (cycles*(2Pi)/w)}, 
AxesLabel -> {"x (rad)","v (rad/s)"},
AspectRatio->1/2
]


(* Plot the phase diagram (using the reduced angle). *)
ListPlot[
newpoints,
PlotRange -> {{-3,3},{-3,3}},
AxesLabel -> {"x (rad)","v (rad/s)"},
AspectRatio->1/2
]


(* Make the poincare plot. *)
poincare = Table[newpoints[[n]],{n, 1+25steps, Length[points], steps}];
Length[poincare];
ListPlot[
poincare,
PlotRange -> {{-3,3},{-3,3}},
AxesLabel -> {"x (rad)","v (rad/s)"},
PlotStyle->PointSize[0.015],
AspectRatio->1/2
]