; Notes: http://integrals.wolfram.com/index.jsp ; Integral 1/√(a+bx+x^2) = log(2√(a+bx+x^2) + b + 2x) (define a 3.6) (define b -0.93) (define (h x) (sqrt (+ a (* b x) (* x x)))) (define (f x) (/ (h x))) (define (I x) (log (+ (* 2 (h x)) b x x))) (define (D f) (let ((ep 0.0000001)) (lambda (x) (/ (- (f (+ x ep)) (f (- x ep))) (* 2 ep))))) (define g (D I)) ; functions f and g should be near. (define (sq x) (* x x)) (define (p1 x y z) (let* ((r2 (+ (sq x) (sq y) (sq z)))) (log (/ (+ (sqrt (+ r2 (* -2 x) 1)) (- x) 1) (- (sqrt r2) x))))) (define (p0 x y z) (/ (sqrt (+ (sq x) (sq y) (sq z))))) (define ((poi p) x y z) (let ((e 0.00001)) (/ (- (* 6 (p x y z)) (+ (p (+ x e) y z) (p (- x e) y z) (p x (+ y e) z) (p x (- y e) z) (p x y (+ z e)) (p x y (- z e)))) (* e e)))) ; Two sticks from [0, 1/2] & [1/2, 1] (define (pv x y z) (+ (p1 (* 2 x) (* 2 y) (* 2 z)) (p1 (- (* 2 x) 1) (* 2 y) (* 2 z)))) ; pv=p1 ; Two sticks together from (-1, 0, 0) to (1, 0, 0) (define (pu x y z) (+ (p1 x y z) (p1 (- x) y z))) (define (px x y z) (p1 (/ (+ x 1) 2) (/ y 2) (/ z 2))) (define (p3 L x y z) (let ((r (/ L))) (px (* x r) (* y r) (* z r)))) ; px=pu (define (pd x y z) (log (/ (+ (sqrt (+ (/ (sq x) 4) (* -0.5 x) 1/4 (/ (sq y) 4) (/ (sq z) 4))) (* -0.5 x) 0.5) (+ (sqrt (+ (/ (sq x) 4) (* 0.5 x) 1/4 (/ (sq y) 4) (/ (sq z) 4))) (* -0.5 x) -0.5)))) (define (pd x y z) (log (/ (- (sqrt (+ (sq (- x 1)) (sq y) (sq z))) (- x 1)) (- (sqrt (+ (sq (+ x 1)) (sq y) (sq z))) (+ x 1))))) (define (pd x y z) (let ((f (lambda (x) (- (sqrt (+ (sq x) (sq y) (sq z))) x)))) (log (/ (f (- x 1)) (f (+ x 1)))))) (define (pd x y z) (let* ((yz (+ (sq y) (sq z))) (f (lambda (x) (- (sqrt (+ (sq x) yz)) x)))) (log (/ (f (- x 1)) (f (+ x 1)))))) ; px=pu=pd ; The last definition of pd is a bit quicker than the previous, but not as pretty. (define ((pg u d) x y z) (let* ((yz (+ (sq y) (sq z))) (f (lambda (x) (- (sqrt (+ (sq x) yz)) x)))) (log (/ (f (- x d)) (f (- x u)))))) (define ni (let* ((n 100) (del (+ 0.0 (/ n)))) (lambda (f l p) (let ((dx (* del (- l f)))) (* dx (let sum ((s 0)(x (+ f (/ dx 2)))(k n)) (if (zero? k) s (sum (+ s (p x))(+ x dx)(- k 1))))))))) (ni 2 5 (lambda (x) 2)) ; => 6.0 (ni 0 1 sq) ; => 1/3 (0.33332500000000054) ; following is the field of a brick with corners (±L, ±B, ±D). ; https://arxiv.org/pdf/1206.3857.pdf (define ((brick L B D) X Y Z) (let ((I (lambda (x y z) (let ((r (sqrt (+ (* x x) (* y y) (* z z)))))(+ (* y z (log (+ x r))) (* -1/2 (* x x) (atan (/ (* y z)(* x r)))) (* x z (log (+ y r))) (* -1/2 (* y y) (atan (/ (* x z)(* y r)))) (* x y (log (+ z r))) (* -1/2 (* z z) (atan (/ (* x y)(* z r))))))))) (let ((P (- D Z)) (p (- (+ D Z))) (Q (- B Y)) (q (- (+ B Y))) (R (- L X)) (r (- (+ L X)))) (- (- (- (I P Q R) (I P Q r)) (- (I P q R) (I P q r))) (- (- (I p Q R) (I p Q r)) (- (I p q R) (I p q r))))))) ; br1 = brick (define ((br1 L B D) X Y Z) (let ((I (lambda (x y z) (let ((xx (* x x)) (xy (* x y)) (xz (* x z)) (yy (* y y)) (yz (* y z)) (zz (* z z))) (let ((r (sqrt (+ xx yy zz)))) (- (+ (* yz (log (+ x r))) (* xz (log (+ y r))) (* xy (log (+ z r)))) (* 1/2 (+ (* xx (atan (/ yz (* x r)))) (* yy (atan (/ xz (* y r)))) (* zz (atan (/ xy (* z r)))))))))))) (let ((P (- D Z)) (p (- (+ D Z))) (Q (- B Y)) (q (- (+ B Y))) (R (- L X)) (r (- (+ L X)))) (- (- (- (I P Q R) (I P Q r)) (- (I P q R) (I P q r))) (- (- (I p Q R) (I p Q r)) (- (I p q R) (I p q r))))))) ; This version works on the brick faces. (define ((br2 L B D) X Y Z) (let ((I (lambda (x y z) (let* ((r (sqrt (+ (* x x) (* y y) (* z z)))) (J (lambda (x y z) (+ (* y z (log (+ x r))) (if (zero? x) 0 (* -1/2 (* x x) (atan (/ (* y z)(* x r))))))))) (+ (J x y z) (J y z x) (J z x y)))))) (let ((P (- D Z)) (p (- (+ D Z))) (Q (- B Y)) (q (- (+ B Y))) (R (- L X)) (r (- (+ L X)))) (- (- (- (I P Q R) (I P Q r)) (- (I P q R) (I P q r))) (- (- (I p Q R) (I p Q r)) (- (I p q R) (I p q r))))))) ; potential at center of square face of long skinny stick: L by B by B. ; good optimization upon inlining. (define (stke L B) (let ((I (lambda (x y z) (let* ((r (sqrt (+ (* x x) (* y y) (* z z)))) (J (lambda (x y z) (+ (* y z (log (+ x r))) (if (zero? x) 0 (* -1/2 (* x x) (atan (/ (* y z)(* x r))))))))) (+ (J x y z) (J y z x) (J z x y)))))) (let ( (p (- B)) (r (* -2 L))) (- (- (- (I B B 0) (I B B r)) (- (I B p 0) (I B p r))) (- (- (I p B 0) (I p B r)) (- (I p p 0) (I p p r))))))) ; Faster; less pretty: (define (stke L B) (let* ((B2 (* B B)) (r (* 2 L)) (r2 (* r r)) (R (sqrt (+ (* 2 B2) r2))) (Br (* B r))) (* 4 (- (+ (* B2 (log (* B (sqrt 2)))) (* Br (log (+ B R)))) (+ (* Br (log (- R B))) (* B2 (log (- R r))) (* 1/2 r2 (atan (/ B2 (* r R)))) (* B2 (atan (/ r R)))))))) (define (f x) (let ((b (brick 1/2 1/2 1/2))) (* 1/2 (+ (b 0 0 (+ 1/2 x)) (b 0 0 (- 1/2 x)))))) (f 0.001) => 1.7928098404628092 ; tran moves field P by vector (p q r). New field is ((tran p q r) P) (define ((tran p q r) P) (lambda (x y z) (P (- x p) (- y q) (- z r)))) ; siz scales field linearly by d. Mass unchanged. (define ((siz d) P) (let ((r (/ d))) (lambda (x y z) (* r (P (* r x) (* r y) (* r z)))))) ; mm multiplies mass by d. (define ((mm d) P) (lambda (x y z) (* d (P x y z)))) Two candidate creators of field for triangle. (define ((p2c xp yp) x y z) (ni 0 1 (lambda (λ) (let ((b (- 1 λ))) (((tran (* λ xp) (* λ yp) 0) ((mm b) ((siz b) p1))) x y z))))) (define ((pqc xp yp) X Y Z) (ni 0 1 (lambda (λ) (let* ((b (- 1 λ))(r (/ b))) ; (((tran (* λ xp) (* λ yp) 0) ((mm b) ((siz b) p1))) X Y Z) ; (((mm b) ((siz b) p1)) (- X (* λ xp)) (- Y (* λ yp)) Z) ; (* b (((siz b) p1) (- X (* λ xp)) (- Y (* λ yp)) Z)) ; (* b (* r (p1 (* r (- X (* λ xp))) (* r (- Y (* λ yp))) (* r Z)))) (p1 (* r (- X (* λ xp))) (* r (- Y (* λ yp))) (* r Z)) )))) (define (ET pxc) (list ; test case: triangle with verts: (0 0 0), (1 0 0), (0 1 0). (let ((pt (pxc 0 1))) (cons (pt 1 0 1) (pt 0 1 1))) ; test case: triangle with verts: (0 0 0), (1 0 0), (1 1 0). (let ((pt (pxc 1 1))) (cons (pt 0 0 1) (pt 1 1 -1))) ; Equilateral triangle. (let* ((t (/ (sqrt 3) 2))(te (pxc 1/2 t))) (list (te 1/2 (- (/ t 3) 1/2) 0) (te (/ (- 1 t) 2) (+ (/ t 3) 1/4) 0) (te (/ (+ t 1) 2) (+ (/ t 3) 1/4) 0))) )) (ET p2c); => ((0.3966777995670832 . 0.3966825245760029) ; (0.3966777995670833 . 0.39668252457600284) ; (1.0135429731716272 1.0135867669194305 1.013586766919481)) (ET pqc); => ((0.3966777995670832 . 0.39668252457600284) ; (0.3966777995670833 . 0.39668252457600284) ; (1.0135429731716272 1.0135867669194303 1.013586766919481)) ;(p1 (* r (- X (* λ xp))) (* r (- Y (* λ yp))) (* r Z)) ;(let* ((r2 (+ (sq (* r (- X (* λ xp)))) (sq (* r (- Y (* λ yp)))) (sq (* r Z))))) ; (log (/ (+ (sqrt (+ r2 (* -2 (* r (- X (* λ xp)))) 1)) (- (* r (- X (* λ xp)))) 1) (- (sqrt r2) (* r (- X (* λ xp))))))) ;(let* ((r2 (* (sq r) (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z))))) ; (log (/ (+ (sqrt (+ r2 (* -2 (* r (- X (* λ xp)))) 1)) (- (* r (- X (* λ xp)))) 1) (- (sqrt r2) (* r (- X (* λ xp))))))) ;(let* ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) ; (log (/ (+ (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1)) (- (* r (- X (* λ xp)))) 1) (- (* r (sqrt s2)) (* r (- X (* λ xp))))))) ;(let* ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) ; (log (/ (+ (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1)) (- (* r (- X (* λ xp)))) 1) (* r (- (sqrt s2) (- X (* λ xp))))))) ;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) ; (log (/ (* b (+ (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1)) (- (* r (- X (* λ xp)))) 1)) (* b (* r (- (sqrt s2) (- X (* λ xp)))))))) ;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) ; (log (/ (+ (* b (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1))) (- (- X (* λ xp))) b) (* b (* r (- (sqrt s2) (- X (* λ xp)))))))) ;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) ; (log (/ (+ (* b (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))) ;(let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) ; (log (/ (+ (* b (sqrt (+ (* s2 (sq r)) (* r -2 (- X (* λ xp))) 1))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))) (define ((prc xp yp) X Y Z) (ni 0 1 (lambda (λ) (let* ((b (- 1 λ))) (let ((s2 (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)))) (log (/ (+ (sqrt (+ s2 (* b -2 (- X (* λ xp))) (sq b))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))) )))) ; (+ (sq (- X (* λ xp))) (sq (- Y (* λ yp))) (sq Z)) ; (+ (* λ λ (+ (sq xp) (sq yp))) (* λ (- (+ (* 2 X xp) (* 2 Y yp)))) (+ (sq X) (sq Y) (sq Z))) (define ((prc xp yp) X Y Z) (ni 0 1 (lambda (λ) (let* ((b (- 1 λ)) (r (/ b))) ; good (let ((s2 (+ (* λ λ (+ (sq xp) (sq yp))) (* λ (- (+ (* 2 X xp) (* 2 Y yp)))) (+ (sq X) (sq Y) (sq Z))) )) (log (/ (+ (sqrt (+ s2 (* b -2 (- X (* λ xp))) (sq b))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))) )))) (define ((prc xp yp) X Y Z) ; good (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (ni 0 1 (lambda (λ) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0))) (log (/ (+ (sqrt (+ s2 (* b -2 (- X (* λ xp))) (sq b))) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp)))))))))) (define ((prc xp yp) X Y Z) ; good (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (ni 0 1 (lambda (λ) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))) (log (/ (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))))))) (define ((prc xp yp) X Y Z) ; good (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))) (ni 0 1 (lambda (λ) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0))) (log (/ (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (- (* λ xp) X) b) (- (sqrt s2) (- X (* λ xp))))))))))) (define ((prc xp yp) X Y Z) ; good (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))) (n1 (- xp 1)) (n0 (- 1 X))) (ni 0 1 (lambda (λ) (log (/ (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0)) (- (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- X (* λ xp)))))))))) (define ((prc xp yp) X Y Z) ; good (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))) (n1 (- xp 1)) (n0 (- 1 X))) (ni 0 1 (lambda (λ) (- (log (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0))) (log (- (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- X (* λ xp)))))))))) (define ((prc xp yp) X Y Z) ; good (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))) (n1 (- xp 1)) (n0 (- 1 X))) (- (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0))))) (ni 0 1 (lambda (λ) (log (- (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- X (* λ xp)))))))))) (define prc (let ((int (lambda (a b c d e) ; good (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e)))))))) (lambda (xp yp) (lambda (X Y Z) (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))) (n1 (- xp 1)) (n0 (- 1 X))) (- ; (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ m2) (* λ m1) m0)) (+ (* λ n1) n0))))) (int m2 m1 m0 n1 n0) ; (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ l2) (* λ l1) l0)) (- (* λ xp) X))))) (int l2 l1 l0 xp (- X)) ))))))) (define prc (let ((int (lambda (a b c d e) ; good (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e)))))))) (lambda (xp yp) (lambda (X Y Z) (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X))) (n1 (- xp 1)) (n0 (- 1 X))) (- (int m2 m1 m0 n1 n0) (int l2 l1 l0 xp (- X))))))))) (define prc (let ((int (lambda (a b c d e) ; good (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e)))))))) (lambda (xp yp) (lambda (X Y Z) (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (- (int (+ l2 1 (* -2 xp)) (+ l1 (* 2 (+ xp X)) -2) (+ l0 1 (* -2 X)) (- xp 1) (- 1 X)) (int l2 l1 l0 xp (- X)))))))) (define prc (let ((int (lambda (a b c d e) ; good (ni 0 1 (lambda (λ) (log (+ (sqrt (+ (* λ λ a) (* λ b) c)) (+ (* λ d) e)))))))) (lambda (xp yp) (lambda (X Y Z) (let ((l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (- (int (+ l2 (* -2 xp) 1) (+ l1 (* 2 (+ xp X)) -2) (+ l0 1 (* -2 X)) (- xp 1) (- 1 X)) (int l2 l1 l0 xp (- X)))))))) ; Log(Sqrt(ax^2 + bx + c) + dx + e) ; Log(Sqrt(-2x^2 + 4x +2) + 2x + 1) (let ((pt (prc 1 1))) (cons (pt 0 0 1) (pt 1 1 -1))) (ET prc) (let ((xp 0.48)(yp 0.5)(X 3/2)(Y 0.6)(Z -0.4)) (let ((λ 1/2)(l2 (+ (sq xp) (sq yp)))(l1 (* -2 (+ (* X xp) (* Y yp))))(l0 (+ (sq X) (sq Y) (sq Z)))) (let ((b (- 1 λ)) (s2 (+ (* λ λ l2) (* λ l1) l0))) (list (+ s2 (* b -2 (- X (* λ xp))) (sq b)) (+ s2 (* b -2 (- X (* λ xp))) (sq (- 1 λ))) (+ s2 (* b -2 (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ))) (+ (+ (* λ λ l2) (* λ l1) l0) (* (- 1 λ) -2 (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ))) (+ (* λ λ l2) (* λ l1) l0 (* (- 1 λ) -2 (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ))) (+ (* λ λ l2) (* λ l1) l0 (* -2 (- 1 λ) (- X (* λ xp))) (+ 1 (* -2 λ) (* λ λ))) (+ (* λ λ l2) (* λ l1) l0 (* -2 (+ X (* λ (- (+ xp X))) (* xp λ λ))) (+ 1 (* -2 λ) (* λ λ))) (+ (* λ λ l2) (* λ l1) l0 (* -2 (+ X (* λ (- (+ xp X))) (* xp λ λ))) 1 (* -2 λ) (* λ λ)) (+ (* λ λ (+ l2 1)) (* λ l1) l0 (* -2 (+ X (* λ (- (+ xp X))) (* xp λ λ))) 1 (* -2 λ)) (+ (* λ λ (+ l2 1)) (* -2 xp λ λ) (* λ l1) (* -2 (+ X (* λ (- (+ xp X))) )) 1 (* -2 λ) l0) (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ l1) (* -2 (+ X (* λ (- (+ xp X))) )) 1 (* -2 λ) l0) (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ (+ l1 (* 2 (+ xp X)))) (* -2 X) 1 (* -2 λ) l0) (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ (+ l1 (* 2 (+ xp X)) -2)) (* -2 X) 1 l0) (+ (* λ λ (+ l2 1 (* -2 xp))) (* λ (+ l1 (* 2 (+ xp X)) -2)) (+ (* -2 X) 1 l0)) (let ((m2 (+ l2 1 (* -2 xp)))(m1 (+ l1 (* 2 (+ xp X)) -2))(m0 (+ l0 1 (* -2 X)))) (+ (* λ λ m2) (* λ m1) m0)))))) ; same