# Continuous energy external source with constant probability for each energy bin

**URL:** <https://openmc.discourse.group/t/continuous-energy-external-source-with-constant-probability-for-each-energy-bin/4960>\
**Category:** User Support\
**Created:** [November 14, 2024, 1:10pm UTC](https://openmc.discourse.group/t/continuous-energy-external-source-with-constant-probability-for-each-energy-bin/4960 "2024-11-14T13:10:04Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jeff](https://avatars.discourse-cdn.com/v4/letter/j/e19b73/32.png) [@Jeff](https://openmc.discourse.group/u/Jeff)\
**Post date:** [November 14, 2024, 1:10pm UTC](https://openmc.discourse.group/t/continuous-energy-external-source-with-constant-probability-for-each-energy-bin/4960/1 "2024-11-14T13:10:04Z")

</div>

Hello, everyone. I am trying to simulate a fixed source problem with openmc. The external neutron source has a continuous energy which are separated into 47 groups. Each group has an identical neutron emission probability. According to the openmc manual, I chose `openmc.stats.Tabular` to set the energy distribution of the neutron source. Here is a part of codes for setting the external source:

```auto
# energy distribution
energy_bins =\
    [1.00000E-05, 1.00000E-01, 4.13990E-01, 8.76430E-01, 1.85540E+00, 5.04350E+00,
     1.06770E+01, 3.72660E+01, 1.01300E+02, 2.14450E+02, 4.54000E+02, 1.58460E+03,
     3.35460E+03, 7.10170E+03, 1.50340E+04, 2.18750E+04, 2.41760E+04, 2.60580E+04,
     3.18280E+04, 4.08680E+04, 6.73790E+04, 1.11090E+05, 1.83160E+05, 2.97210E+05,
     3.68830E+05, 4.97870E+05, 6.08100E+05, 7.42740E+05, 8.20850E+05, 1.00260E+06,
     1.35340E+06, 1.65300E+06, 1.92050E+06, 2.23130E+06, 2.34570E+06, 2.36530E+06,
     2.46600E+06, 2.72530E+06, 3.01190E+06, 3.67880E+06, 4.96590E+06, 6.06530E+06,
     7.40820E+06, 8.60710E+06, 1.00000E+07, 1.22140E+07, 1.41910E+07, 1.73320E+07]
# probabilities for each energy bin
energy_probs =\
    [1.44438E-11, 1.06857E-10, 2.52767E-10, 7.78361E-10, 4.01217E-09, 1.07794E-08,
     8.78665E-08, 3.61095E-07, 9.68270E-07, 2.98238E-06, 2.43405E-05, 5.97813E-05,
     1.83222E-04, 5.64378E-04, 6.27236E-04, 2.35381E-04, 2.00608E-04, 6.59333E-04,
     1.15149E-03, 4.06566E-03, 8.39881E-03, 1.71186E-02, 3.24986E-02, 2.27356E-02,
     4.34652E-02, 3.87843E-02, 4.80123E-02, 2.79514E-02, 6.39272E-02, 1.16353E-01,
     8.90702E-02, 7.06142E-02, 7.15504E-02, 2.36718E-02, 3.90518E-03, 1.93922E-02,
     4.53375E-02, 4.27965E-02, 7.51613E-02, 7.89060E-02, 2.92888E-02, 1.49788E-02,
     5.06871E-03, 2.17995E-03, 8.88027E-04, 1.39089E-04, 3.07600E-05]
   
# create a source for a fuel assembly
source = openmc.Source()
source.space = openmc.stats.Box((x_left, y_down, 0.0), (x_right, y_up, 190.0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Tabular(np.array(energy_bins), np.array(energy_probs), "histogram")                                      
source.strength = strength[count_assem]

```

However, I got an unreasonable result comparing to mcnp. I found that the external source sampled from openmc is different from that I input. Here are the figures of two neutron probability distributions:

 ![image](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/6/6c0f6cc238678278b3322d1f91cca48e53d4ca52.png)  
 ![input_spectrum](https://global.discourse-cdn.com/free1/uploads/openmc/original/2X/9/9be532921b68b37512dd7e0d3c89f8927c76110c.jpeg)

I tried different interpolation schemes of `openmc.stats.Tabular`, but I cannot fix it. Any suggestions? Thank you very much.

---

<div class="post-metadata">

**Author:** ![Jeff](https://avatars.discourse-cdn.com/v4/letter/j/e19b73/32.png) [@Jeff](https://openmc.discourse.group/u/Jeff)\
**Post date:** [November 15, 2024, 2:42am UTC](https://openmc.discourse.group/t/continuous-energy-external-source-with-constant-probability-for-each-energy-bin/4960/2 "2024-11-15T02:42:19Z")

</div>

Here are the codes to plot the energy probability density sampled from openmc:

```auto
# plot external source distribution
openmc.lib.init(output=False)
particles = openmc.lib.sample_external_source(n_samples=50000, prn_seed=1)
openmc.lib.finalize()

e_values = [particle.E for particle in particles]
e_array = np.array(e_values)

# count energy samples for each energy bin
counts, bin_edges = np.histogram(e_array, bins=energy_bins)

# plot probability density
plt.figure(figsize=(10, 6))
plt.bar(bin_edges[:-1], counts/50000, width=np.diff(bin_edges), edgecolor="black", align="edge")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Energy [eV]")
plt.ylabel("Probability")
plt.title("Energy spectrum from OpenMC samples")
plt.savefig("spectrum.jpg", dpi=300)

```

---

<div class="post-metadata">

**Author:** ![paulromano](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/paulromano/32/486_2.png) [@paulromano](https://openmc.discourse.group/u/paulromano)\
**Post date:** [November 16, 2024, 1:55pm UTC](https://openmc.discourse.group/t/continuous-energy-external-source-with-constant-probability-for-each-energy-bin/4960/3 "2024-11-16T13:55:37Z")

</div>

The thing to keep in mind is that the probabilities need to be given “per eV”, so if you haven’t already you should divide by the width of the energy bins. I suspect that is probably the issue here.

---

<div class="post-metadata">

**Author:** ![Jeff](https://avatars.discourse-cdn.com/v4/letter/j/e19b73/32.png) [@Jeff](https://openmc.discourse.group/u/Jeff)\
**Post date:** [November 20, 2024, 6:27am UTC](https://openmc.discourse.group/t/continuous-energy-external-source-with-constant-probability-for-each-energy-bin/4960/4 "2024-11-20T06:27:58Z")

</div>

Thanks for your reply. Dividing the probabbilities by energy bin width fixed this problem.