Merge pull request #377 from bshanahan/sin_alpha-as-F3D

Change sin_alpha to a Field3D in sheath_boundary

Author: bendudson

include/sheath_boundary.hxx

@@ -77,7 +77,7 @@ struct SheathBoundary : public Component {
   void transform(Options &state) override;
 private:
   BoutReal Ge; // Secondary electron emission coefficient
-  BoutReal sin_alpha; // sin of angle between magnetic field and wall.
+  Field3D sin_alpha; // sin of angle between magnetic field and wall.
 
   bool lower_y; // Boundary on lower y?
   bool upper_y; // Boundary on upper y?

src/sheath_boundary.cxx

@@ -61,9 +61,9 @@ SheathBoundary::SheathBoundary(std::string name, Options& alloptions, Solver*) {
   sin_alpha = options["sin_alpha"]
                   .doc("Sin of the angle between magnetic field line and wall surface. "
                        "Should be between 0 and 1")
-                  .withDefault(1.0);
+    .withDefault(Field3D(1.0));
 
-  if ((sin_alpha < 0.0) or (sin_alpha > 1.0)) {
+  if ((min(sin_alpha) < 0.0) or (max(sin_alpha) > 1.0)) {
     throw BoutException("Range of sin_alpha must be between 0 and 1");
   }
 
@@ -202,7 +202,7 @@ void SheathBoundary::transform(Options& state) {
                 (adiabatic * ti + Zi * s_i * te * grad_ne / grad_ni) / Mi, 0., 100.);
 
             // Note: Vzi = C_i * sin(α)
-            ion_sum[i] += s_i * Zi * sin_alpha * sqrt(C_i_sq);
+            ion_sum[i] += s_i * Zi * sin_alpha[ip] * sqrt(C_i_sq);
           }
         }
       }
@@ -238,7 +238,7 @@ void SheathBoundary::transform(Options& state) {
             BoutReal C_i_sq = std::clamp(
                 (adiabatic * ti + Zi * s_i * te * grad_ne / grad_ni) / Mi, 0., 100.);
 
-            ion_sum[i] += s_i * Zi * sin_alpha * sqrt(C_i_sq);
+            ion_sum[i] += s_i * Zi * sin_alpha[im] * sqrt(C_i_sq);
           }
         }
       }