Beam with different divergence in X and Y

Dear FLUKA users,

I have used FLUKA for some time but never needed to create a beam with different divergence in X and Y as primary beam before. By default in FLUKA, you can give the Deltaphi that is then applied to both X and Y.

In order to make it different in X and Y, I use the source.f routine. As far as I understood I should add to FLAIR a card named source that sends the correct numbers to this routine.

The problem is that I did not understand where these values should be inserted.

Wish you can help me!

Best Regards,
Gian Luigi

Hello Gian Luigi,

If your parameters values are fixed, you can define them directly in your custom source.f.

However, if you need to change them on the fly (and want to avoid recompiling every time), you can indeed use the available WHATs in the SOURCE card to and access them from your source.f via WHASOU(1), WHASOU(2), etc.



Thank you for your answer Cesc!

I looked inside the source.f (that it’s the default one that FLUKA comes with) but I did not manage to find the variables linked to divergence. I have only found: momentum, cosines, etc…

Hi Gian Luigi,

You can either

  • hard-code the divergence along X and Y in your source.f
    by adding lines saying DIVX=blabla and DIVY=bublu,
  • or use e.g. the first two WHATs in the SOURCE card and then say
    DIVX=WHASOU(1) and DIVY=WHASOU(2) in your source.f.

… both followed by whichever logics you need to implement. From your
sentence “I use the source.f routine” I understand this is the process you
are in now. In the end you populate the direction cosines in the elements

If your problem is how to actually sample with different divergence along
X and Y, say so.


Hi Cesc,

Thank you again for your suggestions.

I tried to input the divergence in the source.f file but this did not generate a change in the divergence of the primary beam. So maybe I cannot reach what I need modifying source.f.

What I need is a primary beam with a rectangular shape in X and Y and with 11 mrad divergence in X and 0 mrad in Y.

Gian Luigi

Dear Gian Luigi,

your beam definitely can be created with the source.f user routine.

But bear in mind, the source routine does not do the divergence sampling by default. You need to code it, and feed the sampled direction cosines to the TXFLK(NPFLKA), TYFLK(NPFLKA), and TZFLK(NPFLKA) variables as Cesc mentioned.

Moreover, the sampling for the rectangular shape has to be implemented as well. The starting point of the beam particles has to be stored in XFLK (NPFLKA), YFLK (NPFLKA), XFLK (NPFLKA).

If you still need help, please attach the source routine file as well, so we could take a look.


Dear Gian Luigi,

What I need is a primary beam with a rectangular shape in X and Y and with 11 mrad divergence in X and 0 mrad in Y.

Now, THAT should have been the opening statement in the original question.

I understand you are shooting particles starting along the z axis with a certain polar angle theta from 0 to 0.5*DIVX, and azimuthal angle either 0 or pi, such the particle’s initial direction is in the XZ plane.

What you need to do (adapt/optimize as needed) is sample COS(THETA) homogeneously between 1 and 1-COS(0.5D0*DIVX). The sampling is analytical (double check in case I slipped here) and should give you something along the following ~FORTRAN lines:

XI  = FLRNDM(XI)                        ! Homogeneously sampled number in [0,1)
CTH = 1.0D0 - XI * (1.0D0 - COS(0.5D0*DIVX))   ! Sampling cos(theta), note DIVX must be in radians.
STH = SQRT((1.0D0+CTH) * (1.0D0-CTH))   ! sin(theta)

XI2 = FLRNDM(XI2)       ! One more homogeneously sampled number in [0,1)
IF (XI2 .LT. 0.5D0) THEN
  TXFLK(NPFLKA) = STH   ! 50% of the times along +x
  TXFLK(NPFLKA) = -STH  ! 50% of the times along -x

Regarding the rectangular (spatial) shape, you can define XMIN,XMAX,YMIN,YMAX in your source.f and then do

XI3 = FLRNDM(XI3)  ! One more homogeneously sampled number in [0,1)
XI4 = FLRNDM(XI4)  ! One more homogeneously sampled number in [0,1)

Adapt/correct as needed, and test profusely.



PS: I multiply all divergences by 0.5 to remain in line with the way they are used everywhere else in FLUKA.

Dear David,

Thank you for pointing that out! You are right indeed.

Gian Luigi

Dear Cesc,

Sorry for my phrasing, it’s the first time I send a message to the forum and I see that was not particularly clear.

Thank you very much for you help, especially for writing down some lines of code that are very useful!

I am trying to implement your suggestions now! Thank you again!

Gian Luigi