Dear Fluka users,
Some questions about the use of the user written routine source:
1) I would like to use the function subprogram
DGAUSS(fu,emin,emax,0.001) - CERNLIB (MATHLIB) - in my source subroutine
(see the attached file source_prot_min_new.f). I declared the function
"fu" as external (it is reported at the end of the source subroutine). I
used the lflukac script in order to link Fluka with CERN Library but at
the end of the procedure I obtained the following 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 error?
2) In source.f the sum
(TXFLK(NPFLKA))**2+(TYFLK(NPFLKA))**2+(TXFLK(NPFLKA))**2 must be exactly
equal 1 (double precision). In my source_prot_min_new.f I set
TXFLK(NPFLKA)=costx and TYFLK(NPFLKA)=costy (where costx and costy were
read from an external data file) and I obtained TZFLK(NPFLKA) from the
algorithm (which was commented on the original source.f in usermvax)
reported at line 228 TZFLK(NPFLKA) = SQRT(ONEONE - TXFLK(NPFLKA)**2-
TYFLK(NPFLKA)**2 ) but sometime it happen that the argument of SQRT is
negative , moreover by using this algorithm TZFLK(NPFLKA) is always
positive while I need also negative values in order to simulate an
isotropic particle flux. How can I solve these problems?
Thanks in advance,
noemi
-- Address: Dott.ssa Noemi Finetti c/o Dipartimento di Fisica dell'Universita' di L'Aquila Via Vetoio - 67010 Coppito - L'Aquila - Italy Phone(office):+39-0862-433051; Fax(department):+39-0862-433033 "NON SI NASCONDE FUORI DEL MONDO CHI LO SALVA E NON LO SA. E' UNO COME NOI, NON DEI MIGLIORI." EUGENIO MONTALE
*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
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
c OUTPUT FILE:check.dat
c
c Modified 04 DEC 2007: COMMON /NOEMI/ to be shared with usrmed*
c with the values of
c kkpart,ijpart,atomicN,massN,enkintot,pmom,fluxtime.
c Modified 04 DEC 2007: in order to read from Prm.dat the value
c of fluxtime the format 50 has been changed!
c Modified 04 DEC 2007: in order to write in check.dat also the fluxtime
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
c COMMON /NOEMI/ kkpart,ijpart,atomicN,massN,enkintot,pmom,fluxtime
COMMON /NOEMI/ enkintot,pmom,fluxtime,kkpart,ijpart,atomicN,massN
*
SAVE LFIRST
DATA LFIRST / .TRUE. /
*
*======================================================================*
* *
* BASIC VERSION *
* *
*======================================================================*
NOMORE = 0
* +-------------------------------------------------------------------*
* | First call initializations:
IF ( LFIRST ) THEN
* | *** The following 3 cards are mandatory ***
TKESUM = ZERZER
LFIRST = .FALSE.
LUSSRC = .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=
+'/home/finetti/flukawork/lisa/prot_min/new25mar08/'//
+'Prm.dat',
+status='old')
OPEN(22,
+FILE=
+'/home/finetti/flukawork/lisa/prot_min/new25mar08/'//
+'check.dat',
+status='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=7.1 !sphere radius in cm
pi=3.14159265358979323846
GF=4*(pi**2)*((radius/100.)**2) ! sphere geometric factor in m**2 sr
c------------------
c Proton at solar minimum:
emin=0.1 ! minimum energy value in GeV/A
emax=100. ! maximum energy value in GeV/A
c------------------
par1=18000.
par2=1.09
par3=3.66
par4=0.87
c------------------
ab=DGAUSS(fu,emin,emax,0.001)
fluxtime=kk/ab
c------------------
kkpart=kk
ijpart=IJBEAM !Incident particle type
atomicN=1 !Atomic Number Z (integer) 1=protons
massN=1 !Atomic Mass Number A (integer) 1=protons
enkintot=enkin*massN !Total kinetic energy
WRITE(22,*)
+' Event N., IJBEAM, Atomic N., Mass N., En.Kin./A, En.Kin.Tot.,
+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 =0
NPFLKA = NPFLKA + 1
* Wt is the weight of the particle
WTFLK (NPFLKA) = ONEONE
WEIPRI = WEIPRI + WTFLK (NPFLKA)
* Particle type (1=proton.....). Ijbeam is the type set by the BEAM
* card
* +-------------------------------------------------------------------*
* | (Radioactive) isotope:
IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
IARES = IPROA
IZRES = IPROZ
IISRES = IPROM
CALL STISBM ( IARES, IZRES, IISRES )
IJHION = IPROZ * 1000 + IPROA
IJHION = IJHION * 100 + KXHEAV
IONID = IJHION
CALL DCDION ( IONID )
CALL SETION ( IONID )
* |
* +-------------------------------------------------------------------*
* | Heavy ion:
ELSE IF ( IJBEAM .EQ. -2 ) THEN
IJHION = IPROZ * 1000 + IPROA
IJHION = IJHION * 100 + KXHEAV
IONID = IJHION
CALL DCDION ( IONID )
CALL SETION ( IONID )
ILOFLK (NPFLKA) = IJHION
* | Flag this is prompt radiation
LRADDC (NPFLKA) = .FALSE.
* |
* +-------------------------------------------------------------------*
* | Normal hadron:
ELSE
IONID = IJBEAM
ILOFLK (NPFLKA) = IJBEAM
* | Flag this is prompt radiation
LRADDC (NPFLKA) = .FALSE.
END IF
* |
* +-------------------------------------------------------------------*
* From this point .....
* Particle generation (1 for primaries)
LOFLK (NPFLKA) = 1
* User dependent flag:
LOUSE (NPFLKA) = 0
* User dependent spare variables:
DO 100 ISPR = 1, MKBMX1
SPAREK (ISPR,NPFLKA) = ZERZER
100 CONTINUE
* User dependent spare flags:
DO 200 ISPR = 1, MKBMX2
ISPARK (ISPR,NPFLKA) = 0
200 CONTINUE
* Save the track number of the stack particle:
ISPARK (MKBMX2,NPFLKA) = NPFLKA
NPARMA = NPARMA + 1
NUMPAR (NPFLKA) = NPARMA
NEVENT (NPFLKA) = 0
DFNEAR (NPFLKA) = +ZERZER
* ... to this point: don't change anything
* Particle age (s)
AGESTK (NPFLKA) = +ZERZER
AKNSHR (NPFLKA) = -TWOTWO
* Group number for "low" energy neutrons, set to 0 anyway
IGROUP (NPFLKA) = 0
* Kinetic energy of the particle (GeV)
C TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
TKEFLK (NPFLKA) = enkintot
* Particle momentum
C PMOFLK (NPFLKA) = PBEAM
PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
& + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
pmom=PMOFLK (NPFLKA)
write(22,*)' Particle mass (in SOURCE)=',AM (ILOFLK(NPFLKA))
write(22,*)' Particle kin. en. (in SOURCE)=',TKEFLK (NPFLKA)
write(22,*)' Particle momentum (in SOURCE)=',PMOFLK (NPFLKA)
* Cosines (tx,ty,tz)
C TXFLK (NPFLKA) = UBEAM
C TYFLK (NPFLKA) = VBEAM
C TZFLK (NPFLKA) = WBEAM
c* TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
c* & - TYFLK (NPFLKA)**2 )
TXFLK (NPFLKA) = costx
TYFLK (NPFLKA) = costy
TZFLK (NPFLKA) = costz
write(22,*)' '
write(22,*)' Particle costx (in SOURCE)=',TXFLK (NPFLKA)
write(22,*)' Particle costy (in SOURCE)=',TYFLK (NPFLKA)
write(22,*)' Particle costz (in SOURCE)=',TZFLK (NPFLKA)
IF(costz.ne.0.) THEN
TZFLK (NPFLKA) = (costz/abs(costz))* !!!! Noemi Finetti 19mar08
+ SQRT (abs( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 ) )
ELSE
TZFLK (NPFLKA) = !!!! Noemi Finetti 25mar08
+ SQRT (abs( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 ) )
ENDIF
write(22,*)' '
write(22,*)' Particle costz (in SOURCE)=',TZFLK (NPFLKA)
write(22,*)' costx**2+costy**2+costz**2=',
+(TXFLK (NPFLKA))**2+(TYFLK (NPFLKA))**2+(TZFLK (NPFLKA))**2
write(22,*)' '
* Polarization cosines:
TXPOL (NPFLKA) = -TWOTWO
TYPOL (NPFLKA) = +ZERZER
TZPOL (NPFLKA) = +ZERZER
* Particle coordinates
C XFLK (NPFLKA) = XBEAM
C YFLK (NPFLKA) = YBEAM
C ZFLK (NPFLKA) = ZBEAM
XFLK (NPFLKA) = xi
YFLK (NPFLKA) = yi
ZFLK (NPFLKA) = zi
* Calculate the total kinetic energy of the primaries: don't change
IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
& THEN
TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
ELSE IF ( ILOFLK (NPFLKA) .NE. 0 ) THEN
TKESUM = TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) )
& * WTFLK (NPFLKA)
ELSE
TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
END IF
RADDLY (NPFLKA) = ZERZER
* Here we ask for the region number of the hitting point.
* NREG (NPFLKA) = ...
* 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) = MLATTC
CMPATH (NPFLKA) = ZERZER
CALL SOEVSV
RETURN
*=== End of subroutine Source =========================================*
END
******* External function fu ******************************************
FUNCTION fu(X)
REAL*8 GF,par1,par2,par3,par4
COMMON /NOEMIFU/ GF,par1,par2,par3,par4
c
fu=GF*(par1*(x+par2)**(-par3)*x**par4)
c print*,'func in fuction',fu,' Energia=',x,' GeV/A'
return
end
***********************************************************************
Received on Thu Mar 27 00:46:16 2008
This archive was generated by hypermail 2.1.8 : Thu Mar 27 2008 - 00:46:16 CET