(no subject)

040910
X-Accept-Language: fr, en
MIME-Version: 1.0
To: Finetti Noemi <noemi.finetti@aquila.infn.it>
CC: fluka-discuss@fluka.org
Subject: Re: use of CERNLIB with Fluka (source)
References: <200803262343.m2QNhJui013817@smtp1.mi.infn.it>
In-Reply-To: <200803262343.m2QNhJui013817@smtp1.mi.infn.it>
Content-Type: text/plain; charset=3DISO-8859-1; format=3Dflowed
Content-Transfer-Encoding: 8bit
Sender: owner-fluka-discuss@mi.infn.it

Hello Noemi,

About the second point of your mail : isotropy sampling can be done with=20
the use of RACO routine. It returns three randoms numbers Txx, Tyy and=20
Tzz such as
Txx=B2+Tyy=B2+Tzz=B2 =3D 1.D+00 (please see following link :=20
http://www.fluka.org/frequentlyAQ/Functions,_subroutines_and_Fortran_progra=
mming.html=20
)
You could also find a solution for your first problem.

Hope it helps.
Regards.
Sebastien.

Finetti Noemi wrote :

> Dear Fluka users,
> Some questions about the use of the user written routine source:
> 1) I would like to use the function subprogram=20
> DGAUSS(fu,emin,emax,0.001) - CERNLIB (MATHLIB) - in my source=20
> subroutine (see the attached file source_prot_min_new.f). I declared=20
> the function "fu" as external (it is reported at the end of the source=20
> subroutine). I used the lflukac script in order to link Fluka with=20
> CERN Library but at the end of the procedure I obtained the following=20
> error message:
> "source_prot_min_new.f:130: undefined reference to `dgauss_'
> collect2: ld returned 1 exit status". Can you find out where is the=20
> error?
> 2) In source.f the sum=20
> (TXFLK(NPFLKA))**2+(TYFLK(NPFLKA))**2+(TXFLK(NPFLKA))**2 must be=20
> exactly equal 1 (double precision). In my source_prot_min_new.f I set=20
> TXFLK(NPFLKA)=3Dcostx and TYFLK(NPFLKA)=3Dcosty (where costx and costy=20
> were read from an external data file) and I obtained TZFLK(NPFLKA)=20
> from the algorithm (which was commented on the original source.f in=20
> usermvax) reported at line 228 TZFLK(NPFLKA) =3D SQRT(ONEONE -=20
> TXFLK(NPFLKA)**2- TYFLK(NPFLKA)**2 ) but sometime it happen that the=20
> argument of SQRT is negative , moreover by using this algorithm=20
> TZFLK(NPFLKA) is always positive while I need also negative values in=20
> order to simulate an isotropic particle flux. How can I solve these=20
> problems?
> Thanks in advance,
> noemi
>
>------------------------------------------------------------------------
>
>*$ CREATE SOURCE.FOR
>*COPY SOURCE
>*
>*=3D=3D=3D source =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D*
>*
> SUBROUTINE SOURCE ( NOMORE )
>
> INCLUDE '(DBLPRC)'
> INCLUDE '(DIMPAR)'
> INCLUDE '(IOUNIT)'
>*
>*----------------------------------------------------------------------*
>* *
>* Copyright (C) 1990-2006 by Alfredo Ferrari & Paola Sala *
>* All Rights Reserved. *
>* *
>* *
>* New source for FLUKA9x-FLUKA200x: *
>* *
>* Created on 07 january 1990 by Alfredo Ferrari & Paola Sala *
>* Infn - Milan *
>* *
>* Last change on 03-mar-06 by Alfredo Ferrari *
>* *
>* This is just an example of a possible user written source routine. *
>* note that the beam card still has some meaning - in the scoring the *
>* maximum momentum used in deciding the binning is taken from the *
>* beam momentum. Other beam card parameters are obsolete. *
>* *
>*----------------------------------------------------------------------*
>*
>c INPUT FILE: pernoemi.dat
>c=20
>c OUTPUT FILE:check.dat
>c
>c Modified 04 DEC 2007: COMMON /NOEMI/ to be shared with usrmed*=20
>c with the values of
>c kkpart,ijpart,atomicN,massN,enkintot,pmom,fluxtime=
=2E
>c Modified 04 DEC 2007: in order to read from Prm.dat the value=20
>c of fluxtime the format 50 has been changed!
>c Modified 04 DEC 2007: in order to write in check.dat also the fluxtime=
=20
>c value format 51 has been changed.
>c----------------------------------------------------------------------c
>c
> INCLUDE '(BEAMCM)'
> INCLUDE '(FHEAVY)'
> INCLUDE '(FLKSTK)'
> INCLUDE '(IOIOCM)'
> INCLUDE '(LTCLCM)'
> INCLUDE '(PAPROP)'
> INCLUDE '(SOURCM)'
> INCLUDE '(SUMCOU)'
>*
> EXTERNAL fu
> COMMON/PAWC/h(500000)
> REAL h
> COMMON /NOEMIFU/ GF,par1,par2,par3,par4
>*
> LOGICAL LFIRST
>*
> INTEGER kk,atomicN,massN
>c REAL*8 xi,yi,zi,costx,costy,costz,enkin,enkintot
>*
> INTEGER kkpart,ijpart
>c REAL*8 pmom,fluxtime
> =20
>c COMMON /NOEMI/ kkpart,ijpart,atomicN,massN,enkintot,pmom,fluxtime
> COMMON /NOEMI/ enkintot,pmom,fluxtime,kkpart,ijpart,atomicN,massN
>*
> SAVE LFIRST
> DATA LFIRST / .TRUE. /
>*
>*=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D*
>* *
>* BASIC VERSION *
>* *
>*=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D*
> NOMORE =3D 0
>* +-------------------------------------------------------------------*
>* | First call initializations:
> IF ( LFIRST ) THEN
>* | *** The following 3 cards are mandatory ***
> TKESUM =3D ZERZER
> LFIRST =3D .FALSE.
> LUSSRC =3D .TRUE.
>* | *** User initialization ***
>* return message from first call
> WRITE(LUNOUT,*)' '
> WRITE(LUNOUT,*)'*** PROT_MIN Version of Routine SOURCE called ...'
> WRITE(LUNOUT,*)' '
> END IF
>* |
>* +-------------------------------------------------------------------*
> OPEN(21,
> +FILE=3D
> +'/home/finetti/flukawork/lisa/prot_min/new25mar08/'//
> +'Prm.dat',
> +status=3D'old')
> OPEN(22,
> +FILE=3D
> +'/home/finetti/flukawork/lisa/prot_min/new25mar08/'//
> +'check.dat',
> +status=3D'unknown')
>
> READ(21,50) kk,xi,yi,zi,costx,costy,costz,enkin,fulltime
> 50 format(i10,6f7.3,f10.3,f10.3)
>
>Checking ...
> WRITE(22,50)kk,xi,yi,zi,costx,costy,costz,enkin,fulltime
>c------------------
>c Geometric Factor:
> radius=3D7.1 !sphere radius in cm
> pi=3D3.14159265358979323846
> GF=3D4*(pi**2)*((radius/100.)**2) ! sphere geometric factor in m**2 =
sr
>c------------------
>c Proton at solar minimum:
> emin=3D0.1 ! minimum energy value in GeV/A
> emax=3D100. ! maximum energy value in GeV/A
>c------------------
> par1=3D18000.
> par2=3D1.09
> par3=3D3.66
> par4=3D0.87
>c------------------
> ab=3DDGAUSS(fu,emin,emax,0.001)
> fluxtime=3Dkk/ab
>c------------------
> kkpart=3Dkk
> ijpart=3DIJBEAM !Incident particle type
> atomicN=3D1 !Atomic Number Z (integer) 1=3Dprotons
> massN=3D1 !Atomic Mass Number A (integer) 1=3Dprotons
> enkintot=3Denkin*massN !Total kinetic energy
>
> WRITE(22,*)
> +' Event N., IJBEAM, Atomic N., Mass N., En.Kin./A, En.Kin.Tot.,=20
> +fluxtime'
> WRITE(22,51) kk,IJBEAM,atomicN,massN,enkin,enkintot,fluxtime
> 51 FORMAT(i8,1x,i8,1x,i8,1x,i8,1x,1f10.3,1x,1f10.3,1x,f10.3)
>* +-------------------------------------------------------------------*
>* Push one source particle to the stack. Note that you could as well
>* push many but this way we reserve a maximum amount of space in the
>* stack for the secondaries to be generated
>* Npflka is the stack counter: of course any time source is called it
>* must be =3D0
> NPFLKA =3D NPFLKA + 1
>* Wt is the weight of the particle
> WTFLK (NPFLKA) =3D ONEONE
> WEIPRI =3D WEIPRI + WTFLK (NPFLKA)
>* Particle type (1=3Dproton.....). Ijbeam is the type set by the BEAM
>* card
>* +-------------------------------------------------------------------*
>* | (Radioactive) isotope:
> IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
> IARES =3D IPROA
> IZRES =3D IPROZ
> IISRES =3D IPROM
> CALL STISBM ( IARES, IZRES, IISRES )
> IJHION =3D IPROZ * 1000 + IPROA
> IJHION =3D IJHION * 100 + KXHEAV
> IONID =3D IJHION
> CALL DCDION ( IONID )
> CALL SETION ( IONID )
>* |
>* +-------------------------------------------------------------------*
>* | Heavy ion:
> ELSE IF ( IJBEAM .EQ. -2 ) THEN
> IJHION =3D IPROZ * 1000 + IPROA
> IJHION =3D IJHION * 100 + KXHEAV
> IONID =3D IJHION
> CALL DCDION ( IONID )
> CALL SETION ( IONID )
> ILOFLK (NPFLKA) =3D IJHION
>* | Flag this is prompt radiation
> LRADDC (NPFLKA) =3D .FALSE.
>* |
>* +-------------------------------------------------------------------*
>* | Normal hadron:
> ELSE
> IONID =3D IJBEAM
> ILOFLK (NPFLKA) =3D IJBEAM
>* | Flag this is prompt radiation
> LRADDC (NPFLKA) =3D .FALSE.
> END IF
>* |
>* +-------------------------------------------------------------------*
>* From this point .....
>* Particle generation (1 for primaries)
> LOFLK (NPFLKA) =3D 1
>* User dependent flag:
> LOUSE (NPFLKA) =3D 0
>* User dependent spare variables:
> DO 100 ISPR =3D 1, MKBMX1
> SPAREK (ISPR,NPFLKA) =3D ZERZER
> 100 CONTINUE
>* User dependent spare flags:
> DO 200 ISPR =3D 1, MKBMX2
> ISPARK (ISPR,NPFLKA) =3D 0
> 200 CONTINUE
>* Save the track number of the stack particle:
> ISPARK (MKBMX2,NPFLKA) =3D NPFLKA
> NPARMA =3D NPARMA + 1
> NUMPAR (NPFLKA) =3D NPARMA
> NEVENT (NPFLKA) =3D 0
> DFNEAR (NPFLKA) =3D +ZERZER
>* ... to this point: don't change anything
>* Particle age (s)
> AGESTK (NPFLKA) =3D +ZERZER
> AKNSHR (NPFLKA) =3D -TWOTWO
>* Group number for "low" energy neutrons, set to 0 anyway
> IGROUP (NPFLKA) =3D 0
>* Kinetic energy of the particle (GeV)
>C TKEFLK (NPFLKA) =3D SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
> TKEFLK (NPFLKA) =3D enkintot
>* Particle momentum
>C PMOFLK (NPFLKA) =3D PBEAM
> PMOFLK (NPFLKA) =3D SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
> & + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
> pmom=3DPMOFLK (NPFLKA)
>
> write(22,*)' Particle mass (in SOURCE)=3D',AM (ILOFLK(NPFLKA))
> write(22,*)' Particle kin. en. (in SOURCE)=3D',TKEFLK (NPFLKA)
> write(22,*)' Particle momentum (in SOURCE)=3D',PMOFLK (NPFLKA)
>
>* Cosines (tx,ty,tz)
>C TXFLK (NPFLKA) =3D UBEAM
>C TYFLK (NPFLKA) =3D VBEAM
>C TZFLK (NPFLKA) =3D WBEAM
>c* TZFLK (NPFLKA) =3D SQRT ( ONEONE - TXFLK (NPFLKA)**2
>c* & - TYFLK (NPFLKA)**2 )
>
> TXFLK (NPFLKA) =3D costx
> TYFLK (NPFLKA) =3D costy
> TZFLK (NPFLKA) =3D costz
> write(22,*)' '
> write(22,*)' Particle costx (in SOURCE)=3D',TXFLK (NPFLKA)=20
> write(22,*)' Particle costy (in SOURCE)=3D',TYFLK (NPFLKA)
> write(22,*)' Particle costz (in SOURCE)=3D',TZFLK (NPFLKA)
>
> IF(costz.ne.0.) THEN
> TZFLK (NPFLKA) =3D (costz/abs(costz))* !!!! Noemi Finetti 19ma=
r08
> + SQRT (abs( ONEONE - TXFLK (NPFLKA)**2
> & - TYFLK (NPFLKA)**2 ) )
> ELSE
> TZFLK (NPFLKA) =3D !!!! Noemi Finetti 25ma=
r08
> + SQRT (abs( ONEONE - TXFLK (NPFLKA)**2
> & - TYFLK (NPFLKA)**2 ) )
> ENDIF
>
> write(22,*)' '
> write(22,*)' Particle costz (in SOURCE)=3D',TZFLK (NPFLKA)
> write(22,*)' costx**2+costy**2+costz**2=3D',
> +(TXFLK (NPFLKA))**2+(TYFLK (NPFLKA))**2+(TZFLK (NPFLKA))**2
> write(22,*)' '
>
>* Polarization cosines:
> TXPOL (NPFLKA) =3D -TWOTWO
> TYPOL (NPFLKA) =3D +ZERZER
> TZPOL (NPFLKA) =3D +ZERZER
>* Particle coordinates
>C XFLK (NPFLKA) =3D XBEAM
>C YFLK (NPFLKA) =3D YBEAM
>C ZFLK (NPFLKA) =3D ZBEAM
> XFLK (NPFLKA) =3D xi
> YFLK (NPFLKA) =3D yi
> ZFLK (NPFLKA) =3D zi
>
>* Calculate the total kinetic energy of the primaries: don't change
> IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
> & THEN
> TKESUM =3D TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
> ELSE IF ( ILOFLK (NPFLKA) .NE. 0 ) THEN
> TKESUM =3D TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) )
> & * WTFLK (NPFLKA)
> ELSE
> TKESUM =3D TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
> END IF
> RADDLY (NPFLKA) =3D ZERZER
>* Here we ask for the region number of the hitting point.
>* NREG (NPFLKA) =3D ...
>* The following line makes the starting region search much more
>* robust if particles are starting very close to a boundary:
> CALL GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
> CALL GEOREG ( XFLK (NPFLKA), YFLK (NPFLKA), ZFLK (NPFLKA),
> & NRGFLK(NPFLKA), IDISC )
>* Do not change these cards:
> CALL GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
> NLATTC (NPFLKA) =3D MLATTC
> CMPATH (NPFLKA) =3D ZERZER
> CALL SOEVSV
> RETURN
>*=3D=3D=3D End of subroutine Source =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D*
> END
>******* External function fu ******************************************
> FUNCTION fu(X)
> REAL*8 GF,par1,par2,par3,par4
> COMMON /NOEMIFU/ GF,par1,par2,par3,par4
>c
> fu=3DGF*(par1*(x+par2)**(-par3)*x**par4)
>c print*,'func in fuction',fu,' Energia=3D',x,' GeV/A'
> return
> end
>***********************************************************************
>
>
> =20
>
Received on Thu Mar 27 15:02:27 2008