From c9aca677d707f164c7a0d9db695ca57011e8e12b Mon Sep 17 00:00:00 2001
From: GaugeAndGravity <phywyj@gmail.com>
Date: Fri, 10 Nov 2023 13:41:01 -0500
Subject: [PATCH] use uneven sample points

---
 rainbow.py | 33 +++++++++++++++++----------------
 1 file changed, 17 insertions(+), 16 deletions(-)

diff --git a/rainbow.py b/rainbow.py
index a6f5f2f..321cc4e 100644
--- a/rainbow.py
+++ b/rainbow.py
@@ -16,8 +16,7 @@ center = [0,0]
 r = 1
 # n = 1.3
 # disk = Disk(center, r, n)
-dx = 0.00001
-N = int(2*r /dx - 1)
+N = 10000
 min_intensity = 0.00001
 max_ray = 1000000
 
@@ -25,11 +24,11 @@ stack = []
 result = []
 points = []
 
-def init():
-	stack = []
-	for x in np.linspace(-r + dx, r-dx, N):
-		stack.append(Ray([x, 2*r], [0,-1], 10*np.abs(x)))
-	return stack
+# def init():
+# 	stack = []
+# 	for x in np.linspace(-r + dx, r-dx, N):
+# 		stack.append(Ray([x, 2*r], [0,-1], 10*np.abs(x)))
+# 	return stack
 
 def reflection_and_refraction(ray:Ray, intersection_point, normal, n):
 	# print("reflection/refraction at the point:", intersection_point)
@@ -110,8 +109,8 @@ def modified_trace(wavelength, temp, center, r, dx, N, n_theta, d_theta, min_int
 	n = water_refraction_index(temp, wavelength)
 	disk = Disk(center, r, n)
 	stack = []
-	for x in np.linspace(-r + dx, r-dx, N):
-		stack.append(Ray([x, 2*r], [0,-1], 10*np.abs(x)))
+	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))
@@ -163,8 +162,8 @@ def rainbow(n_theta, max_theta, temp):
 
 	maxRGB = np.max(angle_sRGB)
 	angle_sRGB/= maxRGB
+	angle_sRGB = np.clip(angle_sRGB, 0, 1)
 	colorspace.vectorized_gamma_correct(angle_sRGB)
-	np.clip(angle_sRGB, 0, 1)
 	return angles,angle_sRGB
 
 def bin_find(x, list, start):
@@ -221,10 +220,12 @@ def take_picture(angles, angle_sRGB, w, h, distance, filename="image.png"):
 	image.save(filename)
 
 if __name__=="__main__":
-	w = 7680
-	h = 4320
-	dis = 2400
+	w = 7680 *2
+	h = 4320 *2
+	dis = 2400 *2
 	max_angle = 1.1*np.arctan(np.sqrt(1+(h*h+w*w/4)/(dis*dis)))
-	angles,sRGB=rainbow(10000, max_angle, 10)
-	np.savez_compressed("saved.npz", a=angles, b=sRGB)
-	take_picture(angles,sRGB, 7680, 4320, 2400, "image.png")
+	# angles,sRGB=rainbow(10000, max_angle, 10)
+	# np.savez_compressed("saved.npz", a=angles, b=sRGB)
+	loaded = np.load("saved.npz")
+	angles, sRGB = loaded['a'], loaded['b']
+	take_picture(angles,sRGB, w, h, dis, "image16k.png")