naca.py

# Sage code to generate NACA airfoil equations
import operator

var('x, y, xc, m, p')

def run(y_chord):
    g_chord = vector([1, derivative(y_chord(xc), xc)])
    pos_pt = vector([x, y])
    pos_chord = vector([xc, y_chord(xc)])
    delta = pos_pt - pos_chord
    print delta.dot_product(g_chord) == 0
    sols = solve(delta.dot_product(g_chord) == 0, xc)
    return expand(sols[1].op[1])

def to_prefix(expr):
    if expr.operator() == None:
        if str(expr) == 'I':
            return '0'
        return str(expr)
    children = " ".join(map(to_prefix, expr.op))
    f = {operator.add: "+",
         operator.sub: "-",
         operator.pow: "expt",
         operator.mul: "*"}[expr.operator()]
    return "(%s %s)" % (f, children)

def to_func(y_chord):
    return """  (define (yc x)
    %s)
  (define (xc x y)
    %s)
  (lambda (x y)
    (let ((x (xc x y)))
      (values x (yc x))))""" % (to_prefix(y_chord(x)),
                                to_prefix(run(y_chord)))


y_chord_a = lambda x: m * (2*x/p -(x/p)^2)
y_chord_b = lambda x: m / (1-p)^2 * (1-2*p + 2*p*x - x^2)