Add particle BC: angle_threshold (thermalize vs reflective by incidence angle)

[PLASMAsDJ]

Summary
This PR adds a new particle boundary condition, angle_threshold, which switches between thermalize and reflective depending on the particle incidence angle relative to the outward normal of the boundary.

Motivation / reference
This implementation is inspired by the boundary-treatment approach described in:
[Chen, H., et al, Geophysical research letters, 52, e2025GL116478 (2025)]

Definition
Let n be the outward normal vector of the boundary. The particle velocity is decomposed into components parallel and perpendicular to n, and the incidence angle is defined as:

α = tan⁻¹(v⊥ / v∥)

If α ≤ α₀ (user-defined threshold), the particle is handled with thermalize; otherwise, reflective is applied.

Implementation
In the updated implementation:

• the threshold angle is parsed once during initialization,
• the parsed value is stored as a species parameter,
• and the runtime boundary-condition routine reuses the stored value instead of parsing the boundary-condition string repeatedly.

Modified files

• src/ParticleBC/BoundaryConditionType.cpp
• src/ParticleBC/BoundaryConditionType.h
• src/ParticleBC/PartBoundCond.cpp
• src/Python/pyinit.py
• src/Species/Species.h
• src/Species/SpeciesFactory.h

Tests
I rechecked the basic behavior with 1D and 2D test runs after the refactor.
If needed, I can provide a minimal namelist and a small test to reproduce the results.

Example namelists:

• angle_threshold_Test_1D.py
  angle_threshold_Test_1D.py

• angle_threshold_Test_2D.py
  angle_threshold_Test_2D.py

Note
This PR is mainly intended to get feedback on the implementation approach and code structure.

Namelist usage
1D:
boundary_conditions = [
["angle_threshold", "angle_threshold"],
]
threshold_angle = [
[a_xmin, a_xmax],
]

2D (x and y both use angle_threshold):
boundary_conditions = [
["angle_threshold", "angle_threshold"],
["angle_threshold", "angle_threshold"],
]
threshold_angle = [
[a_xmin, a_xmax],
[a_ymin, a_ymax],
]

2D (only y uses angle_threshold):
boundary_conditions = [
["reflective", "reflective"],
["angle_threshold", "angle_threshold"],
]
threshold_angle = [
[a_ymin, a_ymax],
]

[PLASMAsDJ]

This pull request is opened primarily for review/evaluation and to provide an easy-to-read diff of the proposed changes.
There is no expectation that it must be merged; feedback is very welcome.

[mccoys]

I didn't check fully your code but there are some comments

[mccoys]

Inline will likely not be used here

[mccoys]

Same

[mccoys]

You can use M_PI, i believe, for pi. To be checked

[mccoys]

I suggest not giving the choice of degrees. This makes the input too difficult to benchmark when there are too many possibilities. Users can easily switch from the python script

[mccoys]

Again I am not sure you can defined inline functions here. If needed on gpu, this might cost something

[mccoys]

These tests should not be done here. Parse the string art the initialization, otherwise this will be to slow

[PLASMAsDJ]

@mccoys Thanks for the comment. I was able to address the other comments rather straightforwardly, but this one seemed a bit more involved, so I just wanted to check if I understood your suggestion correctly.

Would you recommend parsing the boundary-condition string only once during initialization, then storing the parsed angle-threshold value as a parameter to be reused at runtime, instead of parsing it inside BoundaryConditionType.cpp?

If that is what you mean, I was thinking it might require a small refactor on the initialization side, for example by storing the value in Species and setting it from SpeciesFactory. Would that be the right direction?

[mccoys]

Yes that sounds exactly right

[PLASMAsDJ]

@mccoys Thanks again for the suggestion. I updated the implementation accordingly, and there have been several revisions to the code. I would be very grateful if you could review it again when you have time.

I agree that optimization is important. At the same time, my main focus here is to implement the intended particle boundary-condition mechanism, where the particle behavior at the boundary is determined by the incidence angle with respect to the boundary normal. I would also appreciate your feedback on whether this mechanism is represented appropriately in the current implementation.

The threshold angle is now parsed once during initialization, stored as a species parameter, and reused at runtime instead of being parsed inside BoundaryConditionType.cpp.

I also updated the PR description and rechecked the behavior with 1D and 2D test runs.