(************************************************************
 *   proto.ml		Created      : Wed Mar 13 01:44:21 2002
 *                      Last modified: Wed Mar 20 12:14:18 2002
 *  Compile: make #
 ************************************************************)
(* プロトタイプ
 *  原点にある球
 *
 * 球と直線の交点 => 二次方程式の解
 *  (p - s)*(p - s) = r*r
 *  p = e + vt を解きtの小さい方をとる (負ならば裏)
 *
 *)

module Vector3D =
  struct 
    type t = Vector of float * float * float
    let length (Vector(x,y,z)) = sqrt (x*.x +. y*.y +. z*.z)
    let normalize ((Vector(x,y,z)) as v) = 
      let len = length v in
      Vector(x/.len,y/.len,z/.len)
    let (-^) (Vector(a1,a2,a3)) (Vector(b1,b2,b3)) =
      Vector(a1-.b1,a2-.b2,a3-.b3)
    let (+^) (Vector(a1,a2,a3)) (Vector(b1,b2,b3)) =
      Vector(a1+.b1,a2+.b2,a3+.b3)
    let scale t (Vector(x,y,z)) = Vector(t*.x,t*.y,t*.z)
    let inner_product (Vector(a1,a2,a3)) (Vector(b1,b2,b3)) =
      a1*.b1 +. a2*.b2 +. a3*.b3
  end

open Vector3D

let origin = Vector (0.0,0.0,0.0)
(*let object_color = Graphics.color 0x99 0x00 0x99*)
let object_red = 0x99
let object_green = 0x0
let object_blue = 0x33
let object_center_position = origin
let radius = 80.0
let screen_center_position = Vector(0.0,0.0,-100.0)
let screen_width  = 200
let screen_height = 200
let snx = Vector(1.0,0.0,0.0)
let sny = Vector(0.0,1.0,0.0)
    
let object_red_float = float_of_int object_red
let object_green_float = float_of_int object_green
let object_blue_float = float_of_int object_blue

let screen_width_half = (float_of_int screen_width)/. 2.0
let screen_height_half = (float_of_int screen_height) /. 2.0
let eye_position = Vector(0.0,0.0,-200.0)
let eye_to_screen_vector = screen_center_position -^ eye_position
let eye_to_object_vector = eye_position -^ object_center_position
let c = (inner_product eye_to_object_vector eye_to_object_vector) -. radius*.radius
    
(*let d4 nvvec =*)
(*  let t = (eye_position -^ object_center_position) in*)
(*  (inner_product (eye_position -^ object_center_position) v)*)
(*    -. (inner_product t t) +. radius*.radius*)
    
(* 外積 *)
let cross_product (Vector(a1,a2,a3)) (Vector(b1,b2,b3)) =
  let r1 = a2*.b3 -. a3*.b2
  and r2 = a3*.b1 -. a1*.b3
  and r3 = a1*.b2 -. a2*.b1 in
  Vector(r1,r2,r3)
    
(* 内積 *)
let inner_product (Vector(a1,a2,a3)) (Vector(b1,b2,b3)) =
  a1*.b1 +. a2*.b2 +. a3*.b3
    
let plot ~sx ~sy ~color =
  let rsx = int_of_float(sx+.screen_width_half) and
      rsy = int_of_float(sy+.screen_height_half) in
  begin
    Graphics.set_color color;
    Graphics.plots [| (rsx,rsy) |];
  end
    
(* あるやもしれんしないかもしれんよ
 * (sx,sy)はスクリーン上の座標??
 *
 *)
let delta = 1.0
    
let rec draw_line ~sx ~sy ~ray_vector = 
  if sx <= screen_width_half then
    let vvec = (scale sx snx) +^ (scale sy sny) +^ eye_to_screen_vector in(* 絶対 *)
    let nvvec = normalize vvec in (* 視線単位 *)
    let b = inner_product eye_to_object_vector nvvec in
    let d = b*.b -. c in
    begin
      if d > 0.0 then 
	let t = -.b -. (sqrt d) in
	if t > 0.0 then
	  let p = eye_position +^ (scale t nvvec) in             (* 交点絶対 *)
	  let npvec = normalize (p -^ object_center_position) in (* 表面の法線 *)
	  let tt = inner_product nvvec npvec in
	  let nrefvec = (scale (-2.0*.tt) npvec) +^ nvvec in
	  let ray_ratio = inner_product nrefvec ray_vector in
	  if ray_ratio > 0.0 then
	    let red = int_of_float(ray_ratio *. object_red_float) 
	    and green = int_of_float(ray_ratio *. object_green_float) 
	    and blue = int_of_float(ray_ratio *. object_blue_float) in
	    let c = Graphics.rgb red green blue in
(*	    if(red >= 0xff || green >= 0xff || blue >= 0xff) then*)
(*	      Printf.printf "(%3.0f,%3.0f) %f%% len %f/ref %f/np %f/nvvec %f/tt %f 0x%06X\n"*)
(*	      sx sy ray_ratio*)
(*	      (length ray_vector) (length nrefvec) (length npvec) (length nvvec) tt c;*)
	    plot ~sx:sx ~sy:sy ~color:c
	  else
	    plot ~sx:sx ~sy:sy ~color:Graphics.black (* ここらへんが... *)
	else
	  plot ~sx:sx ~sy:sy ~color:Graphics.black (* ここらへんが... *)
      else
	plot ~sx:sx ~sy:sy ~color:Graphics.black; (* ここらへんが... *)
      draw_line ~sx:(sx+.delta) ~sy:sy ~ray_vector:ray_vector;
    end
  else ()
  
(* 反射ベクトル *)
let rec draw_lines ~sy ~ray_vector =
  if sy < screen_height_half then
    begin
      draw_line ~sx:(-.screen_width_half) ~sy:sy ~ray_vector:ray_vector;
      draw_lines ~sy:(sy +. delta) ~ray_vector:ray_vector
    end
  else ()
      
let draw_all ray_vector =
  draw_lines ~sy:(-.screen_height_half) ~ray_vector:ray_vector
 
let _ =
  let lexbuf = Lexing.from_channel stdin in
  begin
    Graphics.open_graph " 200x200+0-0";
    while true do
      begin
	print_string "ray vector: ";
	flush stdout;
	let ix = Num_lexer.num lexbuf in
	let iy = Num_lexer.num lexbuf in
	let iz = Num_lexer.num lexbuf in
	draw_all (normalize (Vector(ix,iy,iz)));
      end
    done;
  end