diff --git a/rainbow.py b/rainbow.py index 94ad879..5fdb1da 100644 --- a/rainbow.py +++ b/rainbow.py @@ -18,7 +18,7 @@ r = 1 # disk = Disk(center, r, n) N = 100000 min_intensity = 0.00001 -max_ray = 100*N +max_ray = 100 stack = [] result = [] @@ -109,31 +109,31 @@ def modified_trace(wavelength, temp, center, r, N, n_theta, d_theta, min_intensi n = water_refraction_index(temp, wavelength) disk = Disk(center, r, n) stack = [] - for u in np.linspace(0, 1, N, endpoint=False): - stack.append(Ray([r*np.sqrt(u), 2*r], [0,-1], 1)) - ray_count = 0 XYZ = colorspace.wavelength2XYZ(wavelength)*sun_spectral(wavelength) own_angle_XYZ = np.zeros((n_theta, 3)) - while stack and ray_count < max_ray: - ray = stack.pop() - ray_count += 1 - if ray is None: - continue - if isinstance(ray, Ray): - if ray.intensity < min_intensity: + for u in np.linspace(0, 1, N, endpoint=False): + stack=[Ray([r*np.sqrt(u), 2*r], [0,-1], 1)] + ray_count = 0 + while stack and ray_count < max_ray: + ray = stack.pop() + ray_count += 1 + if ray is None: continue - direction = ray.direction - t,intersection_point = disk.find_intersection(ray) - if intersection_point is not None: - points.append(np.concatenate((ray.origin, intersection_point,[ray.intensity]))) - normal = disk.get_normal(intersection_point) - stack.extend(reflection_and_refraction(ray, intersection_point, normal, disk.refractive_index)) - else: - points.append(np.concatenate((ray.origin, ray.origin+ray.direction,[ray.intensity]))) - if direction[1] > 0: - angle = np.arccos(direction[1]) - if angle < max_angle: - own_angle_XYZ[int(angle/d_theta)] += XYZ*ray.intensity*direction[1] + if isinstance(ray, Ray): + if ray.intensity < min_intensity: + continue + direction = ray.direction + t,intersection_point = disk.find_intersection(ray) + if intersection_point is not None: + points.append(np.concatenate((ray.origin, intersection_point,[ray.intensity]))) + normal = disk.get_normal(intersection_point) + stack.extend(reflection_and_refraction(ray, intersection_point, normal, disk.refractive_index)) + else: + points.append(np.concatenate((ray.origin, ray.origin+ray.direction,[ray.intensity]))) + if direction[1] > 0: + angle = np.arccos(direction[1]) + if angle < max_angle: + own_angle_XYZ[int(angle/d_theta)] += XYZ*ray.intensity*direction[1] return own_angle_XYZ def rainbow(n_theta, max_theta, temp):