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

from math import sin
import operator # pour la définition de la classe vector
import sys # pour avoir les arguments passés en ligne de commande

# 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)
   # définitions pour la multiplication (noter le produit scalaire)
   def __mul__(self, other):
     if isinstance(other,(int, long, float, complex)):
       return vector([x*other for x in self])# résultat: vector
     else:
       return sum(map(operator.mul, self, other))# résultat: scalaire
   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])])

def integration_rk4(etat,dt):
    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)
    return etat + dt*(1.0/6)*(k1 + 2*k2 + 2*k3 + k4) # nouvel etat

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

# conditions initiales
try:# tenter de lire les conditions initiales dans les arguments
  theta0 = float(sys.argv[1])
  omega0 = float(sys.argv[2])
except ValueError:# pepin: se rabattre sur les valeurs par défaut
  theta0=0.0
  omega0=1.99
  print "usage: python pendule_rk4_full.py theta0 omega0 > pendule.data"
  print "Erreur à la lecture des arguments; défaut: (theta0, omega0) =",\
      (theta0, omega0)

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

# 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
    etat = integration_rk4(etat,dt)

    n+=1