Re: Source.f Subroutine for 2 Sources (2nd attempt)

From: Alberto Fasso' <fasso_at_mail.cern.ch>
Date: Thu, 14 Apr 2011 15:53:59 +0200

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