Fipy: Beer-Lambert law with multiple reflections

multipass_beer_lambert.py

# http://matforge.org/wd15/blog
from __future__ import division
from fipy import *

# Grid things
N = 1000
D = 1
dx = D / N
m = Grid1D(nx=N+1, dx=dx) + [[-dx/2.]]
I_right = CellVariable(mesh=m, value=1., name='I_right') # Intensity along the positive direction
I_left  = CellVariable(mesh=m, value=1., name='I_left')  # Intensity along the positive direction

# Constants
A = 0.
B = 1.
C = (1 - A)
alpha = CellVariable(mesh=m, value=5., name="absorption coefficient")

x = m.cellCenters[0]
left_anchor = 1e+10 * (x < dx)
right_anchor = 1e+10 * (x > D - dx)

# Equations to solve
## introduce the coupled boundary conditions implicitly
eq1 =  (ExponentialConvectionTerm(coeff=[[1]], var=I_right) 
        + ImplicitSourceTerm(alpha, var=I_right) 
        + ImplicitSourceTerm(left_anchor, var=I_right)
        - ImplicitSourceTerm(left_anchor * A, var=I_left)
        - C * left_anchor 
        == 0)
eq2 = (-ExponentialConvectionTerm(coeff=[[1]], var=I_left) 
       + ImplicitSourceTerm(alpha, var=I_left)
       + ImplicitSourceTerm(right_anchor, var=I_left) 
       - ImplicitSourceTerm(right_anchor * B, var=I_right)
       == 0)

# Beer-lambert (holds when A and B are small)
beer_lambert  = CellVariable(mesh=m, name='beer_lambert', value=numerix.exp(-alpha * m.x))

# Analytical expression for the differential equation above
analytical = CellVariable(mesh=m, name='analytical', value=-( B * numerix.exp( 2 * alpha * m.x ) + numerix.exp( 2 * alpha * D )  ) * C * numerix.exp( -alpha * m.x ) / ( A * B - numerix.exp( 2 * alpha * D ) ))

eq = eq1 & eq2
step = 0
while step < 10:
    print step
    res= eq.sweep()
    step += 1

vi = Viewer(vars=(beer_lambert, analytical, I_left + I_right))
#vi = Viewer(vars=(I_left + I_right))
vi.plot(filename='light_A=0,B=1.png')
raw_input()