;*******************************************************************************************
;* Scherk-Collins Minimal Surface, by Paul Nylander, bugman123.com, 3/11/09
;*******************************************************************************************

;Basic Functions
(defun ² (x) (* x x))
(defun add (x1 x2) (if (listp x1) (mapcar 'add x1 x2) (+ x1 x2)))
(defun subtract (x1 x2) (if (listp x1) (mapcar 'subtract x1 x2) (- x1 x2)))
(defun × (x1 x2) (cond
 ((listp x1) (mapcar '(lambda (x) (× x x2)) x1))
 ((listp x2) (mapcar '(lambda (x) (× x1 x)) x2))
 (T (* x1 x2))
));d
(defun ÷ (x1 x2) (if (listp x1) (mapcar '(lambda (x) (÷ x x2)) x1) (/ x1 x2)))
(defun sinh (x) (* 0.5 (- (exp x) (exp (- x)))))
(defun cosh (x) (* 0.5 (+ (exp x) (exp (- x)))))

;Complex Functions
(defun complex (x y) (list x y))
(setq CR (complex 1.0 0.0) CI (complex 0.0 1.0))
(defun Re (z) (car z)); (if (numberp z) z (car z))
(defun Im (z) (cadr z)); (if (numberp z) 0.0 (cadr z))
(defun CAbs (z) (sqrt (+ (² (Im z)) (² (Re z)))))
(defun arg (z) (atan (Im z) (Re z)))
(defun CPow (z n / r theta) (setq r (CAbs z))
 (if (= r 0.0) (complex 0.0 0.0) (progn (setq theta (* n (arg z)))
  (CMult (expt r n) (complex (cos theta) (sin theta)))
 ));i
);d
(defun CSqrt (z) (CPow z 0.5))
(defun CMult (z1 z2) (cond
 ((numberp z1) (complex (* z1 (Re z2)) (* z1 (Im z2))))
 ((numberp z2) (complex (* z2 (Re z1)) (* z2 (Im z1))))
 (T (complex (- (* (Re z1) (Re z2)) (* (Im z1) (Im z2))) (+ (* (Im z1) (Re z2)) (* (Re z1) (Im z2)))))
));d
(defun CLog (z) (complex (log (CAbs z)) (arg z)))
(defun CCot (z) (÷ (complex (- (sin (* 2 (Re z)))) (sinh (* 2 (Im z)))) (- (cos (* 2 (Re z))) (cosh (* 2 (Im z))))))
(defun CAtan (z / z1 z2)
 (setq z1 (add CR (CMult CI Z)) z2 (subtract CR (CMult CI z)))
 (÷ (complex (- (arg z1) (arg z2)) (- (log (CAbs z2)) (log (CAbs z1)))) 2.0)
);d

;Converts Polygon into AutoCAD's Triangle Elements, only works for convex polygons
;(TriangulateElem '(1 2 3 4 5 6 7 8 9 10)) -> '((1 2 -3) (-1 3 -4) (-1 4 -5) (-1 5 -6) (-1 6 -7) (-1 7 -8) (-1 8 -9) (-1 9 10))
(defun TriangulateElem (elem / n i i0 i1 i2 elems) (setq n (length elem))
 (if (= n 3) (list elem) (progn
  (setq i0 (car elem) i1 (caddr elem) elems (list (list i0 (cadr elem) (- i1))) i 3)
  (while (< i (1- n)) (setq i2 (nth i elem) elems (append elems (list (list (- i0) i1 (- i2)))) i1 i2 i (1+ i)))
  (append elems (list (list (- i0) i1 (last elem))))
 ));i
);d

(defun mesh (nodes elems1 / elems2 nnode nelem elem)
 (setq elems2 (apply 'append (mapcar 'TriangulateElem elems1))); introduces AutoCAD's negative index notation
 (setq nnode (length nodes) nelem (length elems2))
 (entmake (list '(0 . "POLYLINE") '(70 . 64) (cons 71 nnode) (cons 72 nelem)))
 (foreach p nodes (entmake (list '(0 . "VERTEX") (cons 10 p) '(70 . 192))))
 (foreach elem elems2 (entmake (list '(0 . "VERTEX") '(10 0.0 0.0 0.0) '(70 . 128)
  (cons 71 (car elem)) (cons 72 (cadr elem)) (cons 73 (caddr elem))
  (cons 74 (- (abs (car elem))))
 )));f
 (entmake '((0 . "SEQEND")))
);d

(defun ParametricPlot3D (f u1 u2 du v1 v2 dv / u v i j i1 i2 i3 i4 imax jmax nodes elems)
 (setq imax (fix (+ (/ (- u2 u1) du) 0.5)) jmax (fix (+ (/ (- v2 v1) dv) 0.5)))
 (setq j 0); calculate nodes
 (while (<= j jmax) (setq v (+ v1 (* j dv)) i 0)
  (while (<= i imax) (setq u (+ u1 (* i du)) nodes (append nodes (list (f u v))) i (1+ i)))
  (setq j (1+ j))
 );w
 (setq i 1); assemble quad elements
 (while (<= i imax) (setq j 0)
  (while (< j jmax) (setq
   i1 (+ i (* j (1+ imax))) i2 (1+ i1) i3 (+ i1 (1+ imax) 1) i4 (+ i1 (1+ imax))
   elems (append elems (list (list i1 i2 i3 i4)))
   j (1+ j)
  ));w
  (setq i (1+ i))
 );w
 (mesh nodes elems)
);d

;Scherk-Collins Surface
(setq nstory 10 ntwist 1 wflange 1.5 rwarp 3.75)
(defun Twist (p theta / sint cost) (setq sint (sin theta) cost (cos theta))
 (list (- (* (car p) cost) (* (cadr p) sint)) (+ (* (car p) sint) (* (cadr p) cost)) (caddr p))
);d
(defun Warp (p theta / r) (setq r (+ (car p) rwarp))
 (list (* r (cos theta)) (* r (sin theta)) (cadr p))
);d
(defun Scherk (z / t1 t2 x y z)
 (setq t1 (CSqrt (× 2 (CCot z))) t2 (add (CCot z) CR))
 (setq x (/ (* xsign 0.5 (- (Re (CLog (subtract t1 t2))) (Re (CLog (add t1 t2))))) (sqrt 2.0)))
 (setq y (/ (* ysign (Re (CMult CI (add (subtract (CAtan (subtract CR t1)) (CAtan (add CR t1))) (× pi CR))))) (sqrt 2.0)))
 (setq z (Re z))
 (list x y z)
);d
(defun ScherkCollins (x y / theta) (setq theta (/ (* 4 x) nstory))
 (Warp (Twist (Scherk (complex x y)) (* ntwist theta)) theta)
);d
(setq x1 0.0 x2 (* 0.5 nstory pi) dx (/ pi 8))
(setq y1 1.0e-6 y2 (/ wflange 2) dy (/ (- y2 y1) 6))
(foreach xsign '(1 -1) (foreach ysign '(1 -1) (ParametricPlot3D ScherkCollins x1 x2 dx y1 y2 dy)))

(command "vpoint" (list -1.0 -1.0 1.0)); (command "zoom" "extents")