# 시드니 오페라하우스 (Sydney Opera House) — Blender Python 코드만으로 만든 모델
#
# - 지붕 셸: 각 셸을 '좌우 대칭인 두 장의 곡면'으로 만들었다. 옆에서 본 영역(능선 원호·입구 곡선·바닥선으로
#   둘러싸인 삼각형)을 Coons 패치로 격자화하고, 높이마다 뾰족 아치(고딕 아치) 단면을 대입해 x 를 구한다.
#   격자의 한 방향 선이 받침점(발)에서 능선으로 부채꼴로 퍼지므로, 그 선을 따라 실제 셸의 리브 이음줄을 얹었다.
# - 콘서트홀(서쪽)·조앤 서덜랜드 극장(동쪽)은 북쪽(항구)을 향한 셸 3장 + 남쪽(계단)을 향한 셸 1장,
#   베넬롱 레스토랑은 작은 셸 2장. 셸 입구는 청동빛 유리벽(선직면)과 멀리언 격자로 막았다.
# - 기단(포디엄)은 화강암 다각형 돌출, 남쪽 '모뉴멘털 계단'은 단면 프로필을 폭 방향으로 돌출해 만들었다.
# - 단위 m, 실제 크기 기준(전체 길이 약 185m, 최고점은 해수면 기준 약 67m).
import math
from atelier import *

SHOT = {"azimuth": 155, "elevation": 11, "lens": 50, "light": "day", "background": "#22303f"}

G_H = 0.6            # 앞마당(포어코트) 판 두께
P_H = 12.0           # 기단 높이
Z0 = G_H + P_H       # 기단 윗면(브로드워크) 높이

M_TILE = mat('Shell_Tile_Cream', '#ebe6d8', 0.22)
M_RIB = mat('Shell_Rib_Joint', '#cdc3ad', 0.6)
M_UNDER = mat('Shell_Underside_Concrete', '#b8ad98', 0.85)
M_GLASS = mat('Glass_Bronze_Tint', '#3b3229', 0.07, 0.55)
M_MULL = mat('Mullion_Bronze', '#2b251f', 0.45, 0.7)
M_GRANITE = mat('Podium_Granite', '#8c705d', 0.75)
M_GRANITE_D = mat('Podium_Granite_Joint', '#6e5646', 0.85)
M_PAVE = mat('Broadwalk_Paving', '#b49d86', 0.85)
M_STEPS = mat('Monumental_Steps', '#a58e79', 0.8)
M_FORE = mat('Forecourt_Sandstone', '#c4b192', 0.9)


# ---------------------------------------------------------------- 공용 도우미
def mesh_from_bm(name, bm, materials, loc=(0, 0, 0), rot=(0, 0, 0)):
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me)
    bm.free()
    ob = bpy.data.objects.new(name, me)
    link(ob)
    for m in materials:
        me.materials.append(m)
    ob.location = Vector(loc)
    ob.rotation_euler = Euler(rot)
    return ob


def beam(bm, p0, p1, w, d, up=Vector((0, 0, 1))):
    """두 점 사이 직육면체 막대(멀리언, 난간)."""
    p0, p1 = Vector(p0), Vector(p1)
    t = p1 - p0
    if t.length < 1e-4:
        return
    t.normalize()
    s = t.cross(up)
    if s.length < 1e-3:
        s = t.cross(Vector((1, 0, 0)))
    s.normalize()
    n = s.cross(t).normalized()
    vs = []
    for p in (p0, p1):
        for a, b in ((-1, -1), (1, -1), (1, 1), (-1, 1)):
            vs.append(bm.verts.new(p + s * (a * w / 2) + n * (b * d / 2)))
    f = [(0, 1, 2, 3), (7, 6, 5, 4), (0, 4, 5, 1), (1, 5, 6, 2), (2, 6, 7, 3), (3, 7, 4, 0)]
    for q in f:
        bm.faces.new([vs[i] for i in q])


def profile_yz(name, pts, x0, x1, material):
    """YZ 단면 다각형을 x0~x1 로 돌출."""
    bm = bmesh.new()
    a = [bm.verts.new((x0, y, z)) for y, z in pts]
    b = [bm.verts.new((x1, y, z)) for y, z in pts]
    bm.faces.new(a)
    bm.faces.new(list(reversed(b)))
    n = len(pts)
    for i in range(n):
        bm.faces.new((a[i], a[(i + 1) % n], b[(i + 1) % n], b[i]))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces)
    return mesh_from_bm(name, bm, [material])


def smin(a, b, k=2.5):
    return -k * math.log(math.exp(-a / k) + math.exp(-b / k))


# ---------------------------------------------------------------- 셸
def ridge_arc(yB, hT, sag, n=240):
    """B=(yB,0) 에서 끝점 T=(0,hT) 까지 위로 볼록한 원호(옆에서 본 능선)."""
    by, bz, ty, tz = yB, 0.0, 0.0, hT
    dy, dz = ty - by, tz - bz
    L = math.hypot(dy, dz)
    ny, nz = -dz / L, dy / L
    r = (L * L / 4 + sag * sag) / (2 * sag)
    my, mz = (by + ty) / 2, (bz + tz) / 2
    cy, cz = my - ny * (r - sag), mz - nz * (r - sag)
    aB = math.atan2(bz - cz, by - cy)
    aT = math.atan2(tz - cz, ty - cy)
    da = (aT - aB + math.pi) % (2 * math.pi) - math.pi
    return [(cy + r * math.cos(aB + da * i / n), cz + r * math.sin(aB + da * i / n)) for i in range(n + 1)]


def build_shell(name, yB, hT, yF, W, sag=8.0, bow=3.0, k=0.8, nu=44, nv=30, thick=1.1, rib_step=3,
                loc=(0, 0, 0), facing=1, glass_wall=True):
    """끝점(tip)이 로컬 원점 위 (0,0,hT). 앞(입구)은 +y. facing=-1 이면 남쪽을 본다."""
    rp = ridge_arc(yB, hT, sag)
    ry = [p[0] for p in rp]

    def h_of(y):
        if y <= ry[0]:
            return 0.0
        if y >= ry[-1]:
            return rp[-1][1]
        lo, hi = 0, len(ry) - 1
        while hi - lo > 1:
            mid = (lo + hi) // 2
            if ry[mid] <= y:
                lo = mid
            else:
                hi = mid
        t = (y - ry[lo]) / max(ry[hi] - ry[lo], 1e-9)
        return rp[lo][1] * (1 - t) + rp[hi][1] * t

    def x_of(y, z):
        h = h_of(y)
        if h < 1e-4 or z >= h:
            return 0.0
        w = min(smin(W, k * h), 0.97 * h)
        c = (w * w - h * h) / (2 * w)
        rho = w - c
        return max(c + math.sqrt(max(rho * rho - max(z, 0) ** 2, 0.0)), 0.0)

    def top(u):
        f = u * (len(rp) - 1)
        i = min(int(f), len(rp) - 2)
        t = f - i
        return (rp[i][0] * (1 - t) + rp[i + 1][0] * t, rp[i][1] * (1 - t) + rp[i + 1][1] * t)

    F, B, T = (yF, 0.0), (yB, 0.0), (0.0, hT)
    ctrl = ((yF + 0.0) / 2 - bow, hT / 2)

    def mouth(v):
        a, b, c = (1 - v) ** 2, 2 * v * (1 - v), v * v
        return (a * F[0] + b * ctrl[0] + c * T[0], a * F[1] + b * ctrl[1] + c * T[1])

    def coons(u, v):
        Ly, Lz = yF + (yB - yF) * v, 0.0
        Ry, Rz = mouth(v)
        Ty, Tz = top(u)
        y = (1 - u) * Ly + u * Ry + (1 - v) * F[0] + v * Ty - ((1 - u) * (1 - v) * F[0] + u * (1 - v) * F[0] + (1 - u) * v * B[0] + u * v * T[0])
        z = (1 - u) * Lz + u * Rz + (1 - v) * F[1] + v * Tz - ((1 - u) * (1 - v) * F[1] + u * (1 - v) * F[1] + (1 - u) * v * B[1] + u * v * T[1])
        y = min(max(y, yB), 0.0)
        z = min(max(z, 0.0), h_of(y))
        return y, z

    def P(side, u, v):
        y, z = coons(u, v)
        return Vector((side * x_of(y, z), y, z))

    # --- 셸 곡면 (좌우 두 장)
    bm = bmesh.new()
    for side in (1, -1):
        g = [[bm.verts.new(P(side, i / nu, j / nv)) for j in range(nv + 1)] for i in range(nu + 1)]
        for i in range(nu):
            for j in range(nv):
                q = [g[i][j], g[i + 1][j], g[i + 1][j + 1], g[i][j + 1]]
                if side < 0:
                    q.reverse()
                try:
                    bm.faces.new(q)
                except ValueError:
                    pass
    bmesh.ops.remove_doubles(bm, verts=bm.verts, dist=1e-3)
    bmesh.ops.dissolve_degenerate(bm, edges=bm.edges, dist=1e-4)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces)
    bm.faces.ensure_lookup_table()
    # 오른쪽 장의 법선이 바깥(+x)을 보도록
    sx = sum(f.normal.x for f in bm.faces if f.calc_center_median().x > 1.0)
    if sx < 0:
        for f in bm.faces:
            f.normal_flip()
    for e in bm.edges:
        if all(abs(v.co.x) < 1e-3 for v in e.verts):
            e.smooth = False
    for f in bm.faces:
        f.smooth = True
        f.material_index = 0
    rot = (0, 0, 0 if facing > 0 else math.pi)
    sh = mesh_from_bm(name, bm, [M_TILE, M_UNDER], loc, rot)
    sm = sh.modifiers.new('Thickness', 'SOLIDIFY')
    sm.thickness = thick
    sm.offset = -1
    sm.use_even_offset = True
    sm.material_offset = 1
    sm.material_offset_rim = 1
    apply_all(sh)

    # --- 리브 이음줄: 받침점에서 능선으로 부채꼴
    bm = bmesh.new()
    eps = 1e-3
    for side in (1, -1):
        for i in range(rib_step, nu, rib_step):
            u = i / nu
            pts = []
            for j in range(2, nv + 1):
                v = j / nv
                p = P(side, u, v)
                pu = P(side, min(u + eps, 1), v) - P(side, max(u - eps, 0), v)
                pv = P(side, u, min(v + eps, 1)) - P(side, u, max(v - eps, 0))
                n = pu.cross(pv)
                if n.length < 1e-9:
                    continue
                n.normalize()
                if n.x * side < 0 and abs(p.x) > 0.5:
                    n = -n
                if abs(p.x) < 0.3:
                    n = Vector((0, 0, 1))
                pts.append((p, n))
            for a in range(len(pts) - 1):
                (p0, n0), (p1, n1) = pts[a], pts[a + 1]
                t0 = (p1 - p0).normalized()
                s0 = n0.cross(t0).normalized()
                s1 = n1.cross(t0).normalized()
                v = [bm.verts.new(p0 + n0 * 0.16), bm.verts.new(p0 + s0 * 0.32 - n0 * 0.12), bm.verts.new(p0 - s0 * 0.32 - n0 * 0.12),
                     bm.verts.new(p1 + n1 * 0.16), bm.verts.new(p1 + s1 * 0.32 - n1 * 0.12), bm.verts.new(p1 - s1 * 0.32 - n1 * 0.12)]
                bm.faces.new((v[0], v[1], v[4], v[3]))
                bm.faces.new((v[0], v[3], v[5], v[2]))
    rb = mesh_from_bm(name + '_RibJoints', bm, [M_RIB], loc, rot)
    out = [sh, rb]

    # --- 입구 유리벽 (입구 곡선에서 바닥선까지 선직면) + 멀리언 격자
    if glass_wall:
        nr = 18
        ruls = []
        for side, rng in ((1, range(0, nr + 1)), (-1, range(nr - 1, -1, -1))):
            for j in rng:
                y, z = mouth(j / nr)
                x = x_of(y, z) * side * 0.96
                topp = Vector((x, y - 0.9, max(z - 0.9, 0.0)))
                bot = Vector((x, yF - 0.9, 0.0))
                ruls.append((bot, topp))
        bm = bmesh.new()
        vv = [(bm.verts.new(b), bm.verts.new(t)) for b, t in ruls]
        for a in range(len(vv) - 1):
            try:
                bm.faces.new((vv[a][0], vv[a + 1][0], vv[a + 1][1], vv[a][1]))
            except ValueError:
                pass
        bmesh.ops.remove_doubles(bm, verts=bm.verts, dist=1e-3)
        gl = mesh_from_bm(name + '_GlassWall', bm, [M_GLASS], loc, rot)
        bm = bmesh.new()
        fwd = Vector((0, 0.12, 0))
        for b, t in ruls[1:-1]:
            beam(bm, b + fwd, t + fwd, 0.22, 0.3)
        step = 4.5
        zl = step
        while zl < hT - 1:
            row = [b.lerp(t, zl / t.z) + fwd for b, t in ruls if t.z > zl]
            for a in range(len(row) - 1):
                beam(bm, row[a], row[a + 1], 0.18, 0.26)
            zl += step
        mu = mesh_from_bm(name + '_Mullions', bm, [M_MULL], loc, rot)
        out += [gl, mu]
    return out


def shell_group(prefix, gx, specs, scale=1.0, zbase=Z0 - 0.3):
    objs = []
    for nm, y_tip, hT, Lb, foot, W, facing, sag in specs:
        objs += build_shell(f'{prefix}_{nm}', -Lb * scale, hT * scale, -foot * scale, W * scale, sag=sag * scale,
                            bow=3.0 * scale, loc=(gx, y_tip, zbase), facing=facing)
    return objs


# 이름, 끝점 y(세계), 끝점 높이, 뒤로 뻗는 길이, 발까지 거리, 바닥 반폭, 방향(+1 북/-1 남), 능선 볼록량
CONCERT = [
    ('Shell_A1_NorthFoyer', 66.0, 44.0, 48.0, 13.0, 22.5, 1, 7.0),
    ('Shell_A2_Main', 44.0, 54.0, 56.0, 14.0, 25.5, 1, 8.0),
    ('Shell_A3', 13.0, 41.0, 44.0, 11.0, 21.5, 1, 6.5),
    ('Shell_A4_South', -44.0, 29.0, 40.0, 7.0, 19.5, -1, 5.0),
]
THEATRE = [
    ('Shell_B1_NorthFoyer', 62.0, 44.0, 48.0, 13.0, 22.5, 1, 7.0),
    ('Shell_B2_Main', 40.0, 54.0, 56.0, 14.0, 25.5, 1, 8.0),
    ('Shell_B3', 10.0, 41.0, 44.0, 11.0, 21.5, 1, 6.5),
    ('Shell_B4_South', -42.0, 29.0, 40.0, 7.0, 19.5, -1, 5.0),
]
RESTAURANT = [
    ('Shell_R1_South', -68.0, 19.0, 22.0, 5.0, 10.0, -1, 3.0),
    ('Shell_R2_North', -50.0, 14.0, 18.0, 4.0, 8.0, 1, 2.5),
]

shell_group('ConcertHall', -29.0, CONCERT)
shell_group('Theatre', 27.0, THEATRE, scale=0.9)
shell_group('Bennelong', -45.0, RESTAURANT)

# ---------------------------------------------------------------- 기단 · 계단 · 앞마당
PODIUM = [(-60, -80), (60, -80), (60, 46), (36, 86), (-36, 86), (-60, 52)]
FORE = [(-68, -124), (68, -124), (68, 48), (40, 94), (-40, 94), (-68, 56)]


def offset_poly(pts, d):
    n = len(pts)
    out = []
    for i in range(n):
        p0, p1, p2 = Vector(pts[i - 1]), Vector(pts[i]), Vector(pts[(i + 1) % n])
        e0 = (p1 - p0).normalized()
        e1 = (p2 - p1).normalized()
        n0 = Vector((e0.y, -e0.x))
        n1 = Vector((e1.y, -e1.x))
        bis = (n0 + n1).normalized()
        out.append(tuple(p1 + bis * (d / max(bis.dot(n0), 0.3))))
    return out


prism('Forecourt', FORE, G_H, material=M_FORE)
prism('Podium_Granite', PODIUM, P_H - 0.6, loc=(0, 0, G_H - 0.05), material=M_GRANITE)
prism('Podium_Broadwalk_Edge', offset_poly(PODIUM, 1.2), 1.4, loc=(0, 0, Z0 - 1.4), material=M_PAVE)

# 남쪽 모뉴멘털 계단 (3단 계단참)
rise, tread, land = 0.25, 0.55, 4.0
pts = [(-80.0, G_H - 0.05), (-80.0, Z0 - 0.02)]
cy, cz = -80.0, Z0 - 0.02
for f in range(3):
    for s in range(16):
        cz -= rise
        pts.append((cy, max(cz, G_H)))
        if f == 2 and s == 15:
            break
        cy -= tread + (land if s == 15 else 0.0)
        pts.append((cy, max(cz, G_H)))
pts.append((cy, G_H - 0.05))
pts = [pts[0]] + pts[1:]
STEP_X0, STEP_X1 = -30.0, 56.0
profile_yz('Monumental_Steps', list(reversed(pts)), STEP_X0, STEP_X1, M_STEPS)
for nm, x0, x1 in (('Steps_CheekWall_W', STEP_X0 - 2.2, STEP_X0 + 0.05), ('Steps_CheekWall_E', STEP_X1 - 0.05, STEP_X1 + 2.2)):
    profile_yz(nm, [(-80.5, G_H - 0.05), (-80.5, Z0 + 0.5), (cy, G_H + 0.9), (cy, G_H - 0.05)], x0, x1, M_GRANITE)

# ---------------------------------------------------------------- 기단 디테일
# 화강암 벽면의 수평 줄눈(살짝 튀어나온 띠), 수직 판 이음, 아래층 창 띠, 브로드워크 난간
def podium_details():
    bm_band = bmesh.new()
    bm_win = bmesh.new()
    bm_rail = bmesh.new()
    n = len(PODIUM)
    for i in range(n):
        a, b = Vector(PODIUM[i]), Vector(PODIUM[(i + 1) % n])
        e = b - a
        L = e.length
        t = e.normalized()
        out = Vector((t.y, -t.x))          # 바깥쪽 법선 (반시계 다각형)
        a3 = Vector((a.x, a.y, 0)); t3 = Vector((t.x, t.y, 0)); o3 = Vector((out.x, out.y, 0))
        # 수평 띠 (3줄)
        for zz in (G_H + 3.8, G_H + 6.6, G_H + 8.8):
            beam(bm_band, a3 + o3 * 0.12 + Vector((0, 0, zz)), a3 + t3 * L + o3 * 0.12 + Vector((0, 0, zz)), 0.35, 0.25)
        # 수직 판 이음 / 창 띠
        cnt = int(L // 4.0)
        for k in range(1, cnt):
            p = a3 + t3 * (L * k / cnt) + o3 * 0.1
            beam(bm_band, p + Vector((0, 0, G_H + 3.8)), p + Vector((0, 0, Z0 - 1.4)), 0.22, 0.2)
            if i != 0 and k % 2 == 1 and 1 < k < cnt - 1:
                c = a3 + t3 * (L * (k + 0.5) / cnt) + o3 * 0.06
                q = [c + t3 * sx * 3.0 + Vector((0, 0, G_H + sz)) for sx, sz in ((-1, 0.6), (1, 0.6), (1, 3.2), (-1, 3.2))]
                vs = [bm_win.verts.new(v) for v in q]
                f = bm_win.faces.new(vs)
                if f.normal.dot(o3) < 0:
                    f.normal_flip()
        # 브로드워크 난간 (남쪽 계단 구간 제외)
        if i == 0:
            continue
        ro = offset_poly(PODIUM, 0.6)
        ra, rb = Vector((*ro[i], Z0)), Vector((*ro[(i + 1) % n], Z0))
        beam(bm_rail, ra + Vector((0, 0, 1.05)), rb + Vector((0, 0, 1.05)), 0.14, 0.1)
        beam(bm_rail, ra + Vector((0, 0, 0.5)), rb + Vector((0, 0, 0.5)), 0.06, 0.06)
        cntp = max(int((rb - ra).length // 2.4), 1)
        for k in range(cntp + 1):
            p = ra.lerp(rb, k / cntp)
            beam(bm_rail, p - Vector((0, 0, 0.05)), p + Vector((0, 0, 1.08)), 0.08, 0.08)
    mesh_from_bm('Podium_Granite_Joints', bm_band, [M_GRANITE_D])
    mesh_from_bm('Podium_Windows', bm_win, [M_GLASS])
    mesh_from_bm('Broadwalk_Railing', bm_rail, [M_MULL])


podium_details()


# 계단 손잡이 (경사를 따라 놓인 청동 난간 4줄)
def stair_rails():
    bm = bmesh.new()
    for x in (STEP_X0 + 18.0, STEP_X0 + 38.0, STEP_X0 + 58.0, STEP_X0 + 76.0):
        prev = None
        for y, z in pts[1:-1:4]:
            p = Vector((x, y, z + 1.0))
            beam(bm, Vector((x, y, z - 0.02)), p + Vector((0, 0, 0.04)), 0.07, 0.07)
            if prev is not None:
                beam(bm, prev, p, 0.1, 0.1)
            prev = p
    mesh_from_bm('Steps_Handrails', bm, [M_MULL])


stair_rails()