diff -urN rt_v4/cg.py rt_v5/cg.py
--- rt_v4/cg.py	2018-03-18 17:56:56.000000000 +0900
+++ rt_v5/cg.py	2018-03-19 15:36:06.000000000 +0900
@@ -20,7 +20,7 @@
 	copen = 'copen-090419.jpg'
 	beat = 'beat.jpg'
 
-	rtd = { 'reflect': 0.3 }
+	rtd = { 'reflect': 0.3, 'reflact': 0.3, 'density': 1.1 }
 
 	data = [ {
 		'kind': 'square',
diff -urN rt_v4/mt.py rt_v5/mt.py
--- rt_v4/mt.py	2018-03-16 20:38:43.000000000 +0900
+++ rt_v5/mt.py	2018-03-19 14:49:47.000000000 +0900
@@ -34,6 +34,8 @@
 ###
 
 def quadratic_formula(a, b, c):
+	if a == 0:
+		return [ -c / b ] if b != 0 else []
 	d = b * b - 4 * a * c
 	if d < 0:
 		return []
diff -urN rt_v4/rt.py rt_v5/rt.py
--- rt_v4/rt.py	2018-03-18 18:31:18.000000000 +0900
+++ rt_v5/rt.py	2018-03-19 15:39:17.000000000 +0900
@@ -1,5 +1,6 @@
 #!/usr/bin/env python
 
+import math
 import ut
 import v
 import line
@@ -20,7 +21,8 @@
 		return []
 
 	if ref_len:
-		lst = filter( lambda (d, crs): crs.get('t') > ref_len * 0.000001, lst )
+		ref_rate = 0.000001
+		lst = filter( lambda (d, crs): crs.get('t') > ref_len * ref_rate, lst )
 		if not lst:
 			return []
 
@@ -31,10 +33,12 @@
 	base = rtd.get('base', 0.3)
 
 	crsp = crs.get('p')
+	crst = crs.get('t')
 	nv = v.unit( crs.get('nv') )
 	eyev = v.unit(l_g.v)
 	eyev_nega = v.nega(eyev)
-	ang_nv_eye = v.dot_product( nv, eyev_nega )
+	ang_nv_eyev = v.dot_product( nv, eyev )
+	ang_nv_eyev_nega = -ang_nv_eyev
 
 	def diffusivity():
 		diff = rtd.get('diff', 0.9)	
@@ -48,7 +52,7 @@
 			elif 'p' in ltd:
 				ltv = v.sub_unit( ltd.get('p'), crsp )
 			ang_nv_lt = v.dot_product( nv, ltv )
-			if ang_nv_eye * ang_nv_lt < 0:
+			if ang_nv_eyev_nega * ang_nv_lt < 0:
 				return 0
 			return max( ang_nv_lt, 0 ) * ltd.get('e', 1.0)
 
@@ -61,15 +65,43 @@
 		reflect = rtd.get('reflect', 0)
 		if not reflect:
 			return []
-		ref_v = v.op2( '+', v.op1( '*', nv, 2 * ang_nv_eye ), eyev )
+		ref_v = v.op2( '+', v.op1( '*', nv, 2 * ang_nv_eyev_nega ), eyev )
 		ref_l = line.new( crsp, ref_v )
-		col = get_col( data, lights, ref_l, nest_rate * reflect, crs.get('t') )
+		col = get_col( data, lights, ref_l, nest_rate * reflect, crst )
 		return v.op1( '*', col, reflect )
 
 	col_ = reflect_col()
 	if col_:
 		col = v.op2( '+', col, col_ )
 
+	def reflact_v():
+		density = rtd.get('density', 1)
+		if density == 1 or abs(ang_nv_eyev) == 1 or ang_nv_eyev == 0:
+			return eyev
+		cos_v = v.op1( '*', nv, ang_nv_eyev )
+		sin_v = v.sub( eyev, cos_v )
+		if ang_nv_eyev < 0:
+			density = 1.0 / density
+		r_sin_v = v.op1( '*', sin_v, density )
+		sin_len2 = v.len2( r_sin_v )
+		sin_len2 = min( sin_len2, 1 )
+		cos_len = math.sqrt( 1 - sin_len2 )
+		r_cos_v = v.op1( '*', cos_v, cos_len / v.len1( cos_v ) )
+		return v.unit( v.op2( '+', r_sin_v, r_cos_v ) )
+
+	def reflact_col():
+		reflact = rtd.get('reflact', 0)
+		if not reflact:
+			return []
+		ref_v = reflact_v()
+		ref_l = line.new( crsp, ref_v )
+		col = get_col( data, lights, ref_l, nest_rate * reflact, crst )
+		return v.op1( '*', col, reflact )
+
+	col_ = reflact_col()
+	if col_:
+		col = v.op2( '+', col, col_ )
+
 	return map( lambda v: min( int(v), 255 ), col )
 
 def draw(data, lights, eye2g, wh2eye, sc_sz, video):