```python
import bpy
import math
import random
from mathutils import Vector

# A sculptural study in atelier plaster.  All of the visible anatomy and hair
# is geometry; the very fine plaster grain is the only bump-map detail.
random.seed(18)
bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False)

def material(name, color, roughness=0.7):
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    p = m.node_tree.nodes.get("Principled BSDF")
    p.inputs['Base Color'].default_value = (*color, 1)
    p.inputs['Roughness'].default_value = roughness
    return m

plaster = material("Warm white • fine atelier casting plaster", (0.72, .695, .65), .76)
nt = plaster.node_tree
p = nt.nodes.get("Principled BSDF")
p.inputs["Subsurface Weight"].default_value = .045
p.inputs["Subsurface Radius"].default_value = (.55, .38, .23)
noise = nt.nodes.new("ShaderNodeTexNoise")
noise.inputs["Scale"].default_value = 115
noise.inputs["Detail"].default_value = 2.3
bump = nt.nodes.new("ShaderNodeBump")
bump.inputs["Strength"].default_value = .17
bump.inputs["Distance"].default_value = .009
nt.links.new(noise.outputs["Fac"], bump.inputs["Height"])
nt.links.new(bump.outputs["Normal"], p.inputs["Normal"])
grain = nt.nodes.new("ShaderNodeTexNoise")
grain.inputs["Scale"].default_value = 5.5
grain.inputs["Detail"].default_value = 3
ramp = nt.nodes.new("ShaderNodeValToRGB")
ramp.color_ramp.elements[0].position = .15
ramp.color_ramp.elements[0].color = (.62, .60, .56, 1)
ramp.color_ramp.elements[1].position = .85
ramp.color_ramp.elements[1].color = (.78, .755, .715, 1)
nt.links.new(grain.outputs["Fac"], ramp.inputs[0])
nt.links.new(ramp.outputs["Color"], p.inputs["Base Color"])

eye_mat = material("Blank antique eyes • polished plaster", (.735, .716, .676), .56)
eye_mat.node_tree.nodes.get("Principled BSDF").inputs["Subsurface Weight"].default_value = .035
recess = material("Plaster in deep apertures", (.245, .225, .195), .92)
seam_mat = material("Faint casting parting line", (.65, .625, .58), .83)
ped_mat = material("Pedestal • grey limestone", (.265, .256, .237), .82)
floor_mat = material("Neutral warm-grey studio", (.205, .20, .19), .88)

def mesh_obj(name, vertices, faces, mat=plaster, smooth=True):
    mesh = bpy.data.meshes.new(name)
    mesh.from_pydata(vertices, [], faces)
    mesh.update()
    ob = bpy.data.objects.new(name, mesh)
    bpy.context.collection.objects.link(ob)
    if mat:
        mesh.materials.append(mat)
    if smooth:
        for f in mesh.polygons:
            f.use_smooth = True
    return ob

def interp(z, pairs):
    if z <= pairs[0][0]:
        return pairs[0][1]
    if z >= pairs[-1][0]:
        return pairs[-1][1]
    for i in range(len(pairs)-1):
        a, b = pairs[i], pairs[i+1]
        if a[0] <= z <= b[0]:
            h = b[0]-a[0]
            t = (z-a[0])/h
            prev = pairs[max(0, i-1)]
            nex = pairs[min(len(pairs)-1, i+2)]
            m0 = (b[1]-prev[1])/(b[0]-prev[0])
            m1 = (nex[1]-a[1])/(nex[0]-a[0])
            return ((2*t**3-3*t*t+1)*a[1] +
                    (t**3-2*t*t+t)*h*m0 +
                    (-2*t**3+3*t*t)*b[1] +
                    (t**3-t*t)*h*m1)

def gauss(x, z, cx, cz, sx, sz):
    return math.exp(-((x-cx)/sx)**2-((z-cz)/sz)**2)

widths = [(4.79,.03),(4.85,.28),(4.94,.47),(5.08,.66),
          (5.26,.80),(5.48,.90),(5.73,1.00),(5.96,1.065),
          (6.18,1.055),(6.43,1.035),(6.67,1.025),(6.96,1.04),
          (7.19,1.01),(7.39,.89),(7.58,.62),(7.73,.05)]
depths = [(4.79,.33),(4.85,.54),(4.99,.77),(5.20,.84),
          (5.45,.90),(5.72,.95),(5.96,1.015),(6.2,1.035),
          (6.44,1.01),(6.68,.995),(6.98,.94),
          (7.22,.82),(7.42,.62),(7.6,.36),(7.73,.03)]
backs = [(4.79,.38),(5.05,.64),(5.35,.88),(5.7,1.02),
         (6.1,1.10),(6.55,1.15),(6.95,1.13),(7.25,.98),
         (7.48,.73),(7.65,.38),(7.73,.03)]
nose_projection = [(5.59,0),(5.64,.045),(5.69,.21),(5.75,.44),
                   (5.82,.545),(5.91,.51),(6.10,.375),
                   (6.30,.23),(6.47,.12),(6.62,.015),(6.72,0)]
nose_width = [(5.59,.12),(5.69,.16),(5.80,.205),(5.90,.18),
              (6.10,.145),(6.35,.13),(6.6,.15),(6.72,.16)]

def face_y(x, z):
    w = max(.035, interp(z, widths))
    c = math.sqrt(max(0, 1-(min(.99999, abs(x)/w))**2))
    y = -interp(z, depths)*c**.64
    d = 0.
    # Supraorbital rims, frontal eminences and glabella.
    for s in [-1, 1]:
        d -= .130*gauss(x,z,s*.45,6.435,.34,.115)
        d -= .043*gauss(x,z,s*.44,6.82,.43,.28)
        d += .137*gauss(x,z,s*.47,6.225,.29,.175)
        # Zygomatic plane, lower orbital ledge, and buccal hollow.
        d -= .122*gauss(x,z,s*.735,5.995,.255,.18)
        d -= .040*gauss(x,z,s*.48,6.065,.30,.10)
        d += .048*gauss(x,z,s*.73,5.59,.21,.235)
        d -= .026*gauss(x,z,s*.71,5.24,.19,.18)
    d -= .075*gauss(x,z,0,6.43,.17,.21)
    d -= .049*gauss(x,z,0,5.44,.51,.24)
    d -= .075*gauss(x,z,0,5.015,.38,.135)
    d += .031*gauss(x,z,0,5.155,.30,.052)
    # The straight bridge is part of the head surface, not an attached cone.
    if 5.59 < z < 6.72:
        nw = interp(z, nose_width)
        d -= interp(z, nose_projection)*math.exp(-(abs(x)/nw)**2.35)
    for s in [-1,1]:
        d -= .136*gauss(x,z,s*.225,5.759,.107,.079)
        d += .044*gauss(x,z,s*.209,5.694,.073,.028)
        d += .022*gauss(x,z,s*.335,5.742,.031,.076)
        # Philtral pillars.
        d -= .027*gauss(x,z,s*.068,5.527,.031,.103)
        if 5.30 < z < 5.73:
            crease_x = s*(.32+.27*(5.73-z))
            d += .013*gauss(x,z,crease_x,5.52,.036,.22)
    d += .012*gauss(x,z,0,5.535,.035,.087)
    return y+d*c**.45

# Continuous closed craniofacial volume.
verts, faces = [], []
NZ, NA = 270, 256
for j in range(NZ+1):
    z = 4.79+(7.73-4.79)*j/NZ
    w = max(.018, interp(z, widths))
    for i in range(NA):
        a = 2*math.pi*i/NA
        x = w*math.sin(a)
        c = math.cos(a)
        y = face_y(x,z) if c >= 0 else -interp(z,backs)*c
        verts.append((x,y,z))
for j in range(NZ):
    for i in range(NA):
        a=j*NA+i
        b=j*NA+(i+1)%NA
        faces.append((a,b,b+NA,a+NA))
faces.append(tuple(reversed(range(NA))))
faces.append(tuple(NZ*NA+i for i in range(NA)))
head = mesh_obj("01 • continuous craniofacial sculpture", verts, faces)

# A shoulder bust with an architectural flat termination, and a neck
# with broad integrated sternomastoid and collarbone relief.
neck_w = [(3.06,1.36),(3.18,1.57),(3.38,1.71),(3.56,1.62),
          (3.76,1.30),(3.96,.89),(4.16,.62),(4.43,.555),
          (4.78,.59),(5.08,.66),(5.39,.72)]
neck_d = [(3.06,.65),(3.25,.73),(3.48,.75),(3.72,.66),
          (3.96,.56),(4.23,.49),(4.57,.49),(4.9,.56),(5.39,.64)]
v, f = [], []
NN, NC = 142, 160
for j in range(NN+1):
    z=3.06+(5.39-3.06)*j/NN
    w=interp(z,neck_w)
    dep=interp(z,neck_d)
    for i in range(NC):
        a=2*math.pi*i/NC
        x=w*math.sin(a)
        y=.20-dep*math.cos(a)
        if math.cos(a)>0:
            d=0
            for s in [-1,1]:
                line=s*(.24+.34*(z-3.78))
                d-=.076*math.exp(-((x-line)/.105)**2)*math.exp(-((z-4.5)/.83)**4)
                clav_z=3.57+.15*math.sin(min(1,abs(x)/1.3)*math.pi*.88)
                d-=.055*gauss(x,z,s*.72,clav_z,.62,.055)
            d-=.026*gauss(x,z,0,4.41,.15,.18)
            d+=.052*gauss(x,z,0,3.79,.20,.135)
            y+=d*math.cos(a)**.5
        v.append((x,y,z))
for j in range(NN):
    for i in range(NC):
        a=j*NC+i;b=j*NC+(i+1)%NC
        f.append((a,b,b+NC,a+NC))
f.append(tuple(reversed(range(NC))))
f.append(tuple(NN*NC+i for i in range(NC)))
neck=mesh_obj("02 • neck, clavicles and cut shoulders",v,f)

# Join the two large volumes into a single actual sculpt.
bpy.ops.object.select_all(action='DESELECT')
head.select_set(True)
neck.select_set(True)
bpy.context.view_layer.objects.active=head
bpy.ops.object.join()
rm=head.modifiers.new("Unified plaster volume","REMESH")
rm.mode='VOXEL'
rm.voxel_size=.019
rm.use_smooth_shade=True
bpy.ops.object.modifier_apply(modifier=rm.name)
sm=head.modifiers.new("Gentle sculpt finishing","SMOOTH")
sm.factor=.23
sm.iterations=3
bpy.ops.object.modifier_apply(modifier=sm.name)

def tube(name, points, radii, mat=plaster, sides=10):
    pts=[Vector(p) for p in points]
    vs, fs=[],[]
    old_b=None
    for i,pt in enumerate(pts):
        tangent=(pts[min(i+1,len(pts)-1)]-pts[max(i-1,0)]).normalized()
        ref=Vector((0,1,0))
        if abs(tangent.dot(ref))>.92:
            ref=Vector((1,0,0))
        b=tangent.cross(ref).normalized()
        if old_b is not None and b.dot(old_b)<0:
            b=-b
        n=tangent.cross(b).normalized()
        old_b=b
        r=radii[i] if isinstance(radii,(list,tuple)) else radii
        for k in range(sides):
            a=2*math.pi*k/sides
            vs.append(tuple(pt+r*(math.cos(a)*b+math.sin(a)*n)))
    for i in range(len(pts)-1):
        for k in range(sides):
            a=i*sides+k;b=i*sides+(k+1)%sides
            fs.append((a,b,b+sides,a+sides))
    fs.append(tuple(reversed(range(sides))))
    fs.append(tuple((len(pts)-1)*sides+k for k in range(sides)))
    return mesh_obj(name,vs,fs,mat)

# Eyelids bound almond-shaped blank eyeballs, with broad upper lid surfaces
# and incised folds returning into the orbital plane.
for side in [-1,1]:
    def eye_coords(t):
        x=side*(.183+.600*t)
        z=6.207+.020*t
        arch=math.sin(math.pi*t)**.80
        zu=z+.113*arch
        zl=z-.079*arch
        yin=face_y(side*.183,6.207)-.016
        yout=face_y(side*.783,6.227)-.012
        ym=yin*(1-t)+yout*t-.088*math.sin(math.pi*t)**.8
        return x,zu,zl,ym
    vs,fs=[],[]
    nt_, nv=68,20
    for i in range(nt_+1):
        t=i/nt_
        x,zu,zl,ym=eye_coords(t)
        for j in range(nv+1):
            q=j/nv
            z=zl+(zu-zl)*q
            y=ym-.041*math.sin(math.pi*q)*math.sin(math.pi*t)
            vs.append((x,y,z))
    for i in range(nt_):
        for j in range(nv):
            a=i*(nv+1)+j
            fs.append((a,a+1,a+nv+2,a+nv+1))
    mesh_obj(("Left" if side<0 else "Right")+" blank carved eye",vs,fs,eye_mat)

    for upper in [True,False]:
        vs,fs=[],[]
        NU,NB=68,15
        edge_points=[]
        edge_r=[]
        for i in range(NU+1):
            t=i/NU
            x,zu,zl,ym=eye_coords(t)
            zi=zu if upper else zl
            edge_points.append((x,ym-.007,zi))
            edge_r.append((.017 if upper else .012)*math.sin(math.pi*t)**.38+.003)
            for j in range(NB+1):
                s=j/NB
                xx=x+side*(t-.5)*.095*s
                ext=(.119 if upper else .126)*math.sin(math.pi*t)**.42+.014
                zz=zi+(ext*s if upper else -ext*s)
                outer_y=face_y(xx,zz)-.003
                yy=(1-s)*ym+s*outer_y
                yy-=.028*math.sin(math.pi*s)*math.sin(math.pi*t)**.5
                if upper:
                    yy+=.022*math.exp(-((s-.81)/.10)**2)*math.sin(math.pi*t)
                vs.append((xx,yy,zz))
        for i in range(NU):
            for j in range(NB):
                a=i*(NB+1)+j
                fs.append((a,a+1,a+NB+2,a+NB+1))
        mesh_obj(("Upper" if upper else "Lower")+" sculpted lid "+str(side),vs,fs)
        tube("Soft palpebral edge",edge_points,edge_r,sides=10)

# A narrow recessed nostril opening, not a painted nose.
for side in [-1,1]:
    vs=[]
    cx=side*.210
    cz=5.696
    vs.append((cx,face_y(cx,cz)-.003,cz))
    count=48
    for i in range(count):
        a=2*math.pi*i/count
        x=cx+.069*math.cos(a)
        z=cz+.022*math.sin(a)+side*.075*(x-cx)
        vs.append((x,face_y(x,z)-.006,z))
    fs=[(0,1+i,1+(i+1)%count) for i in range(count)]
    mesh_obj("Recessed nostril "+str(side),vs,fs,recess)
    points=[]
    for i in range(33):
        a=.04+math.pi*.95*i/32
        x=cx+.077*math.cos(a)
        z=cz+.031*math.sin(a)+side*.075*(x-cx)
        points.append((x,face_y(x,z)-.010,z))
    tube("Alar edge "+str(side),points,[.005+.009*math.sin(math.pi*i/32) for i in range(33)],sides=8)

# Strong but quiet classical lips: one continuous bowed upper vermilion,
# a fuller lower plane, and a small shadow slit.
def mouth_z(a):
    return 5.324-.020*abs(a)**1.6+.007*math.exp(-((abs(a)-.30)/.19)**2)

for upper in [True,False]:
    vs,fs=[],[]
    NX,NV=110,22
    for i in range(NX+1):
        a=-1+2*i/NX
        x=.414*a
        weight=max(0,1-a*a)**.65
        line=mouth_z(a)
        if upper:
            height=(.077+.032*math.exp(-((abs(a)-.28)/.18)**2)
                    -.012*math.exp(-(a/.13)**2))*weight
        else:
            height=.111*max(0,1-a*a)**.78
        for j in range(NV+1):
            t=j/NV
            zi=line+(height*t+.004*weight if upper else -height*t-.005*weight)
            base=face_y(x,zi)
            relief=weight*(.046*(1-t)+(.056 if upper else .080)*math.sin(math.pi*t))
            y=base-relief-.002
            vs.append((x,y,zi))
    for i in range(NX):
        for j in range(NV):
            a=i*(NV+1)+j
            fs.append((a,a+1,a+NV+2,a+NV+1))
    mesh_obj("Cupid's bow • upper lip" if upper else "Full lower labial plane",vs,fs)
pts=[]; rr=[]
for i in range(101):
    a=-1+2*i/100
    x=.411*a
    z=mouth_z(a)
    w=max(0,1-a*a)**.65
    pts.append((x,face_y(x,z)-.027*w-.001,z))
    rr.append(.0065*math.sin(math.pi*i/100)**.5+.001)
tube("Quiet recessed mouth aperture",pts,rr,recess,8)

# Anatomical ears, attached by a broad root: cupped pinna, helix,
# bifurcated antihelix, tragus, antitragus and soft lobule.
def ellipsoid(name, loc, scale, mat=plaster, seg=48, rings=32):
    bpy.ops.mesh.primitive_uv_sphere_add(segments=seg, ring_count=rings, location=loc)
    ob=bpy.context.object
    ob.name=name
    ob.scale=scale
    ob.data.materials.append(mat)
    for poly in ob.data.polygons:
        poly.use_smooth=True
    return ob

for side in [-1,1]:
    origin=Vector((side*.965,.035,6.025))
    U=Vector((side*.82,.572,0))
    N=Vector((side*.572,-.82,0))
    Z=Vector((0,0,1))
    def ep(u,v,h):
        return origin+U*u+Z*v+N*h
    ellipsoid("Buried auricular root",tuple(ep(.045,-.01,-.008)),(.14,.165,.375))
    v,f=[],[]
    NR,NT=25,88
    for j in range(NR+1):
        r=j/NR
        for i in range(NT):
            a=2*math.pi*i/NT
            u=.180+.223*r*math.cos(a)*(1+.08*math.sin(a))
            zz=.485*r*math.sin(a)
            # broad relief cup, deepest at the concha
            h=-.012+.080*math.exp(-((r-.84)/.16)**2)-.044*math.exp(-(r/.48)**2)
            v.append(tuple(ep(u,zz,h)))
    for j in range(NR):
        for i in range(NT):
            a=j*NT+i;b=j*NT+(i+1)%NT
            f.append((a,b,b+NT,a+NT))
    # Thickness returns the outer border into the skull.
    base=len(v)
    for i in range(NT):
        a=2*math.pi*i/NT
        v.append(tuple(ep(.180+.223*math.cos(a)*(1+.08*math.sin(a)),.485*math.sin(a),-.075)))
    for i in range(NT):
        a=NR*NT+i;b=NR*NT+(i+1)%NT
        f.append((a,b,base+(i+1)%NT,base+i))
    f.append(tuple(base+i for i in reversed(range(NT))))
    mesh_obj("Sculpted auricle "+str(side),v,f)
    pts=[];rs=[]
    for i in range(95):
        a=-1.17+5.84*i/94
        u=.180+.211*math.cos(a)*(1+.08*math.sin(a))
        zz=.462*math.sin(a)
        pts.append(tuple(ep(u,zz,.082)))
        rs.append(.024+.020*math.sin(math.pi*i/94)**.45)
    tube("Rolled helix "+str(side),pts,rs,sides=12)
    # Antihelix stem and its superior fork.
    for branch in [0,1]:
        pts=[];rs=[]
        for i in range(43):
            t=i/42
            if branch==0:
                u=.215+.055*math.sin(math.pi*t)-.063*t
                zz=-.285+.63*t
            else:
                u=.235-.17*t
                zz=.09+.205*t
            pts.append(tuple(ep(u,zz,.036+.012*math.sin(math.pi*t))))
            rs.append(.009+.023*math.sin(math.pi*t)**.5)
        tube("Antihelix fork "+str(side),pts,rs,sides=10)
    ellipsoid("Tragus "+str(side),tuple(ep(.040,-.070,.055)),(.058,.061,.100),seg=32,rings=20)
    ellipsoid("Antitragus "+str(side),tuple(ep(.150,-.245,.061)),(.060,.049,.060),seg=32,rings=20)
    ellipsoid("Lobule "+str(side),tuple(ep(.172,-.411,.005)),(.111,.089,.119),seg=36,rings=24)

# Classical wavy hair.  A real scalp volume carries interlocking, tapered
# locks; each lock has broad carved flutes rather than fine strand hairs.
def theta_edge(phi):
    front=math.cos(phi)
    if front>0:
        # curved temple recession, lower sides, open central forehead
        return 1.10+.54*(1-front)**.70
    return 1.64+.30*(-front)

def scalp(theta,phi):
    return Vector((1.087*math.sin(theta)*math.sin(phi),
                   .105-1.125*math.sin(theta)*math.cos(phi),
                   6.48+1.295*math.cos(theta)))

def scalp_normal(theta,phi):
    v=scalp(theta,phi)-Vector((0,.105,6.48))
    return Vector((v.x/(1.087**2),v.y/(1.125**2),v.z/(1.295**2))).normalized()

v,f=[],[]
NH,NP=66,180
for j in range(NH+1):
    t=j/NH
    for i in range(NP):
        phi=2*math.pi*i/NP
        edge=theta_edge(phi)+.018*math.sin(9*phi)+.012*math.sin(17*phi)
        th=.002+(edge-.002)*t
        p0=scalp(th,phi)
        relief=.024*math.sin(9*phi+th*8)*math.sin(th)**2
        p0+=scalp_normal(th,phi)*relief
        v.append(tuple(p0))
for j in range(NH):
    for i in range(NP):
        a=j*NP+i;b=j*NP+(i+1)%NP
        f.append((a,b,b+NP,a+NP))
# Turn under at the hairline instead of leaving a paper-thin cap.
base=len(v)
for i in range(NP):
    phi=2*math.pi*i/NP
    p0=scalp(theta_edge(phi),phi)-scalp_normal(theta_edge(phi),phi)*.095
    v.append(tuple(p0))
for i in range(NP):
    a=NH*NP+i;b=NH*NP+(i+1)%NP
    f.append((a,b,base+(i+1)%NP,base+i))
mesh_obj("Hair • solid undercut mass",v,f)

def hair_lock(name, params, width, height, phase=0, lift=.02):
    centers=[];normals=[]
    for th,ph in params:
        n=scalp_normal(th,ph)
        centers.append(scalp(th,ph)+n*lift)
        normals.append(n)
    vs,fs=[],[]
    count=len(centers)
    around=24
    for i in range(count):
        t=i/(count-1)
        tangent=(centers[min(i+1,count-1)]-centers[max(0,i-1)]).normalized()
        n=normals[i]
        b=tangent.cross(n).normalized()
        n=b.cross(tangent).normalized()
        if n.dot(normals[i])<0:
            n=-n
        taper=(.17+.83*math.sin(math.pi*t)**.46)
        # Broad asymmetric crest with three quiet longitudinal tool flutes.
        for k in range(around):
            a=2*math.pi*k/around
            lateral=math.cos(a)
            top=math.sin(a)
            h=height*top
            if top>0:
                grooves=0
                for c in [-.52,.03,.54]:
                    grooves+=math.exp(-((lateral-c-.045*math.sin(t*5+phase))/.072)**2)
                h-=.015*grooves*top
            point=centers[i]+b*(width*taper*lateral)+n*(h*taper)
            vs.append(tuple(point))
    for i in range(count-1):
        for k in range(around):
            a=i*around+k;b=i*around+(k+1)%around
            fs.append((a,b,b+around,a+around))
    fs.append(tuple(reversed(range(around))))
    fs.append(tuple((count-1)*around+k for k in range(around)))
    return mesh_obj(name,vs,fs)

# Staggered crown-to-temple flowing locks.
lock_index=0
for row,(th0,num,length) in enumerate([(.15,9,.58),(.55,15,.57),(.94,19,.53),(1.31,20,.43),(1.60,15,.33)]):
    for k in range(num):
        phi=2*math.pi*(k+.46*(row%2))/num+random.uniform(-.055,.055)
        edge=theta_edge(phi)
        start=th0+random.uniform(-.065,.035)
        if start>edge-.16:
            continue
        end=min(start+length+random.uniform(-.07,.07),edge+.024)
        params=[]
        sweep=.17+random.uniform(-.04,.09)
        if math.sin(phi)<0:
            sweep=-sweep
        for i in range(43):
            t=i/42
            th=start+(end-start)*t
            ph=phi+sweep*math.sin(math.pi*t*.83)+.035*math.sin(2*math.pi*t)
            params.append((th,ph))
        hair_lock("Carved wavy lock %03d"%lock_index,params,
                  random.uniform(.135,.196),random.uniform(.085,.133),
                  phase=random.random()*6,lift=random.uniform(.03,.055))
        lock_index+=1

# More freely curled, individually modelled locks at the face and temples.
# These overturn the hairline, breaking the silhouette with real clumped mass.
for k in range(19):
    phi=-1.72+3.44*k/18+random.uniform(-.025,.025)
    edge=theta_edge(phi)
    params=[]
    curl_direction=1 if phi>=0 else -1
    for i in range(58):
        t=i/57
        # S-curve followed by a tucked-in end, on the scalp's own surface.
        th=edge-.29+.40*t-.115*t**5
        ph=phi+curl_direction*(.100*math.sin(t*math.pi*1.65)+.031*t)
        params.append((th,ph))
    hair_lock("Overturned forehead curl %02d"%k,params,
              random.uniform(.116,.16),random.uniform(.11,.155),
              phase=k*.8,lift=.075)

# Tighter hooked waves around the back of the head and above the ears.
for side in [-1,1]:
    for k in range(7):
        phi=side*(1.53+.215*k)
        e=theta_edge(phi)
        params=[]
        for j in range(48):
            t=j/47
            th=e-.25+.23*t+.072*math.sin(math.pi*t)
            ph=phi+side*.12*math.sin(t*math.pi*1.5)
            params.append((th,ph))
        hair_lock("Nape and temporal curl",params,.133,.113,k,.055)

# A barely raised mold seam on the rear of the cast, where a real two-part
# atelier mold would leave it. It is deliberately not drawn across the face.
points=[]
for i in range(90):
    z=3.19+2.62*i/89
    if z<=5.39:
        w=interp(z,neck_w)
        d=interp(z,neck_d)
        a=2.76
        x=w*math.sin(a)
        y=.20-d*math.cos(a)
    else:
        w=interp(z,widths)
        x=w*.37
        y=interp(z,backs)*math.sqrt(1-.37**2)
    points.append((x,y+.004,z))
tube("Subtle rear mold seam",points,.0045,seam_mat,7)

# Integral oval socle, made as a turned mesh with a few small edge losses.
profile=[(2.64,1.04,.735),(2.68,1.115,.78),(2.755,1.115,.78),
         (2.79,1.03,.725),(2.84,.96,.68),(2.94,.89,.635),
         (3.01,1.11,.635),(3.075,1.36,.65)]
v,f=[],[]
NS=160
for j,(z,rx,ry) in enumerate(profile):
    for i in range(NS):
        a=2*math.pi*i/NS
        nick=0
        if j in [1,2,3]:
            for center,spread,depth in [(3.94,.030,.031),(4.51,.022,.019),(5.39,.018,.025),(1.62,.025,.018)]:
                da=abs((a-center+math.pi)%(2*math.pi)-math.pi)
                nick+=depth*math.exp(-(da/spread)**4)
        v.append(((rx-nick)*math.cos(a),.16+(ry-nick)*math.sin(a),z-.35*nick))
for j in range(len(profile)-1):
    for i in range(NS):
        a=j*NS+i;b=j*NS+(i+1)%NS
        f.append((a,b,b+NS,a+NS))
f.append(tuple(reversed(range(NS))))
f.append(tuple((len(profile)-1)*NS+i for i in range(NS)))
socle=mesh_obj("Integral oval socle • small chipped arrises",v,f)

# Simple, unornamented column pedestal.
def bevel_cube(name, loc, scale, mat, bevel):
    bpy.ops.mesh.primitive_cube_add(size=1,location=loc)
    ob=bpy.context.object
    ob.name=name
    ob.dimensions=scale
    bpy.ops.object.transform_apply(location=False,rotation=False,scale=True)
    ob.data.materials.append(mat)
    mod=ob.modifiers.new("Soft stone arris","BEVEL")
    mod.width=bevel
    mod.segments=3
    ob.modifiers.new("Weighted corner normals","WEIGHTED_NORMAL")
    return ob

bevel_cube("Plain square studio pedestal", (0,.16,1.24),(2.05,1.83,2.48),ped_mat,.045)
bevel_cube("Pedestal top slab",(0,.16,2.555),(2.30,2.06,.17),ped_mat,.025)
bevel_cube("Pedestal foot",(0,.16,.10),(2.19,1.97,.20),ped_mat,.025)

# Neutral seamless studio, soft atelier key and restrained reflected fill.
bpy.ops.mesh.primitive_plane_add(size=200,location=(0,0,-.013))
floor=bpy.context.object
floor.name="Seamless neutral studio ground"
floor.data.materials.append(floor_mat)

world=bpy.data.worlds.new("Studio ambient")
world.use_nodes=True
world.node_tree.nodes.get("Background").inputs["Color"].default_value=(.34,.36,.40,1)
world.node_tree.nodes.get("Background").inputs["Strength"].default_value=.22
bpy.context.scene.world=world

def area(name,loc,power,color,size,target):
    data=bpy.data.lights.new(name,'AREA')
    data.energy=power
    data.color=color
    data.shape='DISK'
    data.size=size
    ob=bpy.data.objects.new(name,data)
    bpy.context.collection.objects.link(ob)
    ob.location=loc
    ob.rotation_euler=(Vector(target)-ob.location).to_track_quat('-Z','Y').to_euler()
    return ob

area("Large soft key • raking the facial planes",(-5.5,-7.0,10.5),1150,(1,.94,.85),3.5,(0,-.35,5.95))
area("Restrained studio bounce",(5.0,-4.0,6.5),160,(.84,.9,1),4.0,(0,0,5.6))
area("Soft high rear separation",(2.0,4.0,10.0),750,(1,.96,.89),3.2,(0,.2,6.15))

# Catalogue framing: eye-height camera, a restrained three-quarter turn,
# pedestal cropped low so the sculpt itself occupies the plate.
bpy.ops.object.camera_add(location=(7.2,-15.0,6.30))
camera=bpy.context.object
camera.name="Catalogue • eye-level three-quarter"
target=Vector((0,-.035,5.26))
camera.rotation_euler=(target-camera.location).to_track_quat('-Z','Y').to_euler()
camera.data.lens=93
camera.data.sensor_width=36
camera.data.clip_start=.1
camera.data.clip_end=250
camera.data.dof.use_dof=False
bpy.context.scene.camera=camera

scene=bpy.context.scene
scene.render.engine='CYCLES'
scene.cycles.samples=64
scene.cycles.use_denoising=True
scene.render.resolution_x=1536
scene.render.resolution_y=1536
scene.render.resolution_percentage=100
scene.render.image_settings.file_format='PNG'
scene.render.film_transparent=False
scene.view_settings.view_transform='AgX'
scene.view_settings.look='AgX - Medium High Contrast'
scene.view_settings.exposure=.05
scene.render.filepath='/mnt/data/render.png'
scene.world.color=(.15,.15,.15)

# Helpful organization in the saved, fully editable scene.
scene["Artwork"]="ATELIER / ideal youth in the antique manner"
scene["Medium"]="White casting plaster; geometric blank eyes; carved clumped waves"
scene["Presentation"]="Single cast, integral oval socle, plain column, raking soft studio light"
bpy.ops.wm.save_as_mainfile(filepath='/mnt/data/scene.blend')
```