import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
import openpyxl


print("")
print("")
print("DISEÑO DE SECCIONES")
print("")
print("")

#-------------------------------------MATERIALES------------------------------------------------
dens=7850       #kg/m3
fy=355          #MPA
gammas=1.05     #Factor seguridad acero
Eyoung=210000        #MPa
nu=0.3          #Coef Poisson

#---

fyd=fy*10**6/gammas         # Pa (resistencia diseño)
Eyoung=Eyoung*10**6                   # Pa (rigidez)
G=Eyoung/(2*(1+nu))              # Pa, Módulo Cizalladura




#----------------------------------------------------CARGAR DATOS DESDE EXCEL (debe estar en la misma carpeta que el archivo python)
archivo_excel = 'C:/Users/Daniel.p/Documents/Automatizaciones/Propiedades seccion/datos_generales.xlsx'
#base_datos=openpyxl.load_workbook('/datos_generales.xlsx')
base_datos=openpyxl.load_workbook(archivo_excel)



#obtener los nombres de las secciones.
hojas =base_datos.sheetnames
print("Secciones Disponibles")
nsec=np.size(hojas)

#Enumera todas las hojas menos la primera con el nombre de la sección
for i in range(nsec-1):
    k=str(i+1)
    print("     - "+k+". "+hojas[i+1])


print("")
print("")
sec=input("Elige una sección: ")
sec=int(sec)
nombreseccion=hojas[sec]

#Obtener el nº de filas de la sección
filas=base_datos[hojas[sec]].max_row
#print("filas = "+filas)

npuntos=filas-1


#importamos los valores de los puntos
puntos=np.zeros([npuntos,2])
for i in range(npuntos):
    for j in range(2):
        puntos[i,j]=base_datos[hojas[sec]].cell(row=i+2,column=j+1).value

px=np.zeros([npuntos,1])
py=np.zeros([npuntos,1])
for i in range(npuntos):
    px[i]=puntos[i,0]
    py[i]=puntos[i,1]


#--------------------------CÁLCULO DEL ANCHO y CANTO MÁXIMO (segun ejes de referencia)--------------------------------
bmax=np.amax(px)-np.amin(px)
hmax=np.amax(py)-np.amin(py)




#------------------------CALCULO DEL PERÍMETRO----------------------------- ## ERROR DE PERIMETRO EN PIEZAS COMPUESTAS
long_i=np.zeros([npuntos-1,1])
for i in range(npuntos-1):
       long_i[i]=((puntos[i+1,0]-puntos[i,0])**2+(puntos[i+1,1]-puntos[i,1])**2)**(1/2)
p=abs(sum(long_i))
perimetro=p[0]



#--------------------------------CÁLCULO DEL ÁREA-------------------------------------------------------
area_i=np.zeros([npuntos-1,1])
for i in range(npuntos-1):
       area_i[i]=(puntos[i+1,0]-puntos[i,0])*(puntos[i+1,1]+puntos[i,1])/2
a=(sum(area_i))
area=abs(a[0])


#--------------------------PESO POR METRO DE SECCIÓN-----------------------------------------------
peso=area*dens 


#-----------------------------------------CÁLCULO DEL CENTRO DE GRAVEDAD ----------------------------------------------------------
cdg_i=np.zeros([npuntos,2])
for i in range(npuntos-1):
    h1=puntos[i,1]
    h2=puntos[i+1,1]
    b=puntos[i+1,0]-puntos[i,0]
    d=puntos[i,0]

    #Centro de gravedad y_i
    if h1+h2==0:
        cdg_i[i,1]=0
    else:
        cdg_i[i,1]=1/3*(h1*h1+h1*h2+h2*h2)/(h1+h2)

    #centro de gravedad x_i
    if h1+h2==0:
        cdg_i[i,0]=d+b/2
    else:
        cdg_i[i,0]=d+b/3*(h1+2*h2)/(h1+h2)
       

#Cálculo del momento estático 
statico_i=np.zeros([npuntos,2])
for i in range(npuntos-1):
        #momento_estático respecto eje x (para cdg,y)
        statico_i[i,1]=area_i[i,0]*cdg_i[i,1]    
         #momento_estático respecto eje y (para cdg,x)
        statico_i[i,0]=area_i[i,0]*cdg_i[i,0]   


#cdg=sum(statico_i)/sum(area_i)
cdg=sum(statico_i)/a
xg=(cdg[0])
yg=(cdg[1])


#Cálculo de las fibras más alejadas desde el CDG:
v1y=np.amax(py)-yg
v2y=np.amin(py)-yg

v1x=np.amax(px)-xg
v2x=np.amin(px)-xg




# ----------------- CÁLCULO MOMENTOS DE INERCIA ----------------------------------

inercia_i=np.zeros([npuntos,3])
for i in range(npuntos-1):
    h1=puntos[i,1]
    h2=puntos[i+1,1]
    b=puntos[i+1,0]-puntos[i,0]
    d=puntos[i,0]
    xgi=cdg_i[i,0]
    ygi=cdg_i[i,1]
    ai=area_i[i,0] 

    #momento de inercia Ixg
    if h2>=h1 :
        ixcuad_G_local=1/12*b*(h1*h1*h1)+b*h1*(h1/2-ygi)*(h1/2-ygi)
        ixtriang_G_loc=1/36*b*(h2-h1)*(h2-h1)*(h2-h1)+1/2*b*(h2-h1)*((2*h1+h2)/3-ygi)*((2*h1+h2)/3-ygi)
    else:
        ixcuad_G_local=1/12*b*(h2*h2*h2)+b*h2*(h2/2-ygi)*(h2/2-ygi)
        ixtriang_G_loc=1/36*b*(h1-h2)*(h1-h2)*(h1-h2)+1/2*b*(h1-h2)*((2*h2+h1)/3-ygi)*((2*h2+h1)/3-ygi)
    inercia_i[i,0]=ixcuad_G_local+ixtriang_G_loc+ai*(yg-ygi)*(yg-ygi)
   
    #momento de inercia Iyg
    if h2>=h1:
        iycuad=1/12*h1*b*b*b+h1*b*(b/2+d-xgi)*(b/2+d-xgi)
        iytrian=1/36*(h2-h1)*b*b*b+1/2*b*(h2-h1)*(2/3*b+d-xgi)*(2/3*b+d-xgi)
    else:
        iycuad=1/12*h2*b*b*b+h2*b*(b/2+d-xgi)*(b/2+d-xgi)
        iytrian=1/36*(h1-h2)*b*b*b+1/2*b*(h1-h2)*(1/3*b+d-xgi)*(1/3*b+d-xgi)        
    inercia_i[i,1]=iycuad+iytrian+ai*(xg-xgi)*(xg-xgi)

    #producto de inercia pxyg
    if h2>=h1:
        pxygcuadrado=b*h1*(-h1/2+ygi)*(-d-b/2+xgi)
        pxytriangulo=b*b*(h2-h1)*(h2-h1)/72+b*(h2-h1)/2*(-(h2-h1)/3-h1+ygi)*(-d-2/3*b+xgi)
    else:
        pxygcuadrado=b*h2*(-h2/2+ygi)*(-d-b/2+xgi)
        pxytriangulo=-b*b*(h1-h2)*(h1-h2)/72+b*(h1-h2)/2*(-(h1-h2)/3-h2+ygi)*(-d-1/3*b+xgi)       
    inercia_i[i,2]=pxygcuadrado+pxytriangulo+ai*(-xg+xgi)*(-yg+ygi)


ig=sum(inercia_i)
ixg=abs(ig[0])
iyg=abs(ig[1])
pxyg=ig[2]
if a>=0:
    pxyg=pxyg
else:
    pxyg=-pxyg

#Radios de Giro
rx=(ixg/area)**0.5
ry=(iyg/area)**0.5

#EJES PRINCIPALES DE INERCIA
ic= (ixg+iyg)/2
ir= (((ixg-iyg)/2)**2+pxyg**2)**0.5

imax=ic+ir
imin=ic-ir
rmax=(imax/area)**0.5
rmin=(imin/area)**0.5


#obtención de la orientación de los ejes principales de inercia
rest=ixg-iyg

if abs(rest)<10**-12:
    tetha=45
else:
    tetha=0.5*math.atan((pxyg*2)/(ixg-iyg))
    tetha=abs(tetha)*180/math.pi

#Criterio de signos producto de inercia
if pxyg>0:
    if ixg>iyg:
         tetha=-tetha
    else:
        tetha=tetha
else:
    if ixg>iyg:
         tetha=tetha
    else:
        tetha=-tetha


#--------------------------------CALCULO DE LOS PUNTOS ROTADOS RESPECTO DE LOS EJES PRINCIPALES---------------
pxrot=np.zeros([npuntos,1])
pyrot=np.zeros([npuntos,1])
for i in range(npuntos):
    pxrot[i]=puntos[i,0]
    pyrot[i]=puntos[i,1]

#-----------------------------CÁLCULO DEL MÓDULO RESISTENTE ELÁSTICO---------------------------

##SOLAMENTE ESTÁ EN EJES PRINCIPALES
wel1x=abs(ixg/v1y)
wel2x=abs(ixg/v2y)
wel1y=abs(iyg/v1x)
wel2y=abs(iyg/v2x)









#-----------------AXIL ELÁSTICO
nel=area*fyd

#--------------------MOMENTO ELÁSTICO--------------------
melx=min(wel1x,wel2x)*fyd
mely=min(wel1y,wel2y)*fyd











#--------------------------------------DIBUJO DE LA SECCIÓN (colocar antes de imprimir resultados para no modificar las unidades)-------------------
fig, ax=plt.subplots()
#Nombre de la sección y notación de ejes
plt.title(nombreseccion +" (m)",color="SteelBlue",fontsize=16,fontweight="bold")
plt.xlabel('X',fontsize=14,fontweight="bold")
plt.ylabel('Y',fontsize=14,fontweight="bold")
ax.set_aspect("equal")

#Dibujamos los puntos
ax.plot(px,py)
ax.fill(px,py,facecolor="lightblue")
#Dibujamos el centro de gravedad
ax.plot(xg, yg, 'ro',markersize=4, fillstyle='full', markerfacecolor='red')

#Dibujamos los ejes X,Y
tt=min(bmax,hmax)
ax.arrow(xg,yg,tt/3,0,head_width=tt/30, head_length=tt/30, fc='green', ec='green')
ax.arrow(xg,yg,0,tt/3,head_width=tt/30, head_length=tt/30, fc='green', ec='green')

tet=tetha*(math.pi)/180
#Dibujamos los ejes principales
ax.arrow(xg,yg,tt/3*math.cos(tet),tt/3*math.sin(tet),head_width=tt/30, head_length=tt/30, fc='red', ec='red')
ax.arrow(xg,yg,-tt/3*math.sin(tet),tt/3*math.cos(tet),head_width=tt/30, head_length=tt/30, fc='red', ec='red')



#--------------------------------------------------------------------------IMPRIMIR RESULTADOS (siempre al final)---------------------------------------------------


# UNIDADES REDONDEO
redondeo=2

# ------------------------------------------------CAMBIO DE UNIDADES ------------
# Pasar a cm
pot=2

bmax=bmax*10**pot
hmax=hmax*10**pot
perimetro = perimetro*10**pot
xg=xg*10**pot
yg=yg*10**pot
v1y=v1y*10**pot
v2y=v2y*10**pot
v1x=v1x*10**pot
v2x=v2x*10**pot
rx=rx*10**pot
ry=ry*10**pot
rmax=rmax*10**pot
rmin=rmin*10**pot

area=area*10**(pot*2)

wel1x=wel1x*10**(pot*3)
wel2x=wel2x*10**(pot*3)
wel1y=wel1y*10**(pot*3)
wel2y=wel2y*10**(pot*3)

ixg=ixg*10**(pot*4)
iyg=iyg*10**(pot*4)
pxyg=pxyg*10**(pot*4)
imax=imax*10**(pot*4)
imin=imin*10**(pot*4)


# pasar a KN
nel=nel/1000
melx=melx/1000
mely=mely/1000


Eyoung=Eyoung/(10**6)
G=G/(10**6)




print("")
print("")

print(" ___________________________________________________________________________________________________________")
print("")

print("  MATERIAL DE LA SECCIÓN")
print("   - Densidad = "+str(dens)+ " kg/m3")
print("   - σy = "+str(fy)+ " MPa")
print("   - γm = "+str(gammas)+ " [ ]")
print("   - E = "+str(Eyoung)+ " MPa")
print("   - υ = "+str(nu)+ " [ ]")
print("   - G = "+str(round(G,0))+ " MPa")
print(" ___________________________________________________________________________________________________________")

print("  DATOS GEOMÉTRICOS")
print("   - Área = " + str(round(area,redondeo))+ " cm2")
print("   - Perímetro = " + str(round(perimetro,redondeo))+ " cm")
print("")



print("  DIMENSIONES MÁXIMAS")
print("   - b_max = "+str(round(bmax,redondeo))+" cm")
print("   - h_max = "+str(round(hmax,redondeo))+" cm")
print("")


print("  POSICIÓN DEL CDG (respecto ejes de referencia)")
print("   - CDG = [ "+ str(round(xg,redondeo))+ " ; "+str(round(yg,redondeo))+" ] cm")
print("")

print("  FIBRAS MÁS ALEJADAS DESDE G")
print("   - v1y = "+str(round(v1y,redondeo))+" cm")
print("   - v2y = "+str(round(v2y,redondeo))+" cm")
print("   - v1x = "+str(round(v1x,redondeo))+" cm")
print("   - v2x = "+str(round(v2x,redondeo))+" cm")
print("")

print("  TENSOR DE INERCIA EN G RESPECTO DE LOS EJES Xg E Yg")
print("   - IxG = "+str(round(ixg,redondeo))+" cm4")
print("   - IyG = "+str(round(iyg,redondeo))+" cm4")
print("   - PxyG = "+str(round(pxyg,redondeo))+" cm4")
print("   - rx = "+str(round(rx,redondeo))+" cm")
print("   - ry = "+str(round(ry,redondeo))+" cm")
print("")

print("  TENSOR PRINCIPAL DE INERCIA EN G")
print("   - Imax = "+str(round(imax,redondeo))+" cm4")
print("   - Imin = "+str(round(imin,redondeo))+" cm4")
print("   - ϴ (orientación ejes princiaples respecto Xg e Yg) = "+str(round(tetha,redondeo))+" º")
print("   - rmax = "+str(round(rmax,redondeo))+" cm")
print("   - rmin = "+str(round(rmin,redondeo))+" cm")
print("")


print("  MÓDULO RESISTENTE ELÁSTICO")
print("   - wel1x = "+str(round(wel1x,redondeo))+" cm3")
print("   - wel2x = "+str(round(wel2x,redondeo))+" cm3")
print("   - wel1y = "+str(round(wel1y,redondeo))+" cm3")
print("   - wel2y = "+str(round(wel2y,redondeo))+" cm3")
print("")
print(" ___________________________________________________________________________________________________________")



print("   PROPIEDADES MECÁNICAS")
print("   - Peso_unitario = "+str(round(peso,redondeo))+" kg/m")
print("   - Nel = "+str(round(nel,redondeo))+" kN")
print("   - Melx = "+str(round(melx,redondeo))+" kN·m")
print("   - Mely = "+str(round(mely,redondeo))+" kN·m")
print("")
plt.show()