naca.scm

(use-modules (ice-9 receive))
(use-modules (ao shapes))

(define (naca-thickness t)
  "naca-thickness t
  Returns the half-thickness for a symmetric airfoil
  of chord 1 and thickness t"
  (lambda (x)
    (* 5 t
    (+ (*  0.2969 (sqrt x))
       (* -0.1260 x)
       (* -0.3516 x x)
       (*  0.2843 x x x)
       (* -0.1015 x x x x)))))

(define (naca-symmetric t)
  "naca-symmetric t
  Returns a shape defining a symmetric airfoil
  of chord 1 and thickness t"
  (let* ((thickness (naca-thickness t))
         (side (lambda (x y z) (max (- y) (- y (thickness x))))))
    (union side (reflect-y side))))

(define (remap-xy shape transform)
  "remap-xy shape transform

  Remaps the shape with the given transform
  The transform should be a function of the form
    (f x y) -> x y"
  (lambda (x y z)
    (receive (cx cy) (transform x y)
      (shape cx (sqrt (+ (square (- x cx))
                         (square (- y cy)))) z))))

(define (camber-first m p)
  "camber m p
  For a given airfoil camber, return a solver that finds the
  perpendicular point on the first half of the camber

  m is the maximum camber
  p is the location of max camber (0 - 1)"
  ;; Autogenerated from naca.py
  (define (yc x)
    (* m (+ (* (expt p -1) x 2) (* (expt p -2) (expt x 2) -1))))
  (define (xc x y)
    (+ p (expt (+ (* (expt (+ (* (expt m 6) -16/3) (* (expt m 4) (expt p 2) 8) (* (expt m 2) (expt p 4) 5) (* (expt p 6) 2/3) (* (expt m 2) (expt p 3) x -18) (* (expt m 2) (expt p 2) (expt x 2) 9) (* (expt m 3) (expt y 3) 16/3) (* (+ (* (expt m 4) 2) (* (expt m 2) (expt p 2) -1)) (expt y 2) -8) (* (+ (* (expt m 5) 4) (* (expt m 3) (expt p 2) -4) (* m (expt p 4))) y 4)) 1/2) (expt m -3) (expt p 3) 1/12) (* (+ (expt p 5) (* (expt p 4) x -1)) (expt m -2) -1/4)) 1/3) (* (+ (* (expt m 2) (expt p 2) 2) (* (expt p 4) -1) (* m (expt p 2) y -2)) (expt m -2) (expt (+ (* (expt (+ (* (expt m 6) -16/3) (* (expt m 4) (expt p 2) 8) (* (expt m 2) (expt p 4) 5) (* (expt p 6) 2/3) (* (expt m 2) (expt p 3) x -18) (* (expt m 2) (expt p 2) (expt x 2) 9) (* (expt m 3) (expt y 3) 16/3) (* (+ (* (expt m 4) 2) (* (expt m 2) (expt p 2) -1)) (expt y 2) -8) (* (+ (* (expt m 5) 4) (* (expt m 3) (expt p 2) -4) (* m (expt p 4))) y 4)) 1/2) (expt m -3) (expt p 3) 1/12) (* (+ (expt p 5) (* (expt p 4) x -1)) (expt m -2) -1/4)) -1/3) 1/6)))
  (lambda (x y)
    (let ((x (xc x y)))
      (values x (yc x)))))

(define (camber-second m p)
  "camber-second m p
  For a given airfoil camber, return a solver that finds the
  perpendicular point on the second half of the camber

  m is the maximum camber
  p is the location of max camber (0 - 1)"
  ;; Autogenerated from naca.py
  (define (yc x)
    (* (+ (* p x 2) (* (expt x 2) -1) (* p -2) 1) m (expt (+ p -1) -2)))
  (define (xc x y)
    (+ (* (+ (* (expt 3 1/2) -I) 1) (expt (+ (* (expt (+ (* (expt m 6) -16/3) (* (expt p 6) 2/3) (* (expt m 3) (expt y 3) 16/3) (* (+ (expt m 2) 2) (expt p 4) 5) (* (expt p 5) -4) (* (expt m 4) 8) (* (+ (* (expt m 2) 3) 20) (expt p 3) -2/3) (* (+ (* (expt m 4) 8) (* (expt m 2) -15) 10) (expt p 2)) (* (+ (* (expt m 2) (expt p 2)) (* (expt m 2) p -2) (expt m 2)) (expt x 2) 9) (* (+ (* (expt m 4) 2) (* (expt m 2) (expt p 2) -1) (* (expt m 2) p 2) (* (expt m 2) -1)) (expt y 2) -8) (* (expt m 2) -4) (* (+ (* (expt m 4) 4) (* (expt m 2) -4) 1) p -4) (* (+ (* (expt m 2) (expt p 3)) (* (expt m 2) (expt p 2) -2) (* (expt m 2) p)) x -18) (* (+ (* (expt m 5) 4) (* m (expt p 4)) (* m (expt p 3) -4) (* (expt m 3) -4) (* (+ (* (expt m 3) 2) (* m -3)) (expt p 2) -2) (* (+ (* (expt m 3) 2) (* m -1)) p 4) m) y 4) 2/3) 1/2) (expt m -3) (expt (+ p -1) 3) 1/12) (* (+ (expt p 5) (* (expt p 4) (+ x 4) -1) (* (expt p 3) (+ (* x 2) 3) 2) (* (expt p 2) (+ (* x 3) 2) -2) (* p (+ (* x 4) 1)) (* x -1)) (expt m -2) -1/4)) 1/3) -1/2) p (* (+ (expt p 4) (* (+ (expt p 2) (* p -2) 1) (expt m 2) -2) (* (expt p 3) -4) (* (+ (* (expt p 2) y) (* p y -2) y) m 2) (* (expt p 2) 6) (* p -4) 1) (expt m -2) (+ (* (expt 3 1/2) I) 1) (expt (+ (* (expt (+ (* (expt m 6) -16/3) (* (expt p 6) 2/3) (* (expt m 3) (expt y 3) 16/3) (* (+ (expt m 2) 2) (expt p 4) 5) (* (expt p 5) -4) (* (expt m 4) 8) (* (+ (* (expt m 2) 3) 20) (expt p 3) -2/3) (* (+ (* (expt m 4) 8) (* (expt m 2) -15) 10) (expt p 2)) (* (+ (* (expt m 2) (expt p 2)) (* (expt m 2) p -2) (expt m 2)) (expt x 2) 9) (* (+ (* (expt m 4) 2) (* (expt m 2) (expt p 2) -1) (* (expt m 2) p 2) (* (expt m 2) -1)) (expt y 2) -8) (* (expt m 2) -4) (* (+ (* (expt m 4) 4) (* (expt m 2) -4) 1) p -4) (* (+ (* (expt m 2) (expt p 3)) (* (expt m 2) (expt p 2) -2) (* (expt m 2) p)) x -18) (* (+ (* (expt m 5) 4) (* m (expt p 4)) (* m (expt p 3) -4) (* (expt m 3) -4) (* (+ (* (expt m 3) 2) (* m -3)) (expt p 2) -2) (* (+ (* (expt m 3) 2) (* m -1)) p 4) m) y 4) 2/3) 1/2) (expt m -3) (expt (+ p -1) 3) 1/12) (* (+ (expt p 5) (* (expt p 4) (+ x 4) -1) (* (expt p 3) (+ (* x 2) 3) 2) (* (expt p 2) (+ (* x 3) 2) -2) (* p (+ (* x 4) 1)) (* x -1)) (expt m -2) -1/4)) -1/3) 1/12)))
  (lambda (x y)
    (let ((x (xc x y)))
      (values x (yc x)))))

;(ao-show "test" (remap-xy (naca-symmetric 0.2)
;                          (camber-first 0.05 0.3)))
;(ao-show "b" (remap-xy (naca-symmetric 0.2)
;                       (camber-second 0.05 0.3)))
(define m 0.05)
(define p 0.3)
  (define (yc x)
    (* (+ (* p x 2) (* (expt x 2) -1) (* p -2) 1) m (expt (+ p -1) -2)))
(ao-show "boop" (lambda (x y z) (- y (yc x))))