Dear Ihor Melnyk,
I’d mainly like to comment on one aspect of your question, the implementation of the angular distribution.
Due to the correlation between the three direction cosines (in order to define a direction in three dimensions, you only need two variables, for example the polar and azimuthal angle), the approach you are using at the moment will not give the expected angular distribution.
In the second step, this code almost certainly samples a value for direction_cosz that is incompatible with direction_cosx and direction_cosy. FLUKA re-normalizes the improper direction vector and produces primaries, but their angular distribution is somewhat in between the intended distribution and an isotropic source.
As an extreme example, consider the following modification of your source routine:
call sample_isotropic_direction( direction_cosx, direction_cosy, direction_cosz )
direction_cosz = 1.0D0
One might expect that all particles move along the z axis. However, direction_cosx and direction_cosymay be assigned values with an absolute value of up to 1.0 by the isotropic sampler, which clearly conflict with the second instruction. This may produce direction vectors like (1.0, 0.0, 1.0) which would be interpreted by FLUKA as (1/sqrt(2), 0.0, 1/sqrt(2)), i.e. a primary traveling in the x-z plane at a 45-degree angle to the z axis. Indeed, the resulting angular distribution is focused along the positive z axis, but within a 45-degree cone due to the correlation between the direction cosines.
Seeing what an expected pencil beam is turned into, it is no surprise that your distribution, which has a moderate dependence on theta, becomes almost (as mentioned above, what your USRBIN shows is some combination of an isotropic and your expected distribution) isotropic.
In the second forum thread you are referring to (Implementation of an anisotropic source from MCNP in source_newgen.f), @giorgia.fossati’s second proposed method is similar to your source routine, but @horvathd deems it incorrect. The first proposed method there correctly takes into account the correlation between the direction cosines.
Now, on top of sampling from a complex angular distribution, you would also like to rotate that distribution. Personally, I would sample the direction cosines in the standard coordinate system and then apply a rotation matrix to the resulting vector. All of that can be hard coded in source_newgen.f. However, I have a feeling that there is a less error-prone and more efficient way to do it in FLUKA, and I’ll let the experts answer that part of the question.
Cheers,
Udo