```python
import bpy, bmesh, random, math, os

random.seed(7)

# ---------- CLEANUP ----------
bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False)
for coll in list(bpy.data.collections):
    bpy.data.collections.remove(coll)

def apply_scale(obj):
    bpy.context.view_layer.objects.active = obj
    bpy.ops.object.transform_apply(location=False, rotation=False, scale=True)

def add_sphere(name, loc, scale, segs=20, rings=12):
    bpy.ops.mesh.primitive_uv_sphere_add(location=loc, segments=segs, ring_count=rings)
    o = bpy.context.object
    o.name = name
    o.scale = scale
    apply_scale(o)
    return o

def add_cube(name, loc, scale, rot=(0,0,0)):
    bpy.ops.mesh.primitive_cube_add(location=loc)
    o = bpy.context.object
    o.name = name
    o.rotation_euler = rot
    o.scale = scale
    bpy.context.view_layer.objects.active = o
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=True)
    return o

def add_cone(name, loc, scale, rot=(0,0,0)):
    bpy.ops.mesh.primitive_cone_add(location=loc, radius1=1.0, radius2=0.0, depth=1.0, vertices=16)
    o = bpy.context.object
    o.name = name
    o.rotation_euler = rot
    o.scale = scale
    bpy.context.view_layer.objects.active = o
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=True)
    return o

def add_cylinder(name, loc, radius, depth, rot=(0,0,0), verts=32):
    bpy.ops.mesh.primitive_cylinder_add(location=loc, radius=radius, depth=depth, vertices=verts)
    o = bpy.context.object
    o.name = name
    o.rotation_euler = rot
    bpy.context.view_layer.objects.active = o
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
    return o

def boolean_op(target, operand, op='UNION', solver='EXACT'):
    mod = target.modifiers.new(name="bool", type='BOOLEAN')
    mod.operation = op
    mod.object = operand
    mod.solver = solver
    bpy.context.view_layer.objects.active = target
    bpy.ops.object.modifier_apply(modifier=mod.name)
    bpy.data.objects.remove(operand, do_unlink=True)

# ---------- HEAD BASE ----------
skull = add_sphere("Skull", (0, 0, 0.38), (0.42, 0.40, 0.48))

jaw = add_cube("Jaw", (0, 0.06, 0.02), (0.30, 0.28, 0.16))
boolean_op(skull, jaw, 'UNION')

brow = add_cube("Brow", (0, 0.33, 0.34), (0.30, 0.08, 0.05))
boolean_op(skull, brow, 'UNION')

nose = add_cone("Nose", (0, 0.40, 0.18), (0.09, 0.24, 0.09), rot=(math.radians(100), 0, 0))
boolean_op(skull, nose, 'UNION')

cheekL = add_sphere("CheekL", (0.22, 0.28, 0.16), (0.15, 0.14, 0.12))
boolean_op(skull, cheekL, 'UNION')
cheekR = add_sphere("CheekR", (-0.22, 0.28, 0.16), (0.15, 0.14, 0.12))
boolean_op(skull, cheekR, 'UNION')

chin = add_sphere("Chin", (0, 0.28, -0.19), (0.14, 0.12, 0.12))
boolean_op(skull, chin, 'UNION')

# eye sockets - subtract (blank antique eye convention)
eyeL = add_sphere("EyeL", (0.14, 0.35, 0.24), (0.09, 0.065, 0.08))
boolean_op(skull, eyeL, 'DIFFERENCE')
eyeR = add_sphere("EyeR", (-0.14, 0.35, 0.24), (0.09, 0.065, 0.08))
boolean_op(skull, eyeR, 'DIFFERENCE')

# eyeballs (blank convex forms sitting back in socket, subtly bulging)
eyeballL = add_sphere("EyeballL", (0.14, 0.33, 0.24), (0.075, 0.06, 0.075))
boolean_op(skull, eyeballL, 'UNION')
eyeballR = add_sphere("EyeballR", (-0.14, 0.33, 0.24), (0.075, 0.06, 0.075))
boolean_op(skull, eyeballR, 'UNION')

# mouth groove
mouth = add_cube("Mouth", (0, 0.37, -0.06), (0.11, 0.025, 0.02))
boolean_op(skull, mouth, 'DIFFERENCE')

# philtrum groove
philtrum = add_cylinder("Philtrum", (0, 0.40, 0.02), 0.012, 0.14, rot=(math.radians(90), 0, 0), verts=8)
boolean_op(skull, philtrum, 'DIFFERENCE')

# ears
earL = add_sphere("EarL", (0.40, 0.0, 0.04), (0.06, 0.12, 0.18))
boolean_op(skull, earL, 'UNION')
earR = add_sphere("EarR", (-0.40, 0.0, 0.04), (0.06, 0.12, 0.18))
boolean_op(skull, earR, 'UNION')

# ---------- VOXEL REMESH to fuse organic base ----------
rem = skull.modifiers.new("Remesh", 'REMESH')
rem.mode = 'VOXEL'
rem.voxel_size = 0.008
rem.use_smooth_shade = True
bpy.context.view_layer.objects.active = skull
bpy.ops.object.modifier_apply(modifier=rem.name)

# ---------- HAIR CLUMPS (sculpted curls, not smooth cap) ----------
hair_positions = []
for row in range(6):
    z = 0.55 + row * 0.055
    n = 10 - row
    for i in range(n):
        a = (i / n) * math.pi * 1.85 - math.pi * 0.05
        r = 0.36 - row * 0.01
        x = math.cos(a) * r
        y = math.sin(a) * r - 0.05
        if y > 0.22:  # keep off the face
            continue
        hair_positions.append((x, y, z))

for idx, (x, y, z) in enumerate(hair_positions):
    jx = random.uniform(-0.015, 0.015)
    jy = random.uniform(-0.015, 0.015)
    jz = random.uniform(-0.01, 0.01)
    s = random.uniform(0.06, 0.085)
    curl = add_sphere(f"Curl{idx}", (x+jx, y+jy, z+jz), (s, s*0.85, s*0.7), segs=10, rings=6)
    curl.rotation_euler = (random.uniform(0, math.pi), random.uniform(0, math.pi), random.uniform(0, math.pi))
    bpy.context.view_layer.objects.active = curl
    bpy.ops.object.transform_apply(location=False, rotation=True, scale=False)
    boolean_op(skull, curl, 'UNION', solver='FAST')

# ---------- NECK / BUST TERMINATION ----------
neck = add_cylinder("Neck", (0, 0.02, -0.28), 0.16, 0.34, verts=28)
boolean_op(skull, neck, 'UNION')

# broken-edge classic bust cut at the bottom (irregular termination)
cut = add_cylinder("Cut", (0, 0.02, -0.50), 0.30, 0.20, verts=28)
boolean_op(skull, cut, 'DIFFERENCE')

head = skull
head.name = "ClassicalBust"

# smooth shading + subsurf for final polish
mod_sub = head.modifiers.new("Subsurf", 'SUBSURF')
mod_sub.levels = 1
mod_sub.render_levels = 2
bpy.context.view_layer.objects.active = head
bpy.ops.object.modifier_apply(modifier=mod_sub.name)
bpy.ops.object.shade_smooth()

# ---------- SOCLE / PEDESTAL ----------
socle = add_cylinder("Socle", (0, 0.02, -0.60), 0.28, 0.10, verts=32)
# small chips at socle edge
for i in range(5):
    ang = random.uniform(0, math.pi*2)
    r = 0.28
    chip = add_sphere(f"Chip{i}", (math.cos(ang)*r, math.sin(ang)*r+0.02, -0.55),
                       (0.03, 0.03, 0.02), segs=8, rings=6)
    boolean_op(socle, chip, 'DIFFERENCE', solver='FAST')

column = add_cylinder("Column", (0, 0.02, -1.0), 0.22, 0.75, verts=32)
plinth = add_cylinder("Plinth", (0, 0.02, -1.42), 0.34, 0.10, verts=32)

for o in (socle, column, plinth):
    bpy.ops.object.shade_smooth()

# ---------- MATERIALS ----------
def build_plaster_material(name, base_rgb):
    mat = bpy.data.materials.new(name)
    mat.use_nodes = True
    nt = mat.node_tree
    nodes, links = nt.nodes, nt.links
    for n in list(nodes):
        nodes.remove(n)

    out = nodes.new("ShaderNodeOutputMaterial")
    out.location = (600, 0)
    bsdf = nodes.new("ShaderNodeBsdfPrincipled")
    bsdf.location = (300, 0)
    bsdf.inputs["Base Color"].default_value = (*base_rgb, 1.0)
    bsdf.inputs["Roughness"].default_value = 0.55
    bsdf.inputs["Subsurface Weight"].default_value = 0.18
    bsdf.inputs["Subsurface Radius"].default_value = (0.012, 0.009, 0.006)
    bsdf.inputs["Subsurface Scale"].default_value = 0.3

    geo = nodes.new("ShaderNodeNewGeometry")
    geo.location = (-600, -200)
    ramp = nodes.new("ShaderNodeValToRGB")
    ramp.location = (-400, -200)
    ramp.color_ramp.elements[0].position = 0.35
    ramp.color_ramp.elements[0].color = (0.55, 0.53, 0.5, 1)
    ramp.color_ramp.elements[1].position = 0.6
    ramp.color_ramp.elements[1].color = (1, 1, 1, 1)
    links.new(geo.outputs["Pointiness"], ramp.inputs["Fac"])

    mix = nodes.new("ShaderNodeMixRGB")
    mix.location = (-150, -100)
    mix.blend_type = 'MULTIPLY'
    mix.inputs[0].default_value = 1.0
    mix.inputs[1].default_value = (*base_rgb, 1.0)
    links.new(ramp.outputs["Color"], mix.inputs[2])
    links.new(mix.outputs["Color"], bsdf.inputs["Base Color"])

    # subtle mold seam lines + micro noise bump
    texcoord = nodes.new("ShaderNodeTexCoord")
    texcoord.location = (-800, 200)
    wave = nodes.new("ShaderNodeTexWave")
    wave.location = (-600, 200)
    wave.inputs["Scale"].default_value = 2.0
    wave.wave_type = 'BANDS'
    wave.bands_direction = 'Z'
    noise = nodes.new("ShaderNodeTexNoise")
    noise.location = (-600, 350)
    noise.inputs["Scale"].default_value = 80.0
    combine = nodes.new("ShaderNodeMixRGB")
    combine.location = (-350, 280)
    combine.blend_type = 'ADD'
    combine.inputs[0].default_value = 0.15
    links.new(wave.outputs["Fac"], combine.inputs[1])
    links.new(noise.outputs["Fac"], combine.inputs[2])

    bump = nodes.new("ShaderNodeBump")
    bump.location = (-100, 200)
    bump.inputs["Strength"].default_value = 0.035
    links.new(combine.outputs["Color"], bump.inputs["Height"])
    links.new(bump.outputs["Normal"], bsdf.inputs["Normal"])

    links.new(bsdf.outputs["BSDF"], out.inputs["Surface"])
    return mat

plaster_head = build_plaster_material("PlasterBust", (0.92, 0.90, 0.855))
plaster_pedestal = build_plaster_material("PlasterPedestal", (0.78, 0.77, 0.73))

head.data.materials.append(plaster_head)
for o in (socle, column, plinth):
    o.data.materials.append(plaster_pedestal)

# ---------- BACKDROP ----------
bpy.ops.mesh.primitive_plane_add(size=8, location=(0, -1.6, 0.5))
backdrop = bpy.context.object
backdrop.name = "Backdrop"
backdrop.rotation_euler = (math.radians(90), 0, 0)
bpy.ops.mesh.primitive_plane_add(size=8, location=(0, 0, -1.55))
floor = bpy.context.object
floor.name = "Floor"

bmat = bpy.data.materials.new("StudioGrey")
bmat.use_nodes = True
bmat.node_tree.nodes["Principled BSDF"].inputs["Base Color"].default_value = (0.16, 0.16, 0.17, 1)
bmat.node_tree.nodes["Principled BSDF"].inputs["Roughness"].default_value = 0.9
backdrop.data.materials.append(bmat)
floor.data.materials.append(bmat)

# ---------- WORLD ----------
world = bpy.data.worlds.new("World")
bpy.context.scene.world = world
world.use_nodes = True
world.node_tree.nodes["Background"].inputs[0].default_value = (0.02, 0.02, 0.02, 1)
world.node_tree.nodes["Background"].inputs[1].default_value = 0.3

# ---------- LIGHTING (classical atelier raking light) ----------
def add_area(name, loc, rot, size, power, color=(1,1,1)):
    bpy.ops.object.light_add(type='AREA', location=loc, rotation=rot)
    l = bpy.context.object
    l.name = name
    l.data.size = size
    l.data.energy = power
    l.data.color = color
    return l

key = add_area("Key", (1.6, -1.0, 1.1), (math.radians(65), 0, math.radians(58)), 1.2, 900)
fill = add_area("Fill", (-1.5, -0.8, 0.6), (math.radians(70), 0, math.radians(-55)), 1.8, 180)
rim = add_area("Rim", (0.3, 1.6, 1.3), (math.radians(-100), 0, math.radians(10)), 1.0, 250)

# ---------- CAMERA ----------
cam_loc = (1.05, -1.55, 0.62)
target = (0, 0.05, 0.42)
bpy.ops.object.camera_add(location=cam_loc)
cam = bpy.context.object
cam.name = "Camera"
direction = (target[0]-cam_loc[0], target[1]-cam_loc[1], target[2]-cam_loc[2])
rot_quat = None
import mathutils
vec = mathutils.Vector(direction)
rot_quat = vec.to_track_quat('-Z', 'Y')
cam.rotation_euler = rot_quat.to_euler()
cam.data.lens = 85
bpy.context.scene.camera = cam

# ---------- RENDER SETTINGS ----------
scene = bpy.context.scene
scene.render.engine = 'CYCLES'
scene.cycles.samples = 64
scene.render.resolution_x = 1536
scene.render.resolution_y = 1536
scene.view_settings.view_transform = 'Filmic'
try:
    scene.cycles.device = 'GPU'
except Exception:
    pass

# ---------- SAVE + EXPORT ----------
out_dir = os.path.dirname(bpy.data.filepath) or "/tmp"
blend_path = os.path.join(out_dir, "classical_bust.blend")
glb_path = os.path.join(out_dir, "classical_bust.glb")

bpy.ops.wm.save_as_mainfile(filepath=blend_path)
try:
    bpy.ops.export_scene.gltf(filepath=glb_path, export_format='GLB', use_selection=False)
except Exception as e:
    print("GLB export failed:", e)
```