(* Pendulo Fisico *)
(* La ecuacion diferencial del pendulo fisico. Nota que la igualdad se escribe con == y la derivada es D[f,{x,n}] donde n es el grado de derivación, x es la variable de derivacion y f es la funcion a derivar. Los ultimos elementos son las condiciones iniciales. *)
ecdif={D[θ[t],{t,2}]==-(g/l)Sin[θ[t]]+(η^2)((Cos[θ[t]])/((Sin[θ[t]])^3)),((Sin[θ[t]])^2)ϕ'[t]==η,θ[0]==θ0,θ'[0]==ω0,ϕ[0]==ϕ0}

(* Aqui le damos valores numericos a todos los parametros. Como aqui tomamos η=0, estamos estudiando el caso del pendulo normal. En este caso recuerda que el movimiento ocurre en un plano. *)
ecdifn=ecdif/.{g->9.81,l->1.0,η->0,θ0->0.5,ω0->0.0,ϕ0->π/2}

(*Aqui resolvemos el sistema de ecuaciones diferenciales y guardamos la solucion en la variable sol. Como las soluciones son numericas se guardan en polinomios de interpolacion *)
sol=NDSolve[ecdifn,{θ[t],ϕ[t]},{t,0,6.0}]

(* Aqui aplicamos la primera y unica solucion sol[[1]] a las variables y las graficamos. Nota que ϕ se mantiene constante mientras que θ oscila como uno esperaría. *)
Plot[{θ[t],ϕ[t]}/.sol[[1]],{t,0,6.0}]

(* El vector de posición en coordnadas polares es *)
r={l Sin[θ[t]]Cos[ϕ[t]],l Sin[θ[t]]Sin[ϕ[t]],l Cos[θ[t]]};
(* Ahora trazamos la trayectoria de la particula *)
Show[{Graphics3D[{Opacity[0.5],Sphere[{0,0,0},1]}],ParametricPlot3D[r/.sol/.conds,{t,0,6.0}]}]

(* Ahora vemos el caso en el que θ es contante. Para esto la fuerza centrífuga, la tensión de la cuerda y el peso de la partícula se tienen que compenzar. Del diagrama de cuerpo libre se puede ver que Tan[θ]=(((((ϕ'[t]))^2)l Sin[θ[t]])/g) y por lo tanto η=(Sqrt[(g/(l Cos[θ[t]]))])((Sin[θ[t]])^2). Nota que en la grafica θ se mantiene constante mientras que ϕ'=cte entonces ϕ=ϕ't*)
conds={g->9.81,l->1.0,θ0->0.5,ω0->0.0,ϕ0->0,η-> (Sqrt[(9.81/(1.0 Cos[0.5]))])((Sin[0.5])^2)}
ecdifn=ecdif/.{g->9.81,l->1.0,θ0->0.5,ω0->0.0,ϕ0->0,η-> (Sqrt[(9.81/(1.0 Cos[0.5]))])((Sin[0.5])^2)}
sol=NDSolve[ecdifn,{θ[t],ϕ[t]},{t,0,6.0}]
Plot[{θ[t],ϕ[t]}/.sol[[1]],{t,0,6.0}]

(* El vector de posición en coordnadas polares es *)
r={l Sin[θ[t]]Cos[ϕ[t]],l Sin[θ[t]]Sin[ϕ[t]],l Cos[θ[t]]};
(* Ahora trazamos la trayectoria de la particula *)
Show[{Graphics3D[{Opacity[0.5],Sphere[{0,0,0},1]}],ParametricPlot3D[r/.sol/.conds,{t,0,6.0}]}]

(* Ahora alejamos a la particula de esta condicion.*)
conds={g->9.81,l->1.0,θ0->0.5,ω0->0.0,ϕ0->0,η-> (Sqrt[(9.81/(1.0 Cos[0.5]))])((Sin[0.5])^2)}
ecdifn=ecdif/.{g->9.81,l->1.0,θ0->0.5,ω0->0.0,ϕ0->0,η-> (Sqrt[(9.81/(1.0 Cos[0.5]))])((Sin[0.5])^2)+0.2}
sol=NDSolve[ecdifn,{θ[t],ϕ[t]},{t,0,6.0}]
Plot[{θ[t],ϕ[t]}/.sol[[1]],{t,0,6.0}]
(* El vector de posición en coordnadas polares es *)
r={l Sin[θ[t]]Cos[ϕ[t]],l Sin[θ[t]]Sin[ϕ[t]],l Cos[θ[t]]};
(* Ahora trazamos la trayectoria de la particula *)
Show[{Graphics3D[{Opacity[0.5],Sphere[{0,0,0},1]}],ParametricPlot3D[r/.sol/.conds,{t,0,6.0}]}]