# -*- coding: utf-8 -*-
"""
Created on Wed Jan  1 14:30:33 2015
ok
@author: joseorellana
"""

# import basic modules
import numpy as np
import matplotlib.pyplot as plt
from numpy import linalg as la
import math


def cm0(alpha1):
    return -.3*alpha1+.001
    
def cmg(alpha1):
    return -.1*alpha1
    
def cxg(beta0):
    return .2*beta0
    
def cz0(alpha1):
    return 10.*alpha1-.05*alpha1**2  
    
def czg(alpha1):
    return 1.*alpha1-.01*alpha1**2  
         
def dcm0(alpha1):
    return -.3
    
def dcmg(alpha1):
    return -.1
        
def dcz0(alpha1):
    return 10.-.1*alpha1**2

    






m=1.
ro=1.2
J0=1.
k=9.
c=7.
a=-0.02
S=1.
L=1
V=1.
T=20.
NPT=6000
dt=T/NPT
z0=0.
dz0=.1
alpha0=0.09
dalpha0=0.
beta0=0.1

temp=np.linspace(1,T+1,NPT)-1
beta=np.linspace(1,2,NPT)*0.
z=np.linspace(1,2,NPT)*0.
alpha=np.linspace(1,2,NPT)*0.

M=[[m, a*math.cos(alpha0)],[a*math.cos(alpha0),J0]]
K=[[k ,-ro*S*V**2*dcz0(alpha0)/2],[0., c-ro*S*V**2*L*dcm0(alpha0)/2]]

R=la.solve(M,K)

F0=[ro*S*V**2*cz0(alpha0)/2,ro*S*L*V**2*cm0(alpha0)/2]
F=la.solve(M,np.transpose(F0))

z[0]=z0
z[1]=z0+dt*dz0
alpha[0]=alpha0
alpha[1]=alpha[0]+dt*dalpha0
X1=np.transpose([z[0],alpha[0]])
X2=np.transpose([z[1],alpha[1]])


for i in range(2,NPT):
	X3=2*X2-X1+(dt**2)*(F-np.dot(R,X2))
	z[i]=X3[0]
	alpha[i]=X3[1]
	X1=X2
	X2=X3

plt.figure
plt.plot(temp,z,temp,alpha)
plt.title('Mouvements avec ou sans flutter (fonction du choix de V) ni controle')
plt.xlabel(' Temps')
plt.ylabel('Variations de z et de alpha')
plt.show()