Dear Lewis,
I have several remarks about your source routine.
1) you set the relative sampling between the two sources in the IF-block
of the first call initialization
IF ( LFIRST ) THEN
......
* Set up relative sampling frequency between 2 sources
SMPLEP = FLRNDM(XDUMMY)
END IF
and then, outside that IF-block, you use the value of SMPLEP to decide
which of the two sources will be used in the present call. But SMPLEP will
never change, because the program enters that IF-block only once: the
value of SMPLEP is sampled randomly once and then is fixed for all the rest
of the run.
2) You define the two source particles using the WHASOUs in the input SOURCE
card:
* First Particle Type
PRTONE = WHASOU(3)
* Second particle Type
PRTTWO = WHASOU(4)
I don't know how did you input WHASOU(3) and WHASOU(4) (you didn't send
the input file). Did you use the particle names or their numbers?
If you used names, the PRTONE and PRTTWO variables should have been declared
as CHARACTER*8 at the beginning of the routine: but anyway particle names
cannot be used in a user routine: user routines understand only numbers,
unless you call a special name-number conversion routine.
On the other hand, if you used numbers, they should be integers: but
PRTONE and PRTTWO are not integers: by default variables with names not
beginning with I,J,K,L,M,N are double precision reals. You must declare them
as integers at the beginning of the subroutine.
3) The IF-block starting with
IF (SMPLEP.LE.0.08693) THEN
does not have the corresponding END IF statement.
You said that you "have managed to get the subroutine to compile and build
ok", but I am sure that the compiler should not accept a missing END IF.
4) You sample the energy from a Gaussian as follows:
ENRGY = (1/(SQRT(PI*FWHM*FWHM)/(4.0d0*LOG(2.0d0))))
& * EXP(-(((FLRNDM(XDUMMY)*EPRTONE)-EPRTONE)
& *4.0d0*LOG(2.0d0))/(FWHM*FWHM))
but PI has not been defined (it is defined only later in the routine,
as PI = 3.14159 - with only 6 significant figures, in single precision).
You should use instead the FLUKA constant PIPIPI, predefined in the INCLUDE
file (DBLPRC) as PIPIPI = 3.141592653589793238462643383279D+00
Also, you don't need to write your own gaussian sampling (I have not checked
that is correct). You should use the FLUKA subroutine that samples
from a Gaussian distribution with unit variance:
CALL FLNRRN (RGAUSS)
5) This is not an error, but you have two successive lines:
PMOFLK (NPFLKA) = PBEAM
PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
& + TWOTWO * AM (IONID) ) )
The second lines overrides the first, so the first is useless. You
better comment it out.
6) In the lines
COSTHETA = 1-(1 - COSTMAX) * FLRNDM(XDUMMY)
SINTHETA = SQRT(1 - (COSTHETA)**2)
PI = 3.14159
PHI = 2.0d0 * PI * FLRNDM(XDUMMY)
you have again PI instead of PIPIPI. And 1 repeatedly written as an
integer, instead of double precision 1.D0 (or predefined FLUKA constant
ONEONE).
7) In the lines
* Send electrons in -Z direction
IF PTYPE = PRTONE THEN
TXFLK (NPFLKA) = SINTHETA * COS(PHI)
TYFLK (NPFLKA) = SINTHETA * SIN(PHI)
TZFLK (NPFLKA) = -SQRT( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 )
*
* Send photons in +Z direction
ELSE IF PTYPE = PRTTWO THEN
TXFLK (NPFLKA) = SINTHETA * COS(PHI)
TYFLK (NPFLKA) = SINTHETA * SIN(PHI)
TZFLK (NPFLKA) = SQRT( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 )
END IF
the statement IF PTYPE = PRTONE THEN should be written
IF (PTYPE .EQ. PRTONE) THEN
and ELSE IF PTYPE = PRTTWO THEN should be written
ELSE IF (PTYPE .EQ. PRTTWO) THEN
(with parentheses, and .EQ. instead of =). Again, I cannot believe that
the compiler didn't complain.
Notice also (this is not an error) that you repeat twice, in the IF
and in the ELSE IF:
TXFLK (NPFLKA) = SINTHETA * COS(PHI)
TYFLK (NPFLKA) = SINTHETA * SIN(PHI)
You could take those two lines outside, before the IF-ELSEIF block.
Kind regards,
Alberto
On Thu, 14 Apr 2011, Macfarlane, Lewis wrote:
> Dear Experts,
>
> I have attempted to write a source.f subroutine to send beams of 2 different
> particles, of different energies in different directions. I have managed to
> get the subroutine to compile and build ok but when I used it in a FLUKA case,
> I get an error as no seeds are started. If you could take a look at the
> subroutine and let me know where it's going wrong, that would be great!
> Apologies, it is a little messy at the moment as my programming experience is
> minimal and I just want to get it working. Routine is attached.
>
> Best regards,
> Lewis.
Received on Thu Apr 14 2011 - 16:37:46 CEST
This archive was generated by hypermail 2.2.0 : Thu Apr 14 2011 - 16:37:46 CEST