# -*- coding: utf-8 -*-
# intégration du mouvement du pendule par Runge-Kutta 4
# version minimaliste avec vecteur
# executer comme suit:
#   python pendule_rk4_mini.py > pendule.data
# 2008 Adrian Daerr, Univ Paris Diderot
# domaine public

from math import sin
import operator # pour la définition de la classe vector

# Definir une classe vecteur.
#
# Nous étendons la classe list, de laquelle nous héritons presque tout
# ce que nous voulons. Il reste juste à définir comment additioner
# deux vecteurs ou un vecteur et un scalaire, et de même pour la
# multiplication.
class vector(list):
   def __add__(self, other):# est appelée pour l'opération vecteur+other
     if isinstance(other,(int, long, float, complex)):# vecteur+scalaire
       return vector([x+other for x in self])
     else:# tenter d'appliquer élément par élément
       return vector(map(operator.add, self, other))
   def __radd__(self, other):# X+vecteur
       return (self+other)# simplement égal à vecteur+X (voir ci-dessus)
   # mêmes définitions pour la multiplication
   def __mul__(self, other):
     if isinstance(other,(int, long, float, complex)):
       return vector([x*other for x in self])
     else:
       return vector(map(operator.mul, self, other))
   def __rmul__(self, other):# X*vector
       return (self*other)   # = vector*X

# equations du mouvement du pendule rigide
def deriv(etat):
  return vector([etat[1], -sin(etat[0])])

# --- début du programme ---

# état du système; ordre par convention des variables d'état: [theta, omega]
etat = vector([0.0, 1.99])# conditions initiales

# boucle principale
t, tfinal = 0.0, 50.0; # intervalle de simulation
dt=0.05 # pas de temps
n=0
while t<tfinal:
    t = n*dt

    # afficher l'état du système
    print t, etat[0], etat[1]

    # integrer d'un pas de temps
    k1 = deriv(etat)
    etat_int = etat + 0.5*dt*k1
    k2 = deriv(etat_int)
    etat_int = etat + 0.5*dt*k2
    k3 = deriv(etat_int)
    etat_int = etat +     dt*k3
    k4 = deriv(etat_int)
    etat = etat + dt*(1.0/6)*(k1 + 2*k2 + 2*k3 + k4)

    n+=1