*) from Alex/
*) how to include the divergence in the source.f / photons
---------------------------------------
Dear experts,
I am using the attached source.f to sample the designed photon
spectrum, however, I noticed that BEAM WHAT(3) does not respond
to the values if source routine is activated – just as stated in the manual.
Is there any “simple” way to make my source beam with divergence
200 mrad, or only additional piece of code must be included in the source.f?
Looking through old forum archives I found similar code for that purpose:
Divbm_Sigma = Divbm*0.2
CALL FLNRR2 (RGAUSSDIVX,RGAUSSDIVY)
DIV_VECTORX = RGAUSSDIVX*Divbm_Sigma
DIV_VECTORY = RGAUSSDIVY*Divbm_Sigma
TXHLP = TAN( DIV_VECTORX )
TYHLP = TAN( DIV_VECTORY )
THELP = SQRT( TXHLP*TXHLP + TYHLP*TYHLP + ONEONE )
TXFLK (NPFLKA) = TXHLP / THELP
TYFLK (NPFLKA) = TYHLP / THELP
TZFLK (NPFLKA) = SQRT ( ONEONE - LTXFLK (NPFLKA)**2 - TYFLK (NPFLKA)**2 )
TXFLK (NPFLKA) = UBEAM
TYFLK (NPFLKA) = VBEAM
TZFLK (NPFLKA) = WBEAM
TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 )
But simply applying it the code is not compiled, and my lack of
experience with user routines is a big obstacle to move forward from
this point. In addition, my BEAMPOS is defined with WHAT(1)-(WHAT(3)
and also rotation value for WHAT(5) is used.
Also the beam is directed to positive Z direction.
Could anyone share a similar working source.f file which would directly
use BEAM WhHAT(3) value, or maybe guide me through this coding problem?
Both source.f routines – the original and the one I tried to modify –
are attached.
Thank you very much,
Rgds
Alex
source_circ.f source_div.f
----------------
==> Dear Alex,
Please note following if it can help you.
cos 200 mrad=0.98007
So we have to randomize direction cosine along beam axis
from 0.98007 to 1 (i.e. cos (0 degree) to cos (200 mrad))
you can try following code as
* Cosines (tx,ty,tz)
COSTH = ONEONE-(1-0.98007)*FLRNDM(TTT)
IF(COSTH.GT.ONEONE) COSTH=ONEONE
IF(COSTH.LT.-ONEONE) COSTH=-ONEONE
PHI = TWOPIP*FLRNDM(SSS)
UBEAM=SQRT(ONEONE-COSTH*COSTH)*COS(PHI)
VBEAM=SQRT(ONEONE-COSTH*COSTH)*SIN(PHI)
WBEAM=COSTH
TXFLK (NPFLKA) = UBEAM
TYFLK (NPFLKA) = VBEAM
* TZFLK (NPFLKA) = WBEAM
TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 )
I think it will work.
Regards,
--
SOUGATA RAKSHIT / सौगत रक्षित
-------------------
Dear Sougata Rakshit,
Wonderful, it works indeed!
But may I ask two more moments here:
1) I need to set up the certain values ( 0.2 ) for BEAMPOS WHAT(4) & (5) -
direction cosine of the beam with respect to the x & y-axis.
It does not take values from the card now when I use this modified
routine. What additional line should be included for this?
2) also this code works for angular flat beam only, but not for
rectangular flat beam - does the code differs much in such case?
If you or anybody else could find an opportunity to help with this,
I will greatly appreciate such help,
Regards
Alex
-------------------------
Dear Alex,
For your query (1).
Any parameter from source card may be passed to source.f file
through whasou(#).
e.g
Coding for 200 mrad divergence beam where WHAT(3)=0.98007 im source card.
* Cosines (tx,ty,tz)
COSTH = ONEONE-(1-whasou(3))*FLRNDM(TTT)
IF(COSTH.GT.ONEONE) COSTH=ONEONE
IF(COSTH.LT.-ONEONE) COSTH=-ONEONE
PHI = TWOPIP*FLRNDM(SSS)
UBEAM=SQRT(ONEONE-COSTH*COSTH)*COS(PHI)
VBEAM=SQRT(ONEONE-COSTH*COSTH)*SIN(PHI)
WBEAM=COSTH
TXFLK (NPFLKA) = UBEAM
TYFLK (NPFLKA) = VBEAM
* TZFLK (NPFLKA) = WBEAM
TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 )
Please check.
Regarding your query (2).
I am not so experienced. Better consult some expert person.
Regards,
--
SOUGATA RAKSHIT / सौगत रक्षित
-----------------
*) [fluka-discuss]: how to include the divergence in the source.f / photons
Thank You, Sougata,
I meant BEAMPOS WHAT(4) & (5) , i.e. the ones which gives
the direction cosines of the beam itself.
The divergence works itself great when beam is straight,
with BEAMPOS WHAT(4)&(5) = 0...
In case anyone can help with this question, also with the coding
for the rectangular beam divergence, kindly please,
Regards
Alex
---------------
Hey Alex,
If I get you right, you want to sample beam particle initial directions
with a rectangular divergence (i.e. cosine uniformly sampled in a given
interval) around an arbitrary direction theta_f and phi_f in your
modified source.f.
On paper, what one would ordinarily do is first sample the initial
direction with the desired rectangular divergence around the z axis and,
then, apply a rotation onto the final direction. So you're almost there.
- I believe for a rectangular beam the divergence is to be understood as
a full width. Thus, to sample the (cosine of the) polar angle one
would rather do
CDH = COS ( HLFHLF * DIV )
COSTH = ONEONE - ( ONEONE - CDH ) * FLRNDM ( TTT )
alternatively
COSTH = CDH + ( ONEONE - CDH ) * FLRNDM ( TTT )
where DIV is the divergence in radians, and the azimuthal angle you
sample uniformly from 0 to 2*pi like you already do. By the way, the
DIVBM variable should still have the value you passed in WHAT(3) of
your BEAM card, already converted to radians, but double check in your
particular case.
- You then construct T{X,Y,Z}FLK and ensure it is unit normalized. If
you chose to do smth a la TZFLK(NPFLKA) =SQRT(...) keep an eye on its
sign (for big DIV).
- The resulting T{X,Y,Z}FLK should be fine, but they are still sampled
around the z axis. Let d_i=(tx,ty,tz). To get the "final" direction,
d_f=(tx',ty',tz'), one would apply on d_i a rotation of angle theta_f
around the y axis, R(theta_f,y), and, next, a rotation of angle phi_f
around the z axis, R(phi_f,z).
To be explicit, see the (hopefully not too far off) attached recipe,
without any particular effort to be numerically robust (which is a bad
idea... and exactly the reason why director cosines are preferrable
instead of angles...), but it will get you started. At the very least,
make sure that the final T{X,Y,Z}FLK components are normalized (such
that the sum of their squares is ONEONE).
Cesc
4591_001.pdf -> how to rotate !
----------------
Alex,
As a corollary, you may want to consider calling the slightly more
general FLUKA subroutine TRANS, e.g.
CALL TRANS ( XO, YO, ZO, DE, SFE, CFE, X, Y, Z )
See attached figure. This subroutine is intended for updating a given
original particle direction d_O with director cosines {XO,YO,ZO}, by a
polar scattering angle DE, and an azimuthal scattering angle phi_e with
sine and cosine given by SFE and CFE. The final direction cosines are
{X,Y,Z}.
At first glance I think you can call TRANS passing as {XO,YO,ZO} the
unit vector with theta_f and phi_f of the previous mail, DE is the
sampled polar angle (within the prescribed "divergence") and SFE and CFE
should be the sine and cosine of the uniformly sampled phi_e in the
previous mail. You may want to double check this in case of faux pas on
my side...
The result should be essentially the same as for the previous recipe,
the advantage being that with the subroutine you shouldn't have worries
regarding numerical robustness.
Cesc
-> 4595_001.pdf
//=================================