import numpy as np
import trimesh
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import CubicSpline

# Dados da Tabela de Cotas
# baliza: x (mm), z_quilha (mm), z_borda (mm), y_25, y_50, y_75, y_borda
data = [
    # baliza, x, z_quilha, z_borda, y_25, y_50, y_75, y_borda
    (1, 0.0, 845.7, 2008.6, 1465.2, 1640.0, 1686.5, 1691.4),
    (2, 2467.0, 191.6, 1779.2, 1879.5, 2095.2, 2151.0, 2156.6),
    (3, 4933.0, 14.3, 1632.0, 2064.8, 2285.4, 2339.6, 2344.5),
    (4, 7400.0, 0.0, 1585.7, 2171.3, 2394.6, 2447.9, 2452.6),
    (5, 9867.0, 90.9, 1696.8, 1917.4, 2161.0, 2228.8, 2236.4),
    (6, 12333.0, 267.7, 2014.7, 1105.9, 1322.5, 1403.6, 1418.9),
    (7, 14800.0, 2537.1, 2537.1, 0.0, 0.0, 0.0, 0.0)
]

# Converter para metros para modelagem naval padrão
stations_data = []
for row in data:
    b, x, z_k, z_d, y25, y50, y75, y_d = row
    stations_data.append({
        'b': b,
        'x': x / 1000.0,
        'z_k': z_k / 1000.0,
        'z_d': z_d / 1000.0,
        'y_25': y25 / 1000.0,
        'y_50': y50 / 1000.0,
        'y_75': y75 / 1000.0,
        'y_d': y_d / 1000.0
    })

# Gerar malha de curvas das balizas (half-breadths e alturas)
# Para cada baliza, temos 5 pontos verticais fundamentais na seção transversal:
# Ponto 0: Quilha (y = 0, z = z_k)
# Ponto 1: 25% da altura útil (y = y_25, z = z_k + 0.25*(z_d - z_k))
# Ponto 2: 50% da altura útil (y = y_50, z = z_k + 0.50*(z_d - z_k))
# Ponto 3: 75% da altura útil (y = y_75, z = z_k + 0.75*(z_d - z_k))
# Ponto 4: Borda / Convés (y = y_d, z = z_d)

num_z_pts = 15 # interpolação suave ao longo da altura
num_stations = 25 # interpolação suave ao longo do comprimento (x)

# Criar splines para cada parâmetro ao longo de x
x_orig = np.array([s['x'] for s in stations_data])
z_k_orig = np.array([s['z_k'] for s in stations_data])
z_d_orig = np.array([s['z_d'] for s in stations_data])
y_25_orig = np.array([s['y_25'] for s in stations_data])
y_50_orig = np.array([s['y_50'] for s in stations_data])
y_75_orig = np.array([s['y_75'] for s in stations_data])
y_d_orig = np.array([s['y_d'] for s in stations_data])

cs_zk = CubicSpline(x_orig, z_k_orig)
cs_zd = CubicSpline(x_orig, z_d_orig)
cs_y25 = CubicSpline(x_orig, y_25_orig)
cs_y50 = CubicSpline(x_orig, y_50_orig)
cs_y75 = CubicSpline(x_orig, y_75_orig)
cs_yd = CubicSpline(x_orig, y_d_orig)

x_dense = np.linspace(0, 14.800, num_stations)

# Construir grade de pontos 3D para o semi-casco de Boreste (y >= 0) e Bombordo (y <= 0)
grid_x = np.zeros((num_stations, num_z_pts))
grid_y_be = np.zeros((num_stations, num_z_pts))
grid_z = np.zeros((num_stations, num_z_pts))

for i, x in enumerate(x_dense):
    if x >= 14.800:
        zk = z_k_orig[-1]
        zd = z_d_orig[-1]
        y25, y50, y75, yd = 0.0, 0.0, 0.0, 0.0
    else:
        zk = max(0.0, float(cs_zk(x)))
        zd = float(cs_zd(x))
        y25 = max(0.0, float(cs_y25(x)))
        y50 = max(0.0, float(cs_y50(x)))
        y75 = max(0.0, float(cs_y75(x)))
        yd = max(0.0, float(cs_yd(x)))
    
    # Pontos de controle na seção
    z_ctrl = np.array([zk, zk + 0.25*(zd - zk), zk + 0.50*(zd - zk), zk + 0.75*(zd - zk), zd])
    y_ctrl = np.array([0.0, y25, y50, y75, yd])
    
    if zd - zk > 0.001 and (yd > 0 or y50 > 0):
        # Parametrização suave do contorno da baliza
        # usamos spline cúbica monotônica ou interpolação com tangentes naturais
        z_vals = np.linspace(zk, zd, num_z_pts)
        # Interpolação y(z)
        cs_section = CubicSpline(z_ctrl, y_ctrl)
        y_vals = np.maximum(0.0, cs_section(z_vals))
        y_vals[0] = 0.0 # quilha estrita na linha de centro
    else:
        z_vals = np.linspace(zk, zd, num_z_pts)
        y_vals = np.zeros(num_z_pts)

    grid_x[i, :] = x
    grid_y_be[i, :] = y_vals
    grid_z[i, :] = z_vals

# Gerar malha de triângulos para o casco completo (BE + BB + Espelho de Popa)
vertices = []
faces = []

# Indexação da grade BE: station i, height j -> index = i * num_z_pts + j
# Boreste (BE, y >= 0)
v_be_offset = len(vertices)
for i in range(num_stations):
    for j in range(num_z_pts):
        vertices.append([grid_x[i, j], grid_y_be[i, j], grid_z[i, j]])

for i in range(num_stations - 1):
    for j in range(num_z_pts - 1):
        p1 = v_be_offset + i * num_z_pts + j
        p2 = v_be_offset + (i + 1) * num_z_pts + j
        p3 = v_be_offset + (i + 1) * num_z_pts + (j + 1)
        p4 = v_be_offset + i * num_z_pts + (j + 1)
        faces.append([p1, p2, p3])
        faces.append([p1, p3, p4])

# Bombordo (BB, y <= 0)
v_bb_offset = len(vertices)
for i in range(num_stations):
    for j in range(num_z_pts):
        vertices.append([grid_x[i, j], -grid_y_be[i, j], grid_z[i, j]])

for i in range(num_stations - 1):
    for j in range(num_z_pts - 1):
        p1 = v_bb_offset + i * num_z_pts + j
        p2 = v_bb_offset + (i + 1) * num_z_pts + j
        p3 = v_bb_offset + (i + 1) * num_z_pts + (j + 1)
        p4 = v_bb_offset + i * num_z_pts + (j + 1)
        # Orientação normal invertida para BB
        faces.append([p1, p3, p2])
        faces.append([p1, p4, p3])

# Espelho de popa (Transom na estação 0, x = 0)
# Triangulamos a seção em x = 0 fechando a popa
v_transom_center = len(vertices)
vertices.append([0.0, 0.0, grid_z[0, 0]]) # centro inferior da popa
for j in range(num_z_pts - 1):
    # BE transom
    p_be1 = v_be_offset + j
    p_be2 = v_be_offset + (j + 1)
    # BB transom
    p_bb1 = v_bb_offset + j
    p_bb2 = v_bb_offset + (j + 1)
    
    faces.append([v_transom_center, p_be2, p_be1])
    faces.append([v_transom_center, p_bb1, p_bb2])

mesh = trimesh.Trimesh(vertices=np.array(vertices), faces=np.array(faces))

# Exportar STL e OBJ para uso em CAD / 3D
stl_path = '/home/magdalena/eee/grupos/grupo-3/Casco_3D_Saveiro_14_80m.stl'
obj_path = '/home/magdalena/eee/grupos/grupo-3/Casco_3D_Saveiro_14_80m.obj'
mesh.export(stl_path)
mesh.export(obj_path)

# Também copiar para pasta do Rodolfo
mesh.export('/home/magdalena/eee/alunos/rodolfo-g-giordano/Casco_3D_Saveiro_14_80m.stl')
mesh.export('/home/magdalena/eee/alunos/rodolfo-g-giordano/Casco_3D_Saveiro_14_80m.obj')

print("Modelos 3D STL e OBJ exportados com sucesso.")