# -*- coding: utf-8 -*-
"""
Created on Tue Feb 25 21:38:12 2014

@author: Quentin Renaud
"""
# Read the data file
import numpy as np
from numpy import pi
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
import os

# Input file [Angle_deg Cl Cd Cm]
folder = "."
os.chdir(folder)
filename = "DU21_A17.dat"

coefs = np.genfromtxt(filename, dtype=float, skip_header=14)
# Convert first column to radians
coefs[:, 0] = pi/180*coefs[:, 0]
angles  = coefs[:, 0]
CL      = coefs[:, 1]
CD      = coefs[:, 2]

filename=filename.replace('.dat', '.npy')
np.save(filename, coefs)

def poly_trigo(a, theta, N=5):
    p = a[0]    
    for i in range(1, N+1):
        p = p + a[2*i-1]*np.cos(i*theta)
    for i in range(1, N+1):
        p = p + a[2*i]*np.sin(i*theta)
    return p
    

def residual(a, y, theta, N):
    err = y - poly_trigo(a, theta, N=N)
    return err

N = 15
a = np.ones(2*N+1)
CL_poly = leastsq(residual, a, args=(CL, angles, N))
CL_poly = CL_poly[0]
filename = filename.replace('.npy', '_CL.npy')
np.save(filename, CL_poly)
print(CL_poly)

CD_poly = leastsq(residual, a, args=(CD, angles, N))
CD_poly = CD_poly[0]
filename = filename.replace('CL', 'CD')
np.save(filename, CD_poly)
print(CD_poly)


plt.plot(angles*180/pi, CL, 'r')
plt.plot(angles*180/pi, poly_trigo(CL_poly, angles, N), 'r+')
plt.plot(angles*180/pi, CD, 'b')
plt.plot(angles*180/pi, poly_trigo(CD_poly, angles, N), 'b+')
plt.show()

@@ -1272,16 +1272,22 @@
 y)%0D%0A%0D%0A%0D%0A
+p1, = 
 plt.plot
@@ -1308,20 +1308,27 @@
 , CL, 'r
++
 ')%0D%0A
+p2, = 
 plt.plot
@@ -1377,21 +1377,26 @@
 , N), 'r
-+
 ')%0D%0A
+p3, = 
 plt.plot
@@ -1417,20 +1417,27 @@
 , CD, 'b
++
 ')%0D%0A
+p4, = 
 plt.plot
@@ -1482,26 +1482,155 @@
 gles, N), 'b
-+
 '
+)%0D%0Aplt.title(%22Coefficients de l'aile en fonction de AoA%22)%0D%0Aplt.legend(%5Bp1,p2,p3,p4%5D,%5B%22Clift%22, %22Poly_Clift%22, %22Cdrag%22, %22Poly_Cdrag%22%5D
 )%0D%0Aplt.show(

Add legends