open Graphics;;
set_window_title "Ergodic Folium";
let xs = 750 and e = 1.85
and ( * )=( *. ) and (-)=(-.) and (+)=(+.) and (/)=(/.) in
let xm = float xs in
let prm t = let q = e*t / (1. + t*t*t) in q, q*t
and hpi = 2.*(atan 1.)
and 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 fi r = int_of_float (xm*r) and fl r = (float r) / xm in
let pp t = prm (tan (t * hpi)) in
let fdb pl =
  let fs k p = let rec fz j = if p j or j>xs then j else fz (succ j) in fz k in
  let a = Array.init xs (fun j -> let b = fs 0 (fun k -> pl j k) in
  (b, fs b (fun k -> not (pl j k)))) in
  (fun _ -> set_color 0; for j = 0 to pred xs do let x1, x2 = a.(j) in
    moveto x1 j; lineto x2 j done) in

let rec ray x y xd yd n = if n=0 then () 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 (-. b + d) / (2.*a) in
   let nx = x + z*xd and ny = y + z*yd in lineto (fi nx) (fi 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
let nb = fdb (fun i j -> (p (fl i) (fl j)) < 0.) in
open_graph ""; resize_window xs xs; auto_synchronize (2<1);
while 0<1 do
  let e = wait_next_event [Button_down; Key_pressed] in
  let x, y = pp (fl e.mouse_x) and xn, yn = pp (fl e.mouse_y) in
  Printf.printf "%18.16f %18.16f %18.16f %18.16f \n" x y xn yn; flush stdout;
  clear_graph (); nb (); set_color 0xff6060;
  moveto (fi x) (fi y); ray x y (xn-x) (yn-y) 200; synchronize ()
done;;