# ============================================================================
# CC0 "No zmights zmeserved"
# To the extent possible under law, the author(s) have dedicated all copyright
# and related and neighboring rights to this software to the public domain
# worldwide. See http://creativecommons.org/publicdomain/zero/1.0/.
# ============================================================================

# Warps a reference plane to a surface and keeps off-planar nodes perpendicular to this surface.
# Select nodes before using.


import math

# Lr, Ls:   dimensions of reference plane
# zm:       height of hypar

Lr = 10.0
Ls = 10.0
zm = 2.0


# fx calculates the x-coordinate of the surface
def fx(r, s, Lr, Ls, zm):
    return r


# fx calculates the y-coordinate of the surface
def fy(r, s, Lr, Ls, zm):
    return s


# fx calculates the z-coordinate of the surface
def fz(r, s, Lr, Ls, zm):
    return 4.0 * zm * r * s / Lr / Ls


# dfxdr calculates the r-derivative of the x-coordinate of the surface at point (r, s)
def dfxdr(r, s, Lr, Ls, zm):
    return 1.0


# dfxds calculates the s-derivative of the x-coordinate of the surface at point (r, s)
def dfxds(r, s, Lr, Ls, zm):
    return 0.0


# dfydr calculates the r-derivative of the y-coordinate of the surface at point (r, s)
def dfydr(r, s, Lr, Ls, zm):
    return 0.0


# dfyds calculates the s-derivative of the y-coordinate of the surface at point (r, s)
def dfyds(r, s, Lr, Ls, zm):
    return 1.0


# dfzdr calculates the r-derivative of the z-coordinate of the surface at point (r, s)
def dfzdr(r, s, Lr, Ls, zm):
    return 4.0 * zm * s / Lr / Ls


# dfzds calculates the s-derivative of the z-coordinate of the surface at point (r, s)
def dfzds(r, s, Lr, Ls, zm):
    return 4.0 * zm * r / Lr / Ls


for node_id in mw.selected_nodes():
    coord = mw.node(node_id)
    r = coord.x
    s = coord.y
    t = coord.z

    # v_mid = position vector of surface
    v_mid = [fx(r, s, Lr, Ls, zm), fy(r, s, Lr, Ls, zm), fz(r, s, Lr, Ls, zm)]

    # dvdr = tangent vector to surface = d(v_mid)/dr
    # dvds = tangent vector to surface = d(v_mid)/ds
    dvdr = [dfxdr(r, s, Lr, Ls, zm), dfydr(r, s, Lr, Ls, zm), dfzdr(r, s, Lr, Ls, zm)]
    dvds = [dfxds(r, s, Lr, Ls, zm), dfyds(r, s, Lr, Ls, zm), dfzds(r, s, Lr, Ls, zm)]

    # n = vector normal to surface at (r, s) = dvdr x dvds
    n = [dvdr[1] * dvds[2] - dvdr[2] * dvds[1],
         dvdr[2] * dvds[0] - dvdr[0] * dvds[2],
         dvdr[0] * dvds[1] - dvdr[1] * dvds[0]]

    # norm = length of normal vector
    norm = math.sqrt(n[0] ** 2 + n[1] ** 2 + n[2] ** 2)

    # nn = unit normal vector at (r, s)
    nn = [entry / norm for entry in n]

    mw.set_node_x(node_id, v_mid[0] + t * nn[0])
    mw.set_node_y(node_id, v_mid[1] + t * nn[1])
    mw.set_node_z(node_id, v_mid[2] + t * nn[2])