#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
#%matplotlib notebook


# In[2]:


import numpy as np
import matplotlib.pyplot as plt


# In[ ]:





# In[3]:


#LLT inviscid implementation for rectangular and trapezoidal wings with and without twist


# In[4]:


#Aerodynamics for engineers. J. Bertin, R. Cummings. 6ed. Example 7.2

rho = 1.225
v_inf = 100

target_weight = 4000

#Number of setions
n = 40       #D(4)

#Manual input of AR and TR
AR = 9     #D(9)
TR = 0.4   #D(0.4)

c_root = 0.726  #D(0.726)
c_tip = 0.290   #D(0.290)

b      = 2*2.286  #semi-span - one wing
b_half = 2.286   #semi-span - one wing - D(2.286)

#Area when using wing-span b
#area = ((b)**2)/AR

#Area when using wing semi-span b/2
#This is the area for one wing.
#area = ((2*b_half)**2)/AR
area = (0.5*(c_root + c_tip))*b_half

#Computed AR using compute area
#Attention is area computed using one b/2 you must multiply area by 2
AR_computed = b**2/(2*area)

#Cmputed TR
TR_computed = c_tip/c_root

#Airfoil cl slope
a0 = 2*np.pi

#Alpha zero-lift
alpha_0 = -1.2  #D(-1.2)

#i_w = 2            # wing setting angle (deg)
#alpha_twist = -1;  # Twist angle (deg)

AOA = 2;            # Wing angle of attack (deg) - D(2)
i_w = 2;             # wing setting angle (deg) - This is the angle at the root - D(2)
alpha_twist = 0;    # Twist angle (deg), for washout the value is negative - D(0)

alpha_effective = AOA + i_w    #Actual AOA at root (deg)

#alpha = i_w+alpha_twist:-alpha_twist/(N-1):i_w;
#alpha = np.linspace(i_w+alpha_twist,i_w,n)  
alpha = np.linspace(alpha_effective+alpha_twist,alpha_effective,n)   #(deg)


# In[5]:


print("Computed area for semi wingspan = ", 1*area)
print("Computed area for wingspan = ", 2*area)
print("Computed AR = ", AR_computed)
print("Computed TR = ", TR_computed)


# In[6]:


print("Angle of attack AOA (deg) = ", AOA)
print("Wing setting angle (deg) = ", i_w )
print("Effective angle of attack at the root (AOA + wing setting angle) = ", alpha_effective)
print("Geometrical twist in deg (negative means washout) = ", alpha_twist)
#print("Twist distribution (in degrees) = ", alpha)
print("Twist distribution tip-to-root (deg) = ", alpha)


# In[ ]:





# In[7]:


#theta = np.linspace(np.pi/(2*n),np.pi/2,n)

if n == 4:
    #theta_force = 22.5    #Force theta so we have same setup as in bertin
    #theta = np.linspace(0,90,n-1) + theta_force
    theta = np.array([22.5,45,67.5,90])
else:
    eps = 1e-8     #Do not start from zero in will divide by zero
    theta = np.linspace(eps,90,n)
    

theta_r = theta*np.pi/180
theta_deg = theta


# In[8]:


print("Theta angle (rad) = ", theta_r)


# In[9]:


print("Theta angle (deg) = ", theta_deg)


# In[ ]:





# In[10]:


#Compute spanwise stations using angular values

#When using wing-span b
z = (b/2)*np.cos(theta_r)

#When using wing semi-span b/2
#z = 0.5*(2*b_half)*np.cos(theta_r)


# In[11]:


print("Span location (+s - center line) = ", z)


# In[12]:


print("Twist distribution (in degrees) = ", alpha)


# In[ ]:





# In[13]:


#Compute local chord

#c1 = 1 + (TR - 1) * np.cos(theta_r) * c_root
c = c_root *(1-(1-TR)*np.cos(theta_r))


# In[14]:


print("Local chord = ", c)


# In[ ]:





# In[15]:


#mu = ( a0/(2*(1 + TR)*AR) ) * ( 1+(TR-1)*np.cos(theta_r) );
mu = (c*a0)/(2*AR*c_root*(1+TR))


# In[16]:


print("Local mu = ", mu)


# In[ ]:





# In[17]:


#Assemble matrices

#Ax=b
#x=A^-1*b
#In our case the matrices are Bx=LHS

#Angles alpha and alpha_0 must be in RAD
LHS = mu * (alpha - alpha_0)/57.3

#Coefficients of monoplane equation
B = np.zeros((n,n))

for i in range(0,n):
    for j in range(0,n):
        B[i,j] = np.sin((2*(j+1)-1) * theta_r[i]) * ( 1 + ((mu[i] * (2*(j+1)-1)) / np.sin(theta_r[i])) )


# In[ ]:





# In[18]:


#Solve linear system for x - Assign x to A for consistency in notation

x = np.linalg.solve(B,LHS)
A = x


#Compute lift coefficient of the wing
CL_wing = np.pi * AR * A[0]


# In[19]:


print("An coefficients = ", A)


# In[ ]:





# In[20]:


#Compute delta or induced drag factor

delta=0
for i in range(1,n):
    delta = delta + (2*(i+1)-1)*A[i]**2/A[0]**2

#Attention in the previous equation i=1 is A[3] in the equation of llt

    
#Compute oswald span efficiency
oswald = 1/(1+delta)

#Compute induced drag
CD_ind = (CL_wing**2/(np.pi*AR)) * (1+delta)


# In[21]:


print("Induced drag factor (delta) = ", delta)
print("Oswald span efficiency (e) = ", oswald)
print("Induced drag coefficient CD_i = ", CD_ind)


# In[ ]:





# In[22]:


#For lift coefficient distribution

sum1=np.zeros((n))

for i in range(0,n):
    for j in range(0,n):
    
#For lift coefficient distribution
        sum1[i] = sum1[i] + A[j]*np.sin((2*(j+1)-1)*theta_r[i])

    
#Atention b is +s 
#C_L distribution is equal to 2*gamma/u.c
#gamma is equal to 4*s*u*(sum)
gamma = 4*b*sum1
#gamma = 4*b*v_inf*sum1
#gamma0=gamma[n-1]

C_L = 2*4*b*sum1 / c;

#lift distribution per unit span
Ls = C_L * 0.5 * rho * v_inf * v_inf * c


# In[23]:


print("Lift coefficient CL distribution = ", C_L)


# In[24]:


print("Lift distribution = ", Ls)


# In[ ]:





# In[25]:


#Print lift and drag coefficients

print("CL = ",CL_wing)
print("CD_ind = ",CD_ind)

#CLalpha = CL_wing/( (i_w*np.pi/180) - (alpha_0*np.pi/180))
CLalpha = CL_wing/( (AOA*np.pi/180) - (alpha_0*np.pi/180))

print("CL_alpha (slope lift curve wing per rad) = ",CLalpha)


# In[ ]:





# In[26]:


#Compute lift and darg forces

Lift = CL_wing * 0.5 * rho * v_inf * v_inf * area
Drag_ind = CD_ind * 0.5 * rho * v_inf * v_inf * area


print("Reference velocity (m/s) = ", v_inf)
print("Reference density (kg/m^3) = ", rho)
print("Reference area (m^2) = ", area)


# In[27]:


print("Lift (N) = ", Lift)
print("Induced drag (N) = ", Drag_ind)


# In[28]:


v_target = np.sqrt(target_weight/(0.5*rho*area*CL_wing))
print("Required lift = ", target_weight)
print("Target velocity to generate required lift = ", v_target)


# In[ ]:





# In[29]:


#Induced angle and downwash computation

sum2=np.zeros((n))

for i in range(0,n):
    for j in range(0,n):
        sum2[i] = sum2[i] + ( (2*(j+1)-1) * A[j] * np.sin((2*(j+1)-1) * theta_r[i]) ) / np.sin(theta_r[i])

#for i in range(0,n):
#    sum2[i] = sum2[i] + ( (2*(i+1)-1) * A[i] * np.sin((2*(i+1)-1) * theta_r[i]) ) / np.sin(theta_r[i])

#Downwash
dw = -1*v_inf * sum2

#Induced angle
ind_a = -1*sum2


# In[30]:


print("Downwash (m/s) = ", dw)


# In[31]:


print("Induced angle (rad) = ", ind_a)


# In[32]:


print("Induced angle (deg) = ", ind_a*180/np.pi)


# In[ ]:





# In[33]:


plt.figure(figsize=(12,8))

#plt.plot(z,C_L)
#plt.plot(z,C_L/CL_wing)
#plt.plot(z,C_L/(2*CL_wing))
#plt.plot(z,C_L/np.max(C_L))

#Normalized +s
#plt.plot(z/np.max(z),C_L)
#plt.plot(z/np.max(z),C_L/(CL_wing))
plt.plot(z/np.max(z),C_L/(2*CL_wing))
#plt.plot(z/np.max(z),C_L/np.max(C_L))


#plt.hlines(y=np.max(C_L), xmin=0.0, xmax=b/1, color='b', label="cl_max")
#plt.hlines(y=C_L[n-1], xmin=0.0, xmax=b/1, color='r', label="cl_root")

plt.hlines(y=np.max(C_L/(2*CL_wing)), xmin=0.0, xmax=z/np.max(z),linestyles='dashdot',color='b',label="$C_{l_{max}}$")
plt.hlines(y=C_L[n-1]/(2*CL_wing), xmin=0.0, xmax=z/np.max(z), linestyles='dashdot',color='r', label="$C_{l_{root}}$")


#plt.xlabel('+s (m)')
plt.xlabel('$\dfrac{+s}{+s_{max}}$',fontsize=14)
#plt.ylabel(r'$c_l$',fontsize=14)
plt.ylabel(r'$\dfrac{C_l}{2*C_L}$',fontsize=14)
plt.title('Lift coefficient distribution',fontsize=20)
plt.grid()

plt.plot(0,np.max(C_L), 1,0.4)

plt.legend(loc=0,fontsize=18)


# In[ ]:





# In[34]:


plt.figure(figsize=(12,8))

#plt.plot(z,ind_a*180/np.pi)
plt.plot(z/np.max(z),ind_a*180/np.pi)

#For downwash
#plt.plot(z,dw)

#plt.xlabel('+s (m)')
plt.xlabel('$\dfrac{+s}{+s_{max}}$',fontsize=14)
plt.ylabel('Induced angle $(^{\circ})$',fontsize=14)
plt.title('Induced angle distribution',fontsize=20)
plt.grid()

#plt.legend(loc=1)


# In[ ]:





# In[35]:


plt.figure(figsize=(12,8))

#plt.plot(z/np.max(z),Ls)
#plt.plot(z/np.max(z),Ls/np.max(Ls))
plt.plot(z,Ls,label="Computed lift distribution")
#plt.plot(z/np.max(z),Ls)

#elliptical lift distribution
el=np.sqrt( (1 - (z/(b/2))**2)*np.max(Ls)**2)
plt.plot(z,el,'-.',label="Ideal elliptical lift distribution")

plt.xlabel('$+s$ (m)',fontsize=14)
#plt.xlabel('$\dfrac{+s}{+s_{max}}$',fontsize=14)
plt.ylabel('Lift (N)',fontsize=14)
plt.title('Lift distribution',fontsize=20)
plt.grid()

plt.legend(loc=0,fontsize=14)


# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]:





# In[ ]: