# paneles

import matplotlib.pyplot as plt
import math as mt
import numpy as np
nombre='NACA0012_perfil.txt'    #define nombre de archivo que contiene los datos

def imprimevector(x):
    n=len(x)-1
    i=0
    while i<=n:
        print ("%0 6.4f" % x[i],end="\n")
        i=i+1
        
def ang_beta(a,b,c,ii,jj):
    #print(ii,jj,a,b,c)
    if ii==jj:
        beta=mt.pi
    else:
        beta= mt.acos((a**2+b**2-c**2)/(2*a*b))
    #if beta>mt.pi/2 and beta!=mt.pi:
        #print ("beta",beta,ii,jj)
    return beta

def largo(a,b,c,d):
  # devuelve el largo del panel
  # a coor x del nodo inicial
  # b coord y del nodo inicial
        return ((a-c)**2+(b-d)**2)**0.5

def ang_th(a,b,c,d):
  # devuelve el angulo th del panel
  # a coor x del nodo inicial
  # b coord y del nodo inicial
    ang=mt.atan2((d-b),(c-a))
    #if ang<0:
        #ang=2*mt.pi+ang
        #print("angulo girado",a,b,c,d)
        
    
    return ang

def imprime(x):
    # imprime una matriz en formato tabular
    # tomando en forma dinamica filas y columnas
    # segun forma de la matriz
    for i in range(len(x)):
        for j in range(len(x[0])):
            print ("%0 6.2f" % x[i][j],end=" ")
        #print (end="\n")
        print (end="\n")
    print("-----------------",end="\n")
        
def graba_csv(x):
    # graba una matriz en formato csv
    # tomando en forma dinamica filas y columnas
    # segun forma de la matriz
    f=open("perfil.txt","w")
    for i in range(len(x)):
        for j in range(len(x[0])):
            #print ("%0 6.2f" % x[i][j],end=";")
            f.write(str(int(x[i][j]*100)/100)+";")
        print (end="\n")
        f.write(chr(10)+chr(13))
    f.close()    

def radio(x1,y1,x2,y2):
    return ((x2-x1)**2+(y2-y1)**2)**0.5

class cpunto2D:
    # esto es solo la definicion de una clase, que en este caso tiene tres atributos
    def __init__(self, x, y):
        self.x=x                  # Coordenadas del punto en unidades de pantalla
        self.y=y              # Coordenadas del punto en unidades de pantalla 3D
        

class cpanel:
    def __init__(self, b, e,l):
        self.b=b              # numero del punto de inicio
        self.e=e              # numero del punto de fin
        self.l=l
        self.th=0
        self.pm_x=0
        self.pm_y=0

panel=[]           # creo una lista vacia las paneles

nodo=[]         # creo una lista vacia con los nodo 

cuerda=100.

f=open(nombre,'r')      # el archivo debe estar en la misma carpeta que el programa
for linea in f.readlines():     #este for desgrana una a una la lista de lineas que lee readlines
    if not linea.find("/")==0:  # me fijo si empieza con /
        dumm=linea.split("     ")                      # separo la linea en una lista de nombre dumm buscando las ocurrencias de "  "
        nodo.append(cpunto2D(float(dumm[0])*cuerda,float(dumm[1])*cuerda))            # al ser cu y top definidas en el cuerpo principalson tomadas como globales
        #nodo.y.append(float(dumm[1])*cuerda)
        #print (x[len(x)-1],y[len(x)-1])
#print (len(x))
f.close()                                   # cierro el archivo
nodo=list(reversed(nodo))                         # doy vuelta las listas porque en el txt viene primero el extrados y luego el intra



# cargo los paneles
print("id    X(i) Y(i)   X(i+1)  Y(i+1) L     FI     COS(FI) SIN(FI)")
for i in range(len(nodo)):
  if i<len(nodo)-1:
    panel.append(cpanel(i,i+1,0))
  else:
    panel.append(cpanel(i,0,0))
  
  panel[i].l=largo(nodo[panel[i].b].x,nodo[panel[i].b].y,nodo[panel[i].e].x,nodo[panel[i].e].y)
  panel[i].th=ang_th(nodo[panel[i].b].x,nodo[panel[i].b].y,nodo[panel[i].e].x,nodo[panel[i].e].y)
  panel[i].pm_x=(nodo[panel[i].e].x+nodo[panel[i].b].x)/2
  panel[i].pm_y=(nodo[panel[i].e].y+nodo[panel[i].b].y)/2
  #print(nodo[panel[i].b].x, panel[i].b].y, nodo[panel[i].e].x, nodo[panel[i].e].y, panel[i].l, panel[i].th)
  print("%2d %6.2f %6.2f %6.2f %6.2f %6.2f %6.3f  %6.3f %6.3f "% (i,nodo[panel[i].b].x, nodo[panel[i].b].y, nodo[panel[i].e].x, nodo[panel[i].e].y, panel[i].l, panel[i].th,mt.cos(panel[i].th),mt.sin(panel[i].th)))

#exit()

# dibujo el perfil
x=[]
y=[]
pm_x=[]
pm_y=[]
for i in range(len(panel)):
  x.append(nodo[panel[i].b].x)
  y.append(nodo[panel[i].b].y)
  pm_x.append(panel[i].pm_x)
  pm_y.append(panel[i].pm_y)
fig = plt.figure()                      # aqui se crea o instancia un objeto contenedor
plt.gcf().canvas.set_window_title('Metodo de los paneles') # defino el titulo de la ventana
plt.xlabel('Eje X')
plt.ylabel('Eje Y')
fig.suptitle("perfil")
ax = fig.add_subplot(111)
ax.plot(x, y, color='lightblue', linewidth=1, label='Curva')
#ax.scatter(x, y, color='darkgreen', marker='^',label='puntos')
ax.scatter(pm_x, pm_y, color='red', marker='+',label='puntos')
leg=ax.legend()
plt.ylim(-50.,50.)
plt.show()
#exit()

# armo una estructura vacia para la matriz A
a=[]
for i in range(len(panel)+1):
    a.append([])
    for j in range(len(panel)+1):
        a[i].append(0.0)
#imprime(a)

# calculo los A para las fuentes- ec (I)
for i in range(len(panel)):
    for j in range(len(panel)):   
        dum1=(radio(panel[i].pm_x,panel[i].pm_y,nodo[panel[j].e].x,nodo[panel[j].e].y))
        dum2=(radio(panel[i].pm_x,panel[i].pm_y,nodo[panel[j].b].x,nodo[panel[j].b].y))
        beta=ang_beta(dum1,dum2,panel[j].l,i,j)
        a[i][j]=1/(2*mt.pi)*(mt.sin(panel[i].th-panel[j].th)*mt.log(dum1/dum2)+mt.cos(panel[i].th-panel[j].th)*beta)

#imprime(a)

#exit()
#calculo las A para los vortices (columna N+1) - ec (II)
#print(N)
ult=len(panel)
for i in range(len(panel)):
    tacho=0.
    for j in range(len(panel)):
        dum1=(radio(panel[i].pm_x,panel[i].pm_y,nodo[panel[j].e].x,nodo[panel[j].e].y))
        dum2=(radio(panel[i].pm_x,panel[i].pm_y,nodo[panel[j].b].x,nodo[panel[j].b].y))
        beta=ang_beta(dum1,dum2,panel[j].l,i,j)
        tacho=tacho+(mt.cos(panel[i].th-panel[j].th)*mt.log(dum1/dum2)-mt.sin(panel[i].th-panel[j].th)*beta)
    a[i][ult]=1/(2*mt.pi)*tacho
  
#imprime(a)
#exit()
# calculo la ultima fila de la matriz - ec(III)
ult=len(panel)
#for k in [0,ult-1]:
for j in range(len(panel)):
    tacho=0.
    for k in [0,ult-1]:
    #for j in range(len(panel)):
       dum1=(radio(panel[k].pm_x,panel[k].pm_y,nodo[panel[j].e].x,nodo[panel[j].e].y))
       dum2=(radio(panel[k].pm_x,panel[k].pm_y,nodo[panel[j].b].x,nodo[panel[j].b].y))
       beta=ang_beta(dum1,dum2,panel[j].l,k,j)
       #print(j,k)
       tacho=tacho+(mt.sin(panel[k].th-panel[j].th)*beta-mt.cos(panel[k].th-panel[j].th)*mt.log(dum1/dum2))
    a[ult][j]=1/(2*mt.pi)*tacho
    
#imprime(a)
#exit()


# calculo la ultima celda (fila y columna) - ec (IV)
tacho=0.
for j in range(len(panel)):
#for k in [0,ult-1]:
    #a[N][N]=0.
    #for j in range(len(panel)):
    for k in [0,ult-1]:    
      dum1=(radio(panel[k].pm_x,panel[k].pm_y,nodo[panel[j].e].x,nodo[panel[j].e].y))
      dum2=(radio(panel[k].pm_x,panel[k].pm_y,nodo[panel[j].b].x,nodo[panel[j].b].y))
      beta=ang_beta(dum1,dum2,panel[j].l,k,j)
      #print(j,k)
      tacho=tacho+(mt.sin(panel[k].th-panel[j].th)*mt.log(dum1/dum2)+mt.cos(panel[k].th-panel[j].th)*beta) 

a[ult][ult]= 1/(2*mt.pi)*tacho


imprime(a)
#exit()


alfa=3./57.
vinf=100.



# aca calculo el vector b y lo adjunto a l amatriz a - ec (V)
for i in range(len(panel)):
    #print(i,a[i][0])
    a[i].append(vinf*mt.sin(panel[i].th-alfa))    # esto es el b[] como columna mas de mariz ampliada


# este es el ultimo elemento de b -ec (VI)    
a[ult].append(-vinf*mt.cos(panel[0].th-alfa)-vinf*mt.cos(panel[ult-1].th-alfa) )   #este es el ultimo elemento de b    
    
#imprime(a)
graba_csv(a)
#exit()


#######
#  este bloque lo puse para verificar que la solucion
#  del sistema por la rutina propuesta coincide con la que da numpy
#  para probarlo poner el perfil NACA2412_perfil:simple
##########
a_np=np.asarray(a)
aa_np=a_np[:ult+1,:ult+1]
b_np=a_np[:,ult+1:ult+2]
print("A")
print(a_np)
print("ARED")
print(aa_np)
print("B")
print(b_np)
c_np=np.linalg.inv(aa_np).dot(b_np)
print("C")
print(c_np)

d=[]
#paso el c_np a lista
for i in range(len(c_np)):
        d.append(c_np[i])


    
# el vector d contiene en los primeros n elementos los sigma y en el n+1 el gamma

print("Solucion del sistema")
imprimevector(d)




# calculo las Vt
# la rotacion es el ultimo elemento de la solucion

gamma=d[ult]
print("gamma",gamma,"ult",ult,"len(panel)",len(panel)) 
v_t=[]
cp=[]
for i in range(len(panel)):
    # la circulacion la saco de la fila correspondiente
    sigma=d[i]
    primer_t=0.0
    segundo_t=0.0
    for j in range(len(panel)):

      dum1=(radio(panel[i].pm_x,panel[i].pm_y,nodo[panel[j].b].x,nodo[panel[j].b].y))
      dum2=(radio(panel[i].pm_x,panel[i].pm_y,nodo[panel[j].e].x,nodo[panel[j].e].y))
      beta=ang_beta(dum1,dum2,panel[j].l,i,j)
      primer_t=primer_t+(sigma/(2*mt.pi))*(mt.sin(panel[i].th-panel[j].th)*beta-mt.cos(panel[i].th-panel[j].th)*mt.log(dum1/dum2))   # ec VII 
        
      segundo_t=segundo_t+(mt.sin(panel[i].th-panel[j].th)*mt.log(dum1/dum2)+mt.cos(panel[i].th-panel[j].th)*beta)        # ec VIII

    if nodo[panel[i].b].x>nodo[panel[i].e].x:
        v_t.append(vinf*mt.cos(panel[i].th-alfa)+primer_t+gamma/(2*mt.pi)*segundo_t) #-
    else:    
        v_t.append(vinf*mt.cos(panel[i].th-alfa)+primer_t+gamma/(2*mt.pi)*segundo_t)
    cp.append(1-v_t[i]**2/vinf**2)
    print(i,"vel",v_t[i],"cp",cp[i])

# dibujo la dist d epresiones 
# para eso genero una nueva lista (pres) en el que sumo coordenada y del perfil

pres=[]
xpm=[]
x_ext=[]
x_int=[]
cp_ext=[]
cp_int=[]
for i in range(len(panel)):
    pres.append(cp[i])
    xpm.append(panel[i].pm_x)
    if nodo[panel[i].b].x>nodo[panel[i].e].x:
        x_int.append(pm_x[i])
        cp_int.append(cp[i])
    else:
        x_ext.append(pm_x[i])
        cp_ext.append(cp[i])
fig = plt.figure()                      # aqui se crea o instancia un objeto contenedor
plt.gcf().canvas.set_window_title('Metodo de los paneles') # defino el titulo de la ventana
plt.xlabel('Eje X')
plt.ylabel('Eje Y')
fig.suptitle("Dist. de presiones")
ax = fig.add_subplot(111)

#ax.plot(xpm, pres, color='red', linewidth=1, label='Presion')
ax.plot(x_int, cp_int, color='blue', linewidth=1, label='Intrados')
ax.plot(x_ext, cp_ext, color='green', linewidth=1, label='Extrados')
leg=ax.legend()
plt.ylim(-2.,2.)
plt.show()

densi=1.29 # kg/m3
cuerda_m=cuerda/100 # en m


sust=densi*cuerda_m*gamma*vinf
print("sustentacion;",sust)