Forgot to subtract surface binding energy from isotropic case.

Author: drobnyjt

src/material.rs

@@ -326,17 +326,13 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
             let dy = y2 - y;
             let dz = z2 - z;
             let mag = (dx*dx + dy*dy + dz*dz).sqrt();
-
             let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag;
-            //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.);
                     particle::surface_refraction(particle_1, Vector::new(dx/mag, dy/mag, dz/mag), Es);
                 }
-                _ => ()
+                _ => particle.E += Es;
             }
         }
     }
@@ -358,16 +354,11 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
         if costheta < 0. {
             if leaving_energy > 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.);
-                        //println!("costheta: {}", costheta);
                         particle::surface_refraction(particle_1, Vector::new(-dx/mag, -dy/mag, -dz/mag), -Es);
                     }
-                    _ => ()
+                    _ => particle.E += -Es;
                 }
 
             } else {

src/particle.rs

@@ -211,22 +211,17 @@ pub fn particle_advance(particle_1: &mut particle::Particle, mfp: f64, asymptoti
 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;
 }