photon_walksΒΆ

../_images/photon_walks.png
#!/usr/bin/env python-mpacts
from mpacts.commands.force.body import DirectedForceCommand
from mpacts.commands.geometry.orientation import RandomWalkInPlaneCommand
from mpacts.commands.monitors.progress import ProgressIndicator
from mpacts.commands.monitors.checkpredicate import StopWhenPcEmptyCommand
from mpacts.commands.monitors.progress import DefaultProgressPrinter
from mpacts.commands.monitors.progress import PrinterChain
from mpacts.commands.monitors.progress import ManagerAlivePrinter
from mpacts.commands.onparticles.create_remove import RemoveParticlesCommand
from mpacts.core.arrays import create_array
from mpacts.core.simulation import Simulation
from mpacts.particles.particlecontainer import ParticleContainer
from mpacts.predicates.predicates import Predicate_Random
from mpacts.predicates.predicates import Predicate_InBox
import mpacts.tools.random_seed as rand
import mpacts.particles as prt
import mpacts.geometrygenerators.quadgeometries as quadgen
from mpacts.predicates.shortened_predicates import *
import mpacts.commands.monitors.progress as pp


N = 20000
r = 0.005  # Only for visualization
gamma = 1.  # So force = velocity
v = 0.1
scatter_probability = 0.5
box_size = (0.3, 0.3, 0.3)
box_max = [0.5 * el for el in box_size]
box_min = [-0.5 * el for el in box_size]
init_position = (0, box_max[1], 0)
output_interval = 0.1
absorption_probability = 0.1

mysim = Simulation("simulation", timestep=1e-3)

rand.set_random_seed()  # Comment this out if you want the same results every time!

# photons:
xs = N * [init_position]
photons = prt.ParticleContainer("photons", prt.Sphere, mysim)
direction = create_array("Vector", "direction", photons)
photons.add_and_set_particles(x=xs, r=r, m=0., direction=(0, -1, 0))
photons.TimeIntegration_ForwardEuler_UncoupledOverdamped_Cmd(gamma=gamma)
inboxpred = Predicate_InBox( "Predicate_in_box",
                           array=photons['x'],
                           min=tuple(box_min),
                           max=tuple(box_max))
ranpred = Predicate_Random( "Predictate_Random",
                          p=absorption_probability,
                          sim=mysim )

orcmd = RandomWalkInPlaneCommand( "randomwalk",
                                mysim,
                                pc=photons,
                                p=scatter_probability,
                                direction=direction,
                                predicate=inboxpred )

df = DirectedForceCommand( "directedforce",
                         mysim,
                         pc=photons,
                         magnitude=v,
                         direction=direction )

wrt = photons.VTKWriterCmd(executeInterval=output_interval)
wrt.select_all(True)

RemoveParticlesCommand("Remove_particles_box",
                      mysim('loop_cmds/advance_time_cmds'),
                      pc=photons,
                      predicate=OR(NOT(inboxpred),ranpred))

# box:
(points, quads) = quadgen.create_box(box_size[0], box_size[1], box_size[2])

rigid = prt.ParticleContainer(
                             "box", prt.RigidBody.compose(
                                                         prt.RigidQuad), mysim)
rigid.VTKWriterCmd(executeOnce=True)
box = rigid.add_particle()
box.controlPoints.add_and_set_particles(xBF=points, x=points)
box.RigidQuads.add_and_set_particles(vertexIndices=quads)
box.x = (0, 0, 0)
box.q = (1, 0, 0, 0)

printer_list = [ pp.DefaultProgressPrinter(mysim),
                 pp.ManagerAlivePrinter(photons) ]

printer = PrinterChain(printer_list)

pp.ProgressIndicator("PrintProgress", mysim, printer, print_interval=4)

StopWhenPcEmptyCommand("DoneCommand", mysim, pc=photons)

mysim.run_until_done()