#include "vde/brep/step_import.h"
#include "vde/curves/nurbs_curve.h"
#include "vde/curves/nurbs_surface.h"
#include <fstream>
#include <sstream>
#include <optional>
#include <unordered_map>
#include <unordered_set>
#include <variant>
#include <cmath>
#include <cctype>
#include <algorithm>
#include <stdexcept>

namespace vde::brep {

using core::Point3D;
using core::Vector3D;

// ─────────────────────────────────────────────────────────────
//  Error state (thread-safe enough for single-threaded use)
// ─────────────────────────────────────────────────────────────
static StepError g_last_error = StepError::Ok;
static std::string g_last_error_msg;

StepError step_last_error() { return g_last_error; }
const std::string& step_last_error_message() { return g_last_error_msg; }

static void set_error(StepError e, const std::string& msg) {
    g_last_error = e;
    g_last_error_msg = msg;
}
static void clear_error() { g_last_error = StepError::Ok; g_last_error_msg.clear(); }

// ─────────────────────────────────────────────────────────────
//  Tokenizer
// ─────────────────────────────────────────────────────────────
enum class TokenKind : uint8_t {
    Eof, Id, Keyword, String, Number, Integer,
    LParen, RParen, Comma, Semicolon, Equals,
    Star, Enumeration, Omitted,
};

struct Token {
    TokenKind kind = TokenKind::Eof;
    std::string text;  // raw text
    int entity_id = 0; // valid for Id tokens
};

class StepTokenizer {
public:
    explicit StepTokenizer(std::string data) : data_(std::move(data)) {}

    Token peek() {
        if (!peeked_)
            peeked_ = next_token();
        return *peeked_;
    }

    Token consume() {
        Token t = peek();
        peeked_.reset();
        return t;
    }

private:
    std::string data_;
    size_t pos_ = 0;
    std::optional<Token> peeked_;

    void skip_ws() {
        while (pos_ < data_.size() && (data_[pos_] == ' ' || data_[pos_] == '\t'
                                       || data_[pos_] == '\r' || data_[pos_] == '\n'))
            ++pos_;
        // Skip comments: /* ... */
        if (pos_ < data_.size() && data_[pos_] == '/'
            && pos_+1 < data_.size() && data_[pos_+1] == '*') {
            pos_ += 2;
            while (pos_+1 < data_.size() && !(data_[pos_]=='*' && data_[pos_+1]=='/'))
                ++pos_;
            if (pos_+1 < data_.size()) pos_ += 2;
            skip_ws();
        }
    }

    Token next_token() {
        skip_ws();
        if (pos_ >= data_.size()) return Token{TokenKind::Eof};

        char c = data_[pos_];

        // Entity ID: # followed by digits
        if (c == '#') {
            size_t start = pos_ + 1;
            ++pos_;
            while (pos_ < data_.size() && std::isdigit(static_cast<unsigned char>(data_[pos_])))
                ++pos_;
            Token t{TokenKind::Id};
            t.text = data_.substr(start, pos_ - start);
            t.entity_id = std::stoi(t.text);
            return t;
        }

        // Single-quoted string
        if (c == '\'') {
            ++pos_; // skip opening quote
            std::string val;
            while (pos_ < data_.size()) {
                if (data_[pos_] == '\'') {
                    if (pos_+1 < data_.size() && data_[pos_+1] == '\'') {
                        val += '\'';
                        pos_ += 2; // escaped single quote
                    } else {
                        ++pos_; // closing quote
                        break;
                    }
                } else {
                    val += data_[pos_];
                    ++pos_;
                }
            }
            Token t{TokenKind::String};
            t.text = val;
            return t;
        }

        // Operators / punctuation
        switch (c) {
            case '(': ++pos_; return Token{TokenKind::LParen, "("};
            case ')': ++pos_; return Token{TokenKind::RParen, ")"};
            case ',': ++pos_; return Token{TokenKind::Comma, ","};
            case ';': ++pos_; return Token{TokenKind::Semicolon, ";"};
            case '=': ++pos_; return Token{TokenKind::Equals, "="};
            case '*': ++pos_; return Token{TokenKind::Star, "*"};
        }

        // Enumeration: .T. or .F. or .U. etc.
        if (c == '.') {
            size_t start = pos_;
            ++pos_;
            while (pos_ < data_.size() && data_[pos_] != '.') ++pos_;
            if (pos_ < data_.size()) ++pos_; // skip closing dot
            Token t{TokenKind::Enumeration};
            t.text = data_.substr(start, pos_ - start);
            return t;
        }

        // Omitted: $
        if (c == '$') { ++pos_; return Token{TokenKind::Omitted, "$"}; }

        // Number or keyword
        {
            size_t start = pos_;
            bool is_number = (c == '-' || c == '+' || c == '.' || std::isdigit(static_cast<unsigned char>(c)));
            while (pos_ < data_.size()
                   && data_[pos_] != ' ' && data_[pos_] != '\t' && data_[pos_] != '\r' && data_[pos_] != '\n'
                   && data_[pos_] != '(' && data_[pos_] != ')' && data_[pos_] != ',' && data_[pos_] != ';'
                   && data_[pos_] != '=' && data_[pos_] != '\'' && data_[pos_] != '$') {
                ++pos_;
            }
            std::string tok = data_.substr(start, pos_ - start);

            if (is_number && !tok.empty()) {
                // Check if it's a pure integer (no decimal point, no exponent)
                bool is_int = true;
                for (char ch : tok) {
                    if (ch == '.' || ch == 'e' || ch == 'E') { is_int = false; break; }
                }
                if (is_int && tok != "-" && tok != "+") {
                    try {
                        Token t{TokenKind::Integer};
                        t.text = tok;
                        t.entity_id = std::stoi(tok);
                        return t;
                    } catch (...) {}
                }
                Token t{TokenKind::Number};
                t.text = tok;
                return t;
            }

            Token t{TokenKind::Keyword};
            t.text = tok;
            return t;
        }
    }
};

// ─────────────────────────────────────────────────────────────
//  STEP Value type for parsed arguments
// ─────────────────────────────────────────────────────────────
struct StepValue;
using StepList = std::vector<StepValue>;

struct StepValue {
    enum class Type : uint8_t { Omitted, Integer, Float, String, Enum, Ref, List };
    Type type = Type::Omitted;
    int ival = 0;
    double dval = 0.0;
    std::string sval;
    StepList list;
};

// ─────────────────────────────────────────────────────────────
//  Entity record
// ─────────────────────────────────────────────────────────────
struct StepEntity {
    std::string type;
    StepList args;
};

// ─────────────────────────────────────────────────────────────
//  Parsed STEP model (raw entity map)
// ─────────────────────────────────────────────────────────────
class StepParser {
public:
    explicit StepParser(std::string data) : tokenizer_(std::move(data)) {}

    bool parse_data_section() {
        // Find DATA; section
        while (true) {
            Token t = tokenizer_.consume();
            if (t.kind == TokenKind::Eof) {
                set_error(StepError::ParseError, "DATA section not found");
                return false;
            }
            if (t.kind == TokenKind::Keyword && t.text == "DATA") {
                Token s = tokenizer_.consume();
                if (s.kind == TokenKind::Semicolon) break;
            }
        }

        // Parse entities until ENDSEC;
        while (true) {
            Token t = tokenizer_.peek();
            if (t.kind == TokenKind::Keyword && t.text == "ENDSEC") {
                tokenizer_.consume();
                Token s = tokenizer_.consume();
                if (s.kind == TokenKind::Semicolon) break;
                set_error(StepError::ParseError, "Expected ; after ENDSEC");
                return false;
            }
            if (t.kind == TokenKind::Eof) {
                set_error(StepError::ParseError, "Unexpected EOF in DATA section");
                return false;
            }

            if (!parse_entity()) return false;
        }
        return true;
    }

    const std::unordered_map<int, StepEntity>& entities() const { return entities_; }

private:
    StepTokenizer tokenizer_;
    std::unordered_map<int, StepEntity> entities_;

    bool parse_entity() {
        // #ID = TYPE ( args ) ;
        Token id_tok = tokenizer_.consume();
        if (id_tok.kind != TokenKind::Id) {
            set_error(StepError::ParseError, "Expected entity ID, got: " + id_tok.text);
            return false;
        }

        Token eq = tokenizer_.consume();
        if (eq.kind != TokenKind::Equals) {
            set_error(StepError::ParseError, "Expected = after entity ID");
            return false;
        }

        Token type_tok = tokenizer_.consume();
        if (type_tok.kind != TokenKind::Keyword) {
            set_error(StepError::ParseError, "Expected entity type keyword");
            return false;
        }

        Token lp = tokenizer_.consume();
        if (lp.kind != TokenKind::LParen) {
            set_error(StepError::ParseError, "Expected ( after entity type");
            return false;
        }

        StepList args = parse_arg_list();
        // parse_arg_list already consumes the closing )

        Token semi = tokenizer_.consume();
        if (semi.kind != TokenKind::Semicolon) {
            set_error(StepError::ParseError, "Expected ; at end of entity");
            return false;
        }

        entities_[id_tok.entity_id] = {type_tok.text, std::move(args)};
        return true;
    }

    StepList parse_arg_list() {
        StepList args;
        while (true) {
            Token t = tokenizer_.peek();
            if (t.kind == TokenKind::RParen) {
                tokenizer_.consume();
                return args;
            }
            if (t.kind == TokenKind::Comma) {
                tokenizer_.consume();
                continue;
            }
            if (t.kind == TokenKind::LParen) {
                tokenizer_.consume();
                StepValue sv;
                sv.type = StepValue::Type::List;
                sv.list = parse_arg_list();
                args.push_back(std::move(sv));
                continue;
            }

            StepValue sv = parse_value();
            args.push_back(std::move(sv));
        }
    }

    StepValue parse_value() {
        Token t = tokenizer_.consume();
        StepValue sv;
        switch (t.kind) {
            case TokenKind::Omitted:
                sv.type = StepValue::Type::Omitted;
                break;
            case TokenKind::Id:
                sv.type = StepValue::Type::Ref;
                sv.ival = t.entity_id;
                break;
            case TokenKind::Integer:
                sv.type = StepValue::Type::Integer;
                sv.ival = std::stoi(t.text);
                sv.dval = sv.ival;
                break;
            case TokenKind::Number:
                sv.type = StepValue::Type::Float;
                sv.dval = std::stod(t.text);
                break;
            case TokenKind::String:
                sv.type = StepValue::Type::String;
                sv.sval = t.text;
                break;
            case TokenKind::Enumeration:
                sv.type = StepValue::Type::Enum;
                sv.sval = t.text;
                break;
            case TokenKind::Star:
                // * used as a placeholder in some STEP files — treat as Omitted-deriv
                sv.type = StepValue::Type::Omitted;
                break;
            default:
                sv.type = StepValue::Type::Omitted;
                break;
        }
        return sv;
    }
};

// ─────────────────────────────────────────────────────────────
//  Helper: extract values from parsed args
// ─────────────────────────────────────────────────────────────
static double to_double(const StepValue& v) {
    switch (v.type) {
        case StepValue::Type::Integer: return static_cast<double>(v.ival);
        case StepValue::Type::Float:   return v.dval;
        default: return 0.0;
    }
}

static std::string to_string(const StepValue& v) {
    if (v.type == StepValue::Type::String) return v.sval;
    return "";
}

static Point3D to_point(const StepList& lst) {
    if (lst.size() >= 3)
        return Point3D(to_double(lst[0]), to_double(lst[1]), to_double(lst[2]));
    return Point3D(0, 0, 0);
}

static Vector3D to_vector(const StepList& lst) {
    if (lst.size() >= 3)
        return Vector3D(to_double(lst[0]), to_double(lst[1]), to_double(lst[2]));
    return Vector3D(0, 0, 1);
}

static int to_ref(const StepValue& v) { return (v.type == StepValue::Type::Ref) ? v.ival : -1; }

static bool is_omitted(const StepValue& v) { return v.type == StepValue::Type::Omitted; }

// ─────────────────────────────────────────────────────────────
//  AXIS2_PLACEMENT_3D → coordinate frame
// ─────────────────────────────────────────────────────────────
struct Placement3D {
    Point3D origin = Point3D(0, 0, 0);
    Vector3D axis = Vector3D(0, 0, 1);   // Z-axis (primary)
    Vector3D ref_dir = Vector3D(1, 0, 0); // X-axis (secondary, optional)
};

static Placement3D decode_axis2_placement_3d(const StepEntity& ent) {
    Placement3D p;
    // Args: name, location (CARTESIAN_POINT ref), axis (DIRECTION ref, optional), ref_direction (DIRECTION ref, optional)
    if (ent.args.size() >= 2) {
        p.origin = to_point({}); // will be resolved later via reference
        // Store refs for later resolution
    }
    return p;
}

// For internal use: we store placement's refs
struct ResolvedPlacement {
    Point3D origin = Point3D(0, 0, 0);
    Vector3D z_axis = Vector3D(0, 0, 1);
    Vector3D x_axis = Vector3D(1, 0, 0);
    Vector3D y_axis = Vector3D(0, 1, 0);
};

// ─────────────────────────────────────────────────────────────
//  Main converter: entity map → BrepModel(s)
// ─────────────────────────────────────────────────────────────
class StepToBrep {
public:
    StepToBrep(const std::unordered_map<int, StepEntity>& entities) : entities_(entities) {}

    std::vector<BrepModel> convert() {
        std::vector<BrepModel> result;
        resolve_all_points();
        resolve_all_directions();
        resolve_all_placements();

        // Find top-level products or directly find MANIFOLD_SOLID_BREP
        std::vector<int> root_entities;

        // Try finding PRODUCT entities first
        for (const auto& [id, ent] : entities_) {
            if (ent.type == "PRODUCT" || ent.type == "PRODUCT_DEFINITION") {
                root_entities.push_back(id);
            }
        }

        // If no PRODUCT, find SHAPE_DEFINITION_REPRESENTATION or directly MANIFOLD_SOLID_BREP
        if (root_entities.empty()) {
            for (const auto& [id, ent] : entities_) {
                if (ent.type == "SHAPE_DEFINITION_REPRESENTATION") {
                    root_entities.push_back(id);
                }
            }
        }

        // Fallback: find all MANIFOLD_SOLID_BREP and SHELL_BASED_SURFACE_MODEL
        std::vector<int> solid_ids;
        for (const auto& [id, ent] : entities_) {
            if (ent.type == "MANIFOLD_SOLID_BREP" || ent.type == "SHELL_BASED_SURFACE_MODEL") {
                solid_ids.push_back(id);
            }
        }

        if (!solid_ids.empty()) {
            for (int sid : solid_ids) {
                BrepModel body = convert_solid(sid);
                if (body.num_faces() > 0) {
                    // Try to find product name
                    std::string name = find_product_name(sid);
                    result.push_back(std::move(body));
                }
            }
        } else if (!root_entities.empty()) {
            // Walk PRODUCT → PRODUCT_DEFINITION_FORMATION → ...
            for (int rid : root_entities) {
                auto bodies = walk_product_tree(rid);
                for (auto& b : bodies)
                    result.push_back(std::move(b));
            }
        }

        return result;
    }

private:
    const std::unordered_map<int, StepEntity>& entities_;

    // Caches
    std::unordered_map<int, Point3D> points_;
    std::unordered_map<int, Vector3D> directions_;
    std::unordered_map<int, ResolvedPlacement> placements_;
    std::unordered_map<int, Point3D> vectors_; // VECTOR entities
    std::unordered_map<int, curves::NurbsCurve> curves_;
    std::unordered_map<int, curves::NurbsSurface> surfaces_;

    // ── Resolution helpers ──

    void resolve_all_points() {
        for (const auto& [id, ent] : entities_) {
            if (ent.type == "CARTESIAN_POINT" && ent.args.size() >= 1) {
                // Args: name (string), coordinates (list of 3)
                if (ent.args.size() >= 2 && ent.args[1].type == StepValue::Type::List) {
                    points_[id] = to_point(ent.args[1].list);
                }
            }
        }
    }

    void resolve_all_directions() {
        for (const auto& [id, ent] : entities_) {
            if (ent.type == "DIRECTION" && ent.args.size() >= 1) {
                // Args: name, direction_ratios (list of 2 or 3)
                // The list might be args[0] or args[1] depending on schema
                const StepList* ratios = nullptr;
                if (ent.args[0].type == StepValue::Type::List) ratios = &ent.args[0].list;
                else if (ent.args.size() >= 2 && ent.args[1].type == StepValue::Type::List) ratios = &ent.args[1].list;

                if (ratios && ratios->size() >= 3) {
                    directions_[id] = to_vector(*ratios);
                } else if (ratios && ratios->size() == 2) {
                    directions_[id] = Vector3D(to_double((*ratios)[0]), to_double((*ratios)[1]), 0.0);
                }
            }
        }
    }

    void resolve_all_placements() {
        for (const auto& [id, ent] : entities_) {
            if (ent.type == "AXIS2_PLACEMENT_3D") {
                ResolvedPlacement rp;
                // Args sequence depends on schema version, but typically:
                // name, location_ref, axis_ref (optional), ref_direction_ref (optional)
                if (ent.args.size() >= 2) {
                    int loc_ref = to_ref(ent.args[1]);
                    if (loc_ref > 0 && points_.count(loc_ref))
                        rp.origin = points_[loc_ref];
                }
                if (ent.args.size() >= 3 && !is_omitted(ent.args[2])) {
                    int ax_ref = to_ref(ent.args[2]);
                    if (ax_ref > 0 && directions_.count(ax_ref))
                        rp.z_axis = directions_[ax_ref].normalized();
                }
                if (ent.args.size() >= 4 && !is_omitted(ent.args[3])) {
                    int rd_ref = to_ref(ent.args[3]);
                    if (rd_ref > 0 && directions_.count(rd_ref))
                        rp.x_axis = directions_[rd_ref].normalized();
                }
                // Compute Y = Z × X
                // If X is not perpendicular to Z, project it
                if (rp.x_axis.dot(rp.z_axis) > 1e-12) {
                    rp.x_axis = (rp.x_axis - rp.x_axis.dot(rp.z_axis) * rp.z_axis).normalized();
                }
                rp.y_axis = rp.z_axis.cross(rp.x_axis).normalized();
                placements_[id] = rp;
            }
        }
    }

    const ResolvedPlacement* get_placement(int ref) {
        if (ref > 0 && placements_.count(ref))
            return &placements_[ref];
        return &default_placement_;
    }

    static ResolvedPlacement default_placement_;

    // ── Curve construction ──

    curves::NurbsCurve build_line_curve(int id) {
        const auto& ent = entities_.at(id);
        // LINE args: name, pnt (CARTESIAN_POINT ref), dir (VECTOR ref)
        if (ent.args.size() < 3) throw std::runtime_error("LINE missing args");

        int pnt_ref = to_ref(ent.args[1]);
        int vec_ref = to_ref(ent.args[2]);

        Point3D p = (points_.count(pnt_ref)) ? points_[pnt_ref] : Point3D(0,0,0);
        // Eagerly resolve VECTOR entity if not yet cached
        Vector3D v = resolve_vector(vec_ref);

        std::vector<Point3D> cps = {p, p + v};
        return curves::NurbsCurve(cps, {0,0,1,1}, {1,1}, 1);
    }

    Vector3D resolve_vector(int vec_ref) {
        if (vectors_.count(vec_ref)) return vectors_[vec_ref];
        if (!entities_.count(vec_ref)) return Vector3D(1, 0, 0);
        const auto& vent = entities_.at(vec_ref);
        // VECTOR: name, orientation (DIRECTION ref), magnitude
        if (vent.args.size() >= 3) {
            int dir_ref = to_ref(vent.args[1]);
            double mag = to_double(vent.args[2]);
            Vector3D d = (directions_.count(dir_ref))
                ? directions_[dir_ref].normalized()
                : Vector3D(1, 0, 0);
            Vector3D v = d * mag;
            vectors_[vec_ref] = v;
            return v;
        }
        return Vector3D(1, 0, 0);
    }

    curves::NurbsCurve build_circle_curve(int id) {
        const auto& ent = entities_.at(id);
        // CIRCLE args: name, placement (AXIS2_PLACEMENT_3D ref), radius
        if (ent.args.size() < 3) throw std::runtime_error("CIRCLE missing args");

        int pl_ref = to_ref(ent.args[1]);
        double R = to_double(ent.args[2]);
        const auto* pl = get_placement(pl_ref);

        return make_circle_nurbs(pl->origin, pl->x_axis, pl->y_axis, R);
    }

    curves::NurbsCurve build_ellipse_curve(int id) {
        const auto& ent = entities_.at(id);
        // ELLIPSE args: name, placement, semi_axis1, semi_axis2
        if (ent.args.size() < 4) throw std::runtime_error("ELLIPSE missing args");

        int pl_ref = to_ref(ent.args[1]);
        double r1 = to_double(ent.args[2]);
        double r2 = to_double(ent.args[3]);
        const auto* pl = get_placement(pl_ref);

        return make_ellipse_nurbs(pl->origin, pl->x_axis, pl->y_axis, r1, r2);
    }

    curves::NurbsCurve build_bspline_curve(int id) {
        const auto& ent = entities_.at(id);
        // B_SPLINE_CURVE_WITH_KNOTS args:
        // degree, control_points_list, curve_form, closed_curve, self_intersect,
        // knot_multiplicities, knots, knot_spec
        // Optional: weights_data (last arg)
        if (ent.args.size() < 7) throw std::runtime_error("B_SPLINE_CURVE_WITH_KNOTS missing args");

        int degree = static_cast<int>(to_double(ent.args[0]));

        // Parse control points (ent.args[1] is a list of CARTESIAN_POINT refs)
        std::vector<Point3D> cps;
        if (ent.args[1].type == StepValue::Type::List) {
            for (const auto& v : ent.args[1].list) {
                int ref = to_ref(v);
                if (ref > 0 && points_.count(ref))
                    cps.push_back(points_[ref]);
            }
        }

        // Parse knot multiplicities (ent.args[5])
        std::vector<int> mults;
        if (ent.args[5].type == StepValue::Type::List) {
            for (const auto& v : ent.args[5].list)
                mults.push_back(static_cast<int>(to_double(v)));
        }

        // Parse knots (ent.args[6])
        std::vector<double> uknots;
        if (ent.args[6].type == StepValue::Type::List) {
            for (const auto& v : ent.args[6].list)
                uknots.push_back(to_double(v));
        }

        // Expand knot multiplicities into full knot vector
        std::vector<double> knots;
        for (size_t i = 0; i < uknots.size() && i < mults.size(); ++i) {
            for (int j = 0; j < mults[i]; ++j)
                knots.push_back(uknots[i]);
        }

        // Parse weights if provided (last argument)
        std::vector<double> weights;
        if (ent.args.size() > 7 && ent.args.back().type == StepValue::Type::List
            && !is_omitted(ent.args.back())) {
            for (const auto& v : ent.args.back().list)
                weights.push_back(to_double(v));
        }

        if (weights.empty())
            weights.assign(cps.size(), 1.0);

        // The STEP B-spline may store control points in a different convention
        // Enforce degree consistency
        if (knots.empty()) {
            int n = static_cast<int>(cps.size());
            for (int i = 0; i < n + degree + 1; ++i)
                knots.push_back(static_cast<double>(i));
        }

        return curves::NurbsCurve(cps, knots, weights, degree);
    }

    curves::NurbsCurve build_vector_entity(int id) {
        const auto& ent = entities_.at(id);
        // VECTOR: name, orientation (DIRECTION ref), magnitude
        if (ent.args.size() < 3) throw std::runtime_error("VECTOR missing args");
        int dir_ref = to_ref(ent.args[1]);
        double mag = to_double(ent.args[2]);

        Vector3D d = (directions_.count(dir_ref)) ? directions_[dir_ref] : Vector3D(1,0,0);
        d = d.normalized() * mag;

        Point3D o(0,0,0);
        vectors_[id] = d;
        return curves::NurbsCurve({o, o+d}, {0,0,1,1}, {1,1}, 1);
    }

    // ── Surface construction ──

    curves::NurbsSurface build_plane_surface(int id) {
        const auto& ent = entities_.at(id);
        // PLANE: name, position (AXIS2_PLACEMENT_3D ref)
        if (ent.args.size() < 2) throw std::runtime_error("PLANE missing args");

        int pl_ref = to_ref(ent.args[1]);
        const auto* pl = get_placement(pl_ref);

        // Create a large plane as a degree-1×1 NURBS surface
        double s = 1e6; // large extent
        std::vector<std::vector<Point3D>> grid(2);
        grid[0].resize(2); grid[1].resize(2);
        grid[0][0] = pl->origin - s*pl->x_axis - s*pl->y_axis;
        grid[0][1] = pl->origin - s*pl->x_axis + s*pl->y_axis;
        grid[1][0] = pl->origin + s*pl->x_axis - s*pl->y_axis;
        grid[1][1] = pl->origin + s*pl->x_axis + s*pl->y_axis;

        return curves::NurbsSurface(grid, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
    }

    curves::NurbsSurface build_cylindrical_surface(int id) {
        const auto& ent = entities_.at(id);
        // CYLINDRICAL_SURFACE: name, position (AXIS2_PLACEMENT_3D ref), radius
        if (ent.args.size() < 3) throw std::runtime_error("CYLINDRICAL_SURFACE missing args");

        int pl_ref = to_ref(ent.args[1]);
        double R = to_double(ent.args[2]);
        const auto* pl = get_placement(pl_ref);

        return make_cylinder_surface(pl->origin, pl->z_axis, pl->x_axis, pl->y_axis, R);
    }

    curves::NurbsSurface build_conical_surface(int id) {
        const auto& ent = entities_.at(id);
        // CONICAL_SURFACE: name, position, radius, semi_angle
        if (ent.args.size() < 4) throw std::runtime_error("CONICAL_SURFACE missing args");

        int pl_ref = to_ref(ent.args[1]);
        double R = to_double(ent.args[2]);
        double angle = to_double(ent.args[3]); // semi-angle in radians
        const auto* pl = get_placement(pl_ref);

        return make_cone_surface(pl->origin, pl->z_axis, pl->x_axis, pl->y_axis, R, angle);
    }

    curves::NurbsSurface build_spherical_surface(int id) {
        const auto& ent = entities_.at(id);
        // SPHERICAL_SURFACE: name, position (AXIS2_PLACEMENT_3D ref), radius
        if (ent.args.size() < 3) throw std::runtime_error("SPHERICAL_SURFACE missing args");

        int pl_ref = to_ref(ent.args[1]);
        double R = to_double(ent.args[2]);
        const auto* pl = get_placement(pl_ref);

        return make_sphere_surface(pl->origin, pl->z_axis, pl->x_axis, pl->y_axis, R);
    }

    curves::NurbsSurface build_toroidal_surface(int id) {
        const auto& ent = entities_.at(id);
        // TOROIDAL_SURFACE: name, position, major_radius, minor_radius
        if (ent.args.size() < 4) throw std::runtime_error("TOROIDAL_SURFACE missing args");

        int pl_ref = to_ref(ent.args[1]);
        double R = to_double(ent.args[2]);  // major radius
        double r = to_double(ent.args[3]);  // minor radius
        const auto* pl = get_placement(pl_ref);

        return make_torus_surface(pl->origin, pl->z_axis, pl->x_axis, pl->y_axis, R, r);
    }

    curves::NurbsSurface build_bspline_surface(int id) {
        const auto& ent = entities_.at(id);
        // B_SPLINE_SURFACE_WITH_KNOTS args:
        // u_degree, v_degree, control_points_list (list of lists), surface_form,
        // u_closed, v_closed, self_intersect,
        // u_mults, v_mults (?), u_knots, v_knots, knot_spec
        // Arguments: 0=udeg, 1=vdeg, 2=cps, 3=form, 4=uclosed, 5=vclosed, 6=selfint,
        //            7=umults, 8=vmults, 9=uknots, 10=vknots, 11=spec
        if (ent.args.size() < 11) throw std::runtime_error("B_SPLINE_SURFACE_WITH_KNOTS missing args");

        int du = static_cast<int>(to_double(ent.args[0]));
        int dv = static_cast<int>(to_double(ent.args[1]));

        // Control points: list of lists of CARTESIAN_POINT refs
        std::vector<std::vector<Point3D>> grid;
        if (ent.args[2].type == StepValue::Type::List) {
            for (const auto& row_val : ent.args[2].list) {
                if (row_val.type == StepValue::Type::List) {
                    std::vector<Point3D> row;
                    for (const auto& pv : row_val.list) {
                        int ref = to_ref(pv);
                        row.push_back((ref > 0 && points_.count(ref)) ? points_[ref] : Point3D(0,0,0));
                    }
                    grid.push_back(std::move(row));
                }
            }
        }

        // Knots
        auto parse_knot_list = [&](const StepValue& mults_v, const StepValue& knots_v) -> std::vector<double> {
            std::vector<int> mults;
            std::vector<double> uknots;
            if (mults_v.type == StepValue::Type::List)
                for (const auto& v : mults_v.list) mults.push_back(static_cast<int>(to_double(v)));
            if (knots_v.type == StepValue::Type::List)
                for (const auto& v : knots_v.list) uknots.push_back(to_double(v));
            std::vector<double> result;
            for (size_t i = 0; i < uknots.size() && i < mults.size(); ++i)
                for (int j = 0; j < mults[i]; ++j) result.push_back(uknots[i]);
            return result;
        };

        std::vector<double> ku = parse_knot_list(ent.args[7], ent.args[9]);
        std::vector<double> kv = parse_knot_list(ent.args[8], ent.args[10]);

        // Weights (optional)
        std::vector<std::vector<double>> wg;
        if (ent.args.size() > 12 && ent.args[12].type == StepValue::Type::List) {
            for (const auto& row_val : ent.args[12].list) {
                if (row_val.type == StepValue::Type::List) {
                    std::vector<double> wrow;
                    for (const auto& v : row_val.list) wrow.push_back(to_double(v));
                    wg.push_back(std::move(wrow));
                }
            }
        }

        return curves::NurbsSurface(grid, ku, kv, wg, du, dv);
    }

    // ── NURBS geometry constructors ──

    static curves::NurbsCurve make_circle_nurbs(const Point3D& C,
                                                  const Vector3D& X,
                                                  const Vector3D& Y, double R) {
        double w = std::sqrt(2.0) / 2.0; // cos(π/4)
        std::vector<Point3D> cps = {
            C + R*X,
            C + R*X + R*Y,
            C + R*Y,
            C - R*X + R*Y,
            C - R*X,
            C - R*X - R*Y,
            C - R*Y,
            C + R*X - R*Y,
            C + R*X, // closed
        };
        std::vector<double> wgs = {1, w, 1, w, 1, w, 1, w, 1};
        std::vector<double> kts;
        for (int i = 0; i < 12; ++i) kts.push_back(i / 4.0); // [0,0,0, 0.25,0.25, 0.5,0.5, 0.75,0.75, 1,1,1] — 12 entries
        return curves::NurbsCurve(cps, kts, wgs, 2);
    }

    static curves::NurbsCurve make_ellipse_nurbs(const Point3D& C,
                                                   const Vector3D& X,
                                                   const Vector3D& Y, double r1, double r2) {
        double w = std::sqrt(2.0) / 2.0;
        std::vector<Point3D> cps = {
            C + r1*X,
            C + r1*X + r2*Y,
            C + r2*Y,
            C - r1*X + r2*Y,
            C - r1*X,
            C - r1*X - r2*Y,
            C - r2*Y,
            C + r1*X - r2*Y,
            C + r1*X,
        };
        std::vector<double> kts;
        for (int i = 0; i < 12; ++i) kts.push_back(i / 4.0);
        return curves::NurbsCurve(cps, kts, {1,w,1,w,1,w,1,w,1}, 2);
    }

    static curves::NurbsSurface make_cylinder_surface(const Point3D& C,
                                                        const Vector3D& Z,
                                                        const Vector3D& X,
                                                        const Vector3D& Y, double R) {
        // Build the circular profile in the (X,Y) plane
        double w = std::sqrt(2.0) / 2.0;
        std::vector<Point3D> circ_cps = {
            C + R*X,        C + R*X + R*Y,   C + R*Y,
            C - R*X + R*Y,  C - R*X,         C - R*X - R*Y,
            C - R*Y,        C + R*X - R*Y,   C + R*X,
        };

        // Build control grid: 2 rows in v (Z direction), 9 cols in u (circle)
        // Bottom row: circ_cps at z=0, Top row: circ_cps + large extent in Z
        double h = 1e6;
        std::vector<std::vector<Point3D>> grid(2);
        grid[0] = circ_cps;
        for (const auto& p : circ_cps) grid[1].push_back(p + h * Z);

        std::vector<double> ku;
        for (int i = 0; i < 12; ++i) ku.push_back(i / 4.0);
        std::vector<double> kv = {0, 0, 1, 1};

        std::vector<std::vector<double>> weights(2);
        const std::vector<double> cw = {1, w, 1, w, 1, w, 1, w, 1};
        weights[0] = cw;
        weights[1] = cw;

        return curves::NurbsSurface(grid, ku, kv, weights, 2, 1);
    }

    static curves::NurbsSurface make_cone_surface(const Point3D& C,
                                                    const Vector3D& Z,
                                                    const Vector3D& X,
                                                    const Vector3D& Y,
                                                    double R, double angle) {
        double h = 1e6;
        double R_top = R + h * std::tan(angle);
        double w = std::sqrt(2.0) / 2.0;

        auto circle_at = [&](const Point3D& ctr, double r) -> std::vector<Point3D> {
            return {
                ctr + r*X,   ctr + r*X + r*Y, ctr + r*Y,
                ctr - r*X + r*Y, ctr - r*X,    ctr - r*X - r*Y,
                ctr - r*Y,   ctr + r*X - r*Y,  ctr + r*X,
            };
        };

        std::vector<std::vector<Point3D>> grid = {
            circle_at(C, R),
            circle_at(C + h * Z, R_top),
        };

        std::vector<double> ku;
        for (int i = 0; i < 12; ++i) ku.push_back(i / 4.0);
        std::vector<double> kv = {0, 0, 1, 1};
        std::vector<double> cw = {1, w, 1, w, 1, w, 1, w, 1};
        return curves::NurbsSurface(grid, ku, kv, {cw, cw}, 2, 1);
    }

    static curves::NurbsSurface make_sphere_surface(const Point3D& C,
                                                      const Vector3D& Z,
                                                      const Vector3D& X,
                                                      const Vector3D& Y,
                                                      double R) {
        // A sphere as a NURBS surface of revolution about Z
        // Profile: a semi-circle (degree 2, 3 control points) in the XZ plane
        double w = std::sqrt(2.0) / 2.0;
        // Semi-circle profile points: from +Z pole, through equator (+X), to -Z pole
        // Actually for sphere of revolution, profile is a half-circle in XZ plane
        // The half circle: start at (+R, 0, 0), arc through (R, 0, R)*w to (0,0,R)
        // More standard: a rational Bezier half-circle of degree 2:
        // P0=(0,0,R) w=1, P1=(R,0,R)*w(?) -- wait let me use the right construction

        // Sphere as surface of revolution of half-circle profile about Z:
        // We construct 3 rows of control points with the equator as the middle row
        // and degenerate poles at top/bottom.
        std::vector<std::vector<Point3D>> grid(3);
        for (int u = 0; u < 9; ++u) {
            double theta = u * (M_PI / 4.0); // 0, π/4, π/2, ..., 2π
            double cx = std::cos(theta), sy = std::sin(theta);
            Vector3D rad_dir = cx * X + sy * Y;

            // Top row (v=0): degenerate at north pole
            grid[0].push_back(C + R * Z);
            // Middle row (v=1): equator ring
            grid[1].push_back(C + R * rad_dir);
            // Bottom row (v=2): degenerate at south pole
            grid[2].push_back(C - R * Z);
        }

        // Weights: for proper sphere, each row uses circle arc weights in u
        std::vector<std::vector<double>> wg = {
            {1, w, 1, w, 1, w, 1, w, 1},
            {1, w, 1, w, 1, w, 1, w, 1},
            {1, w, 1, w, 1, w, 1, w, 1},
        };

        // Knot vectors
        std::vector<double> ku;
        for (int i = 0; i < 12; ++i) ku.push_back(i / 4.0);
        std::vector<double> kv = {0, 0, 0, 1, 1, 1};

        return curves::NurbsSurface(grid, ku, kv, wg, 2, 2);
    }

    static curves::NurbsSurface make_torus_surface(const Point3D& C,
                                                     const Vector3D& Z,
                                                     const Vector3D& X,
                                                     const Vector3D& Y,
                                                     double R, double r) {
        // Torus: circle of radius r swept around circle of radius R
        // Control grid: 9×9 for proper torus (degree 2×2)
        double w = std::sqrt(2.0) / 2.0;
        std::vector<std::vector<Point3D>> grid(9);
        std::vector<std::vector<double>> wg(9);

        // Sweep minor circle around major circle
        for (int u = 0; u < 9; ++u) {
            grid[u].resize(9);
            wg[u].resize(9);
            double theta = u * (M_PI / 4.0); // 0 to 2π
            double ct = std::cos(theta), st = std::sin(theta);
            Point3D mc = C + R * (ct * X + st * Y); // major circle center

            for (int v = 0; v < 9; ++v) {
                double phi = v * (M_PI / 4.0);
                double cp = std::cos(phi), sp = std::sin(phi);
                Vector3D rad = ct * X + st * Y;
                grid[u][v] = mc + r * (cp * rad + sp * Z);

                // Weight = product of circle arc weights for u and v directions
                double wu = (u % 2 == 0) ? 1.0 : w;
                double wv = (v % 2 == 0) ? 1.0 : w;
                wg[u][v] = wu * wv;
            }
        }

        std::vector<double> ku, kv;
        for (int i = 0; i < 12; ++i) { ku.push_back(i / 4.0); kv.push_back(i / 4.0); }

        return curves::NurbsSurface(grid, ku, kv, wg, 2, 2);
    }

    // ── Topological entity conversion ──

    BrepModel convert_solid(int solid_id) {
        BrepModel model;

        const auto& solid_ent = entities_.at(solid_id);
        // MANIFOLD_SOLID_BREP args: name, outer (CLOSED_SHELL ref)
        // SHELL_BASED_SURFACE_MODEL args: name, sbsm_boundary (SET of SHELL refs)
        if (solid_ent.args.size() < 2) return model;

        std::string solid_name;
        if (!is_omitted(solid_ent.args[0]))
            solid_name = to_string(solid_ent.args[0]);

        // Look for shell refs — might be a single ref or a list of refs
        std::vector<int> shell_refs;
        if (solid_ent.args[1].type == StepValue::Type::Ref) {
            shell_refs.push_back(to_ref(solid_ent.args[1]));
        } else if (solid_ent.args[1].type == StepValue::Type::List) {
            for (const auto& sv : solid_ent.args[1].list) {
                int r = to_ref(sv);
                if (r > 0) shell_refs.push_back(r);
            }
        }

        // Also scan all args for Ref-type shells as fallback
        if (shell_refs.empty()) {
            for (const auto& arg : solid_ent.args) {
                int r = to_ref(arg);
                if (r > 0 && entities_.count(r)) {
                    const auto& e = entities_.at(r);
                    if (e.type == "CLOSED_SHELL" || e.type == "OPEN_SHELL")
                        shell_refs.push_back(r);
                }
            }
        }

        std::vector<int> all_face_ids;
        for (int sh_ref : shell_refs) {
            auto face_ids = convert_shell(sh_ref, model);
            all_face_ids.insert(all_face_ids.end(), face_ids.begin(), face_ids.end());
        }

        if (!all_face_ids.empty()) {
            int shell_id = model.add_shell(all_face_ids, true);
            model.add_body({shell_id}, solid_name);
        }

        return model;
    }

    std::vector<int> convert_shell(int shell_id, BrepModel& model) {
        std::vector<int> face_ids;
        if (!entities_.count(shell_id)) return face_ids;

        const auto& shell_ent = entities_.at(shell_id);
        // CLOSED_SHELL args: name, cfs_faces (SET of ADVANCED_FACE refs)
        if (shell_ent.args.size() < 2) return face_ids;

        // The faces set is args[1] (or args[0] if no name)
        size_t faces_idx = (shell_ent.args[0].type == StepValue::Type::String) ? 1 : 0;
        if (faces_idx >= shell_ent.args.size()) return face_ids;

        const StepValue& faces_val = shell_ent.args[faces_idx];
        if (faces_val.type != StepValue::Type::List) return face_ids;

        for (const auto& fv : faces_val.list) {
            int face_ref = to_ref(fv);
            if (face_ref > 0) {
                int face_id = convert_face(face_ref, model);
                if (face_id >= 0) face_ids.push_back(face_id);
            }
        }
        return face_ids;
    }

    int convert_face(int face_id, BrepModel& model) {
        if (!entities_.count(face_id)) return -1;
        const auto& face_ent = entities_.at(face_id);

        // ADVANCED_FACE args:
        // 0: name
        // 1: bounds (SET of FACE_BOUND/FACE_OUTER_BOUND refs)
        // 2: face_geometry (SURFACE ref)
        // 3: same_sense (BOOLEAN)

        if (face_ent.args.size() < 3) return -1;

        // Resolve surface
        int surf_ref = to_ref(face_ent.args[2]);
        int surf_id = convert_or_get_surface(surf_ref, model);
        if (surf_id < 0) return -1;

        // Resolve bounds
        std::vector<int> loop_ids;
        size_t bounds_idx = 1;
        if (face_ent.args[bounds_idx].type == StepValue::Type::List) {
            for (const auto& bv : face_ent.args[bounds_idx].list) {
                int bound_ref = to_ref(bv);
                if (bound_ref > 0) {
                    int loop_id = convert_face_bound(bound_ref, model);
                    if (loop_id >= 0) loop_ids.push_back(loop_id);
                }
            }
        }

        if (loop_ids.empty()) return -1;
        return model.add_face(surf_id, loop_ids);
    }

    int convert_face_bound(int bound_id, BrepModel& model) {
        if (!entities_.count(bound_id)) return -1;
        const auto& ent = entities_.at(bound_id);

        // FACE_BOUND / FACE_OUTER_BOUND args:
        // 0: name
        // 1: bound (EDGE_LOOP ref)
        // 2: orientation (BOOLEAN)

        bool is_outer = (ent.type == "FACE_OUTER_BOUND");

        int loop_ref = -1;
        for (const auto& arg : ent.args) {
            int r = to_ref(arg);
            if (r > 0 && entities_.count(r) && entities_.at(r).type == "EDGE_LOOP") {
                loop_ref = r;
                break;
            }
        }
        if (loop_ref < 0) return -1;

        return convert_edge_loop(loop_ref, model, is_outer);
    }

    int convert_edge_loop(int loop_id, BrepModel& model, bool is_outer) {
        if (!entities_.count(loop_id)) return -1;
        const auto& ent = entities_.at(loop_id);

        // EDGE_LOOP args: name, edge_list (SET of ORIENTED_EDGE refs)
        const StepValue* edge_list = nullptr;
        for (const auto& arg : ent.args) {
            if (arg.type == StepValue::Type::List) {
                edge_list = &arg;
                break;
            }
        }
        if (!edge_list) return -1;

        std::vector<int> edge_indices;
        for (const auto& ev : edge_list->list) {
            int oe_ref = to_ref(ev);
            if (oe_ref > 0) {
                int ei = convert_oriented_edge(oe_ref, model);
                if (ei >= 0) edge_indices.push_back(ei);
            }
        }

        if (edge_indices.empty()) return -1;
        return model.add_loop(edge_indices, is_outer);
    }

    int convert_oriented_edge(int oe_id, BrepModel& model) {
        if (!entities_.count(oe_id)) return -1;
        const auto& ent = entities_.at(oe_id);

        // ORIENTED_EDGE args:
        // 0: name, 1: edge_element (EDGE_CURVE ref), 2: orientation
        int edge_ref = -1;
        bool reversed = false;
        for (const auto& arg : ent.args) {
            int r = to_ref(arg);
            if (r > 0 && entities_.count(r) && (entities_.at(r).type == "EDGE_CURVE" || entities_.at(r).type == "EDGE_CURVE_IMPLICIT")) {
                edge_ref = r;
            }
        }
        // Check orientation: arg[2] is T/F boolean
        if (ent.args.size() >= 3 && ent.args[2].type == StepValue::Type::Enum) {
            reversed = (ent.args[2].sval == ".F.");
        }

        if (edge_ref < 0) return -1;
        return convert_edge_curve(edge_ref, model, reversed);
    }

    int convert_edge_curve(int edge_id, BrepModel& model, bool reversed) {
        if (!entities_.count(edge_id)) return -1;
        const auto& ent = entities_.at(edge_id);

        // EDGE_CURVE args:
        // 0: name
        // 1: edge_start (VERTEX_POINT ref)
        // 2: edge_end (VERTEX_POINT ref)
        // 3: edge_geometry (CURVE ref)
        // 4: same_sense (BOOLEAN)

        if (ent.args.size() < 5) return -1;

        int v_s = convert_vertex_point(to_ref(ent.args[1]), model);
        int v_e = convert_vertex_point(to_ref(ent.args[2]), model);

        int curve_ref = to_ref(ent.args[3]);
        curves::NurbsCurve nc = get_or_build_curve(curve_ref);

        return model.add_edge(v_s, v_e, nc);
    }

    int convert_vertex_point(int vp_id, BrepModel& model) {
        if (vp_cache_.count(vp_id)) return vp_cache_[vp_id];
        if (!entities_.count(vp_id)) return -1;

        const auto& ent = entities_.at(vp_id);
        // VERTEX_POINT args: name, vertex_geometry (CARTESIAN_POINT ref)
        int pt_ref = -1;
        for (const auto& arg : ent.args) {
            pt_ref = to_ref(arg);
            if (pt_ref > 0 && points_.count(pt_ref)) break;
        }

        Point3D pt = (pt_ref > 0 && points_.count(pt_ref)) ? points_[pt_ref] : Point3D(0,0,0);
        int vi = model.add_vertex(pt);
        vp_cache_[vp_id] = vi;
        return vi;
    }

    int convert_or_get_surface(int surf_ref, BrepModel& model) {
        if (surf_cache_.count(surf_ref)) return surf_cache_[surf_ref];

        curves::NurbsSurface surf = get_or_build_surface(surf_ref);
        int si = model.add_surface(surf);
        surf_cache_[surf_ref] = si;
        return si;
    }

    curves::NurbsCurve get_or_build_curve(int ref) {
        auto cit = curve_cache_.find(ref);
        if (cit != curve_cache_.end()) return cit->second;
        if (!entities_.count(ref)) {
            // Return a dummy curve
            return curves::NurbsCurve({Point3D(0,0,0), Point3D(1,0,0)}, {0,0,1,1}, {1,1}, 1);
        }

        const auto& ent = entities_.at(ref);
        // Default-construct with dummy line (NurbsCurve has no default ctor)
        curves::NurbsCurve nc({Point3D(0,0,0), Point3D(1,0,0)}, {0,0,1,1}, {1,1}, 1);
        try {
            if (ent.type == "LINE")
                nc = build_line_curve(ref);
            else if (ent.type == "CIRCLE")
                nc = build_circle_curve(ref);
            else if (ent.type == "ELLIPSE")
                nc = build_ellipse_curve(ref);
            else if (ent.type == "B_SPLINE_CURVE_WITH_KNOTS")
                nc = build_bspline_curve(ref);
            else if (ent.type == "VECTOR")
                nc = build_vector_entity(ref);
            else
                nc = curves::NurbsCurve({Point3D(0,0,0), Point3D(1,0,0)}, {0,0,1,1}, {1,1}, 1);
        } catch (const std::exception& e) {
            set_error(StepError::EntityTypeError, "Error building curve #" + std::to_string(ref) + ": " + e.what());
            nc = curves::NurbsCurve({Point3D(0,0,0), Point3D(1,0,0)}, {0,0,1,1}, {1,1}, 1);
        }
        curve_cache_.insert_or_assign(ref, nc);
        return nc;
    }

    curves::NurbsSurface get_or_build_surface(int ref) {
        auto sit = surface_cache_.find(ref);
        if (sit != surface_cache_.end()) return sit->second;
        if (!entities_.count(ref)) {
            // Return a dummy surface
            std::vector<std::vector<Point3D>> g = {{Point3D(0,0,0), Point3D(1,0,0)},
                                                    {Point3D(0,1,0), Point3D(1,1,0)}};
            return curves::NurbsSurface(g, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
        }

        const auto& ent = entities_.at(ref);
        // Default-construct with a dummy surface (NurbsSurface has no default ctor)
        std::vector<std::vector<Point3D>> g0 = {{Point3D(0,0,0), Point3D(1,0,0)},
                                                 {Point3D(0,1,0), Point3D(1,1,0)}};
        curves::NurbsSurface ns(g0, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
        try {
            if (ent.type == "PLANE")
                ns = build_plane_surface(ref);
            else if (ent.type == "CYLINDRICAL_SURFACE")
                ns = build_cylindrical_surface(ref);
            else if (ent.type == "CONICAL_SURFACE")
                ns = build_conical_surface(ref);
            else if (ent.type == "SPHERICAL_SURFACE")
                ns = build_spherical_surface(ref);
            else if (ent.type == "TOROIDAL_SURFACE")
                ns = build_toroidal_surface(ref);
            else if (ent.type == "B_SPLINE_SURFACE_WITH_KNOTS")
                ns = build_bspline_surface(ref);
            else {
                std::vector<std::vector<Point3D>> g = {{Point3D(0,0,0), Point3D(1,0,0)},
                                                        {Point3D(0,1,0), Point3D(1,1,0)}};
                ns = curves::NurbsSurface(g, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
            }
        } catch (const std::exception& e) {
            set_error(StepError::EntityTypeError, "Error building surface #" + std::to_string(ref) + ": " + e.what());
            std::vector<std::vector<Point3D>> g = {{Point3D(0,0,0), Point3D(1,0,0)},
                                                    {Point3D(0,1,0), Point3D(1,1,0)}};
            ns = curves::NurbsSurface(g, {0,0,1,1}, {0,0,1,1}, {}, 1, 1);
        }
        surface_cache_.insert_or_assign(ref, ns);
        return ns;
    }

    // ── Product tree walk (for assembly) ──

    std::string find_product_name(int solid_id) {
        // Walk back through the entity graph to find PRODUCT name
        // This is a heuristic
        for (const auto& [id, ent] : entities_) {
            if (ent.type == "PRODUCT") {
                // PRODUCT args: name, description, id, ...
                if (ent.args.size() >= 1 && ent.args[0].type == StepValue::Type::String)
                    return ent.args[0].sval;
            }
        }
        (void)solid_id;
        return "";
    }

    std::vector<BrepModel> walk_product_tree(int root_id) {
        std::vector<BrepModel> result;
        if (!entities_.count(root_id)) return result;

        const auto& ent = entities_.at(root_id);
        // SHAPE_DEFINITION_REPRESENTATION → find MANIFOLD_SOLID_BREP refs
        std::vector<int> solid_refs;
        find_nested_refs(ent, "MANIFOLD_SOLID_BREP", solid_refs);
        find_nested_refs(ent, "SHELL_BASED_SURFACE_MODEL", solid_refs);

        if (solid_refs.empty()) {
            // Try just iterating all args
            collect_all_refs(ent, solid_refs);
            // Filter to solids
            std::vector<int> filtered;
            for (int r : solid_refs) {
                if (entities_.count(r)) {
                    const auto& e = entities_.at(r);
                    if (e.type == "MANIFOLD_SOLID_BREP" || e.type == "SHELL_BASED_SURFACE_MODEL")
                        filtered.push_back(r);
                }
            }
            solid_refs = filtered;
        }

        for (int sid : solid_refs) {
            auto body = convert_solid(sid);
            if (body.num_faces() > 0) result.push_back(std::move(body));
        }
        return result;
    }

    void find_nested_refs(const StepEntity& ent, const std::string& target_type, std::vector<int>& out) {
        for (const auto& arg : ent.args) {
            if (arg.type == StepValue::Type::Ref) {
                int r = arg.ival;
                if (entities_.count(r)) {
                    if (entities_.at(r).type == target_type)
                        out.push_back(r);
                    else
                        find_nested_refs(entities_.at(r), target_type, out);
                }
            } else if (arg.type == StepValue::Type::List) {
                for (const auto& v : arg.list) {
                    StepEntity dummy;
                    dummy.args = {v};
                    find_nested_refs(dummy, target_type, out);
                }
            }
        }
    }

    void collect_all_refs(const StepEntity& ent, std::vector<int>& out) {
        for (const auto& arg : ent.args) {
            if (arg.type == StepValue::Type::Ref)
                out.push_back(arg.ival);
            else if (arg.type == StepValue::Type::List)
                for (const auto& v : arg.list)
                    if (v.type == StepValue::Type::Ref) out.push_back(v.ival);
        }
    }

    // Cache maps
    std::unordered_map<int, int> vp_cache_;    // vertex_point_id → BrepModel vertex index
    std::unordered_map<int, int> surf_cache_;  // surface entity id → BrepModel surface index
    std::unordered_map<int, curves::NurbsCurve> curve_cache_;
    std::unordered_map<int, curves::NurbsSurface> surface_cache_;
};

ResolvedPlacement StepToBrep::default_placement_;

// ─────────────────────────────────────────────────────────────
//  Public API
// ─────────────────────────────────────────────────────────────

std::vector<BrepModel> import_step_from_string(const std::string& step_data) {
    clear_error();

    if (step_data.empty()) {
        set_error(StepError::ParseError, "Empty STEP data");
        return {};
    }

    StepParser parser(step_data);
    if (!parser.parse_data_section()) {
        return {};
    }

    StepToBrep converter(parser.entities());

    try {
        return converter.convert();
    } catch (const std::exception& e) {
        set_error(StepError::ParseError, std::string("Conversion error: ") + e.what());
        return {};
    }
}

std::vector<BrepModel> import_step(const std::string& filepath) {
    clear_error();

    std::ifstream file(filepath);
    if (!file.is_open()) {
        set_error(StepError::FileNotFound, "Cannot open file: " + filepath);
        return {};
    }

    std::ostringstream buf;
    buf << file.rdbuf();
    return import_step_from_string(buf.str());
}

}  // namespace vde::brep