From 386db2d6b84c09c403b2e44c9763ac44cf43942b Mon Sep 17 00:00:00 2001 From: GaugeAndGravity Date: Mon, 30 Oct 2023 00:41:13 -0400 Subject: [PATCH] init --- __pycache__/raytrace_2D.cpython-311.pyc | Bin 0 -> 6336 bytes points.csv | 380 ++++++++++++++++++++++++ rainbow.py | 78 +++++ raytrace_2D.py | 78 +++++ result.csv | 5 + 5 files changed, 541 insertions(+) create mode 100644 __pycache__/raytrace_2D.cpython-311.pyc create mode 100644 points.csv create mode 100644 rainbow.py create mode 100644 raytrace_2D.py create mode 100644 result.csv diff --git a/__pycache__/raytrace_2D.cpython-311.pyc b/__pycache__/raytrace_2D.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..55b06964a0e61b4af556ebab1b1be940d10f49c1 GIT binary patch literal 6336 zcmb_gTWk|o8lL-%FA2F30u%`3(%K<}K+6)gMGa+JS{7?*TB)0D*K#})aB!UTj6)mO zMnh$(R^>{eN<~&eT0bB*0jl!2(rPcSEA4{^S;9Ljv`DKy^erenRMkH0|DUnPGjUSw zN_*`2=gj%fZT`#m|7U*J(&A^J{Po|bwLkk9=3n?wFJi5-aR(|l8I@7lBs0mvolA1F z{3JgsObRT+F>f#`pJ!Cz4#zO};Ihl4NI4>KB#pnr0}rkmPog{-c*>?xGUfPy;MYGlKJu zNW*;fGUKdY!quBfcGK8QfMK#sw5~nS+XS=5MnaOsrRFU$Mo4VYlIBR_QX*wZ$wVra zykrTfIWlYc)CAGu>BL+L-d)hO>jLi4;F-C;PLLgw|Yt;bNAuO?!#Y``-|UPyL;{T33KFJW#pW>`yBLXEx>OBhJ&sJ zTq6~m)uK_$AC1n=sS8Q02cyvs7h*}9L)xM58bdnZM%s`_b%IjWUjSk_zjcli-(=UB zt*Jc;j*r9Tl)9caPC}E*pk(k4VUpJbRnWv~P8IV4del8<;H0EVV01DyI>nMEF309I z-I9owULdJR+%X;)fpuJ{X0=%sRP9WJswN6NhFGhOXwzu!cxbx+vnuW-A{^iN;j;z7%veW9__eW6$oENVL;-SAvr z`Z=h(m$urcVQa9xG#Hd^c$uEsMw&Q594rGdm?xZM^wW*}_j_2gU*8hv=U_e0F8r6Y zbQB#Rmdv!*mTxqn&^GLKtIoPG+iqx&gcQ;VcY+u~b|OJ(5&;P{u>Dv&gaifQGQS>t zL;|qO=A=Z{e4y5M0oI@t^3avqb#ePY=- z48&!!ifyv00A?wuq9y`_OVffXgCW~y>H{;4U)E+dTMizObZqo3I;+=4vm98g3mQw< zyQz2`gl*Up0F$KK{vm)3>XuT+n-Fi?H+nJ&0jVRF$I9uR9Vp{rRD4 zYX{sypy|(sKOMe4k{>bX&$c#a-!KH#OD+CSP`L?^I|FySoMk?EC&#Mn`@qj}IX))< z*BJ}ka8AsMTfh0{;PF1lCnx3PEO%YWO4ox~Da&1<%Mc0ba#qXc%5nx;?Rg<+N{zKfn3;vykF5&Fnb8=yY8_M{LMXPl(Wfmo zZFRW4I7*$&3QQ-6o{my|S)zU8(xX;qbKg-NaSP%LEw;+S&;r0@tE{f(UC^li76_=U zrM>veb;iddR3kOw)fY~h{D3i9mU>o(s%>4%fpvx#kdM?@sO~%TXy1zu_r3UF!rV7j z**9hidyF$>c~4bVjK!s?rPS)N2khMwt0$`A!{sC6e|c%;d|BCJT>QQ&$ws0o1fK}X zBjM0P;ZW)LYM}hm8>VosBAhGJ|9WR*FSf^_bjVR49tS~Fpr*4WSn|XfMUFsEuN|YO zW_s&F*-kqYar&UEy^VgTZ39&(WE-j|BH{F#p=uD@kYMm^f-0Wt0O@*&;^$E0F4~RA zb)4eS>U0B^wyN1X#^vI0Ii(I08z{}sTwqmTroq+n3MOru0Pz*0tw^CuDli|jlQI~% zXbEQ%`jvPc@Ninx0M)|hQ2C0z$+&Z~JJ8mM&?U!tq1h?mPSD=S^H?VBWL=OC;@q@s zc4+`WXCD*}ol#EMgltXJ7v?o$#}kYZ#1F)h;u=OvNDg9ZVnGjywilkn4r53#yj#9# zG@gvRd@JZ#eJl?aS!vmmJ`5I_* zsxYygop6xLgRB~7DECV4I$cN1fH1F6%OnS~Lclw#;Yz!T(`2 zXv;e~3;C`4-8_1h-SVgMfYozKR!KWr0Z+l~ik1=vfu@pGHoXDvT7l4~Oye<4id?D< zdG$M-l^|5CagoTSKm760-~XBY$15S9eU71Cg1%(&>Rg&0YpBz}S?5ENCYWK879h}{ zrph$7adUBtjalrJ#m23mn$Y7UF$<@rRGJLq01DAKDS8<#wx?nI%yd&vYszKAz|+H^ zItXI2Kc87TyL{f14;Y6YOG`JbNg53uA`=O6b7|-EREjiT>5aO{f53WdSl_nUgL08 zY+o9GDDEtaJFEV#Lik&{Ty2IG;Z~re-o9e?j6V38*)#S~8Y@d{&doCcyA1DX+m;L+8{C)~D9(+MvQ+!E51@Q$u z1Zs9cpM+-zx2dIre|hZALFkgp8tXH__$^M>Ihdn<$P!5db}2$P45*fi9)!s($v#?m z!}uZKWOh*NXdE?OhDQAu5HO2iJ2i^nfl}A)y=LIBF<1?>e|2X0&6|-gBc&5&SGdv@ zHUnX3_w6u-sd(15K8+w05a%FbwbX(`TuUaS(KVKw!evoXuiBRT681+jfQ%_15E<&f zb<DVZd)9N;`J|PbYBOO!|2^y}+35i5)|5(KuXAscUEhrQ(kIs!ZcPpASOZ7iipv zb2Y?$JV`g=zGC?=$CA_0`D83rKdQKOjdf~681Mf-IRHKGAOz6@7Uv)_Af2;z2_i&u z*~;P|TX1G*3lpqX8sISQK;>~McUMo$TyM`;fhRwENOcs75OK>#sH(11*gEhPtqDVuc zuFT3CEzvhs_6qGuIQ*L8&gfjTx%2=tJDR9rM&0HZ@9bHQVhEq*-DBvJv$uMLqj8Ck z&`3mHffox&wt`MW-Htm{A5=3C%@{a&IT44@P8eHs6b@?16;YmZDXA&0;H3Jmfj}mw z_={3OFXmR%TUSgaTv5VCUsY5JVl?{Ah}jY<`$J_uM4i=x#Vo8z{Rlr9c@M*Nv}CDGl3wWG z{`x{cj%0Sa-l~Cz7qBk?IlNBQ*HYYD=DKX0=!Z7)29h_i3$9+nuJhQ41Z`Hw#9$M4 z9YGx0nkW(%MgW7`X41-XJMplxJ+U{2D~J4+m|B>fU!*Fx_}EmOpln=Ug?jLD5DCRv zx@qtCp&yi&?GKr!$wlagZ>al#uJbI*R+)}+{a 0: + n1,n2 = n,1 + else: + n1, n2 = 1,n + kp = ray.direction - k_n*normal + k2p = n1*kp/n2 + if np.dot(k2p, k2p) >= 1: + ray1 = Ray(intersection_point, ray.direction-2*k_n*normal, ray.intensity) + return ray1,None + else: + k2n = np.sqrt(1 - np.dot(k2p, k2p))*np.sign(k_n) + Rs = (n1*k_n - n2*k2n)*(n1*k_n - n2*k2n)/((n1*k_n + n2*k2n)*(n1*k_n + n2*k2n)) + Rp = (n1*k2n - n2*k_n)*(n1*k2n - n2*k_n)/((n1*k2n + n2*k_n)*(n1*k2n + n2*k_n)) + R = (Rs+Rp)/2 + ray1 = Ray(intersection_point, ray.direction-2*k_n*normal, ray.intensity*R) + ray2 = Ray(intersection_point, k2n*normal+k2p, ray.intensity*(1-R)) + return ray1,ray2 + +def trace(disk:Disk, stack:list): + while stack: + ray = stack.pop() + if ray is None: + continue + 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: + result.append([np.arccos(direction[1]), ray.intensity*direction[1]]) + +# init() +stack.append(Ray([0.8, 2*r], [0,-1], 1)) +trace(disk, stack) + +with open("result.csv", 'w', newline='') as f: + writer = csv.writer(f) + writer.writerows(result) + +with open("points.csv", 'w', newline='') as f: + writer = csv.writer(f) + writer.writerows(points) diff --git a/raytrace_2D.py b/raytrace_2D.py new file mode 100644 index 0000000..d7e37a5 --- /dev/null +++ b/raytrace_2D.py @@ -0,0 +1,78 @@ +import numpy as np +from abc import ABC, abstractmethod + +class Ray: + def __init__(self, origin, direction): + self.origin = np.array(origin) + self.direction = np.array(direction) / np.linalg.norm(direction) + +class Shapes(ABC): + def __init__(self,refractive_index) -> None: + self.refractive_index = refractive_index + @abstractmethod + def find_intersection(self, ray:Ray): + pass + def get_normal(self, point): + pass + +class Scheme: + def __init__(self, shapes, rays) -> None: + self.shapes = shapes + self.rays = rays + self.raytrace() + def raytrace(self): + for ray in self.rays: + t, intersection_point, first_shape = None, None, None + for shape in self.shapes: + t1,intersection_point1 = shape.find_intersection(ray) + if t is None or ((t1 is not None) and (t1 < t)): + t, intersection_point, first_shape = t1, intersection_point1, shape + if t is None: + continue + else: + self.reflection_and_refraction(ray, first_shape, intersection_point) + def reflection_and_refraction(ray:Ray, shape:Shapes, intersection_point): + pass + +class Disk(Shapes): + def __init__(self, center, radius, refractive_index): + super().__init__(refractive_index) + self.center = np.array(center) + self.radius = radius + def find_intersection(self, ray:Ray): + oc = ray.origin - self.center + a = np.dot(ray.direction, ray.direction) + b = 2.0 * np.dot(oc, ray.direction) + c = np.dot(oc, oc) - self.radius * self.radius + discriminant = b*b - 4*a*c + if discriminant < 0: + return None, None # 没有交点 + else: + if np.dot(oc, oc) <= self.radius * self.radius + 1e-10: + t = (-b + np.sqrt(discriminant)) / (2.0 * a) + else: + t = (-b - np.sqrt(discriminant)) / (2.0 * a) + if t < 0: + return None, None + intersection_point = ray.origin + t * ray.direction + intersection_point = intersection_point/np.linalg.norm(intersection_point) + return t,intersection_point + def get_normal(self, point): + normal = (point - self.center)/self.radius + normal = normal / np.linalg.norm(normal) + return normal + +class half_plane(Shapes): + def __init__(self, point, normal, refractive_index) -> None: + super().__init__(refractive_index) + self.point = point + self.normal = normal/np.linalg.norm(normal) + def find_intersection(self, ray:Ray): + dot = np.dot(ray.direction, self.normal) + if dot == 0: + return None,None + else: + t = np.dot(self.point - ray.origin, self.normal) / dot + return t, ray.origin + t*ray.direction + def get_normal(self, point): + return self.normal diff --git a/result.csv b/result.csv new file mode 100644 index 0000000..22349c0 --- /dev/null +++ b/result.csv @@ -0,0 +1,5 @@ +0.7993010857076976,0.0007048815056688454 +1.3928938872949743,0.003657341850518903 +1.2482567878728648,0.0003454071966921608 +0.8743957180247236,0.0006520310605573929 +0.1162362891411113,0.0010096197790698755