# -*- coding: utf-8 -*-
"""
Created on Tue Jan 20 17:42:59 2015

Listing 11.4 : mouvement du pendule par l'algorithme de
prédiction-correction d'ordre 2
z = y'

@author: grivet
"""

from pylab import *

def smy(t,y,z):
    return z 

def smz(t,y,z):
    return -sin(y)

def pred(t,y,z):
    yf = y + 2*pas* smy(t,y,z)
    zf = z + 2*pas*smz(t,y,z)
    return yf,zf

def corr(t,yi,zi,yp,zp):
    yf = yi+ 0.5*pas*(smy(t,yi,zi)+smy(t,yp,zp));
    zf = zi + 0.5*pas*(smz(t,yi,zi)+smz(t,yp,zp));
    return yf, zf

def pc(tnm1,ynm1,znm1,tn,yn,zn):
    yp,zp = pred(tnm1,ynm1,znm1)
    ynp1,znp1 = corr(tn,yn,zn,yp,zp)
    tnp1 = tn + pas
    return ynp1,znp1,tnp1
    
np = 2000
t = zeros((np,1))
y = zeros_like(t); z = zeros_like(t)
y0 = float(input("position initiale: "))
z0 = float(input("vitesse initiale: "))
pas = float(input("longueur du pas: "))
t[0] = 0; y[0] = y0; z[0] = z0
t[1] = pas; y[1] = y0+pas*z0+0.5*pas*pas*z0;
z[1] = z0+pas*z0;
for i in range(1,np-1):
    y[i+1],z[i+1],t[i+1] = pc(t[i-1],y[i-1],z[i-1],t[i],y[i],z[i])

plot(t,y)
title("pendule simple")
xlabel("temps")
ylabel("angle(rd)")