

In [20]:

    
%pylab inline
import pylab as pl
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

## coefficients. There are differrent coefficients for different r/c 9=(percentage)

# r/c = 0.0
# number of sampling points 
i_Nx = 200
i_x = 0.15
# define coefficients
a = 0.06666666666667
d = 0.08833866691267
b = 0.21501600073802
c = -0.25668266740469
# define the initials coordinates of the leading edge 
x_0 = 0.0
x_le = x_0
x_l = 0.15

# define the number of sampling points to be considered 
x = np.zeros(i_Nx+1)
x[0] = 0.0
for i in range(1,i_Nx+1):
    x[i] = x[i-1] + i_x/i_Nx

#print ('the vector x = '), x
# calculate psi values 
psi = np.zeros(i_Nx+1)

psi[0] = 0.0 
for i in range(1,i_Nx+1):
    psi[i] = ( x[i]- x_0 ) / ( x_l )
#print ('The sampling points for the quatinty psi are  =  '), psi
#calculate the psi function for the upper surface 
phi_of_psi_u = np.zeros(i_Nx+1)

for i in range(0,i_Nx+1):
    phi_of_psi_u[i] = x_l * (a * math.sqrt(psi[i]) + b * psi[i] + c * psi[i]*psi[i] + d * psi[i] * psi[i] * psi[i] )
#print ('the function psi is  = '),phi_of_psi_u

# mirror the upper surface 
phi_of_psi_l = -1.0 * phi_of_psi_u   

#print ('The lower surface phi of psi is  = '),phi_of_psi_l




#________________________________________________________


# CREATE THE FLAT PLATE
#________________________________________________________
x_0_plate = 0.15
x_l_plate = 0.9 - x_0_plate

i_x_flat = x_l_plate
x_plate = np.zeros(i_Nx)

x_plate[0] = x[i_Nx] + i_x_flat/i_Nx
for i in range(1,i_Nx):
    x_plate[i] = x_plate[i-1] + i_x_flat/i_Nx
#print x_plate
omega_psi_plate_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
    omega_psi_plate_u[i] = (d + (a * 3/8)) * x_l

omega_psi_plate_l = -1.0 * omega_psi_plate_u    
#print omega_psi_plate_u
#print omega_psi_plate_l
#_______________________________________________________


#CONSTRUCT THE TRAILING  EDGE CLOSURE REGION
#_______________________________________________________
#coefficients fom table A2
a_tr_edge = 0.0
d_tr_edge = 0.17000800036901
b_tr_edge = 3.0 * d_tr_edge
c_tr_edge = -1.0 * b_tr_edge

x_0_tr_edge = 0.9
x_l_tr_edge = 1.0 - x_0_tr_edge
i_x_tr_edge = x_l_tr_edge
   
x_tr_edge = np.zeros(i_Nx)
x_tr_edge[0] = x_plate[i_Nx-1] + i_x_tr_edge/i_Nx
for i in range(1,i_Nx):
    x_tr_edge[i] = x_tr_edge[i-1] + i_x_tr_edge/i_Nx
#print x_tr_edge

psi_tr_edge = np.zeros(i_Nx)
psi_tr_edge[0] = 0.0
for i in range(1,i_Nx):
    psi_tr_edge[i] = (x_tr_edge[i] - x_0_tr_edge)/x_l_tr_edge
#print ('the function psi for the trailing edge is  = ') ,psi_tr_edge
phi_of_psi_tr_edge_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
    phi_of_psi_tr_edge_u[i] = x_l_tr_edge * (a_tr_edge * math.sqrt(psi_tr_edge[i_Nx-1-i]) + b_tr_edge * psi_tr_edge[i_Nx-1-i] + c_tr_edge * psi_tr_edge[i_Nx-1-i]*psi_tr_edge[i_Nx-1-i] + d_tr_edge * psi_tr_edge[i_Nx-1-i] * psi_tr_edge[i_Nx-1-i] * psi_tr_edge[i_Nx-1-i] )

#print ('The function phi of psi for the trailing edge upper surface is  = '), phi_of_psi_tr_edge_u

phi_of_psi_tr_edge_l = -1.0 * phi_of_psi_tr_edge_u

#print ('The function phi of psi for the trailing edge lower surface is  = '), phi_of_psi_tr_edge_l

#_______________________________________________________


#CONSTRUCT THE STING FAIRING 

#coefficients for the stig fairing.....table A2

a_sting = 0.10040234847327
b_sting = 0.33279822819157
c_sting = -0.39554969598736
d_sting = 0.13603332984884

x_0_sting = 0.61057051
x_l_sting = 0.9797 - x_0_sting
i_x_sting = x_l_sting

x_sting = np.zeros(i_Nx+1)

x_sting[0] = x_0_sting #+ i_x_sting_f/i_Nx

for i in range(1,i_Nx+1):
        x_sting[i] = x_sting[i-1] + i_x_sting/i_Nx
#print ('The coordinate x for the sting is  = '), x_sting

psi_sting = np.zeros(i_Nx + 1)
for i in range(0,i_Nx+1):
    psi_sting[i] = (x_sting[i]-x_0_sting)/x_l_sting
#print ('the function psi for the sting fairing is  = '), psi_sting 

phi_of_psi_sting_u = np.zeros(i_Nx+1)
phi_of_psi_sting_u[0] = 0.0
phi_of_psi_sting_u[i_Nx] = 0.06412
for i in range(1,i_Nx):
    phi_of_psi_sting_u[i] = x_l_sting * ( a_sting * math.sqrt(psi_sting[i]) + b_sting * psi_sting[i] + c_sting * psi_sting[i] * psi_sting[i] + d_sting * psi_sting[i] * psi_sting[i] * psi_sting[i] )
phi_of_psi_sting_l = (-1.0)* phi_of_psi_sting_u

#print ('the function phi_of_psi_sting is  = '), phi_of_psi_sting_u
#_______________________________________________________

#CONSTRUCT THE FORE FAIRING 

## the profile for the fore fairing is the same at psi=0.06412 unitl the end of the sting from x/c = 0.9797 to x/c =1.758
## therefore it can be assumed that the fore fairing is a straight line. !!!!! Ideal treatment for  hangover 


x_0_fore = 0.9797
x_l_fore = 1.758 - x_0_fore
i_x_fore = x_l_fore
i_Nx_fore = 500
x_fore = np.zeros(i_Nx_fore+1)
x_fore[0] = x_0_fore

for i in range(1,i_Nx_fore+1):
    x_fore[i] = x_fore[i-1] + i_x_fore/i_Nx_fore
#print ('the coordinate x for the fore fairing is  = '), x_fore   

phi_of_psi_fore_u = np.zeros(i_Nx_fore + 1)

for i in range(0,i_Nx_fore+1):
    phi_of_psi_fore_u[i] = 0.06412
phi_of_psi_fore_l = (-1.0) * phi_of_psi_fore_u

#print ('the function psi for fore fairing is  = '), phi_of_psi_fore


# CREATE ANOTHER  PROFILE AND COPY IT ALONG THE SPAN

#_______________create the leading edge at x/c = 0.5 _________________________________________
x_0_cp = 0.5
x_le_cp = x_0
x_l_cp = 0.15
x_cp = np.zeros(i_Nx+1)
i_x_cp = x_l_cp 
x_cp[0] = 0.5
for i in range(1,i_Nx+1):
    x_cp[i] = x_cp[i-1] + i_x_cp/i_Nx

#print ('the vector x_cp = '), x_cp
# calculate psi values 

psi_cp = np.zeros(i_Nx+1)
for i in range(1,i_Nx+1):
    psi_cp[i] = ( x_cp[i]- x_0_cp ) / ( x_l_cp )
#print ('The sampling points for the quatinty psi are  =  '), psi_cp
#calculate the psi function for the upper surface 
phi_of_psi_cp_u = np.zeros(i_Nx+1)

for i in range(0,i_Nx+1):
    phi_of_psi_cp_u[i] = x_l_cp * (a * math.sqrt(psi_cp[i]) + b * psi_cp[i] + c * psi_cp[i]*psi_cp[i] + d * psi_cp[i] * psi_cp[i] * psi_cp[i] )
#print ('the function psi is  = '),phi_of_psi_cp_u

# mirror the upper surface 
phi_of_psi_cp_l = -1.0 * phi_of_psi_cp_u  


# ____________create the flat plate __________

x_0_plate_cp = 0.65
x_l_plate_cp = 0.9 - x_0_plate_cp

i_x_flat_cp = x_l_plate_cp
x_plate_cp = np.zeros(i_Nx)

x_plate_cp[0] = x_cp[i_Nx] + i_x_flat_cp/i_Nx
for i in range(1,i_Nx):
    x_plate_cp[i] = x_plate_cp[i-1] + i_x_flat_cp/i_Nx
    
#print ('the coordinate x for the flat plate is '), x_plate_cp
omega_psi_plate_cp_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
    omega_psi_plate_cp_u[i] = (d + (a * 3/8)) * x_l

omega_psi_plate_cp_l = -1.0 * omega_psi_plate_cp_u   
#print ('The function omega of psi for the flat plate is  = '), omega_psi_plate_u
#print omega_psi_plate_l
#___________________________________________________

# ____________ create the trailing edge_____________

#coefficients fom table A2
a_tr_edge_cp = 0.0
d_tr_edge_cp = 0.17000800036901
b_tr_edge_cp = 3.0 * d_tr_edge
c_tr_edge_cp = -1.0 * b_tr_edge

x_0_tr_edge_cp = 0.9
x_l_tr_edge_cp = 1.0 - x_0_tr_edge
i_x_tr_edge_cp = x_l_tr_edge
   
x_tr_edge_cp = np.zeros(i_Nx)
x_tr_edge_cp[0] = x_plate_cp[i_Nx-1] + i_x_tr_edge_cp/i_Nx
for i in range(1,i_Nx):
    x_tr_edge_cp[i] = x_tr_edge_cp[i-1] + i_x_tr_edge_cp/i_Nx
#print ('the x coordinate for the copied is profile is '), x_tr_edge_cp

psi_tr_edge_cp = np.zeros(i_Nx)
psi_tr_edge_cp[0] = 0.0
for i in range(1,i_Nx):
    psi_tr_edge_cp[i] = (x_tr_edge_cp[i] - x_0_tr_edge_cp)/x_l_tr_edge_cp
#print ('the function psi for the trailing edge is  = ') ,psi_tr_edge_cp
phi_of_psi_tr_edge_cp_u = np.zeros(i_Nx)
for i in range(0,i_Nx):
    phi_of_psi_tr_edge_cp_u[i] = x_l_tr_edge_cp * (a_tr_edge_cp * math.sqrt(psi_tr_edge_cp[i_Nx-1-i]) + b_tr_edge_cp * psi_tr_edge_cp[i_Nx-1-i] + c_tr_edge_cp * psi_tr_edge_cp[i_Nx-1-i]*psi_tr_edge_cp[i_Nx-1-i] + d_tr_edge_cp * psi_tr_edge_cp[i_Nx-1-i] * psi_tr_edge_cp[i_Nx-1-i] * psi_tr_edge_cp[i_Nx-1-i] )

#print ('The function phi of psi for the trailing edge upper surface is  = '), phi_of_psi_tr_edge_cp_u

phi_of_psi_tr_edge_cp_l = -1.0 * phi_of_psi_tr_edge_cp_u

#print ('The function phi of psi for the trailing edge lower surface is  = '), phi_of_psi_tr_edge_cp_l


#_______________________________________________________

#      WRITE THE DATA TO FILES

#_______________________________________________________write in matrix format for rhinoceros
#----------symmetry plane

# leading edge 
LE = np.zeros(shape=(i_Nx+1,3))

for i in range(0,i_Nx+1):
        LE[i,0] = x[i] 
        LE[i,1] = phi_of_psi_u[i]
        LE[i,2] = 0.0
#print LE_u
#print x
# FLAT PLATE 
FLAT_PLATE = np.zeros(shape=(i_Nx,3))
for i in range(0,i_Nx):
    FLAT_PLATE[i,0] = x_plate[i]
    FLAT_PLATE[i,1] = omega_psi_plate_u[i]
    FLAT_PLATE[i,2] = 0.0
#print FLAT_PLATE

# TRAILING EDGE

TE = np.zeros(shape=(i_Nx,3))
for i in range(0,i_Nx):
    TE[i,0] = x_tr_edge[i]
    TE[i,1] = phi_of_psi_tr_edge_u[i]
    TE[i,2] = 0.0
#print TE

# STING 
STING = np.zeros(shape=(i_Nx+1,3))
for i in range(0,i_Nx+1):
    STING[i,0] = x_sting[i]
    STING[i,1] = phi_of_psi_sting_u[i]
    STING[i,2] = 0.0
#print STING

# FORE 
FORE = np.zeros(shape=(i_Nx_fore+1,3))
for i in range(0,i_Nx_fore+1):
    FORE[i,0] = x_fore[i]
    FORE[i,1] = phi_of_psi_fore_u[i]
    FORE[i,2] = 0.0
#print FORE


# write in matrix format the profile copied at x/c 0.5  and exposed to y/c 0.233

    # leading edge 
    
LE_CP = np.zeros(shape=(i_Nx+1,3))

for i in range(0,i_Nx+1):
        LE_CP[i,0] = x_cp[i] 
        LE_CP[i,1] = phi_of_psi_cp_u[i]
        LE_CP[i,2] = 0.233 
#print LE_CP
    
    # flat plate 
    
FLAT_PLATE_CP = np.zeros(shape=(i_Nx,3))

for i in range(0,i_Nx):
        FLAT_PLATE_CP[i,0] = x_plate_cp[i] 
        FLAT_PLATE_CP[i,1] = omega_psi_plate_cp_u[i]
        FLAT_PLATE_CP[i,2] = 0.233 
#print FLAT_PLATE_CP  

    # trailing edge 
TE_CP = np.zeros(shape=(i_Nx,3))

for i in range(0,i_Nx):
        TE_CP[i,0] = x_tr_edge_cp[i] 
        TE_CP[i,1] = phi_of_psi_tr_edge_cp_u[i]
        TE_CP[i,2] = 0.233 
#print TE_CP    
    
#-------------------SAVE THE DATA FOR RHINOCEROS----------------------    
np.savetxt('LE_s.txt',LE,delimiter=' ')
np.savetxt('TE_s.txt',TE,delimiter=' ')
np.savetxt('FLAT_PLATE_s.txt',FLAT_PLATE,delimiter=' ')
np.savetxt('FORE_s.txt',FORE,delimiter=' ')
np.savetxt('STING_s.txt',STING,delimiter=' ')
np.savetxt('LE_CP_s.txt',LE_CP,delimiter=' ')
np.savetxt('TE_CP_s.txt',TE_CP,delimiter=' ')
np.savetxt('FLAT_PLATE_CP_s.txt',FLAT_PLATE_CP,delimiter=' ')

#_______________________________________________________

#      PLOT THE POINTS

#_______________________________________________________
#fig = plt.figure()#(figsize=(17, 1.5 , 2),dpi=1200, facecolor='w', edgecolor='k')
fig = plt.figure(figsize=(17, 1.5),dpi=1200, facecolor='w', edgecolor='k')

#plt.plot(x,phi_of_psi_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x,phi_of_psi_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_plate,omega_psi_plate_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_plate,omega_psi_plate_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_tr_edge,phi_of_psi_tr_edge_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_tr_edge,phi_of_psi_tr_edge_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_sting, phi_of_psi_sting_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_sting, phi_of_psi_sting_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_fore,phi_of_psi_fore_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_fore,phi_of_psi_fore_l,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
plt.plot(x_cp,phi_of_psi_cp_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_plate_cp,omega_psi_plate_cp_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)
#plt.plot(x_tr_edge_cp,phi_of_psi_tr_edge_cp_u,'x',linestyle='-',linewidth = '2',marker='o',markersize=5)

plt.show()

#print ('The function phi is ='),phi_of_psi_cp_u
#print ('The function omega is = '),omega_psi_plate_cp_u
#print LE_CP









    



Populating the interactive namespace from numpy and matplotlib