Is there a package or algorithm you could recommend for calculating numerical gradients of a function that takes the form:

f(pos: Vector3<f64>) -> f64

Strictly speaking, I need it to get the intensity gradient (E-field-squared-gradient) for cases of many-beam interference, where an analytical solution is hard to obtain (Although I will think about it some more... for some special cases (like linear polarized Gaussian beams only).

I've seen https://docs.rs/GSL/4.0.0/rgsl/numerical_differentiation/index.html but I don't know how to "install" it.

Obviously, I could just write a 2,3 or 5 -point algorithm myself, and implement it in maths - but I was wondering if there existed an optimized algorithm for that already.

Hi,

is it a known problem that Anti-Virus software (such as Kaspersky Internet Security for instance) think that the

AtomECS-master\target\rls\debug\build\crossbeam-epoch-9c14771e7fefb252\build-script-build.exe

is some sort of malware? It is a bit annoying since it always forces me to "desinfect" (aka delete) files from the project directory. Moreover, because the Anti-Virus software "keeps" the file occupied, cargo run --release fails due to

error: failed to run custom build command for crossbeam-epoch v0.8.2

Caused by:
could not execute process D:...\AtomECS-master\target\release\build\crossbeam-epoch-576fdc92c8f92a03\build-script-build (never executed)

Caused by:
Access is denied. (os error 5)
warning: build failed, waiting for other jobs to finish...
error: build failed

The software gives me the following information on what it thinks about the "malware":

https://threats.kaspersky.com/en/threat/Trojan.Win32.Agent/

Best wishes and thank you!

Maurice

I was wondering if the far detuned dipole beams of a FORT / ODT lead to an A.C. Stark shift of the optical transitions used for a MOT. For the example of a 1064nm dipole trap that is briefly superimposed with a red Sr-MOT, I tried to calculate the light shift according to 7.7 Foot (see PDF attached). As you can see, the result I get is pretty much zero but I'm confused by the sentence "Normally light shifts are most important at large frequency detuning [...]" in 7.7 Foot.

What do you think - is my rough estimation of the magnitude of the light shift correct in principle or have I missed something crucial in my approach? Does someone have practical experience with MOT-->ODT transitions and knows if that effect becomes relevant in real life?

Why I am needing this: I want to simulate MOT --> ODT transitions with AtomECS and I want to know if an a.c. Stark shift is something we'd have to implement for this to work correctly.
Stark-shift_ODT.pdf

bevy has a really nice implementation of 'Plugins', where you can define individual parts of a program to add/remove. A similar plugin model could be used in AtomECS to add different models of laser cooling, different forces, collisions, etc. This would improve the usability and increase program flexibility (possibly also performance by only running required/declared modules).

Related to #3.

To simulate evaporation in optical/magnetic traps, AtomECS requires a way to calculate the effects of short-ranged (eg, s-wave) collisions.
Given that the collisions are short-ranged, all methods used should a spatial binning method for speed.
Beyond that, there are subtle variations in the implementation details, and we could support a number of approaches:

1. Monte-Carlo method: determine a number of collisions to occur within a box, eg based on density and mean velocity. Perform collisions between velocities of random atomic pairs.

   • How to do mixtures with this, eg fermions+bosons? Sr87+88?
   • Will be faster, because it's one loop per box with comparatively few collisions. Scales as O(N), where N is average atom number per box.

2. 'Pair-wise' Monte-Carlo method; Iterate over each pair within a box. Calculate the chance for a collision to occur, then draw a random number to determine if it did.

   • Easier to do mixtures, because you can specify collision lengths between pairs_
   • Will be slower, and scales O(N^2) where N is average atom number per box.

All models should also have the ability to specify some kind of 'macro atom', so that the behaviour of a large number of atoms can be represented by a numerically tractable ensemble. There's a few ways we could do this:

1. Have a simulation resource which specifies a global scaling parameter.
2. Specify per-atom MacroAtom components.

I like the global scaling parameter because it is straight forward. Also, I can see how we might be able to simulate multiple orders of magnitude of atom loss by rescaling the parameter as atoms are lost. For instance, this could allow one to start with 1e4 atoms representing 1e10 atoms and end with 1e4 atoms representing 1e4 atoms.

The QuadrupoleField3D and QuadrupoleField2D have their position associated with an external Position component, while PrecalculatedMagneticFieldGrid has its position defined internally. I was wondering if there is any reason for such inconsistency?

I guess this could be extended to laser beam as well because GaussianBeam also has its "position" (intersection) defined internally.

Hi @MauriceZeuner and @d-garrick ,

I added a WIP implementation of ramps using a generic system here: https://github.com/TeamAtomECS/AtomECS/blob/collisions/src/ramp.rs

I intend to define a macro, rampable, so you can auto derive the required trait impl.

It looks like the scattering force abruptly cuts off at weaker scattering. Running plot_molasses.m with more points illustrates this well:

You can see that for velocities above ~22 the force abruptly cuts to zero.

I expect this is related to the Poisson implementation when the mean is low, but will investigate more.

I have discovered that the length of the integration time-step has a very significant effect on the "numerical" heating in a dipole trap. I attached a picture, showing the simulation results of 3 (otherwise identical) simulations with:

• blue: 0.1 µs
• yellow: 1 µs
• red: 2 µs (only a few atoms left ^^)

as the delta - timestep. They started with 100 atoms and ran for a total real time of 0.2 seconds each.
As you can see from the picture, for ODT one should choose a smaller timestep than in a MOT setup. I just felt this is something to be aware of. I'll try to find the largest time-step with negligible heating over a couple of seconds. Any comments?

I think it would be great to have a very simple way to write not only Position and Velocity to a textfile and analyse it in Python/Matlab but maybe a general WriteComponentToFileSystem that users can configure to write the values they are interested in.

For example: If I wanted to print the contents of the TotalPhotonsScattered component for all atoms over time, I might want to plot that against the power of all GaussianBeam entities. Currently (please correct me if I'm wrong) I write my own writer-system for that every time I need it and there exists no nice and most importantly run-time optimized way of doing that. (See the super-slow xyz-writer, ahem, ahem... ).

But especially if one wants to monitor all sorts of interesting parameters over complicated runs (red Sr-MOT to dipole transition, for example), it would be so cool to get that data printed to a text file in one line of code.

What do you think?

Following a good discussion with Tiffany:

Especially for dipole trapping, it can actually become problematic to pretend the Atoms have only one transition that is relevant. For example, what if you have a dipole beam for Strontium, close to the red transition? You'd have to respect the blue because it has a huge linewidth but the red one also, because you are close.

A far more realistic treatment would be to use polarizabilities instead and also E-fields instead of intensity gradients. I'd propose this as an enhancement for the foreseeable future. What do you think?

There are quite a few branches left over, it needs a tidy up.

@MauriceZeuner is working on a dipole force trap implementation

The GaussianReferenceFrame already exists (and has an "ellipticity" attribute) in the new-dipole-force branch but elliptical Gaussian beams are yet to be implemented. Could be useful not only for dipole traps.

I need to update to the latest specs+rayon to fix the panic issue, so we can get back to nightly toolchain.

The specs API has some breaking changes, so I need to update all examples aswell.

1d MOT is currently broken - it's probably related to the fact that ScatteringFluctuationsOption::Off can give NaNs

There are two problems though the converted power in the current implementation is correct.

First question

Shouldn't the right hand side be multiplying e_radius by sqrt(2)?

Shouldn't the right hand side be like divide the product of peak_intensity, pi e^2_radius by 2 instead of multiplying it by 2?

Given that

(Edmund's optics, RP Photonics)

Rearranging would give

Second question

According to the documentation

making the changes suggested in Q1 would give the same result because following the suggestion above would give

And the current implementation is

which gives the same result.

However, this leads to the problem of a small optimization where the square root and the multiply by two can be simply neglected.

Sorry for spamming so many issues here - let me know when to stop. ^^

Can you debug AtomECS in VS Code? It seems it never stops at my breakpoints in systems... am I getting something wrong? Hast it something to do with the strong parallelism? It would be great to stop at points in the simulation have a look at the entries in all components, is that possible?

Thanks in advance!

Just to let people know that I am working on extending the Gaussian beam description to correct treatment of (at least fully) polarized light with arbitrary ellipticity of the Jones vector.

If you wonder why: Well, I think we do already have 50% of the stuff we need to model Sisyphus-cooling and I suppose that this would definitely be a nice feature in the mid-future. Any thoughts?

The branch will be "alternative-dipole-force" where I also switch to a full complex electric field treatment for dipole traps.

At the moment, I imagine it working like this: Your basic component always is the LinearGaussianEBeam and then you just add a (name provisional) JonesVector component to the beam entity. Since you can separate the rotation of the electric field vector and the other gaussian beam characteristics, you can then in a SampleElectricFieldSystem just (complex) multiply the jones vector (which is basically a rotation-operator) to the other terms of the gaussian E-field and you should get a (for example circularly) polarized gaussian beam.

For Sysiphus cooling I think, you should be able to get a standing wave in polarization (not intensity!) by just adding two counterpropagating beams with orthogonal polarization to the simulation. At positions where the counter-propagating beams have a phase difference of π/2 the polarization is circular, and where there is no phase difference, the polarization is linear. In the intermediate regions, there is a gradient elipticity of the superposed fields.

What, of course would also be required for this to work (will be a separate issue, probably) is a slightly more complicated atom-component which can respect the important M_J substrates of the atom.

We have to redocument the newer classes, as sometimes doc coverage has been lost due to changes (all in the name of progress, which is good :) )

Make sure the final code is consistent with cargo doc --open.

Hey @MauriceZeuner ,

Could you check and fix this unit test? It's been broken since I fixed the scattering rates, so I'm confident it's the test is outdated rather than the implementation (for instance, the integration test for scattering rate is working fine and agrees with the analytic formulae).

Cheers!

On most recent toolchain, running AtomECS gives a long error which ends in:
thread panicked while panicking. aborting.

As a temporary fix, use the stable toolchain from 2020-06-17. Type the following into the console:

rustup install stable-2020-07-16
rustup default stable-2020-07-16

Hey Maurice,

We now have ApplyEmissionForceOption which enables to use of random emission kicks required for recoil limit and part of the Doppler limit.

I say part of Doppler T, because even with this option disabled the simulation will still draw random numbers from the mean dists, and so there will be fluctuations in the absorption forces. That contributes part of the Doppler T.

Suggested course of action, we change ApplyEmissionForceOption to something else (eg, EnableScatteringFluctuations). When this resource is present, the systems which transfer 'mean numbers' to 'actual numbers' use Poisson dist; when it isnt, they just transfer the numbers exactly (so the execution is faster).

Add support for magnetic traps.

For the purposes of this issue, this refers to traps where the potential energy can be written in the form E = g_F \mu_B | B(r) |, and where B is a real magnetic field. This includes static magnetic field traps, such as quadrupole traps, and traps where the magnetic field is slowly varying, e.g. Time-Orbiting Potential (TOP) traps.

This issue does not include support for RF-dressed traps (where the magnetic field is oscillating at a rate comparable to the Zeeman splitting).

I just realized that the current implementation of the ramps does not allow for ramping two components of the same entity at the same time. This is difficult for me since I would like to replace my dipole::transition_switcher::RampMOTBeamsSystem with them to ramp detuning and power of MOT beams. I believe there is no simple way to allow for two Ramps to be associated with an entity, right?

One solution would be to make the wavelength an attribute of GaussianBeam, so only one component needs to be ramped or we need a different approach to ramping in general. Or I leave it roughly like it is currently. (Custom systems to ramp them)

Also, I find it a little unhandy that I have to re-initiate the entire Component for every keyframe, even if I only want to change one attribute of it. This happens, if you, for example, want to ramp the power of a GaussianBeam component. Just takes a lot of code, I think.

Any ideas?

Hi,

@tiffanyharte and myself would like to know what other people are using AtomECS for. For planning future features it would be very helpful to know what kinds of experiments people are interested in simulating - so please, anyone, feel welcome to comment on this issue :)

I'll add information from the Oxford setups below and tag others who have been using it.

There is a known issue related to specs, and it is caused by one of its dependencies, hibitset. This restricts the maximum number of entities to usize^4 (ie. 32^4 for 32-bit systems or 64^4 for 64-bit systems). I don't think there really is an easy solution other than writing a PR for hibitset. However, I am just wondering if you happened to have some workaround in the occasions where there is a large number of atoms? Plus, is there some kind of plan to overcome this in the future?

We are unsure if the total scattering rate in photons_scattered.rs should be calculated with the linewidth in Hertz or angular units. The same question has been raised for the two-level excited state population in twolevel.rs where the Einstein A coefficient is used.

The corresponding section in C. Foot: 7.5.2

Something Tiffany and I would be very interested to see some results of in the mid-future.

Hi all,

sorry for being quiet - I have been desperately trying to find the cause for a heating effect which seems to be still in the code but is most likely a dipole trap specific feature. I apologize in advance, this is super nasty and I'd be glad for new ideas how this could happen.

First observation: If you simply put atoms in a dipole trap, but also have some CoolingLight Entities in the game, however, very far detuned, so they shouldn't do anything - you'll notice that there appears to be a preferred direction for their escape - the x-direction. It can best be observed in the (branch: new-dipole-force) cross_beam_dipole_trap.rs example - where there should be no difference between x and y since we have two identical, horizontal beams in exactly these two directions (and the MOT beams are on the diagonals). However, if you plot the absolute coordinate values of the ensemble (or just have a very sharp eye when plotting the data in 3D in Blender) you'll see that they are more likely to be further away from the centre in (plus or minus) x-direction.

In addition to that, there is still a good amount of numerical heating happening on top of that, this effect, however can be reduced by increasing the detuning of the MOT beams even further:

Before you ask: I checked in a pure symmetrical 3D MOT if there was a difference between x and y and I don't think it is present there. Of course, I also checked the way the dipole force is calculated for any symmetry breaking but couldn't find one. At that point, however, I am running a bit out of ideas where to look for it.

For comparison, here are the same diagrams as above just without any CoolingLight entities in the code:

How did I find this behaviour in the first place? Well, I was running MOT -> DipoleTrap transitions with three stages:

1. MOT only (first 0.2 s)
2. MOT ramp down (power and detuning) (another 0.2s)
3. Dipole only - CoolingLight entities are deleted - (rest 1.4 s)

(see file: example/red_mod_xodt_transition.rs)

which looks like this:

and

I'd be extremely grateful for any, even far-fetched, ideas what could cause this behaviour.

I actually just wanted to see if an error that I encountered on the dipole branch also happens to be in the master branch (good news: it's not) but I noticed that when I ran cargo test, two tests fail, namely:

magnetic::tests::test_field_samplers_are_added
magnetic::tests::test_magnetics_systems

I am happy to look for a solution if no one else feels responsible or has time but I thought it would be nice to make you aware of these fails first.

Full terminal output is attached and also a screenshot of the relevant section.
magnetic_test_fail_full_output.txt

Long-ranged forces are important for 3D MOTs with lots of atoms, as the repulsive forces reduce the phase space density.

Some thoughts:

• May need to use a MacroAtom approach, where relevant quantities like repulsive forces are scaled so that a small number of atoms can simulate the PSD of an (intractable) larger number of atoms.
• Force goes as 1/r^2 which is long ranged and a real pain to parallelise.

Possible implementation:

• Hash spatial position to sort atoms into bins.
• Use a tree structure to allow fine-grained forces at close range, and coarse-grained forces for long range.
  • Atoms are sorted into each leaf of the tree.
  • Each leaf in the tree calculates total 'charge' and 'centre-of-mass'.
  • The tree structure is traversed to aggregate charges into each branch.
  • Tree structure from this point is immutable.
  • Calculate force exerted on each leaf due to tree structure (parallel over each leaf). Only calculate leaves which are occupied (total charge > 0). Iterate branches and leaves, to depth based on distance from calculated leaf.
  • For each atom, determine force by looking up leaf node. Use leaf.force + force calculated pair-wise between other atoms in the leaf.

Building tree will be O(L) O(N).

Hi!
I finally started properly fixing the first of the numerical heating problems. For this, I want to replace the faulty maths::random_direction with rand_distr::UnitSphere since Elliot recommended that

Hmm, I didn't realise we were defining our own function for it! I assumed we were using one from another crate. Can you replace ours with https://docs.rs/rand_distr/0.4.0/rand_distr/struct.UnitSphere.html ?

see #27

So, I was trying that, however, I encountered a problem where the compiler complains about the argument of UnitSphere.sample in force.rs

                                    use rand_distr::{Distribution, UnitSphere};
                                    // explicit random walk implementation
                                    for _i in 0..total {
                                        let v: [f64; 3] = UnitSphere.sample(&mut rng);
                                        force.force = force.force
                                            + force_one_kick * Vector3::new(v[0], v[1], v[2]);
                                    }

saying:

error[E0277]: the trait bound `atom_sources::rand::prelude::ThreadRng: rand_core::RngCore` is not satisfied
   --> src\laser_cooling\force.rs:174:77
    |
174 | ...                   let v: [f64; 3] = UnitSphere.sample(&mut rng);
    |                                                           ^^^^^^^^ the trait `rand_core::RngCore` is not implemented for `atom_sources::rand::prelude::ThreadRng`
    |
    = note: required because of the requirements on the impl of `rand::rng::Rng` for `atom_sources::rand::prelude::ThreadRng`

which makes no sense to me since I don't see any rand being imported from atom_sources nor is it existing there at all. I don't quite know where to start, to be honest, so I'd be glad about suggestions. ^^

Thank you!

I would like to use the uom crate to define quantities in atomecs.
https://github.com/iliekturtles/uom

Advantages:

• Compile-time dimensional analysis
• No runtime overhead
• Convenience to define in whatever unit is appropriate at the time
• No incorrect conversions

I had tried previously but was not able to because of the incompatibility of Vector3<> from nalgebra with uom. There have been a few updates to the package since then and it is worth revisiting.

Performance is slightly better with Legion than specs (Amethyst has now dropped specs and is moving to Legion). Consider moving to Legion.
https://github.com/amethyst/legion

Would anyone object if we added the num = "0.4.0" to the list of project dependencies? It enables the use of complex numbers with:

num::complex

See: https://docs.rs/num-complex/0.4.0/num_complex/index.html

Euler integration is currently used - it's straightforward but has terrible integration error. I'll change it so that AtomECS uses Velocity-Verlet by default. VV has much smaller integration error. This is particularly important for future thermodynamic simulations of atoms in traps, as integration error leads to heating similar to parametric excitation.

I'll admit I thought we were already using the Velocity-Verlet integrator - I confused atomecs with the other ECS code (Unity/C#). My other code uses the follwing implementation:

/// <summary>
/// Update velocity according to velocity verlet.
/// </summary>
[UpdateAfter(typeof(ForceCalculationSystems))]
[UpdateInGroup(typeof(FixedUpdateGroup))]
public class UpdateVelocitySystem : JobComponentSystem
{
    protected override JobHandle OnUpdate(JobHandle inputDependencies)
    {
        float DeltaTime = FixedUpdateGroup.FIXED_TIME_DELTA;
        return Entities.ForEach(
            (
                ref Velocity velocity,
                in Mass mass, 
                in Force force, 
                in PrevForce oldForce
                ) =>
                 velocity.Value = velocity.Value + 0.5f * (force.Value + oldForce.Value) * DeltaTime / mass.Value
            ).Schedule(inputDependencies);
    }
}

/// <summary>
/// Update position according to velocity verlet.
/// </summary>
[UpdateBefore(typeof(ForceCalculationSystems))]
[UpdateInGroup(typeof(FixedUpdateGroup))]
public class UpdatePositionSystem : JobComponentSystem
{
    protected override void OnCreateManager()
    {
        Enabled = false;
    }
    protected override JobHandle OnUpdate(JobHandle inputDependencies)
    {
        float DeltaTime = FixedUpdateGroup.FIXED_TIME_DELTA;
        return Entities.ForEach(
            (
                ref Translation translation,
                in Mass mass,
                in PrevForce force,
                in Velocity velocity
                ) =>
                translation.Value = translation.Value + velocity.Value * DeltaTime + 0.5f * (force.Value / mass.Value) * DeltaTime * DeltaTime
            ).Schedule(inputDependencies);
    }
}

I'll add this to AtomECS when I get time, it is a straightforward addition and worth it.

Hi,

I've been using AtomECS a bit lately and I just had a question on how it is possible to handle linear polarization. In the doc it's mentionned that it should be possible to have a linear polarization by combining sigma_+/- beams, but I don't understand how that can be.

For example I'd like to consider the following situation: a magnetic field along e_z and a cooling light with wavevector along e_x linearly polarized along e_y. This light can induce sigma_+/- transition, but not pi.

However, when I look at the calculation for CalculateRateCoefficientsSystem, the term scatter3 for the pi component will always be non zero whatever the configuration I set for the polarization of the cooling lights.

What would be the correct way to handle this situation?

Thanks

Starting to profile the latest changes to red-mot to get a head start on optimisations.

Interestingly, it seems the vast majority of program time was being taken up by rayon implementations - specifically SwitchToThread and SleepConditionVariableSRW. There was also poor spin time.

This seems likely related to a known issue in rayon, where sleeping threads would burn a lot of cycles/power. For example, see
rayon-rs/rayon#642
rayon-rs/rayon#795
rayon-rs/rfcs#5
It looks like some of these issues were fixed in rayon 1.4.0, although we are currently on rayon 1.3.0. I'll try updating rayon (new version is 1.5.0) and see if it improves things.

I was trying to increase the number of cooling lasers to replicate sidebands.
Every time I add more than 16 lasers to the world, a thread panicks.

Any thoughts on what could be causing the problem ?

Hi here,

I am new to here, and I want to fix compiling and clippy first, if there is any other issues, I may take one.

Hey @d-garrick,

I noticed your recent merge incorrectly added some matlab code to the repository under the +util folder. Can you please delete it? That code should be added to https://github.com/TeamAtomECS/matlab instead.

What is the correct way to generate atoms from an oven with multiple transitions ?

I assumed that I could just add an AtomicTransition component to the oven creation, but it doesn't seem to work.

For example, replacing in 2d_plus_mot_from_oven the line

l156: .with(AtomicTransition::strontium())

either by

l156: .with(AtomicTransition::strontium())
l157: .with(AtomicTransition::strontium_red())

or

l157: .with(AtomicTransition::strontium_red())
l156: .with(AtomicTransition::strontium())

doesn't give the same simulation results in the two cases.

It seems that in the first case, adding the red transition overwrites the blue.

Are you considering publishing this on crates.io? There have been some "bad actors" who claims the crate's name by just publishing an empty crate on crates.io, and you should probably consider claiming this name to avoid any potential issues in the future?

The current ApplyEmissionForceSystem uses a normal distribution when the number of scattered photons N > 5. I don't think this respect recoil limit for reasons we have discussed beforehand.

The code runs faster when you approximate N kicks with a normal distribution, and the recoil limit isn't applicable for some systems. Let's add an optional flag to either always perform individual kicks, or use the distribution. This could be an enum in the ApplyEmissionForceOption struct.

I would like to see a high level API for creating entities in the simulation, using something like factory structs with a fluent interface. e.g.:

BeamBuilder()
   .in_direction(Vector3::x)
   .with_gaussian_profile(waist, power)
   .with_sigma_plus_polarisation()
   .for_cooling(Strontium)
   .with_detuning(-32.0)
   .build(&mut world);

The advantages of this are:

• user doesn't need to know components required.
• guarantees safe use/initialization of components.
• prosaic code.
• reduces redundancy for duplicate variables (e.g. wavevector)
• nicer editor integration, e.g. autocomplete/prompting. Helpful for new atomecs users.

Disadvantages:

• more code to maintain.
• multiple ways to define things (though this can be alleviated by funnelling users towards higher-level functionality in an atomecs::prelude::* mod.

Originally posted by @ElliotB256 in #29 (comment)

Following the nasty heating error in the dipole trap, I'd propose we should think about how to generally clean component data (in each iteration step). Two things come to mind:

1. We could just clear ALL dynamical components every step, like we do with Force. Pro: clean, Con: additional code that runs every step
2. We could make a more dedicated "EntityRemovalSystem" that not only removes (for example MOT beam entities) but also makes sure that corresponding entries in all the vectors are set to 0 or similar. Pro: minimal runtime impact, Con: More prone to errors

I believe this topic becomes more and more important as we move towards simulation files that combine various side-arm stages.

What do you think?

There are currently a number of ways of configuring simulation options which should be standardised in format. For instance, we currently have EnableScatteringFluctuations and ApplyEmissionForceOption, which are dissimilar naming conventions and should be made consistent.

The options also currently favor speed over accuracy. I propose we change it so that accuracy is favored over speed; users disable functionality to improve simulation speed when required, but they only compromise physics when they choose to.

Makes much more sense, this is already half done for the Matlab scripts (some of which are being used to analyse examples for the accompanying paper).

The code seems to assume a Gaussian beam profile, I would like to simulate a MOT geometry using laser beams with a tophat profile or more generally be able to modify the beam profile to be an abritrary function. Adding a tophat module in lasers is not sufficient as the code continues to assume a gaussian beam in other instances.