open Graphics;; open_graph ""; set_window_title "Ergodic Folium"; let xs = 750 and ( * )=( *. ) and (-)=(-.) and (+)=(+.) and (/)=(/.) in let xm = float xs and e = 1.85 in let p x y = x*x*x + y*y*y - e*x*y and g x y = 3.*x*x - e*y, 3.*y*y - e*x and prm t = let q = e*t / (1. + t*t*t) in q, q*t and i r = int_of_float (xm*r) in let dr x y xe ye = moveto (i x) (i y); lineto (i xe) (i ye) in resize_window xs xs; set_color 0; auto_synchronize (2<1); let x, y = prm 0.5 in (let c = p x y in for i = 0 to pred xs do for j = 0 to pred xs do if (p ((float i) / xm) ((float j) / xm)) < c then plot i j done done); synchronize (); set_color 0xff6060; let rec ray x y xd yd n = if n=0 then (synchronize (); print_char (read_key ())) else let a = xd*xd*xd + yd*yd*yd and b = 3.*(x*xd*xd + y*yd*yd) - e*xd*yd and c = (3.*(x*x*xd + y*y*yd) - e*(x*yd + y*xd)) in let d = (sqrt (b*b - 4.*a*c)) in let z = if (abs_float a) < (abs_float c) then 2.*c / (-. (b + d)) else (d - b) / (2.*a) in let nx = x + z*xd and ny = y + z*yd in dr x y nx ny; let px, py = g nx ny in let r = 2.*(xd*px + yd*py) / (px*px + py*py) in ray nx ny (xd - r*px) (yd - r*py) (pred n) in ray x y 1. 0.7 500;;