# -*- coding: utf-8 -*-
"""
Created on Tue Jan 20 15:55:08 2015

Listing 11.2 : pendule par l'algorithme de Runge Kutta
d'ordre 2

@author: grivet
"""

from pylab import *

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

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

def rk2(h,ti,yi,zi):
    yf = yi + 0.5*h*(smy(ti,yi,zi)+smy(ti+h,yi+h*smy(ti,yi,zi),zi+h*smz(ti,yi,zi))) 
    zf = zi + 0.5*h*(smz(ti,yi,zi)+smz(ti+h,yi+h*smy(ti,yi,zi),zi+h*smz(ti,yi,zi)))
    tf = ti + h
    return   yf,zf,tf
  
np = 1200
t = zeros((np,1)) 
y = zeros_like(t); z = zeros_like(t)
y0 = float(input('position initiale(rd) : '))
z0 = float(input('vitesse initiale(rd/s) : '))
pas = float(input('longueur du pas: '))
t[0] = 0; y[0] = y0; z[0] = z0
for i in range (1,np):
    y[i],z[i],t[i] = rk2(pas,t[i-1],y[i-1],z[i-1]);

plot(t,y)

#ep = 1-cos(y);
#ec = 0.5*z*z
#plot(t,[ep,ec,(ep+ec-ep(1)-ec(1))*1000])