#!/usr/bin/env python3
"""
build_rocket-fighter-retro.py
Builds the 1950s pulp science-fiction rocket fighter (atomic-age raygun gothic)
from an empty Blender 5.2 scene, downloads/loads CC0 Poly Haven textures,
bakes PBR textures (base color, normal, ORM) with Cycles, verifies moving pivots,
exports rocket-fighter-retro.glb and nodes.json, and renders final.png, side.png,
clay.png, and paints.png.
"""

import os
import sys
import json
import math
import struct
import urllib.request
import bpy
import bmesh
from mathutils import Vector, Matrix, Euler

# ==============================================================================
# CONFIGURATION & PATHS
# ==============================================================================
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
TEX_DIR = os.path.join(SCRIPT_DIR, "textures")
os.makedirs(TEX_DIR, exist_ok=True)

GLB_PATH = os.path.join(SCRIPT_DIR, "rocket-fighter-retro.glb")
NODES_JSON_PATH = os.path.join(SCRIPT_DIR, "nodes.json")
FINAL_PNG_PATH = os.path.join(SCRIPT_DIR, "final.png")
SIDE_PNG_PATH = os.path.join(SCRIPT_DIR, "side.png")
CLAY_PNG_PATH = os.path.join(SCRIPT_DIR, "clay.png")
PAINTS_PNG_PATH = os.path.join(SCRIPT_DIR, "paints.png")
PIVOT_TEST_PNG_PATH = "/tmp/pivot_test_rocket_fighter.png"

POLYHAVEN_IDS = [
    "metal_plate_02",
    "metal_plate",
    "rusty_metal_03",
    "fabric_leather_02",
]

PAINTS = [
    {"name": "Cherry Red", "linear_rgb": [0.65, 0.05, 0.04], "rgb": [0.65, 0.05, 0.04]},
    {"name": "Atomic Mint", "linear_rgb": [0.14, 0.58, 0.42], "rgb": [0.14, 0.58, 0.42]},
    {"name": "Vintage Cream", "linear_rgb": [0.82, 0.74, 0.56], "rgb": [0.82, 0.74, 0.56]},
    {"name": "Stratosphere Sky Blue", "linear_rgb": [0.10, 0.40, 0.78], "rgb": [0.10, 0.40, 0.78]},
    {"name": "Solar Tangerine", "linear_rgb": [0.84, 0.22, 0.03], "rgb": [0.84, 0.22, 0.03]},
    {"name": "Obsidian Black", "linear_rgb": [0.04, 0.045, 0.055], "rgb": [0.04, 0.045, 0.055]},
]

MOVING_AXES = {
    "Canopy": [1.0, 0.0, 0.0],
    "Turret_1": [0.0, 1.0, 0.0],
    "Gun_1": [1.0, 0.0, 0.0],
    "Turret_2": [0.0, 1.0, 0.0],
    "Gun_2": [1.0, 0.0, 0.0],
}


# ==============================================================================
# 1. CLEAN SCENE & SETUP CYCLES GPU
# ==============================================================================
def reset_scene():
    bpy.ops.wm.read_factory_settings(use_empty=True)
    scene = bpy.context.scene
    scene.unit_settings.system = 'METRIC'
    scene.unit_settings.scale_length = 1.0

    scene.render.engine = 'CYCLES'
    prefs = bpy.context.preferences
    cprefs = prefs.addons['cycles'].preferences
    cprefs.refresh_devices()
    used_gpu = False
    for dev_type in ('OPTIX', 'CUDA'):
        try:
            cprefs.compute_device_type = dev_type
            for d in cprefs.devices:
                d.use = True
            scene.cycles.device = 'GPU'
            used_gpu = True
            print(f"[Cycles] Enabled GPU device type: {dev_type}")
            break
        except Exception:
            pass
    if not used_gpu:
        scene.cycles.device = 'CPU'
        print("[Cycles] Fallback to CPU")

    scene.cycles.samples = 64
    scene.cycles.use_denoising = True
    return scene


# ==============================================================================
# 2. DOWNLOAD / VERIFY POLY HAVEN TEXTURES
# ==============================================================================
def ensure_polyhaven_textures():
    maps = ["diff", "nor_gl", "arm"]
    paths = {}
    for ph_id in POLYHAVEN_IDS:
        paths[ph_id] = {}
        for m in maps:
            fname = f"{ph_id}_{m}_1k.jpg"
            fpath = os.path.join(TEX_DIR, fname)
            if not os.path.exists(fpath) or os.path.getsize(fpath) < 1000:
                url = f"https://dl.polyhaven.org/file/ph-assets/Textures/jpg/1k/{ph_id}/{fname}"
                print(f"[PolyHaven] Downloading {url} -> {fpath}")
                req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"})
                with urllib.request.urlopen(req, timeout=30) as resp, open(fpath, "wb") as out:
                    out.write(resp.read())
            paths[ph_id][m] = fpath
    return paths


# ==============================================================================
# 3. CYCLES MATERIAL BAKING (POLY HAVEN + PROCEDURAL WEATHERING -> IMAGE MAPS)
# ==============================================================================
def load_image(filepath, colorspace='sRGB'):
    img = bpy.data.images.load(filepath, check_existing=True)
    img.colorspace_settings.name = colorspace
    return img


def build_procedural_bake_mat(name, ph_paths, mode):
    """
    Creates a Cycles node material mixing Poly Haven maps + procedural
    panel lines, rivets, grime, scorch, and edge wear for baking onto a 2x2m proxy tile.
    """
    mat = bpy.data.materials.new(name=f"_BakeSrc_{name}")
    nodes = mat.node_tree.nodes
    links = mat.node_tree.links
    nodes.clear()

    out = nodes.new('ShaderNodeOutputMaterial')
    bsdf = nodes.new('ShaderNodeBsdfPrincipled')
    links.new(bsdf.outputs['BSDF'], out.inputs['Surface'])

    tex_coord = nodes.new('ShaderNodeTexCoord')
    mapping = nodes.new('ShaderNodeMapping')
    mapping.inputs['Scale'].default_value = (2.0, 2.0, 2.0)
    links.new(tex_coord.outputs['UV'], mapping.inputs['Vector'])

    if mode in ('Paint', 'Chrome', 'Enamel_White'):
        ph_main = ph_paths['metal_plate_02']
        ph_wear = ph_paths['rusty_metal_03']
    elif mode == 'Scorched_Metal':
        ph_main = ph_paths['metal_plate']
        ph_wear = ph_paths['rusty_metal_03']
    elif mode == 'Cockpit':
        ph_main = ph_paths['fabric_leather_02']
        ph_wear = ph_paths['metal_plate_02']
    else:
        ph_main = ph_paths['metal_plate_02']
        ph_wear = ph_paths['metal_plate']

    n_diff = nodes.new('ShaderNodeTexImage')
    n_diff.image = load_image(ph_main['diff'], 'sRGB')
    links.new(mapping.outputs['Vector'], n_diff.inputs['Vector'])

    n_arm = nodes.new('ShaderNodeTexImage')
    n_arm.image = load_image(ph_main['arm'], 'Non-Color')
    links.new(mapping.outputs['Vector'], n_arm.inputs['Vector'])

    n_nor = nodes.new('ShaderNodeTexImage')
    n_nor.image = load_image(ph_main['nor_gl'], 'Non-Color')
    links.new(mapping.outputs['Vector'], n_nor.inputs['Vector'])

    n_wear_diff = nodes.new('ShaderNodeTexImage')
    n_wear_diff.image = load_image(ph_wear['diff'], 'sRGB')
    links.new(mapping.outputs['Vector'], n_wear_diff.inputs['Vector'])

    brick = nodes.new('ShaderNodeTexBrick')
    brick.offset = 0.5
    brick.inputs['Scale'].default_value = 5.0
    brick.inputs['Mortar Size'].default_value = 0.008
    brick.inputs['Mortar Smooth'].default_value = 0.35
    brick.inputs['Color1'].default_value = (0.98, 0.98, 0.99, 1.0)
    brick.inputs['Color2'].default_value = (0.94, 0.94, 0.96, 1.0)
    brick.inputs['Mortar'].default_value = (0.76, 0.76, 0.78, 1.0)
    links.new(mapping.outputs['Vector'], brick.inputs['Vector'])

    sep_arm = nodes.new('ShaderNodeSeparateColor')
    links.new(n_arm.outputs['Color'], sep_arm.inputs['Color'])

    if mode == 'Paint':
        # PAINT CONTRACT: Neutral near-white texture with panel lines, grey grime,
        # engine scorch, micro-meteorite pitting and chipped edges baked in
        mix_base = nodes.new('ShaderNodeMix')
        mix_base.data_type = 'RGBA'
        mix_base.blend_type = 'MULTIPLY'
        mix_base.inputs['Factor'].default_value = 0.22
        links.new(brick.outputs['Color'], mix_base.inputs['A'])
        links.new(n_diff.outputs['Color'], mix_base.inputs['B'])

        bright = nodes.new('ShaderNodeBrightContrast')
        bright.inputs['Bright'].default_value = 0.20
        bright.inputs['Contrast'].default_value = -0.12
        links.new(mix_base.outputs['Result'], bright.inputs['Color'])

        grime_mix = nodes.new('ShaderNodeMix')
        grime_mix.data_type = 'RGBA'
        grime_mix.blend_type = 'MIX'
        grime_mix.inputs['Factor'].default_value = 0.08
        links.new(bright.outputs['Color'], grime_mix.inputs['A'])
        links.new(n_wear_diff.outputs['Color'], grime_mix.inputs['B'])

        rough_map = nodes.new('ShaderNodeMapRange')
        rough_map.inputs['To Min'].default_value = 0.14
        rough_map.inputs['To Max'].default_value = 0.30
        links.new(sep_arm.outputs['Green'], rough_map.inputs['Value'])

        base_color_socket = grime_mix.outputs['Result']
        rough_socket = rough_map.outputs['Result']
        metal_val = 0.15

    elif mode == 'Chrome':
        mix_chrome = nodes.new('ShaderNodeMix')
        mix_chrome.data_type = 'RGBA'
        mix_chrome.blend_type = 'MULTIPLY'
        mix_chrome.inputs['Factor'].default_value = 0.25
        links.new(brick.outputs['Color'], mix_chrome.inputs['A'])
        links.new(n_diff.outputs['Color'], mix_chrome.inputs['B'])

        bright = nodes.new('ShaderNodeBrightContrast')
        bright.inputs['Bright'].default_value = 0.25
        bright.inputs['Contrast'].default_value = 0.05
        links.new(mix_chrome.outputs['Result'], bright.inputs['Color'])

        rough_map = nodes.new('ShaderNodeMapRange')
        rough_map.inputs['To Min'].default_value = 0.10
        rough_map.inputs['To Max'].default_value = 0.22
        links.new(sep_arm.outputs['Green'], rough_map.inputs['Value'])

        base_color_socket = bright.outputs['Color']
        rough_socket = rough_map.outputs['Result']
        metal_val = 0.96

    elif mode == 'Enamel_White':
        mix_w = nodes.new('ShaderNodeMix')
        mix_w.data_type = 'RGBA'
        mix_w.blend_type = 'MIX'
        mix_w.inputs['Factor'].default_value = 0.88
        links.new(n_diff.outputs['Color'], mix_w.inputs['A'])
        mix_w.inputs['B'].default_value = (0.95, 0.93, 0.89, 1.0)

        rough_map = nodes.new('ShaderNodeMapRange')
        rough_map.inputs['To Min'].default_value = 0.15
        rough_map.inputs['To Max'].default_value = 0.28
        links.new(sep_arm.outputs['Green'], rough_map.inputs['Value'])

        base_color_socket = mix_w.outputs['Result']
        rough_socket = rough_map.outputs['Result']
        metal_val = 0.08

    elif mode == 'Scorched_Metal':
        mix_s = nodes.new('ShaderNodeMix')
        mix_s.data_type = 'RGBA'
        mix_s.blend_type = 'MULTIPLY'
        mix_s.inputs['Factor'].default_value = 0.65
        links.new(n_diff.outputs['Color'], mix_s.inputs['A'])
        links.new(n_wear_diff.outputs['Color'], mix_s.inputs['B'])

        tint = nodes.new('ShaderNodeMix')
        tint.data_type = 'RGBA'
        tint.blend_type = 'MULTIPLY'
        tint.inputs['Factor'].default_value = 0.62
        links.new(mix_s.outputs['Result'], tint.inputs['A'])
        tint.inputs['B'].default_value = (0.48, 0.43, 0.38, 1.0)

        rough_map = nodes.new('ShaderNodeMapRange')
        rough_map.inputs['To Min'].default_value = 0.22
        rough_map.inputs['To Max'].default_value = 0.44
        links.new(sep_arm.outputs['Green'], rough_map.inputs['Value'])

        base_color_socket = tint.outputs['Result']
        rough_socket = rough_map.outputs['Result']
        metal_val = 0.82

    elif mode == 'Cockpit':
        tint = nodes.new('ShaderNodeMix')
        tint.data_type = 'RGBA'
        tint.blend_type = 'MULTIPLY'
        tint.inputs['Factor'].default_value = 0.65
        links.new(n_diff.outputs['Color'], tint.inputs['A'])
        tint.inputs['B'].default_value = (0.24, 0.21, 0.19, 1.0)

        rough_map = nodes.new('ShaderNodeMapRange')
        rough_map.inputs['To Min'].default_value = 0.42
        rough_map.inputs['To Max'].default_value = 0.68
        links.new(sep_arm.outputs['Green'], rough_map.inputs['Value'])

        base_color_socket = tint.outputs['Result']
        rough_socket = rough_map.outputs['Result']
        metal_val = 0.18

    links.new(base_color_socket, bsdf.inputs['Base Color'])
    links.new(rough_socket, bsdf.inputs['Roughness'])
    bsdf.inputs['Metallic'].default_value = metal_val

    nor_map = nodes.new('ShaderNodeNormalMap')
    nor_map.inputs['Strength'].default_value = 0.70
    links.new(n_nor.outputs['Color'], nor_map.inputs['Color'])

    bump = nodes.new('ShaderNodeBump')
    bump.inputs['Strength'].default_value = 0.12
    bump.inputs['Distance'].default_value = 0.004
    links.new(brick.outputs['Fac'], bump.inputs['Height'])
    links.new(nor_map.outputs['Normal'], bump.inputs['Normal'])
    links.new(bump.outputs['Normal'], bsdf.inputs['Normal'])

    comb_orm = nodes.new('ShaderNodeCombineColor')
    links.new(sep_arm.outputs['Red'], comb_orm.inputs['Red'])
    links.new(rough_socket, comb_orm.inputs['Green'])
    comb_orm.inputs['Blue'].default_value = metal_val

    return mat, base_color_socket, comb_orm.outputs['Color'], out


def bake_all_materials(scene, ph_paths):
    print("[Bake] Starting Cycles GPU texture baking...")
    orig_samples = scene.cycles.samples
    orig_denoise = scene.cycles.use_denoising
    scene.cycles.samples = 2
    scene.cycles.use_denoising = False

    bpy.ops.mesh.primitive_plane_add(size=2.0, location=(0, 0, -50.0))
    plane = bpy.context.active_object
    plane.name = "_BakeProxyPlane"

    baked_mats = {}
    modes = ['Paint', 'Chrome', 'Enamel_White', 'Scorched_Metal', 'Cockpit']

    for mode in modes:
        b_path = os.path.join(TEX_DIR, f"{mode}_basecolor.jpg")
        n_path = os.path.join(TEX_DIR, f"{mode}_normal.png")
        o_path = os.path.join(TEX_DIR, f"{mode}_orm.jpg")

        src_mat, base_soc, orm_soc, out_node = build_procedural_bake_mat(mode, ph_paths, mode)
        plane.data.materials.clear()
        plane.data.materials.append(src_mat)
        nodes = src_mat.node_tree.nodes
        links = src_mat.node_tree.links
        bsdf = next(n for n in nodes if n.type == 'BSDF_PRINCIPLED')

        bpy.ops.object.select_all(action='DESELECT')
        plane.select_set(True)
        bpy.context.view_layer.objects.active = plane

        # 1) Bake BaseColor (1024x1024)
        img_base = bpy.data.images.new(f"{mode}_basecolor", width=1024, height=1024, alpha=False)
        img_base.colorspace_settings.name = 'sRGB'
        emit_node = nodes.new('ShaderNodeEmission')
        links.new(base_soc, emit_node.inputs['Color'])
        links.new(emit_node.outputs['Emission'], out_node.inputs['Surface'])

        target_node = nodes.new('ShaderNodeTexImage')
        target_node.image = img_base
        for n in nodes:
            n.select = False
        target_node.select = True
        nodes.active = target_node

        bpy.ops.object.bake(type='EMIT')
        img_base.filepath_raw = b_path
        img_base.file_format = 'JPEG'
        img_base.save()

        # 2) Bake ORM (512x512)
        img_orm = bpy.data.images.new(f"{mode}_orm", width=512, height=512, alpha=False)
        img_orm.colorspace_settings.name = 'Non-Color'
        links.new(orm_soc, emit_node.inputs['Color'])
        target_node.image = img_orm
        for n in nodes:
            n.select = False
        target_node.select = True
        nodes.active = target_node

        bpy.ops.object.bake(type='EMIT')
        img_orm.filepath_raw = o_path
        img_orm.file_format = 'JPEG'
        img_orm.save()

        # 3) Bake Normal (1024x1024)
        img_nor = bpy.data.images.new(f"{mode}_normal", width=1024, height=1024, alpha=False)
        img_nor.colorspace_settings.name = 'Non-Color'
        links.new(bsdf.outputs['BSDF'], out_node.inputs['Surface'])
        target_node.image = img_nor
        for n in nodes:
            n.select = False
        target_node.select = True
        nodes.active = target_node

        scene.render.bake.normal_space = 'TANGENT'
        bpy.ops.object.bake(type='NORMAL')
        img_nor.filepath_raw = n_path
        img_nor.file_format = 'PNG'
        img_nor.save()

        bpy.data.materials.remove(src_mat)

        # Build final glTF-ready PBR material
        final_mat = bpy.data.materials.new(name=mode)
        fnodes = final_mat.node_tree.nodes
        flinks = final_mat.node_tree.links
        fnodes.clear()

        fout = fnodes.new('ShaderNodeOutputMaterial')
        fbsdf = fnodes.new('ShaderNodeBsdfPrincipled')
        flinks.new(fbsdf.outputs['BSDF'], fout.inputs['Surface'])

        tex_b = fnodes.new('ShaderNodeTexImage')
        tex_b.image = load_image(b_path, 'sRGB')
        tex_b.name = "BaseColor_Tex"

        tex_o = fnodes.new('ShaderNodeTexImage')
        tex_o.image = load_image(o_path, 'Non-Color')
        tex_o.name = "ORM_Tex"

        tex_n = fnodes.new('ShaderNodeTexImage')
        tex_n.image = load_image(n_path, 'Non-Color')
        tex_n.name = "Normal_Tex"

        if mode == 'Paint':
            mix_factor = fnodes.new('ShaderNodeMix')
            mix_factor.name = "Paint_Factor_Mix"
            mix_factor.data_type = 'RGBA'
            mix_factor.blend_type = 'MULTIPLY'
            mix_factor.inputs['Factor'].default_value = 1.0
            flinks.new(tex_b.outputs['Color'], mix_factor.inputs['A'])
            r, g, b = PAINTS[0]['linear_rgb']
            mix_factor.inputs['B'].default_value = (r, g, b, 1.0)
            flinks.new(mix_factor.outputs['Result'], fbsdf.inputs['Base Color'])
            if 'Coat Weight' in fbsdf.inputs:
                fbsdf.inputs['Coat Weight'].default_value = 0.45
                fbsdf.inputs['Coat Roughness'].default_value = 0.06
        else:
            flinks.new(tex_b.outputs['Color'], fbsdf.inputs['Base Color'])

        sep_orm = fnodes.new('ShaderNodeSeparateColor')
        flinks.new(tex_o.outputs['Color'], sep_orm.inputs['Color'])
        flinks.new(sep_orm.outputs['Green'], fbsdf.inputs['Roughness'])
        flinks.new(sep_orm.outputs['Blue'], fbsdf.inputs['Metallic'])

        f_normap = fnodes.new('ShaderNodeNormalMap')
        f_normap.inputs['Strength'].default_value = 0.85
        flinks.new(tex_n.outputs['Color'], f_normap.inputs['Color'])
        flinks.new(f_normap.outputs['Normal'], fbsdf.inputs['Normal'])

        baked_mats[mode] = final_mat
        print(f"[Bake] Finished baking material: {mode}")

    bpy.data.objects.remove(plane, do_unlink=True)

    # Glass material (alpha BLEND, tinted)
    mat_glass = bpy.data.materials.new(name="Glass")
    if hasattr(mat_glass, 'blend_method'):
        mat_glass.blend_method = 'BLEND'
    if hasattr(mat_glass, 'surface_render_method'):
        mat_glass.surface_render_method = 'BLENDED'
    gnodes = mat_glass.node_tree.nodes
    gbsdf = next(n for n in gnodes if n.type == 'BSDF_PRINCIPLED')
    gbsdf.inputs['Base Color'].default_value = (0.62, 0.80, 0.90, 1.0)
    gbsdf.inputs['Metallic'].default_value = 0.08
    gbsdf.inputs['Roughness'].default_value = 0.04
    gbsdf.inputs['IOR'].default_value = 1.48
    gbsdf.inputs['Alpha'].default_value = 0.30
    baked_mats['Glass'] = mat_glass

    # Engines material (warm orange-amber emissive nozzle interior)
    mat_eng = bpy.data.materials.new(name="Engines")
    enodes = mat_eng.node_tree.nodes
    ebsdf = next(n for n in enodes if n.type == 'BSDF_PRINCIPLED')
    ebsdf.inputs['Base Color'].default_value = (1.0, 0.40, 0.05, 1.0)
    ebsdf.inputs['Roughness'].default_value = 0.25
    ebsdf.inputs['Emission Color'].default_value = (1.0, 0.44, 0.06, 1.0)
    ebsdf.inputs['Emission Strength'].default_value = 10.0
    baked_mats['Engines'] = mat_eng

    # Lamps material (emissive running lights, ray-gun coils, cockpit dials)
    mat_lamps = bpy.data.materials.new(name="Lamps")
    lnodes = mat_lamps.node_tree.nodes
    lbsdf = next(n for n in lnodes if n.type == 'BSDF_PRINCIPLED')
    lbsdf.inputs['Base Color'].default_value = (0.40, 0.90, 1.0, 1.0)
    lbsdf.inputs['Roughness'].default_value = 0.15
    lbsdf.inputs['Emission Color'].default_value = (0.40, 0.90, 1.0, 1.0)
    lbsdf.inputs['Emission Strength'].default_value = 4.5
    baked_mats['Lamps'] = mat_lamps

    scene.cycles.samples = orig_samples
    scene.cycles.use_denoising = orig_denoise
    return baked_mats


# ==============================================================================
# 4. GEOMETRY HELPERS & FUSELAGE PROFILE MATH
# ==============================================================================
def smooth_interp(z, table):
    if z <= table[0][0]:
        return table[0][1]
    if z >= table[-1][0]:
        return table[-1][1]
    for i in range(len(table) - 1):
        z0, v0 = table[i]
        z1, v1 = table[i + 1]
        if z0 <= z <= z1:
            t = (z - z0) / max(1e-6, (z1 - z0))
            st = 0.5 - 0.5 * math.cos(math.pi * t)
            return v0 + (v1 - v0) * (0.40 * t + 0.60 * st)
    return table[-1][1]


def linear_interp(val, table):
    if val <= table[0][0]:
        return table[0][1]
    if val >= table[-1][0]:
        return table[-1][1]
    for i in range(len(table) - 1):
        x0, v0 = table[i]
        x1, v1 = table[i + 1]
        if x0 <= val <= x1:
            t = (val - x0) / max(1e-6, (x1 - x0))
            return v0 + (v1 - v0) * t
    return table[-1][1]


# Fuselage radius R(z) and vertical center Yc(z) from Z = -4.60 to +4.42
FUSELAGE_R_TABLE = [
    (-4.60, 0.45),
    (-4.30, 0.50),
    (-3.65, 0.59),
    (-2.85, 0.68),
    (-2.05, 0.74),
    (-1.25, 0.78),
    (-0.40, 0.80),
    ( 0.50, 0.79),
    ( 1.35, 0.75),
    ( 2.15, 0.66),
    ( 2.85, 0.52),
    ( 3.45, 0.37),
    ( 3.88, 0.23),
    ( 4.22, 0.10),
    ( 4.38, 0.03),
    ( 4.42, 0.008),
]

FUSELAGE_YC_TABLE = [
    (-4.60,  0.00),
    (-3.65,  0.01),
    (-2.05,  0.02),
    (-0.40,  0.03),
    ( 0.85,  0.02),
    ( 2.15,  0.00),
    ( 3.45, -0.01),
    ( 4.42, -0.01),
]


def get_hull_r_yc(z):
    return smooth_interp(z, FUSELAGE_R_TABLE), smooth_interp(z, FUSELAGE_YC_TABLE)


def hull_surface_point(z, angle_rad, radial_offset=0.0):
    r, yc = get_hull_r_yc(z)
    rx = r * 1.02
    ry = r * 0.98
    x = (rx + radial_offset) * math.cos(angle_rad)
    y = yc + (ry + radial_offset) * math.sin(angle_rad)
    if -2.5 < z < 1.2 and math.sin(angle_rad) < -0.5:
        y *= 0.96
    return Vector((x, y, z))


def ensure_mat_slot(obj, mat):
    for i, m in enumerate(obj.data.materials):
        if m == mat:
            return i
    obj.data.materials.append(mat)
    return len(obj.data.materials) - 1


def create_prim(bm, op_func, mat_idx, smooth=True, **kwargs):
    """
    Runs a bmesh.ops.create_* operator and reliably assigns `mat_idx` and `smooth`
    to all newly created faces (avoiding BMFaceSeq negative-slice issues).
    """
    before = set(bm.faces)
    op_func(bm, **kwargs)
    new_faces = [f for f in bm.faces if f not in before]
    for f in new_faces:
        f.material_index = mat_idx
        f.smooth = smooth
    return new_faces


def add_rivet_dome(bm, center, normal, radius=0.0075, height=0.004, mat_idx=0, segs=6):
    n = normal.normalized()
    up = Vector((0, 0, 1)) if abs(n.z) < 0.9 else Vector((1, 0, 0))
    t1 = n.cross(up).normalized()
    t2 = n.cross(t1).normalized()

    base_verts = []
    mid_verts = []
    for i in range(segs):
        ang = 2.0 * math.pi * i / segs
        ca, sa = math.cos(ang), math.sin(ang)
        dir_v = t1 * ca + t2 * sa
        v_b = bm.verts.new(center - n * 0.001 + dir_v * radius)
        v_m = bm.verts.new(center + n * (height * 0.65) + dir_v * (radius * 0.68))
        base_verts.append(v_b)
        mid_verts.append(v_m)
    v_top = bm.verts.new(center + n * height)

    for i in range(segs):
        ni = (i + 1) % segs
        f1 = bm.faces.new((base_verts[i], base_verts[ni], mid_verts[ni], mid_verts[i]))
        f1.material_index = mat_idx
        f1.smooth = True
        f2 = bm.faces.new((mid_verts[i], mid_verts[ni], v_top))
        f2.material_index = mat_idx
        f2.smooth = True


# ==============================================================================
# 5. BUILD MAIN HULL MESH (FUSELAGE, FAIRINGS, WINGS, FINS, RING FIN, NOZZLE, COCKPIT)
# ==============================================================================
def build_hull_object(mats):
    mesh = bpy.data.meshes.new("Hull_Mesh")
    obj = bpy.data.objects.new("Hull", mesh)
    bpy.context.collection.objects.link(obj)

    idx_paint = ensure_mat_slot(obj, mats['Paint'])
    idx_chrome = ensure_mat_slot(obj, mats['Chrome'])
    idx_white = ensure_mat_slot(obj, mats['Enamel_White'])
    idx_scorch = ensure_mat_slot(obj, mats['Scorched_Metal'])
    idx_cockpit = ensure_mat_slot(obj, mats['Cockpit'])
    idx_engines = ensure_mat_slot(obj, mats['Engines'])
    idx_lamps = ensure_mat_slot(obj, mats['Lamps'])

    bm = bmesh.new()

    # --------------------------------------------------------------------------
    # A. MIRROR-SMOOTH CIGAR FUSELAGE WITH SUPPORT-LOOP PANEL SEAMS
    # --------------------------------------------------------------------------
    # Notice: we put the fuselage panel seam grooves on the bare-chrome body
    # (-4.30 to +2.05) and at +3.85 (chrome nose tip), so the painted/white nose
    # cone (2.05 to 3.85) is a single smooth aerodynamic ogive without groove dips!
    base_stations = [
        -4.60, -4.45, -4.30, -3.98, -3.65, -3.25, -2.85, -2.45, -2.05,
        -1.65, -1.25, -0.82, -0.40, 0.00, 0.45, 0.88, 1.30, 1.70,
        2.05, 2.25, 2.45, 2.65, 2.85, 3.05, 3.25, 3.45, 3.65, 3.85,
        4.10, 4.28, 4.38, 4.42
    ]
    seam_stations = {-4.30, -3.65, -2.85, -2.05, -1.25, -0.40, 0.45, 1.30, 2.05, 3.85}

    z_rings = []
    for z in base_stations:
        if z in seam_stations:
            z_rings.append((z - 0.016, 0.0, False))
            z_rings.append((z - 0.005, 0.0, True))
            z_rings.append((z, -0.0038, True))
            z_rings.append((z + 0.005, 0.0, True))
            z_rings.append((z + 0.016, 0.0, False))
        else:
            z_rings.append((z, 0.0, False))

    def is_inside_cockpit_tub(x, y, z):
        if y <= 0.45 or not (0.28 <= z <= 1.96):
            return False
        # Elliptical cockpit opening strictly inside the red cockpit collar
        dz = (z - 1.12) / 0.86
        dx = x / 0.34
        return (dz * dz + dx * dx) < 1.0

    n_circ = 72
    ring_verts = []
    for z, r_off, is_groove in z_rings:
        ring = []
        for i in range(n_circ):
            ang = 2.0 * math.pi * i / n_circ
            pt = hull_surface_point(z, ang, r_off)
            if is_inside_cockpit_tub(pt.x, pt.y, z):
                pt.y = 0.56
            ring.append(bm.verts.new(pt))
        ring_verts.append((z, is_groove, ring))

    for r_i in range(len(ring_verts) - 1):
        z0, g0, ring0 = ring_verts[r_i]
        z1, g1, ring1 = ring_verts[r_i + 1]
        z_mid = 0.5 * (z0 + z1)
        in_groove = (g0 and g1)

        for i in range(n_circ):
            ni = (i + 1) % n_circ
            f = bm.faces.new((ring0[i], ring0[ni], ring1[ni], ring1[i]))
            c_med = f.calc_center_median()

            mat_id = idx_chrome
            if -4.60 <= z_mid < -4.30:
                mat_id = idx_paint
            elif is_inside_cockpit_tub(c_med.x, c_med.y + 0.15, z_mid) and c_med.y <= 0.62:
                mat_id = idx_cockpit

            f.material_index = mat_id
            f.smooth = not in_groove

    front_ring = ring_verts[-1][2]
    v_nose_tip = bm.verts.new(Vector((0.0, -0.01, 4.435)))
    for i in range(n_circ):
        ni = (i + 1) % n_circ
        f = bm.faces.new((front_ring[i], front_ring[ni], v_nose_tip))
        f.material_index = idx_chrome
        f.smooth = True

    # --------------------------------------------------------------------------
    # B. CURVE-ALIGNED CONFORMAL PAINT & WHITE ENAMEL LIVERY PANELS + CHROME TRIMS
    # --------------------------------------------------------------------------
    def add_trim_strip(pts_normals, width=0.020, height=0.007, mat_id=idx_chrome):
        strip_rings = []
        for pt, nrm in pts_normals:
            n = nrm.normalized()
            tang = Vector((0, 0, 1)) if abs(n.z) < 0.9 else Vector((0, 1, 0))
            side = n.cross(tang).normalized()
            v0 = bm.verts.new(pt - side * (width * 0.5) - n * 0.001)
            v1 = bm.verts.new(pt - side * (width * 0.32) + n * height)
            v2 = bm.verts.new(pt + side * (width * 0.32) + n * height)
            v3 = bm.verts.new(pt + side * (width * 0.5) - n * 0.001)
            strip_rings.append((v0, v1, v2, v3))
        for i in range(len(strip_rings) - 1):
            r0 = strip_rings[i]
            r1 = strip_rings[i + 1]
            for j in range(3):
                f = bm.faces.new((r0[j], r0[j + 1], r1[j + 1], r1[j]))
                f.material_index = mat_id
                f.smooth = True

    def red_bounds_ang(z):
        if z >= 3.45:
            return -0.5 * math.pi, 0.5 * math.pi
        elif z >= 2.15:
            t = (3.45 - z) / (3.45 - 2.15)
            st = t ** 0.78
            top_a = 0.5 * math.pi * (1.0 - st) + 0.17 * st
            bot_a = -0.5 * math.pi * (1.0 - st) - 0.21 * st
            return bot_a, top_a
        else:
            t = max(0.0, min(1.0, (2.15 - z) / (2.15 - (-3.95))))
            half_w = 0.18 * (1.0 - 0.45 * t)
            return -0.04 - half_w, -0.04 + half_w

    def white_upper_top_ang(z):
        if z >= 2.95:
            return 0.5 * math.pi
        t = max(0.0, min(1.0, (2.95 - z) / (2.95 - 1.95)))
        st = t ** 0.8
        _, r_top = red_bounds_ang(z)
        return 0.5 * math.pi * (1.0 - st) + r_top * st

    def white_lower_bot_ang(z):
        if z >= 3.05:
            return -0.5 * math.pi
        t = max(0.0, min(1.0, (3.05 - z) / (3.05 - 2.05)))
        st = t ** 0.8
        r_bot, _ = red_bounds_ang(z)
        return -0.5 * math.pi * (1.0 - st) + r_bot * st

    # Use EXACT matching Z stations (including every 0.10m on nose cone) and +0.0055m offset
    # so the conformal Paint & Enamel_White shells never intersect the hull!
    z_livery_pts = [
        -3.98, -3.65, -3.25, -2.85, -2.45, -2.05, -1.65, -1.25, -0.82, -0.40,
        0.00, 0.45, 0.88, 1.30, 1.70, 2.05,
        2.15, 2.25, 2.35, 2.45, 2.55, 2.65, 2.75, 2.85, 2.95, 3.05, 3.15, 3.25,
        3.35, 3.45, 3.55, 3.65, 3.75, 3.85
    ]
    n_span_red = 20
    for sign_x in (-1.0, 1.0):
        red_grid = []
        upper_trim = []
        lower_trim = []
        for z in z_livery_pts:
            a_bot, a_top = red_bounds_ang(z)
            col = []
            for j in range(n_span_red + 1):
                frac = j / n_span_red
                a_local = a_bot + frac * (a_top - a_bot)
                ang = a_local if sign_x > 0 else (math.pi - a_local)
                pt = hull_surface_point(z, ang, 0.0055)
                col.append(bm.verts.new(pt))
                if j == 0 and z < 3.45:
                    lower_trim.append((pt, Vector((math.cos(ang), math.sin(ang), 0.0))))
                if j == n_span_red and z < 3.45:
                    upper_trim.append((pt, Vector((math.cos(ang), math.sin(ang), 0.0))))
            red_grid.append(col)

        for k in range(len(red_grid) - 1):
            c0, c1 = red_grid[k], red_grid[k + 1]
            for j in range(n_span_red):
                f = bm.faces.new((c0[j], c0[j + 1], c1[j + 1], c1[j]))
                if sign_x < 0:
                    f.normal_flip()
                f.material_index = idx_paint
                f.smooth = True

        add_trim_strip(upper_trim, width=0.016, height=0.006, mat_id=idx_chrome)
        add_trim_strip(lower_trim, width=0.016, height=0.006, mat_id=idx_chrome)

        # Conformal White Enamel Upper & Lower Nose Cheek Bands
        z_white_pts = [1.95, 2.05, 2.15, 2.25, 2.35, 2.45, 2.55, 2.65, 2.75, 2.85, 2.95, 3.05, 3.15, 3.25, 3.35, 3.45]
        n_span_w = 12
        w_up_grid = []
        w_up_trim = []
        for z in z_white_pts:
            _, a_bot = red_bounds_ang(z)
            a_top = white_upper_top_ang(z)
            col = []
            for j in range(n_span_w + 1):
                frac = j / n_span_w
                a_local = a_bot + frac * max(0.005, (a_top - a_bot))
                ang = a_local if sign_x > 0 else (math.pi - a_local)
                pt = hull_surface_point(z, ang, 0.0052)
                col.append(bm.verts.new(pt))
                if j == n_span_w and z < 2.95:
                    w_up_trim.append((pt, Vector((math.cos(ang), math.sin(ang), 0.0))))
            w_up_grid.append(col)
        for k in range(len(w_up_grid) - 1):
            c0, c1 = w_up_grid[k], w_up_grid[k + 1]
            for j in range(n_span_w):
                f = bm.faces.new((c0[j], c0[j + 1], c1[j + 1], c1[j]))
                if sign_x < 0:
                    f.normal_flip()
                f.material_index = idx_white
                f.smooth = True
        add_trim_strip(w_up_trim, width=0.016, height=0.006, mat_id=idx_chrome)

        z_white_low = [2.05, 2.15, 2.25, 2.35, 2.45, 2.55, 2.65, 2.75, 2.85, 2.95, 3.05, 3.15, 3.25, 3.35, 3.45]
        w_low_grid = []
        w_low_trim = []
        for z in z_white_low:
            a_top, _ = red_bounds_ang(z)
            a_bot = white_lower_bot_ang(z)
            col = []
            for j in range(n_span_w + 1):
                frac = j / n_span_w
                a_local = a_bot + frac * max(0.005, (a_top - a_bot))
                ang = a_local if sign_x > 0 else (math.pi - a_local)
                pt = hull_surface_point(z, ang, 0.0052)
                col.append(bm.verts.new(pt))
                if j == 0 and z < 3.05:
                    w_low_trim.append((pt, Vector((math.cos(ang), math.sin(ang), 0.0))))
            w_low_grid.append(col)
        for k in range(len(w_low_grid) - 1):
            c0, c1 = w_low_grid[k], w_low_grid[k + 1]
            for j in range(n_span_w):
                f = bm.faces.new((c0[j], c0[j + 1], c1[j + 1], c1[j]))
                if sign_x < 0:
                    f.normal_flip()
                f.material_index = idx_white
                f.smooth = True
        add_trim_strip(w_low_trim, width=0.016, height=0.006, mat_id=idx_chrome)

    # Circumferential chrome rings at nose tip (z=3.85) and aft collar (z=-4.30, z=-4.58)
    for z_ring_pos, w_strip in [(3.85, 0.022), (-4.30, 0.026), (-4.58, 0.028)]:
        ring_pts = []
        for i in range(n_circ + 1):
            ang = 2.0 * math.pi * (i % n_circ) / n_circ
            pt = hull_surface_point(z_ring_pos, ang, 0.004)
            nrm = Vector((math.cos(ang), math.sin(ang), 0.0))
            ring_pts.append((pt, nrm))
        add_trim_strip(ring_pts, width=w_strip, height=0.009, mat_id=idx_chrome)

    # --------------------------------------------------------------------------
    # B2. SCULPTED FLUSH COCKPIT COMBING COLLAR (RED PAINT + CHROME BORDER)
    # --------------------------------------------------------------------------
    collar_rings = []
    inner_collar_trim = []
    outer_collar_trim = []
    n_col = 48
    for i in range(n_col):
        ang = 2.0 * math.pi * i / n_col
        ca, sa = math.cos(ang), math.sin(ang)
        z_base = 1.14 + 1.08 * ca
        x_base = 0.43 * sa
        if z_base > 1.20:
            x_base *= (1.0 - 0.32 * ((z_base - 1.20) / 1.06))

        def collar_pt(scale_x, add_z, lift):
            zc = 1.14 + (1.08 + add_z) * ca
            xc = x_base * scale_x
            rh, ych = get_hull_r_yc(zc)
            yc = ych + math.sqrt(max(0.04, rh * rh - xc * xc)) + lift
            return Vector((xc, yc, zc))

        p_tub = Vector((x_base * 0.85, 0.54, 1.14 + 1.02 * ca))
        p_in = collar_pt(0.95, 0.00, 0.016)
        p_mid = collar_pt(1.14, 0.09, 0.011)
        p_out = collar_pt(1.32, 0.18, 0.0065)
        collar_rings.append((
            bm.verts.new(p_tub),
            bm.verts.new(p_in),
            bm.verts.new(p_mid),
            bm.verts.new(p_out),
        ))
        inner_collar_trim.append((p_in, Vector((p_in.x, p_in.y, 0.0)).normalized()))
        outer_collar_trim.append((p_out, Vector((p_out.x, p_out.y, 0.0)).normalized()))
    inner_collar_trim.append(inner_collar_trim[0])
    outer_collar_trim.append(outer_collar_trim[0])

    for i in range(n_col):
        ni = (i + 1) % n_col
        r0 = collar_rings[i]
        r1 = collar_rings[ni]
        f0 = bm.faces.new((r0[0], r0[1], r1[1], r1[0]))
        f0.material_index = idx_cockpit
        f0.smooth = True
        f1 = bm.faces.new((r0[1], r0[2], r1[2], r1[1]))
        f1.material_index = idx_paint
        f1.smooth = True
        f2 = bm.faces.new((r0[2], r0[3], r1[3], r1[2]))
        f2.material_index = idx_paint
        f2.smooth = True

    add_trim_strip(inner_collar_trim, width=0.015, height=0.006, mat_id=idx_chrome)
    add_trim_strip(outer_collar_trim, width=0.018, height=0.006, mat_id=idx_chrome)

    # --------------------------------------------------------------------------
    # C. WING-ROOT TORPEDO FAIRINGS, FORWARD INTAKES & SIDE BALL-TIP STRAKES
    # --------------------------------------------------------------------------
    for sign_x in (-1.0, 1.0):
        f_stations = [
            ( 1.02, 0.03, 0.025, 0.72),
            ( 0.58, 0.14, 0.105, 0.78),
            ( 0.00, 0.18, 0.135, 0.80),
            (-1.00, 0.18, 0.135, 0.78),
            (-2.20, 0.15, 0.115, 0.72),
            (-3.30, 0.10, 0.085, 0.64),
            (-3.88, 0.055, 0.055, 0.60),
        ]
        n_f = 20
        f_rings = []
        for z_f, rx_f, ry_f, x_base in f_stations:
            ring = []
            for j in range(n_f):
                a = 2.0 * math.pi * j / n_f
                px = sign_x * (x_base + rx_f * math.cos(a))
                py = -0.04 + ry_f * math.sin(a)
                ring.append(bm.verts.new(Vector((px, py, z_f))))
            f_rings.append(ring)
        for k in range(len(f_rings) - 1):
            r0, r1 = f_rings[k], f_rings[k + 1]
            for j in range(n_f):
                nj = (j + 1) % n_f
                a = 2.0 * math.pi * (j + 0.5) / n_f
                mat_f = idx_paint if abs(math.sin(a)) > 0.36 else idx_chrome
                f = bm.faces.new((r0[j], r0[nj], r1[nj], r1[j]))
                if sign_x < 0:
                    f.normal_flip()
                f.material_index = mat_f
                f.smooth = True

        # Side fairing ball tip at z = -3.96
        create_prim(
            bm, bmesh.ops.create_uvsphere, idx_chrome, smooth=True,
            u_segments=20, v_segments=12, radius=0.095,
            matrix=Matrix.Translation((sign_x * 0.64, -0.04, -3.96))
        )

        # Forward intake scoop recess at z = +0.52
        create_prim(
            bm, bmesh.ops.create_cone, idx_scorch, smooth=True,
            cap_ends=True, cap_tris=False, segments=16,
            radius1=0.055, radius2=0.035, depth=0.18,
            matrix=Matrix.Translation((sign_x * 0.83, -0.04, 0.48)) @ Matrix.Scale(0.65, 4, (0, 1, 0))
        )

    # --------------------------------------------------------------------------
    # D. THREE SWEPT TAIL FINS WITH BALL TIPS (LEFT WING, RIGHT WING, DORSAL FIN)
    # --------------------------------------------------------------------------
    def build_airfoil_section(x, y, z_le, z_te, thickness, up_vec, n_pts=24):
        chord = z_le - z_te
        pts = []
        for i in range(n_pts):
            angle = 2.0 * math.pi * i / n_pts
            s = 0.5 * (1.0 - math.cos(angle))
            sign_side = 1.0 if i < n_pts // 2 else -1.0
            if i == 0 or i == n_pts // 2:
                sign_side = 0.0
            yt = 5.0 * (thickness / max(0.1, chord)) * (
                0.2969 * math.sqrt(max(0.0, s))
                - 0.1260 * s
                - 0.3516 * (s ** 2)
                + 0.2843 * (s ** 3)
                - 0.0950 * (s ** 4)
            ) * chord
            z_pos = z_le - s * chord
            pt = Vector((x, y, z_pos)) + up_vec * (sign_side * yt)
            pts.append((pt, s, sign_side))
        return pts

    # 1) LEFT (+X) AND RIGHT (-X) SWEPT MAIN WINGS (Span = 6.80 m across ball tips!)
    wing_stations = [
        # (abs_x, y, z_le, z_te, thickness)
        (0.65, -0.05,  0.42, -3.45, 0.20),
        (1.10, -0.17, -0.32, -3.62, 0.17),
        (1.70, -0.34, -1.35, -3.86, 0.14),
        (2.30, -0.52, -2.42, -4.12, 0.115),
        (2.82, -0.68, -3.35, -4.46, 0.092),
        (3.26, -0.82, -4.12, -4.76, 0.075),
    ]

    for sign_x in (-1.0, 1.0):
        up_w = Vector((-sign_x * 0.24, 0.97, 0.0)).normalized()
        sec_rings = []
        for ax, wy, zle, zte, thk in wing_stations:
            pts = build_airfoil_section(sign_x * ax, wy, zle, zte, thk, up_w, n_pts=24)
            ring = [bm.verts.new(p[0]) for p in pts]
            sec_rings.append((ax, wy, zle, zte, ring, pts))

        for k in range(len(sec_rings) - 1):
            ax0, wy0, zle0, zte0, r0, pts0 = sec_rings[k]
            ax1, wy1, zle1, zte1, r1, pts1 = sec_rings[k + 1]
            n_r = len(r0)
            for i in range(n_r):
                ni = (i + 1) % n_r
                f = bm.faces.new((r0[i], r0[ni], r1[ni], r1[i]))
                if sign_x < 0:
                    f.normal_flip()
                f.material_index = idx_chrome
                f.smooth = True

        for ax, wy, zle, zte, _ in wing_stations[1:-1]:
            rib_pts = []
            for m_i in range(10):
                t = m_i / 9.0
                z_p = zle - t * (zle - zte)
                y_off = 0.052 * math.sin(math.pi * t)
                rib_pts.append((Vector((sign_x * ax, wy + y_off, z_p)), Vector((0, 1, 0))))
            add_trim_strip(rib_pts, width=0.015, height=0.005, mat_id=idx_chrome)

        # Wingtip horizontal torpedo pod along Z axis
        pod_mat = Matrix.Translation((sign_x * 3.26, -0.82, -4.44))
        create_prim(
            bm, bmesh.ops.create_cone, idx_chrome, smooth=True,
            cap_ends=True, cap_tris=False, segments=20,
            radius1=0.062, radius2=0.045, depth=0.68, matrix=pod_mat
        )

        # Wingtip Spherical Ball Tip at X = ±3.26, radius = 0.14 -> outer X = ±3.40 m (6.80 m span!)
        ball_center = Vector((sign_x * 3.26, -0.82, -4.82))
        ball_faces = create_prim(
            bm, bmesh.ops.create_uvsphere, idx_chrome, smooth=True,
            u_segments=28, v_segments=16, radius=0.14,
            matrix=Matrix.Translation(ball_center)
        )
        for f in ball_faces:
            if (f.calc_center_median() - ball_center).x * sign_x > 0.115:
                f.material_index = idx_lamps

    # 2) TALL SWEPT DORSAL TAIL FIN WITH BALL TIP (Top reaches Y = +2.35 m!)
    dorsal_stations = [
        # (x, y, z_le, z_te, thickness)
        (0.0, 0.48, -0.55, -4.32, 0.18),
        (0.0, 0.78, -1.55, -4.28, 0.15),
        (0.0, 1.05, -2.15, -3.98, 0.135),
        (0.0, 1.48, -2.88, -4.30, 0.115),
        (0.0, 1.84, -3.50, -4.55, 0.096),
        (0.0, 2.06, -3.86, -4.72, 0.084),
        (0.0, 2.21, -4.12, -4.85, 0.076),
    ]
    up_d = Vector((1.0, 0.0, 0.0))
    d_rings = []
    for dx, dy, zle, zte, thk in dorsal_stations:
        pts = build_airfoil_section(dx, dy, zle, zte, thk, up_d, n_pts=24)
        ring = [bm.verts.new(p[0]) for p in pts]
        d_rings.append((dy, zle, zte, ring, pts))

    for k in range(len(d_rings) - 1):
        dy0, zle0, zte0, r0, pts0 = d_rings[k]
        dy1, zle1, zte1, r1, pts1 = d_rings[k + 1]
        n_r = len(r0)
        for i in range(n_r):
            ni = (i + 1) % n_r
            f = bm.faces.new((r0[i], r0[ni], r1[ni], r1[i]))
            f.material_index = idx_chrome
            f.smooth = True

    # Dorsal fin horizontal tip torpedo cap along Z axis (at Y = 2.21 m)
    d_pod_mat = Matrix.Translation((0.0, 2.21, -4.48))
    create_prim(
        bm, bmesh.ops.create_cone, idx_chrome, smooth=True,
        cap_ends=True, cap_tris=False, segments=20,
        radius1=0.062, radius2=0.044, depth=0.78, matrix=d_pod_mat
    )

    # Dorsal Ball Tip at (0.0, 2.21, -4.92), radius = 0.14 m -> max Y = +2.35 m!
    d_ball_center = Vector((0.0, 2.21, -4.92))
    d_ball_faces = create_prim(
        bm, bmesh.ops.create_uvsphere, idx_chrome, smooth=True,
        u_segments=28, v_segments=16, radius=0.14,
        matrix=Matrix.Translation(d_ball_center)
    )
    for f in d_ball_faces:
        if (f.calc_center_median() - d_ball_center).z < -0.115:
            f.material_index = idx_lamps

    # Small retro dorsal radio mast / pitot stub behind canopy (z = -0.25)
    create_prim(
        bm, bmesh.ops.create_cone, idx_chrome, smooth=True,
        cap_ends=True, cap_tris=False, segments=12,
        radius1=0.022, radius2=0.012, depth=0.14,
        matrix=Matrix.Translation((0.0, 0.86, -0.25)) @ Matrix.Rotation(math.radians(-68), 4, 'X')
    )

    # --------------------------------------------------------------------------
    # E. REAR RING FIN (ANNULAR SHROUD) & 4 RADIAL STRUTS
    # --------------------------------------------------------------------------
    ring_profile = [
        (-4.15, 1.15, 1.13),
        (-4.20, 1.215, 1.085),
        (-4.32, 1.250, 1.065),
        (-4.48, 1.242, 1.075),
        (-4.64, 1.205, 1.105),
        (-4.74, 1.160, 1.142),
    ]
    n_rf_circ = 64
    rf_loops = []
    closed_prof = [(z, ro) for (z, ro, ri) in ring_profile] + [(z, ri) for (z, ro, ri) in reversed(ring_profile)]
    for z_p, r_p in closed_prof:
        loop = []
        for i in range(n_rf_circ):
            ang = 2.0 * math.pi * i / n_rf_circ
            loop.append(bm.verts.new(Vector((r_p * math.cos(ang), r_p * math.sin(ang), z_p))))
        rf_loops.append(loop)

    n_prof = len(rf_loops)
    for p_i in range(n_prof):
        np_i = (p_i + 1) % n_prof
        l0, l1 = rf_loops[p_i], rf_loops[np_i]
        for i in range(n_rf_circ):
            ni = (i + 1) % n_rf_circ
            f = bm.faces.new((l0[i], l0[ni], l1[ni], l1[i]))
            f.material_index = idx_chrome
            f.smooth = True

    for strut_ang in (0.0, 0.5 * math.pi, math.pi, 1.5 * math.pi):
        ca, sa = math.cos(strut_ang), math.sin(strut_ang)
        up_s = Vector((-sa, ca, 0.0))
        s_rings = []
        for r_s in (0.44, 0.65, 0.90, 1.10):
            pts = build_airfoil_section(r_s * ca, r_s * sa, -4.22, -4.68, 0.068, up_s, n_pts=16)
            s_rings.append([bm.verts.new(p[0]) for p in pts])
        for k in range(len(s_rings) - 1):
            r0, r1 = s_rings[k], s_rings[k + 1]
            for i in range(len(r0)):
                ni = (i + 1) % len(r0)
                f = bm.faces.new((r0[i], r0[ni], r1[ni], r1[i]))
                f.material_index = idx_chrome
                f.smooth = True

    # --------------------------------------------------------------------------
    # F. CORRUGATED REAR ROCKET NOZZLE (Z = -4.55 to -5.50 m) + GLOWING INTERIOR
    # --------------------------------------------------------------------------
    noz_stations = [
        (-4.55, 0.455, 0.34),
        (-4.75, 0.465, 0.28),
        (-5.00, 0.445, 0.31),
        (-5.28, 0.415, 0.34),
        (-5.50, 0.380, 0.355),
    ]
    n_noz = 48
    outer_noz_rings = []
    inner_noz_rings = []
    for z_n, r_out, r_in in noz_stations:
        o_ring = []
        i_ring = []
        for i in range(n_noz):
            ang = 2.0 * math.pi * i / n_noz
            corrug = 0.012 * math.cos(24.0 * ang)
            ca, sa = math.cos(ang), math.sin(ang)
            o_ring.append(bm.verts.new(Vector(((r_out + corrug) * ca, (r_out + corrug) * sa, z_n))))
            i_ring.append(bm.verts.new(Vector((r_in * ca, r_in * sa, z_n))))
        outer_noz_rings.append(o_ring)
        inner_noz_rings.append(i_ring)

    for k in range(len(noz_stations) - 1):
        or0, or1 = outer_noz_rings[k], outer_noz_rings[k + 1]
        ir0, ir1 = inner_noz_rings[k], inner_noz_rings[k + 1]
        for i in range(n_noz):
            ni = (i + 1) % n_noz
            fo = bm.faces.new((or0[i], or0[ni], or1[ni], or1[i]))
            fo.material_index = idx_scorch
            fo.smooth = True
            fi = bm.faces.new((ir1[i], ir1[ni], ir0[ni], ir0[i]))
            fi.material_index = idx_engines if k < 3 else idx_scorch
            fi.smooth = True

    or_end = outer_noz_rings[-1]
    ir_end = inner_noz_rings[-1]
    for i in range(n_noz):
        ni = (i + 1) % n_noz
        f = bm.faces.new((or_end[i], or_end[ni], ir_end[ni], ir_end[i]))
        f.material_index = idx_scorch
        f.smooth = True

    v_comb_center = bm.verts.new(Vector((0.0, 0.0, -4.62)))
    ir_start = inner_noz_rings[0]
    for i in range(n_noz):
        ni = (i + 1) % n_noz
        f = bm.faces.new((ir_start[ni], ir_start[i], v_comb_center))
        f.material_index = idx_engines
        f.smooth = True

    for z_band, r_band in [(-4.72, 0.478), (-5.12, 0.442)]:
        create_prim(
            bm, bmesh.ops.create_cone, idx_scorch, smooth=True,
            cap_ends=False, cap_tris=False, segments=36,
            radius1=r_band, radius2=r_band, depth=0.045,
            matrix=Matrix.Translation((0.0, 0.0, z_band))
        )

    # --------------------------------------------------------------------------
    # G. BELLY VENT HATCHES & CANTED VENTRAL AUXILIARY THRUSTERS
    # --------------------------------------------------------------------------
    for sign_x in (-1.0, 1.0):
        # Louvered Vent Hatch bay at z = -1.35 to -0.85 (Scorched_Metal slats!)
        for slat_i in range(7):
            z_s = -1.32 + slat_i * 0.072
            slat_mat = (
                Matrix.Translation((sign_x * 0.44, -0.63, z_s))
                @ Matrix.Rotation(math.radians(sign_x * 32), 4, 'Z')
                @ Matrix.Rotation(math.radians(25), 4, 'X')
                @ Matrix.Diagonal((0.18, 0.012, 0.038, 1.0))
            )
            create_prim(bm, bmesh.ops.create_cube, idx_scorch, smooth=False, size=1.0, matrix=slat_mat)

        # Ventral Auxiliary Thruster pod (Red Paint shroud + Scorched nozzle + Orange Engines core)
        thr_pos = Vector((sign_x * 0.28, -0.68, -2.08))
        thr_rot = Matrix.Rotation(math.radians(-22), 4, 'X') @ Matrix.Rotation(math.radians(sign_x * 12), 4, 'Y')
        create_prim(
            bm, bmesh.ops.create_cone, idx_paint, smooth=True,
            cap_ends=True, cap_tris=False, segments=20,
            radius1=0.145, radius2=0.105, depth=0.52,
            matrix=Matrix.Translation(thr_pos) @ thr_rot
        )

        exit_pos = thr_pos + Vector((0.0, -0.08, -0.22))
        noz_faces = create_prim(
            bm, bmesh.ops.create_cone, idx_scorch, smooth=True,
            cap_ends=True, cap_tris=False, segments=20,
            radius1=0.115, radius2=0.095, depth=0.16,
            matrix=Matrix.Translation(exit_pos) @ thr_rot
        )
        for f in noz_faces:
            if f.normal.z < -0.5:
                f.material_index = idx_engines

    # --------------------------------------------------------------------------
    # H. NOSE RAY-GUN FAIRING BLISTERS ON HULL
    # --------------------------------------------------------------------------
    for (tx, ty, z_start) in [(0.16, 0.16, 3.25), (-0.16, -0.16, 3.05)]:
        blister_mat = Matrix.Translation((tx, ty, 0.5 * (z_start + 4.05)))
        create_prim(
            bm, bmesh.ops.create_cone, idx_chrome, smooth=True,
            cap_ends=True, cap_tris=False, segments=18,
            radius1=0.082, radius2=0.068, depth=(4.05 - z_start),
            matrix=blister_mat
        )

    # --------------------------------------------------------------------------
    # I. DETAILED 1950S COCKPIT INTERIOR (SEAT, INSTRUMENT PANEL, RADAR SCOPE, STICK)
    # --------------------------------------------------------------------------
    seat_parts = [
        ((0.0, 0.42, 0.92), (0.40, 0.09, 0.42), -5.0, idx_cockpit),
        ((0.0, 0.65, 0.72), (0.36, 0.46, 0.09), -16.0, idx_cockpit),
        ((0.0, 0.89, 0.64), (0.24, 0.14, 0.09), -10.0, idx_cockpit),
        ((0.21, 0.56, 0.84), (0.04, 0.22, 0.36), -8.0, idx_chrome),
        ((-0.21, 0.56, 0.84), (0.04, 0.22, 0.36), -8.0, idx_chrome),
    ]
    for pos, scl, rx_deg, m_id in seat_parts:
        m = Matrix.Translation(pos) @ Matrix.Rotation(math.radians(rx_deg), 4, 'X') @ Matrix.Diagonal((*scl, 1.0))
        create_prim(bm, bmesh.ops.create_cube, m_id, smooth=False, size=1.0, matrix=m)

    # Curved Instrument Panel & Coaming Shroud inside front windscreen
    create_prim(
        bm, bmesh.ops.create_cone, idx_cockpit, smooth=True,
        cap_ends=True, cap_tris=False, segments=20,
        radius1=0.34, radius2=0.22, depth=0.45,
        matrix=Matrix.Translation((0.0, 0.72, 1.68)) @ Matrix.Scale(0.75, 4, (0, 1, 0))
    )

    for dx, dy, r_dial in [(0.0, 0.82, 0.070), (0.13, 0.78, 0.042), (-0.13, 0.78, 0.042)]:
        d_mat = Matrix.Translation((dx, dy, 1.46))
        dial_faces = create_prim(
            bm, bmesh.ops.create_cone, idx_chrome, smooth=True,
            cap_ends=True, cap_tris=False, segments=16,
            radius1=r_dial, radius2=r_dial * 0.92, depth=0.08, matrix=d_mat
        )
        for f in dial_faces:
            if f.normal.z < -0.5:
                f.material_index = idx_lamps

    stick_mat = Matrix.Translation((0.0, 0.56, 1.26)) @ Matrix.Rotation(math.radians(12), 4, 'X')
    create_prim(
        bm, bmesh.ops.create_cone, idx_scorch, smooth=True,
        cap_ends=True, cap_tris=False, segments=12,
        radius1=0.022, radius2=0.016, depth=0.34, matrix=stick_mat
    )

    # --------------------------------------------------------------------------
    # J. ROWS OF 3D DOME-HEAD RIVETS ALONG PANEL SEAMS, COLLAR & FINS
    # --------------------------------------------------------------------------
    for z_seam in sorted(seam_stations):
        for offset_z in (-0.028, 0.028):
            z_r = z_seam + offset_z
            for i in range(40):
                ang = 2.0 * math.pi * i / 40.0
                if 0.02 < z_r < 2.28 and math.sin(ang) > 0.58:
                    continue
                pt = hull_surface_point(z_r, ang, 0.004)
                nrm = Vector((math.cos(ang), math.sin(ang), 0.0))
                add_rivet_dome(bm, pt, nrm, radius=0.007, height=0.0038, mat_idx=idx_chrome, segs=6)

    for sign_x in (-1.0, 1.0):
        base_ang = 0.0 if sign_x > 0 else math.pi
        for d_ang in (-0.23, 0.19, 0.48, -0.48):
            for k in range(36):
                z_r = -4.20 + k * (7.95 / 35.0)
                ang = base_ang + (d_ang if sign_x > 0 else -d_ang)
                if 0.02 < z_r < 2.28 and math.sin(ang) > 0.58:
                    continue
                pt = hull_surface_point(z_r, ang, 0.006)
                nrm = Vector((math.cos(ang), math.sin(ang), 0.0))
                add_rivet_dome(bm, pt, nrm, radius=0.0068, height=0.0035, mat_idx=idx_chrome, segs=6)

    # Dorsal fin rivets strictly using dorsal_stations interpolation so none float!
    d_zle_tab = [(s[1], s[2]) for s in dorsal_stations]
    d_zte_tab = [(s[1], s[3]) for s in dorsal_stations]
    d_thk_tab = [(s[1], s[4]) for s in dorsal_stations]
    for dy in (0.72, 1.02, 1.38, 1.75, 2.05):
        z_le = linear_interp(dy, d_zle_tab) - 0.06
        z_te = linear_interp(dy, d_zte_tab) + 0.08
        thk = linear_interp(dy, d_thk_tab)
        for m_i in range(10):
            t = (m_i + 0.5) / 10.0
            z_p = z_le - t * (z_le - z_te)
            x_off = 0.48 * thk * math.sin(math.pi * t)
            for sx in (-1.0, 1.0):
                add_rivet_dome(
                    bm, Vector((sx * x_off, dy, z_p)), Vector((sx, 0.0, 0.0)),
                    radius=0.0065, height=0.0035, mat_idx=idx_chrome, segs=6
                )

    # Wing rivets strictly using wing_stations interpolation!
    w_wy_tab = [(s[0], s[1]) for s in wing_stations]
    w_zle_tab = [(s[0], s[2]) for s in wing_stations]
    w_zte_tab = [(s[0], s[3]) for s in wing_stations]
    w_thk_tab = [(s[0], s[4]) for s in wing_stations]
    for sign_x in (-1.0, 1.0):
        for ax in (1.05, 1.55, 2.05, 2.55, 2.98):
            wy = linear_interp(ax, w_wy_tab)
            z_le = linear_interp(ax, w_zle_tab) - 0.06
            z_te = linear_interp(ax, w_zte_tab) + 0.08
            thk = linear_interp(ax, w_thk_tab)
            for m_i in range(10):
                t = (m_i + 0.5) / 10.0
                z_p = z_le - t * (z_le - z_te)
                y_off = 0.48 * thk * math.sin(math.pi * t)
                for sy in (-1.0, 1.0):
                    add_rivet_dome(
                        bm, Vector((sign_x * ax, wy + sy * y_off, z_p)), Vector((0.0, sy, 0.0)),
                        radius=0.0065, height=0.0035, mat_idx=idx_chrome, segs=6
                    )

    bm.normal_update()
    bm.to_mesh(mesh)
    bm.free()

    bpy.ops.object.select_all(action='DESELECT')
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.uv.cube_project(cube_size=2.2)
    bpy.ops.object.mode_set(mode='OBJECT')

    return obj


# ==============================================================================
# 6. BUILD MOVING BUBBLE CANOPY NODE (FLUSH SILL + PIVOT ON REAR HINGE)
# ==============================================================================
def build_canopy_object(mats, hull_obj):
    pivot = Vector((0.0, 0.80, 0.10))
    mesh = bpy.data.meshes.new("Canopy_Mesh")
    obj = bpy.data.objects.new("Canopy", mesh)
    bpy.context.collection.objects.link(obj)
    obj.parent = hull_obj
    obj.location = pivot

    idx_chrome = ensure_mat_slot(obj, mats['Chrome'])
    idx_glass = ensure_mat_slot(obj, mats['Glass'])

    bm = bmesh.new()

    canopy_stations = [
        (0.06, 0.08, 0.03),
        (0.24, 0.25, 0.18),
        (0.55, 0.37, 0.38),
        (0.92, 0.42, 0.50),
        (1.26, 0.40, 0.50),
        (1.58, 0.35, 0.40),
        (1.86, 0.26, 0.26),
        (2.08, 0.15, 0.12),
        (2.22, 0.04, 0.02),
    ]

    n_arc = 32
    outer_rings = []
    inner_rings = []
    for z_w, wx, h_dome in canopy_stations:
        r_h, yc_h = get_hull_r_yc(z_w)
        y_sill = yc_h + math.sqrt(max(0.04, r_h * r_h - wx * wx)) * 0.98 + 0.015

        o_ring = []
        i_ring = []
        for j in range(n_arc + 1):
            ang = math.pi * j / n_arc
            ca = -math.cos(ang)
            sa = math.sin(ang)
            x_w = wx * ca * (1.0 + 0.08 * sa)
            y_w = y_sill + h_dome * (sa ** 0.80)

            p_out = Vector((x_w, y_w, z_w)) - pivot
            nrm = Vector((ca * 0.7, max(0.25, sa), 0.0)).normalized()
            p_in = p_out - nrm * 0.016
            o_ring.append(bm.verts.new(p_out))
            i_ring.append(bm.verts.new(p_in))
        outer_rings.append((z_w, wx, y_sill, h_dome, o_ring))
        inner_rings.append((z_w, i_ring))

    for k in range(len(outer_rings) - 1):
        z0, _, _, _, or0 = outer_rings[k]
        z1, _, _, _, or1 = outer_rings[k + 1]
        _, ir0 = inner_rings[k]
        _, ir1 = inner_rings[k + 1]
        z_mid = 0.5 * (z0 + z1)

        for j in range(n_arc):
            nj = j + 1
            is_sill = (j <= 1 or j >= n_arc - 2)
            is_arch = (1.50 <= z_mid <= 1.65)
            is_front_nose = (z_mid >= 2.05)
            is_rear_cap = (z_mid <= 0.25)
            is_windscreen_center_bar = (z_mid > 1.58 and abs(j - n_arc / 2.0) < 1.2)

            is_frame = is_sill or is_arch or is_front_nose or is_rear_cap or is_windscreen_center_bar
            m_id = idx_chrome if is_frame else idx_glass

            fo = bm.faces.new((or0[j], or0[nj], or1[nj], or1[j]))
            fo.material_index = m_id
            fo.smooth = True

            fi = bm.faces.new((ir1[j], ir1[nj], ir0[nj], ir0[j]))
            fi.material_index = m_id
            fi.smooth = True

            if j == 0:
                fs = bm.faces.new((ir0[0], or0[0], or1[0], ir1[0]))
                fs.material_index = idx_chrome
            if nj == n_arc:
                fs = bm.faces.new((or0[n_arc], ir0[n_arc], ir1[n_arc], or1[n_arc]))
                fs.material_index = idx_chrome

    # Raised external Chrome Frame Arch at windscreen hoop (z = 1.58)
    z_arch = 1.58
    r_h, yc_h = get_hull_r_yc(z_arch)
    wx_a = 0.355
    y_sill_a = yc_h + math.sqrt(max(0.04, r_h * r_h - wx_a * wx_a)) * 0.98 + 0.015
    arch_ring = []
    for j in range(n_arc + 1):
        ang = math.pi * j / n_arc
        ca = -math.cos(ang)
        sa = math.sin(ang)
        x_w = wx_a * ca * (1.0 + 0.08 * sa)
        y_w = y_sill_a + 0.405 * (sa ** 0.80)
        p_base = Vector((x_w, y_w, z_arch)) - pivot
        nrm = Vector((ca * 0.7, max(0.2, sa), 0.0)).normalized()
        v0 = bm.verts.new(p_base + Vector((0, 0, -0.025)))
        v1 = bm.verts.new(p_base + nrm * 0.012 + Vector((0, 0, -0.016)))
        v2 = bm.verts.new(p_base + nrm * 0.012 + Vector((0, 0,  0.016)))
        v3 = bm.verts.new(p_base + Vector((0, 0,  0.025)))
        arch_ring.append((v0, v1, v2, v3))
    for j in range(n_arc):
        r0, r1 = arch_ring[j], arch_ring[j + 1]
        for s_i in range(3):
            f = bm.faces.new((r0[s_i], r0[s_i + 1], r1[s_i + 1], r1[s_i]))
            f.material_index = idx_chrome
            f.smooth = True

    bm.normal_update()
    bm.to_mesh(mesh)
    bm.free()

    bpy.ops.object.select_all(action='DESELECT')
    obj.select_set(True)
    bpy.context.view_layer.objects.active = obj
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.uv.cube_project(cube_size=1.5)
    bpy.ops.object.mode_set(mode='OBJECT')

    return obj


# ==============================================================================
# 7. BUILD ARTICULATED NOSE RAY-GUN CANNONS (TURRET_1->GUN_1, TURRET_2->GUN_2)
# ==============================================================================
def build_raygun_assembly(idx, pivot_world, z_tip_world, coil_z_list, coil_radii, mats, hull_obj):
    t_mesh = bpy.data.meshes.new(f"Turret_{idx}_Mesh")
    t_obj = bpy.data.objects.new(f"Turret_{idx}", t_mesh)
    bpy.context.collection.objects.link(t_obj)
    t_obj.parent = hull_obj
    t_obj.location = pivot_world

    t_chrome = ensure_mat_slot(t_obj, mats['Chrome'])
    t_scorch = ensure_mat_slot(t_obj, mats['Scorched_Metal'])

    bm_t = bmesh.new()
    create_prim(
        bm_t, bmesh.ops.create_cone, t_chrome, smooth=True,
        cap_ends=True, cap_tris=False, segments=20,
        radius1=0.072, radius2=0.062, depth=0.18,
        matrix=Matrix.Translation((0.0, 0.0, 0.02))
    )

    for sx in (-0.055, 0.055):
        create_prim(
            bm_t, bmesh.ops.create_cone, t_scorch, smooth=True,
            cap_ends=True, cap_tris=False, segments=14,
            radius1=0.036, radius2=0.036, depth=0.025,
            matrix=Matrix.Translation((sx, 0.0, 0.08)) @ Matrix.Rotation(math.radians(90), 4, 'Y')
        )

    bm_t.normal_update()
    bm_t.to_mesh(t_mesh)
    bm_t.free()

    bpy.ops.object.select_all(action='DESELECT')
    t_obj.select_set(True)
    bpy.context.view_layer.objects.active = t_obj
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.uv.cube_project(cube_size=0.5)
    bpy.ops.object.mode_set(mode='OBJECT')

    gun_local_pos = Vector((0.0, 0.0, 0.08))
    gun_world_origin_z = pivot_world.z + gun_local_pos.z
    g_mesh = bpy.data.meshes.new(f"Gun_{idx}_Mesh")
    g_obj = bpy.data.objects.new(f"Gun_{idx}", g_mesh)
    bpy.context.collection.objects.link(g_obj)
    g_obj.parent = t_obj
    g_obj.location = gun_local_pos

    g_scorch = ensure_mat_slot(g_obj, mats['Scorched_Metal'])
    g_chrome = ensure_mat_slot(g_obj, mats['Chrome'])
    g_lamps = ensure_mat_slot(g_obj, mats['Lamps'])

    bm_g = bmesh.new()
    L_total = z_tip_world - gun_world_origin_z

    barrel_steps = [
        (-0.04, 0.056, g_scorch),
        ( 0.28, 0.054, g_scorch),
        ( 0.30, 0.044, g_chrome),
        ( 0.52, 0.042, g_chrome),
        ( 0.54, 0.034, g_scorch),
        ( L_total - 0.24, 0.028, g_scorch),
        ( L_total - 0.12, 0.020, g_chrome),
        ( L_total - 0.06, 0.040, g_chrome),
        ( L_total,        0.024, g_scorch),
    ]
    n_b = 20
    b_rings = []
    for zb, rb, m_id in barrel_steps:
        ring = []
        for i in range(n_b):
            ang = 2.0 * math.pi * i / n_b
            flute = 0.0035 * math.cos(8.0 * ang) if zb < 0.52 else 0.0
            ring.append(bm_g.verts.new(Vector(((rb + flute) * math.cos(ang), (rb + flute) * math.sin(ang), zb))))
        b_rings.append((ring, m_id))

    for k in range(len(b_rings) - 1):
        r0, m0 = b_rings[k]
        r1, _ = b_rings[k + 1]
        for i in range(n_b):
            ni = (i + 1) % n_b
            f = bm_g.faces.new((r0[i], r0[ni], r1[ni], r1[i]))
            f.material_index = m0
            f.smooth = True

    v_muz_center = bm_g.verts.new(Vector((0.0, 0.0, L_total - 0.01)))
    end_ring = b_rings[-1][0]
    for i in range(n_b):
        ni = (i + 1) % n_b
        f = bm_g.faces.new((end_ring[i], end_ring[ni], v_muz_center))
        f.material_index = g_lamps
        f.smooth = True

    for z_coil_w, r_coil in zip(coil_z_list, coil_radii):
        z_c = z_coil_w - gun_world_origin_z
        n_maj = 24
        n_min = 10
        r_tube = 0.016
        t_rings = []
        for i in range(n_maj):
            a_maj = 2.0 * math.pi * i / n_maj
            ca, sa = math.cos(a_maj), math.sin(a_maj)
            c_ring = []
            for j in range(n_min):
                a_min = 2.0 * math.pi * j / n_min
                rr = r_coil + r_tube * math.cos(a_min)
                zz = z_c + r_tube * math.sin(a_min) * 0.78
                c_ring.append(bm_g.verts.new(Vector((rr * ca, rr * sa, zz))))
            t_rings.append(c_ring)
        for i in range(n_maj):
            ni = (i + 1) % n_maj
            for j in range(n_min):
                nj = (j + 1) % n_min
                f = bm_g.faces.new((t_rings[i][j], t_rings[ni][j], t_rings[ni][nj], t_rings[i][nj]))
                f.material_index = g_chrome
                f.smooth = True

        create_prim(
            bm_g, bmesh.ops.create_cone, g_scorch, smooth=True,
            cap_ends=True, cap_tris=False, segments=16,
            radius1=r_coil - 0.008, radius2=r_coil - 0.008, depth=0.014,
            matrix=Matrix.Translation((0.0, 0.0, z_c))
        )

    bm_g.normal_update()
    bm_g.to_mesh(g_mesh)
    bm_g.free()

    bpy.ops.object.select_all(action='DESELECT')
    g_obj.select_set(True)
    bpy.context.view_layer.objects.active = g_obj
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.uv.cube_project(cube_size=0.5)
    bpy.ops.object.mode_set(mode='OBJECT')

    muz_empty = bpy.data.objects.new(f"Muzzle_{idx}", None)
    muz_empty.empty_display_type = 'PLAIN_AXES'
    muz_empty.empty_display_size = 0.15
    bpy.context.collection.objects.link(muz_empty)
    muz_empty.parent = g_obj
    muz_empty.location = Vector((0.0, 0.0, L_total))

    return t_obj, g_obj, muz_empty


# ==============================================================================
# 8. CREATE CONTRACT EMPTY NODES
# ==============================================================================
def create_contract_empties(hull_obj):
    empties_spec = [
        ("Thruster_1",  ( 0.00,  0.00, -5.50)),
        ("Thruster_2",  ( 0.28, -0.76, -2.30)),
        ("Thruster_3",  (-0.28, -0.76, -2.30)),
        ("Hardpoint_1", ( 1.65, -0.38, -2.20)),
        ("Hardpoint_2", (-1.65, -0.38, -2.20)),
        ("Seat_Pilot",  ( 0.00,  0.46,  0.92)),
        ("Camera_Chase",( 0.00,  2.80, -11.50)),
        ("Light_Nav_L", ( 3.40, -0.82, -4.82)),
        ("Light_Nav_R", (-3.40, -0.82, -4.82)),
        ("Light_Beacon",( 0.00,  2.35, -4.92)),
    ]
    created = []
    for name, loc in empties_spec:
        e = bpy.data.objects.new(name, None)
        e.empty_display_type = 'PLAIN_AXES'
        e.empty_display_size = 0.20
        bpy.context.collection.objects.link(e)
        e.parent = hull_obj
        e.location = Vector(loc)
        created.append(e)
    return created


def compute_scene_stats():
    depsgraph = bpy.context.evaluated_depsgraph_get()
    min_v = Vector(( 1e9,  1e9,  1e9))
    max_v = Vector((-1e9, -1e9, -1e9))
    total_tris = 0
    used_mats = set()

    for obj in bpy.context.scene.objects:
        if obj.type == 'MESH' and not obj.name.startswith("_"):
            eval_obj = obj.evaluated_get(depsgraph)
            m = eval_obj.to_mesh()
            m.calc_loop_triangles()
            total_tris += len(m.loop_triangles)
            mat_world = obj.matrix_world
            for v in m.vertices:
                wv = mat_world @ v.co
                min_v.x = min(min_v.x, wv.x)
                min_v.y = min(min_v.y, wv.y)
                min_v.z = min(min_v.z, wv.z)
                max_v.x = max(max_v.x, wv.x)
                max_v.y = max(max_v.y, wv.y)
                max_v.z = max(max_v.z, wv.z)
            for slot in obj.material_slots:
                if slot.material:
                    used_mats.add(slot.material.name)
            eval_obj.to_mesh_clear()

    size_v = max_v - min_v
    return min_v, max_v, size_v, total_tris, sorted(used_mats)


def patch_glb_paint_factor(glb_path, factor_rgba):
    with open(glb_path, "rb") as f:
        data = f.read()

    magic, version, length = struct.unpack_from("<III", data, 0)
    if magic != 0x46546C67:
        return

    json_len, json_type = struct.unpack_from("<II", data, 12)
    json_bytes = data[20:20 + json_len]
    gltf = json.loads(json_bytes.decode("utf-8"))

    for mat in gltf.get("materials", []):
        if mat.get("name") == "Paint":
            pbr = mat.setdefault("pbrMetallicRoughness", {})
            pbr["baseColorFactor"] = list(factor_rgba)

    new_json_bytes = json.dumps(gltf, separators=(",", ":")).encode("utf-8")
    pad = (4 - (len(new_json_bytes) % 4)) % 4
    new_json_bytes += b" " * pad

    rest = data[20 + json_len:]
    total_len = 12 + 8 + len(new_json_bytes) + len(rest)

    with open(glb_path, "wb") as f:
        f.write(struct.pack("<III", magic, version, total_len))
        f.write(struct.pack("<II", len(new_json_bytes), json_type))
        f.write(new_json_bytes)
        f.write(rest)


# ==============================================================================
# 9. LIGHTING, CAMERA & RENDER PIPELINE
# ==============================================================================
def setup_space_lighting_and_world(scene):
    world = bpy.data.worlds.new("SpaceWorld")
    scene.world = world
    wnodes = world.node_tree.nodes
    wlinks = world.node_tree.links
    wnodes.clear()

    w_out = wnodes.new('ShaderNodeOutputWorld')
    w_bg = wnodes.new('ShaderNodeBackground')
    w_tc = wnodes.new('ShaderNodeTexCoord')
    w_sep = wnodes.new('ShaderNodeSeparateXYZ')
    w_ramp = wnodes.new('ShaderNodeValToRGB')
    w_ramp.color_ramp.elements[0].position = 0.0
    w_ramp.color_ramp.elements[0].color = (0.004, 0.007, 0.016, 1.0)
    w_ramp.color_ramp.elements[1].position = 1.0
    w_ramp.color_ramp.elements[1].color = (0.018, 0.032, 0.065, 1.0)

    wlinks.new(w_tc.outputs['Window'], w_sep.inputs['Vector'])
    wlinks.new(w_sep.outputs['Y'], w_ramp.inputs['Fac'])
    wlinks.new(w_ramp.outputs['Color'], w_bg.inputs['Color'])
    w_bg.inputs['Strength'].default_value = 1.0
    wlinks.new(w_bg.outputs['Background'], w_out.inputs['Surface'])

    lights = []
    # 1) Strong white Key Light from upper front left
    key_data = bpy.data.lights.new("Key_Light", type='AREA')
    key_data.energy = 2600.0
    key_data.size = 9.5
    key_data.color = (1.0, 0.98, 0.95)
    key_obj = bpy.data.objects.new("Key_Light", key_data)
    bpy.context.collection.objects.link(key_obj)
    key_obj.location = (9.0, 9.5, 10.0)
    lights.append(key_obj)

    # 2) Soft cool Fill Light from opposite side
    fill_data = bpy.data.lights.new("Fill_Light", type='AREA')
    fill_data.energy = 1150.0
    fill_data.size = 11.0
    fill_data.color = (0.74, 0.86, 1.0)
    fill_obj = bpy.data.objects.new("Fill_Light", fill_data)
    bpy.context.collection.objects.link(fill_obj)
    fill_obj.location = (-10.0, 4.5, 6.5)
    lights.append(fill_obj)

    # 3) Wide overhead studio reflection softbox so polished Chrome fuselage shines smoothly
    rim_data = bpy.data.lights.new("Rim_Light", type='AREA')
    rim_data.energy = 1900.0
    rim_data.size = 16.0
    rim_data.color = (0.94, 0.97, 1.0)
    rim_obj = bpy.data.objects.new("Rim_Light", rim_data)
    bpy.context.collection.objects.link(rim_obj)
    rim_obj.location = (-2.5, 11.5, -6.5)
    lights.append(rim_obj)

    # 4) Under-belly soft bounce fill
    bot_data = bpy.data.lights.new("Bottom_Fill", type='AREA')
    bot_data.energy = 600.0
    bot_data.size = 10.0
    bot_data.color = (0.68, 0.80, 0.96)
    bot_obj = bpy.data.objects.new("Bottom_Fill", bot_data)
    bpy.context.collection.objects.link(bot_obj)
    bot_obj.location = (3.5, -8.0, 2.0)
    lights.append(bot_obj)

    for l_obj in lights:
        direction = Vector((0.0, 0.15, 0.0)) - l_obj.location
        z_axis = (-direction).normalized()
        up_ref = Vector((0.0, 1.0, 0.0))
        x_axis = up_ref.cross(z_axis).normalized()
        y_axis = z_axis.cross(x_axis).normalized()
        rot_mat = Matrix((x_axis, y_axis, z_axis)).transposed()
        l_obj.matrix_world = Matrix.Translation(l_obj.location) @ rot_mat.to_4x4()

    return lights, world


def create_camera_y_up(name, loc, target, ortho=False, ortho_scale=12.4, lens=46.0):
    cam_data = bpy.data.cameras.new(name)
    if ortho:
        cam_data.type = 'ORTHO'
        cam_data.ortho_scale = ortho_scale
    else:
        cam_data.type = 'PERSP'
        cam_data.lens = lens
    cam_obj = bpy.data.objects.new(name, cam_data)
    bpy.context.collection.objects.link(cam_obj)

    loc_v = Vector(loc)
    tgt_v = Vector(target)
    z_axis = (loc_v - tgt_v).normalized()
    up_ref = Vector((0.0, 1.0, 0.0))
    x_axis = up_ref.cross(z_axis).normalized()
    y_axis = z_axis.cross(x_axis).normalized()
    rot_mat = Matrix((x_axis, y_axis, z_axis)).transposed()
    cam_obj.matrix_world = Matrix.Translation(loc_v) @ rot_mat.to_4x4()
    return cam_obj


def set_paint_factor(paint_mat, rgb):
    for n in paint_mat.node_tree.nodes:
        if n.name == "Paint_Factor_Mix":
            n.inputs['B'].default_value = (rgb[0], rgb[1], rgb[2], 1.0)


def stitch_grid_3x2(img_paths, out_path, cell_w=640, cell_h=360):
    import numpy as np
    grid = np.zeros((cell_h * 2, cell_w * 3, 4), dtype=np.float32)

    for idx, p in enumerate(img_paths):
        img = bpy.data.images.load(p, check_existing=False)
        px = np.array(img.pixels[:], dtype=np.float32).reshape((cell_h, cell_w, 4))
        bpy.data.images.remove(img)
        row = 1 - (idx // 3)
        col = idx % 3
        grid[row * cell_h:(row + 1) * cell_h, col * cell_w:(col + 1) * cell_w, :] = px

    out_img = bpy.data.images.new("PaintsGrid", width=cell_w * 3, height=cell_h * 2, alpha=True)
    out_img.pixels = grid.flatten().tolist()
    out_img.filepath_raw = out_path
    out_img.file_format = 'PNG'
    out_img.save()
    bpy.data.images.remove(out_img)


# ==============================================================================
# MAIN PIPELINE
# ==============================================================================
def main():
    scene = reset_scene()
    ph_paths = ensure_polyhaven_textures()
    mats = bake_all_materials(scene, ph_paths)

    hull_obj = build_hull_object(mats)
    canopy_obj = build_canopy_object(mats, hull_obj)

    t1_obj, g1_obj, m1_obj = build_raygun_assembly(
        idx=1,
        pivot_world=Vector((0.16, 0.16, 4.02)),
        z_tip_world=5.32,
        coil_z_list=[4.60, 4.71, 4.82, 4.93],
        coil_radii=[0.072, 0.102, 0.094, 0.068],
        mats=mats,
        hull_obj=hull_obj,
    )

    t2_obj, g2_obj, m2_obj = build_raygun_assembly(
        idx=2,
        pivot_world=Vector((-0.16, -0.16, 4.02)),
        z_tip_world=5.50,
        coil_z_list=[4.74, 4.85, 4.96, 5.07],
        coil_radii=[0.060, 0.084, 0.078, 0.056],
        mats=mats,
        hull_obj=hull_obj,
    )

    empties = create_contract_empties(hull_obj)

    bpy.context.view_layer.update()
    min_v, max_v, size_v, total_tris, used_mats = compute_scene_stats()
    print(f"[Stats] BBox Min: {tuple(round(x, 3) for x in min_v)}")
    print(f"[Stats] BBox Max: {tuple(round(x, 3) for x in max_v)}")
    print(f"[Stats] BBox Size (X=Span, Y=Height, Z=Length): {tuple(round(x, 3) for x in size_v)}")
    print(f"[Stats] Total Triangles: {total_tris}, Materials ({len(used_mats)}): {used_mats}")

    # Export clean GLB
    p_nodes = mats['Paint'].node_tree.nodes
    p_links = mats['Paint'].node_tree.links
    p_bsdf = next(n for n in p_nodes if n.type == 'BSDF_PRINCIPLED')
    p_tex = p_nodes["BaseColor_Tex"]
    p_mix = p_nodes["Paint_Factor_Mix"]
    p_links.new(p_tex.outputs['Color'], p_bsdf.inputs['Base Color'])

    bpy.ops.export_scene.gltf(
        filepath=GLB_PATH,
        export_format='GLB',
        export_yup=False,
        export_apply=True,
        export_cameras=False,
        export_lights=False,
        export_image_format='AUTO',
    )
    patch_glb_paint_factor(GLB_PATH, [0.65, 0.05, 0.04, 1.0])
    glb_bytes = os.path.getsize(GLB_PATH)
    print(f"[Export] Exported {GLB_PATH} ({glb_bytes / (1024*1024):.2f} MB)")

    p_links.new(p_mix.outputs['Result'], p_bsdf.inputs['Base Color'])

    # Write nodes.json
    node_list = []
    ordered_objs = [
        hull_obj, canopy_obj,
        t1_obj, g1_obj, m1_obj,
        t2_obj, g2_obj, m2_obj,
    ] + empties

    for o in ordered_objs:
        wpos = o.matrix_world.translation
        node_list.append({
            "name": o.name,
            "type": o.type,
            "parent": o.parent.name if o.parent else None,
            "world_position": [round(wpos.x, 4), round(wpos.y, 4), round(wpos.z, 4)],
            "axis": MOVING_AXES.get(o.name, None),
        })

    bbox_dict = {
        "min": [round(min_v.x, 4), round(min_v.y, 4), round(min_v.z, 4)],
        "max": [round(max_v.x, 4), round(max_v.y, 4), round(max_v.z, 4)],
        "size": [round(size_v.x, 4), round(size_v.y, 4), round(size_v.z, 4)],
    }

    nodes_data = {
        "asset": "rocket-fighter-retro",
        "units": "metres",
        "frame": {
            "forward": "+Z",
            "up": "+Y",
            "left": "+X",
        },
        "axes": MOVING_AXES,
        "bounding_box": bbox_dict,
        "bbox": bbox_dict,
        "triangles": total_tris,
        "materials": used_mats,
        "glb_bytes": glb_bytes,
        "paints": PAINTS,
        "nodes": node_list,
    }
    with open(NODES_JSON_PATH, "w") as f:
        json.dump(nodes_data, f, indent=2)
    print(f"[Export] Wrote {NODES_JSON_PATH}")

    # Render final.png, side.png, clay.png, paints.png, and pivot test
    lights, world = setup_space_lighting_and_world(scene)

    cam_34 = create_camera_y_up(
        "Cam_34",
        loc=(8.9, 4.3, 10.1),
        target=(0.0, 0.15, -0.15),
        ortho=False,
        lens=45.0,
    )
    cam_side = create_camera_y_up(
        "Cam_Side",
        loc=(15.0, 0.45, 0.0),
        target=(0.0, 0.45, 0.0),
        ortho=True,
        ortho_scale=12.2,
    )

    scene.render.image_settings.file_format = 'PNG'
    scene.render.resolution_x = 1920
    scene.render.resolution_y = 1080
    scene.cycles.samples = 64
    scene.cycles.use_denoising = True

    # 1) Render final.png
    set_paint_factor(mats['Paint'], PAINTS[0]['linear_rgb'])
    scene.camera = cam_34
    scene.render.filepath = FINAL_PNG_PATH
    bpy.ops.render.render(write_still=True)
    print(f"[Render] Wrote {FINAL_PNG_PATH}")

    # 2) Render side.png
    scene.camera = cam_side
    scene.render.filepath = SIDE_PNG_PATH
    bpy.ops.render.render(write_still=True)
    print(f"[Render] Wrote {SIDE_PNG_PATH}")

    # 3) Render paints.png (3x2 grid)
    scene.camera = cam_34
    scene.render.resolution_x = 640
    scene.render.resolution_y = 360
    scene.cycles.samples = 32
    cell_paths = []
    for idx_p, p_info in enumerate(PAINTS):
        set_paint_factor(mats['Paint'], p_info['linear_rgb'])
        c_path = f"/tmp/paint_cell_{idx_p}.png"
        scene.render.filepath = c_path
        bpy.ops.render.render(write_still=True)
        cell_paths.append(c_path)
    stitch_grid_3x2(cell_paths, PAINTS_PNG_PATH, 640, 360)
    for cp in cell_paths:
        if os.path.exists(cp):
            os.remove(cp)
    set_paint_factor(mats['Paint'], PAINTS[0]['linear_rgb'])
    print(f"[Render] Wrote {PAINTS_PNG_PATH}")

    # 4) Render pivot test (30 deg rotation on moving parts)
    canopy_obj.rotation_euler = Euler((math.radians(-30.0), 0.0, 0.0), 'XYZ')
    t1_obj.rotation_euler = Euler((0.0, math.radians(30.0), 0.0), 'XYZ')
    g1_obj.rotation_euler = Euler((math.radians(-15.0), 0.0, 0.0), 'XYZ')
    t2_obj.rotation_euler = Euler((0.0, math.radians(-30.0), 0.0), 'XYZ')
    g2_obj.rotation_euler = Euler((math.radians(15.0), 0.0, 0.0), 'XYZ')
    scene.render.resolution_x = 1280
    scene.render.resolution_y = 720
    scene.render.filepath = PIVOT_TEST_PNG_PATH
    bpy.ops.render.render(write_still=True)
    for mo in (canopy_obj, t1_obj, g1_obj, t2_obj, g2_obj):
        mo.rotation_euler = Euler((0.0, 0.0, 0.0), 'XYZ')
    print(f"[Render] Wrote {PIVOT_TEST_PNG_PATH}")

    # 5) Render clay.png
    clay_mat = bpy.data.materials.new(name="Clay_Override")
    c_bsdf = next(n for n in clay_mat.node_tree.nodes if n.type == 'BSDF_PRINCIPLED')
    c_bsdf.inputs['Base Color'].default_value = (0.55, 0.56, 0.58, 1.0)
    c_bsdf.inputs['Metallic'].default_value = 0.0
    c_bsdf.inputs['Roughness'].default_value = 0.65

    scene.render.resolution_x = 1920
    scene.render.resolution_y = 1080
    scene.cycles.samples = 48
    bpy.context.view_layer.material_override = clay_mat
    scene.camera = cam_34
    scene.render.filepath = CLAY_PNG_PATH
    bpy.ops.render.render(write_still=True)
    bpy.context.view_layer.material_override = None
    print(f"[Render] Wrote {CLAY_PNG_PATH}")


if __name__ == "__main__":
    main()