Code

1import os 2import time 3 4import matplotlib 5matplotlib.use('Agg') 6import matplotlib.pyplot as plt 7import matplotlib.font_manager as fm 8import matplotlib.colors as mcolors 9 10import numpy as np 11from numba import njit 12 13@njit(fastmath=True) 14def compute_attractor(xs, ys, a, b, c, d): 15 x, y = 0.1, 0.0 16 for i in range(xs.shape[0]): 17 xn = np.sin(a * y) + c * np.cos(a * x) 18 yn = np.sin(b * x) + d * np.cos(b * y) 19 x, y = xn, yn 20 xs[i] = x 21 ys[i] = y 22 23def main(): 24 libertinus_math_path = os.path.expanduser('~/.local/share/fonts/LibertinusMath-Regular.otf') 25 font_name = 'Libertinus Math' 26 27 if os.path.exists(libertinus_math_path): 28 try: 29 fm.fontManager.addfont(libertinus_math_path) 30 print(f"Registered font: {libertinus_math_path}") 31 matplotlib.rcParams['font.family'] = [font_name, 'serif'] 32 matplotlib.rcParams['mathtext.fontset'] = 'custom' 33 matplotlib.rcParams['mathtext.rm'] = font_name 34 matplotlib.rcParams['mathtext.it'] = font_name 35 matplotlib.rcParams['mathtext.bf'] = font_name 36 print(f"Matplotlib configured to use '{font_name}' for equations.") 37 except Exception as e: 38 print(f"Error registering or configuring {font_name}: {e}") 39 print("Falling back to 'serif' for equation text.") 40 matplotlib.rcParams['font.family'] = ['serif'] 41 matplotlib.rcParams['mathtext.fontset'] = 'cm' 42 else: 43 print(f"Warning: {font_name} font not found at {libertinus_math_path}. Ensure it's installed or check the path.") 44 print("Falling back to 'serif' for equation text.") 45 matplotlib.rcParams['font.family'] = ['serif'] 46 matplotlib.rcParams['mathtext.fontset'] = 'cm' 47 48 49 n_points = 300_000_000 50 burn_in = 10_000 51 a, b, c, d = 1.8, -1.9, 1.0, 1.5 52 53 xs = np.empty(n_points, dtype=np.float64) 54 ys = np.empty(n_points, dtype=np.float64) 55 56 t0 = time.perf_counter() 57 compute_attractor(xs, ys, a, b, c, d) 58 t1 = time.perf_counter() 59 60 xs = xs[burn_in:] 61 ys = ys[burn_in:] 62 63 print(f"Computed {n_points} points in {t1 - t0:.2f} s.") 64 print(f"x ∈ [{xs.min():.3f}, {xs.max():.3f}], y ∈ [{ys.min():.3f}, {ys.max():.3f}]") 65 66 dark_bg = '#111111' 67 fig = plt.figure(figsize=(12, 8), facecolor=dark_bg) 68 ax = fig.add_axes((0, 0, 1, 1), facecolor=dark_bg) 69 x_min, x_max = -2.5, 2.5 70 y_min, y_max = -2.5, 2.5 71 ax.set_xlim(x_min, x_max) 72 ax.set_ylim(y_min, y_max) 73 ax.axis('off') 74 75 H, xedges, yedges = np.histogram2d(xs, ys, bins=2000, range=[[x_min, x_max], [y_min, y_max]]) 76 77 ax.imshow( 78 H.T, 79 origin='lower', 80 extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]], 81 cmap='inferno', 82 norm=mcolors.LogNorm(vmin=1), 83 aspect='auto', 84 interpolation='bilinear' 85 ) 86 87 eq_text = ( 88 r"$x_{t+1} = \sin(%.1f y_t) + %.1f \cos(%.1f x_t)$" "\n" 89 r"$y_{t+1} = \sin(%.1f x_t) + %.1f \cos(%.1f y_t)$" 90 ) % (a, c, a, b, d, b) 91 92 fig.text( 93 0.01, 0.99, eq_text, 94 color='#BF40BF', 95 fontsize=16, 96 va='top', 97 ha='left', 98 alpha=0.8 99 ) 100 101 interval_text = ( 102 r"$x \in [%.3f, %.3f]$" "\n" 103 r"$y \in [%.3f, %.3f]$" 104 ) % (xs.min(), xs.max(), ys.min(), ys.max()) 105 106 fig.text( 107 0.01, 0.92, 108 interval_text, 109 color='#BF40BF', 110 fontsize=14, 111 va='top', 112 ha='left', 113 alpha=0.8 114 ) 115 116 outpath = 'clifford_attractor.png' 117 118 plt.savefig( 119 outpath, 120 dpi=300, 121 bbox_inches='tight', 122 pad_inches=0, 123 facecolor=fig.get_facecolor() 124 ) 125 plt.close(fig) 126 127if __name__ == '__main__': 128 main()

1#version 300 es 2precision highp float; 3 4uniform vec2 u_resolution; 5uniform vec2 u_center; 6uniform float u_zoom; 7uniform vec2 u_c; 8 9out vec4 fragColor; 10 11vec3 palette(float t) { 12 // Attempt at a continuous, smooth coloring. For each color channel, 13 // the expression is: 0.5 + 0.5 * cos(2π * (t + offset)) 14 return 0.5 + 0.5 * cos(6.28318 * (t + vec3(0.0, 0.1, 0.2))); 15} 16 17void main() { 18 // Map fragment coordinates to complex plane 19 vec2 uv = (gl_FragCoord.xy - u_resolution * 0.5) / min(u_resolution.x, u_resolution.y); 20 vec2 z = uv / u_zoom + u_center; 21 22 // Iterate z = z^2 + c 23 float i; 24 for(i = 0.0; i < 256.0; i++) { 25 if(dot(z, z) > 4.0) break; // |z|^2 > 4 means escape 26 z = vec2(z.x*z.x - z.y*z.y, 2.0*z.x*z.y) + u_c; 27 } 28 29 if(i >= 256.0) { 30 // Point is in the Julia set (didn't escape) 31 fragColor = vec4(0, 0, 0, 1); 32 } else { 33 // Smooth iteration count to avoid color banding 34 // smooth_i = i - log2(log2(|z|^2)) + 4 35 float smooth_i = i - log2(log2(dot(z, z))) + 4.0; 36 fragColor = vec4(palette(smooth_i * 0.02), 1); 37 } 38} Mandelbrot Picker 1#version 300 es 2precision highp float; 3 4uniform vec2 u_resolution; 5uniform vec2 u_c; // Current c value to mark 6 7out vec4 fragColor; 8 9void main() { 10 vec2 uv = (gl_FragCoord.xy - u_resolution * 0.5) / min(u_resolution.x, u_resolution.y); 11 12 // Mandelbrot uses c as the varying parameter, z starts at 0 13 vec2 c = uv * 2.5 - vec2(0.5, 0.0); // Shift to show main cardioid 14 vec2 z = vec2(0); 15 16 float i; 17 for(i = 0.0; i < 128.0; i++) { 18 if(dot(z, z) > 4.0) break; 19 z = vec2(z.x*z.x - z.y*z.y, 2.0*z.x*z.y) + c; 20 } 21 22 vec3 col; 23 if(i >= 128.0) { 24 col = vec3(0.1, 0.1, 0.15); // Inside the set 25 } else { 26 col = vec3(0.2, 0.3, 0.4) * (i / 128.0); // Outside 27 } 28 29 // Draw red marker for current c value 30 vec2 cPos = (u_c + vec2(0.5, 0.0)) / 2.5; 31 vec2 pixelC = (cPos + 0.5) * u_resolution; 32 if(length(gl_FragCoord.xy - pixelC) < 4.0) { 33 col = vec3(1, 0.3, 0.3); 34 } 35 36 fragColor = vec4(col, 1); 37}

1float h(vec2 p){p=fract(p*vec2(123.34,456.21));p+=dot(p,p+45.32);return fract(p.x*p.y);} 2vec3 m(float t){ 3 t=clamp(t,0.,1.); 4 return t<.25?mix(vec3(.02,.05,.2),vec3(.1,.2,.5),t/.25):t<.5?mix(vec3(.1,.2,.5),vec3(0.,.2,.8),(t-.25)/.25):t<.85?mix(vec3(0.,.2,.8),vec3(.3,.7,1.),(t-.5)/.35):t<.97?mix(vec3(.3,.7,1.),vec3(.7,.9,1.),(t-.85)/.12):vec3(.7,.9,1.); 5} 6void mainImage(out vec4 o,in vec2 f){ 7 vec2 r=iResolution.xy,u=(floor(f/6.)*6.+3.)/r,c=(u-.5)*2.; 8 c.x*=r.x/r.y; 9 float d=length(c); 10 if(d>.92){o=vec4(0);return;} 11 vec2 w=c; 12 w.x+=.4*sin(w.y*2.+iTime*1.2); 13 w.y+=.4*cos(w.x*2.-iTime); 14 w+=.15*vec2(sin(w.y*5.+iTime*.8),cos(w.x*5.-iTime*.8)); 15 w+=.05*vec2(h(floor(w*30.)),h(floor(w*30.+10.))); 16 vec2 q=u-vec2(.65); 17 float s=length(q),a=atan(q.y,q.x)+5.*s; 18 w*=mix(.85,1.,smoothstep(0.,.35,s)); 19 float v=mix(sin(8.*w.x+10.*w.y+iTime*1.5),sin(5.*a+iTime*1.2+2.*sin(5.*s-iTime*2.)),smoothstep(.5,0.,s)); 20 o=vec4(m(v*.5+.5)*smoothstep(.92,.7,d),1.); 21}