

In [10]:

    
%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 = 50
i_x = 0.15
# define coefficients
a = 0.16329931618554
d = 0.05210142334309
b = 0.03382978289013
c = -0.13589185550609
# 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.500
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)
phi_of_psi_cp_u [0] = 0.0
for i in range(1,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_l.txt',LE,delimiter=' ')
np.savetxt('TE_l.txt',TE,delimiter=' ')
np.savetxt('FLAT_PLATE_l.txt',FLAT_PLATE,delimiter=' ')
np.savetxt('FORE_l.txt',FORE,delimiter=' ')
np.savetxt('STING_l.txt',STING,delimiter=' ')
np.savetxt('LE_CP_l.txt',LE_CP,delimiter=' ')
np.savetxt('TE_CP_l.txt',TE_CP,delimiter=' ')
np.savetxt('FLAT_PLATE_CP_l.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 , 2),dpi=1200, facecolor='w', edgecolor='k')
fig = plt.figure(figsize=(17, 1),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  x_cp
print phi_of_psi_cp_u
print LE_CP
print LE









    



Populating the interactive namespace from numpy and matplotlib






    












    



[ 0.5    0.503  0.506  0.509  0.512  0.515  0.518  0.521  0.524  0.527
  0.53   0.533  0.536  0.539  0.542  0.545  0.548  0.551  0.554  0.557
  0.56   0.563  0.566  0.569  0.572  0.575  0.578  0.581  0.584  0.587
  0.59   0.593  0.596  0.599  0.602  0.605  0.608  0.611  0.614  0.617
  0.62   0.623  0.626  0.629  0.632  0.635  0.638  0.641  0.644  0.647
  0.65 ]
[ 0.          0.0035575   0.00506984  0.00623277  0.00720771  0.00805739
  0.0088142   0.0094975   0.01012006  0.01069085  0.01121652  0.01170215
  0.0121518   0.01256877  0.0129558   0.01331522  0.01364903  0.01395898
  0.01424664  0.01451339  0.01476049  0.0149891   0.01520027  0.015395
  0.01557418  0.0157387   0.01588935  0.01602692  0.01615213  0.01626569
  0.01636827  0.01646053  0.01654309  0.01661657  0.01668154  0.0167386
  0.01678829  0.01683117  0.01686777  0.01689863  0.01692425  0.01694514
  0.01696181  0.01697475  0.01698443  0.01699135  0.01699597  0.01699877
  0.0170002   0.01700073  0.0170008 ]
[[ 0.5         0.          0.233     ]
 [ 0.503       0.0035575   0.233     ]
 [ 0.506       0.00506984  0.233     ]
 [ 0.509       0.00623277  0.233     ]
 [ 0.512       0.00720771  0.233     ]
 [ 0.515       0.00805739  0.233     ]
 [ 0.518       0.0088142   0.233     ]
 [ 0.521       0.0094975   0.233     ]
 [ 0.524       0.01012006  0.233     ]
 [ 0.527       0.01069085  0.233     ]
 [ 0.53        0.01121652  0.233     ]
 [ 0.533       0.01170215  0.233     ]
 [ 0.536       0.0121518   0.233     ]
 [ 0.539       0.01256877  0.233     ]
 [ 0.542       0.0129558   0.233     ]
 [ 0.545       0.01331522  0.233     ]
 [ 0.548       0.01364903  0.233     ]
 [ 0.551       0.01395898  0.233     ]
 [ 0.554       0.01424664  0.233     ]
 [ 0.557       0.01451339  0.233     ]
 [ 0.56        0.01476049  0.233     ]
 [ 0.563       0.0149891   0.233     ]
 [ 0.566       0.01520027  0.233     ]
 [ 0.569       0.015395    0.233     ]
 [ 0.572       0.01557418  0.233     ]
 [ 0.575       0.0157387   0.233     ]
 [ 0.578       0.01588935  0.233     ]
 [ 0.581       0.01602692  0.233     ]
 [ 0.584       0.01615213  0.233     ]
 [ 0.587       0.01626569  0.233     ]
 [ 0.59        0.01636827  0.233     ]
 [ 0.593       0.01646053  0.233     ]
 [ 0.596       0.01654309  0.233     ]
 [ 0.599       0.01661657  0.233     ]
 [ 0.602       0.01668154  0.233     ]
 [ 0.605       0.0167386   0.233     ]
 [ 0.608       0.01678829  0.233     ]
 [ 0.611       0.01683117  0.233     ]
 [ 0.614       0.01686777  0.233     ]
 [ 0.617       0.01689863  0.233     ]
 [ 0.62        0.01692425  0.233     ]
 [ 0.623       0.01694514  0.233     ]
 [ 0.626       0.01696181  0.233     ]
 [ 0.629       0.01697475  0.233     ]
 [ 0.632       0.01698443  0.233     ]
 [ 0.635       0.01699135  0.233     ]
 [ 0.638       0.01699597  0.233     ]
 [ 0.641       0.01699877  0.233     ]
 [ 0.644       0.0170002   0.233     ]
 [ 0.647       0.01700073  0.233     ]
 [ 0.65        0.0170008   0.233     ]]
[[ 0.          0.          0.        ]
 [ 0.003       0.0035575   0.        ]
 [ 0.006       0.00506984  0.        ]
 [ 0.009       0.00623277  0.        ]
 [ 0.012       0.00720771  0.        ]
 [ 0.015       0.00805739  0.        ]
 [ 0.018       0.0088142   0.        ]
 [ 0.021       0.0094975   0.        ]
 [ 0.024       0.01012006  0.        ]
 [ 0.027       0.01069085  0.        ]
 [ 0.03        0.01121652  0.        ]
 [ 0.033       0.01170215  0.        ]
 [ 0.036       0.0121518   0.        ]
 [ 0.039       0.01256877  0.        ]
 [ 0.042       0.0129558   0.        ]
 [ 0.045       0.01331522  0.        ]
 [ 0.048       0.01364903  0.        ]
 [ 0.051       0.01395898  0.        ]
 [ 0.054       0.01424664  0.        ]
 [ 0.057       0.01451339  0.        ]
 [ 0.06        0.01476049  0.        ]
 [ 0.063       0.0149891   0.        ]
 [ 0.066       0.01520027  0.        ]
 [ 0.069       0.015395    0.        ]
 [ 0.072       0.01557418  0.        ]
 [ 0.075       0.0157387   0.        ]
 [ 0.078       0.01588935  0.        ]
 [ 0.081       0.01602692  0.        ]
 [ 0.084       0.01615213  0.        ]
 [ 0.087       0.01626569  0.        ]
 [ 0.09        0.01636827  0.        ]
 [ 0.093       0.01646053  0.        ]
 [ 0.096       0.01654309  0.        ]
 [ 0.099       0.01661657  0.        ]
 [ 0.102       0.01668154  0.        ]
 [ 0.105       0.0167386   0.        ]
 [ 0.108       0.01678829  0.        ]
 [ 0.111       0.01683117  0.        ]
 [ 0.114       0.01686777  0.        ]
 [ 0.117       0.01689863  0.        ]
 [ 0.12        0.01692425  0.        ]
 [ 0.123       0.01694514  0.        ]
 [ 0.126       0.01696181  0.        ]
 [ 0.129       0.01697475  0.        ]
 [ 0.132       0.01698443  0.        ]
 [ 0.135       0.01699135  0.        ]
 [ 0.138       0.01699597  0.        ]
 [ 0.141       0.01699877  0.        ]
 [ 0.144       0.0170002   0.        ]
 [ 0.147       0.01700073  0.        ]
 [ 0.15        0.0170008   0.        ]]