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
ELSE
TXFLK(NPFLKA) = -STH ! 50% of the times along -x
END IF
TYFLK(NPFLKA) = ZERZER ! Zero in FLUKA
TZFLK(NPFLKA) = CTH
```

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)
XFLK(NPFLKA) = XMIN + XI3 * (XMAX - XMIN)
YFLK(NPFLKA) = YMIN + XI4 * (YMAX - YMIN)
```

Adapt/correct as needed, and test profusely.

Cheers,

Cesc

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