calc_b working for Neumann boundary condition ("N") and experimental ("E").

eff_Drop.jl

@@ -49,29 +49,59 @@ return Ptx
 end 
 
 """
-    calc_b_spherical(c, α_i, α_o, Nb)
+    calc_b_spherical(geometry, B, c, α_i, α_o, Δxi, Δxo, x, Nb)
 
     Calculate b in Au=b. 
     # Arguments
-    - `u::Vec`: concentration vector
-    - `α::Float`: defined in solve_D
+    - `geometry::String`: "1D" | "3D"
+    - `B::String`: "N" Neumann boundaries, calc_b inside and 
+        outside the boundary, 
+        | "E" Experimental right boundary, calc_b from L to BP
+    - `c::Vec`: concentration vector
+    - `α_i::Float`: defined in solve_D
+    - `α_o::Float`: defined in solve_D
+    - `Δxi::Float`: x spacing inside
+    - `Δxo::Float`: x spacing outside
+    - `x::Vec`: x grid
     - `Nb::Int`: index of left ghost point"""
-function calc_b(geometry, c, α_i, α_o, Δxi, Δxo, x, Nb)
+function calc_b(geometry, B, c, α_i, α_o, Δxi, Δxo, x, Nb)
     b = zeros(length(c))
-    if geometry=="3D"
-        for i = 2:Nb-1
-            b[i] = α_i*((Δxi/x[i]+1)*c[i+1]+(-Δxi/x[i]+1)*c[i-1])+(1-2*α_i)*c[i]
-        end
-        for i = Nb+2:length(c)-1
-            b[i] = α_o*((Δxo/x[i]+1)*c[i+1]+(-Δxo/x[i]+1)*c[i-1])+(1-2*α_o)*c[i]
+    if B=="N"
+        if geometry=="3D"
+            for i in length(b)
+                if i>=2 && i<=Nb-1
+                    b[i] = α_i*((Δxi/x[i]+1)*c[i+1]+(-Δxi/x[i]+1)*c[i-1])+(1-2*α_i)*c[i]
+                elseif i>=Nb+2 && i<=length(x)-1
+                    b[i] = α_o*((Δxo/x[i]+1)*c[i+1]+(-Δxo/x[i]+1)*c[i-1])+(1-2*α_o)*c[i]
+                end
+            end 
+        elseif geometry=="1D"
+            for i in length(b)
+                if i>=2 && i<=Nb-1
+                    b[i] = α_i*(c[i+1]+c[i-1])+(1-2*α_i)*c[i]
+                elseif i>=Nb+2 && i<=length(x)-1
+                    b[i] = α_o*(c[i+1]+c[i-1])+(1-2*α_o)*c[i]
+                end
+            end 
         end
-    elseif geometry=="1D"
-        for i = 2:Nb-1
-            b[i] = α_i*(c[i+1]+c[i-1])+(1-2*α_i)*c[i]
+    end
+
+    if B=="E"
+        if geometry=="3D"
+            for i in length(b)
+                if i>=2 && i<Nb-1
+                    b[i] = α_i*((Δxi/x[i]+1)*c[i+1]+(-Δxi/x[i]+1)*c[i-1])+(1-2*α_i)*c[i]
+                end
+            end
         end
-        for i = Nb+2:length(c)-1
-            b[i] = α_o*(c[i+1]+c[i-1])+(1-2*α_o)*c[i]
+        if geometry=="1D"
+            for i in length(b)
+                if i>=2 && i<Nb-1
+                    b[i] = α_i*(c[i+1]+c[i-1])+(1-2*α_i)*c[i]
+                end
+            end
         end
+        b[Nb-1]=1
     end
 
     return b
@@ -89,7 +119,9 @@ end
         | "E" Experimental right boundary, make_A from L to BP
     - `α_i::Float`: defined in solve_D.
     - `α_o::Float`: defined in solve_D.
-    - `Δx::Float`: x spacing
+    - `Δxi::Float`: x spacing inside
+    - `Δxo::Float`: x spacing outside
+    - `x::Vec`: x grid
     - `D_i::Float`: diffusion coefficient inside.
     - `D_o::Float`: diffusion coefficient outside.
     - `N::Int`: size of matrix A.