New surface refraction code. Much more mathematically sound than the snell's law form I was using before. It should be guaranteed to work for any direction the particles are moving out of a surface with any direction of surface normal.

Author: drobnyjt

src/material.rs

@@ -328,12 +328,13 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
             let mag = (dx*dx + dy*dy + dz*dz).sqrt();
 
             let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag;
-            particle_1.E += Es;
+            //particle_1.E += Es;
 
             match material.surface_binding_model {
                 SurfaceBindingModel::PLANAR{..} => {
-                    let delta_theta = particle::refraction_angle(costheta, E, E + Es);
-                    particle::rotate_particle(particle_1, delta_theta, 0.);
+                    //let delta_theta = particle::refraction_angle(costheta, E, E + Es);
+                    //particle::rotate_particle(particle_1, delta_theta, 0.);
+                    particle::surface_refraction(particle_1, Vector::new(dx/mag, dy/mag, dz/mag), Es);
                 }
                 _ => ()
             }
@@ -348,17 +349,23 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
         let dz = z2 - z;
         let mag = (dx*dx + dy*dy + dz*dz).sqrt();
         let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag;
-        let leaving_energy = E*costheta*costheta;
+
+        let leaving_energy = match material.surface_binding_model {
+            SurfaceBindingModel::ISOTROPIC{..} => E,
+            SurfaceBindingModel::PLANAR{..} => E*costheta*costheta,
+        };
 
         if costheta < 0. {
             if leaving_energy > Es {
 
-                particle_1.E += -Es;
+                //particle_1.E += -Es;
 
                 match material.surface_binding_model {
                     SurfaceBindingModel::PLANAR{..} => {
-                        let delta_theta = particle::refraction_angle(costheta, E, E + Es);
-                        particle::rotate_particle(particle_1, delta_theta, 0.);
+                        //let delta_theta = particle::refraction_angle(costheta, E, E + Es);
+                        //particle::rotate_particle(particle_1, delta_theta, 0.);
+                        //println!("costheta: {}", costheta);
+                        particle::surface_refraction(particle_1, Vector::new(-dx/mag, -dy/mag, -dz/mag), -Es);
                     }
                     _ => ()
                 }

src/particle.rs

@@ -208,6 +208,28 @@ pub fn particle_advance(particle_1: &mut particle::Particle, mfp: f64, asymptoti
     return distance_traveled;
 }
 
+pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) {
+    let E = particle.E;
+    let costheta = particle.dir.dot(&normal);
+    //println!("costheta: {}", costheta);
+
+    let a = (E/(E + Es)).sqrt();
+    let b = -(E).sqrt()*costheta;
+    let c = (E*costheta.powi(2) + Es).sqrt();
+    //println!("0: {} {} {}", particle.dir.x, particle.dir.y, particle.dir.z);
+    //println!("1: {} {} {}", a, b, c);
+
+    let u1x = (E/(E + Es)).sqrt()*particle.dir.x + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.x;
+    let u1y = (E/(E + Es)).sqrt()*particle.dir.y + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.y;
+    let u1z = (E/(E + Es)).sqrt()*particle.dir.z + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.z;
+    particle.dir.x = u1x;
+    particle.dir.y = u1y;
+    particle.dir.z = u1z;
+    //println!("2: {} {} {}", particle.dir.x, particle.dir.y, particle.dir.z);
+
+    particle.E += Es;
+}
+
 /// Calcualte the refraction angle based on the surface binding energy of the material.
 pub fn refraction_angle(costheta: f64, energy_old: f64, energy_new: f64) -> f64 {
     let costheta = if costheta.abs() > 1. {costheta.signum()} else {costheta};

src/structs.rs

@@ -38,6 +38,10 @@ impl Vector {
     pub fn add(&self, other: &Vector) -> Vector {
         Vector::new(self.x + other.x, self.y + other.y, self.z + other.z)
     }
+
+    pub fn dot(&self, other: &Vector) -> f64 {
+        self.x*other.x + self.y*other.y + self.z*other.z
+    }
 }
 
 /// Vector4 is a trajectory-tracking object that includes x, y, z, and the current energy.

src/tests.rs

@@ -475,8 +475,11 @@ fn test_surface_refraction() {
     let cosy_new = cosy_new/dir_mag;
     let cosz_new = cosz_new/dir_mag;
 
-    let delta_theta = particle::refraction_angle(cosx, E, E + Es);
-    particle::rotate_particle(&mut particle_1, delta_theta, 0.);
+    //let delta_theta = particle::refraction_angle(cosx, E, E + Es);
+    //particle::rotate_particle(&mut particle_1, delta_theta, 0.);
+
+    let normal = Vector::new(1.0, 0.0, 0.0);
+    particle::surface_refraction(&mut particle_1, normal, Es);
 
     if print_output {
         println!("dir_mag: {}", dir_mag);
@@ -502,17 +505,22 @@ fn test_surface_refraction() {
         println!("{} {} {}", particle_1.dir.x, particle_1.dir.y, particle_1.dir.z);
     }
 
-    let delta_theta = particle::refraction_angle(particle_1.dir.x, particle_1.E, particle_1.E - Es);
-    particle::rotate_particle(&mut particle_1, delta_theta, 0.);
+    //let delta_theta = particle::refraction_angle(particle_1.dir.x, particle_1.E, particle_1.E - Es);
+    //particle::rotate_particle(&mut particle_1, delta_theta, 0.);
+
+    let normal = Vector::new(1.0, 0.0, 0.0);
+    particle::surface_refraction(&mut particle_1, normal, -Es);
+
 
     if print_output {
         println!("{} {} {}", cosx_new, cosy_new, cosz_new);
         println!("{} {} {}", particle_1.dir.x, particle_1.dir.y, particle_1.dir.z);
     }
 
-    assert!(approx_eq!(f64, particle_1.dir.x, cosx_new, epsilon=1E-12));
-    assert!(approx_eq!(f64, particle_1.dir.y, cosy_new, epsilon=1E-12));
-    assert!(approx_eq!(f64, particle_1.dir.z, cosz_new, epsilon=1E-12));
+    assert!(approx_eq!(f64, particle_1.dir.x, cosx_new, epsilon=1E-9));
+    assert!(approx_eq!(f64, particle_1.dir.y, cosy_new, epsilon=1E-9));
+    assert!(approx_eq!(f64, particle_1.dir.z, cosz_new, epsilon=1E-9));
+
 }
 
 #[test]