*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
      SUBROUTINE SOURCE ( NOMORE )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2009      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     New source for FLUKA9x-FLUKA20xy:                                *
*                                                                      *
*     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 08-feb-09     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.            *
*                                                                      *
*       Output variables:                                              *
*                                                                      *
*              Nomore = if > 0 the run will be terminated              *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'

********************USER defined variables**************************
      INCLUDE '(CASLIM)'
* Some variables defined in START card
* Here Ncases is used
********************************************************************
*
      LOGICAL LFIRST
*

      SAVE LFIRST
      DATA LFIRST / .TRUE. /

********************USER defined variables******************************
*angle divergence of x and y
      common/angdivs/xpp,ypp
************************************************************************

*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
*      WRITE(*,*) 'Hello world!'
*       OPEN(32,FILE='coor.txt',STATUS='NEW',FORM='FORMATTED')
*      WRITE(*,*) Ncases
*       write(*,*) whasou(1), whasou(2), whasou(3), whasou(4)

* h and (aah, bbh)  are the boundaries
* For beam core of Gaussian distribution, 97% is integrated by r from
-b*sigma to b*sigma
      aac = 0.0
      bbc = WHASOU(5)*2.355/WHASOU(1)

* For beam halo of Gaussian distribution, 3% is integrated by r from 0 to
2.649*sigma
*      aah = 0.0
      aah = WHASOU(11) * 2.355 / WHASOU(7)
      bbh = WHASOU(12) * 2.355 / WHASOU(9)
      write(*,*) 'Beam core parameters for inner and outer
     &            boundaries:', aac, bbc
      write(*,*) 'Beam halo parameters for inner and outer
     &  boundaries:', aah, bbh

      END IF
*  |
*  +-------------------------------------------------------------------*
*  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.
*  |  Group number for "low" energy neutrons, set to 0 anyway
         IGROUP (NPFLKA) = 0
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID = IJBEAM
         ILOFLK (NPFLKA) = IJBEAM
*  |  Flag this is prompt radiation
         LRADDC (NPFLKA) = .FALSE.
*  |  Group number for "low" energy neutrons, set to 0 anyway
         IGROUP (NPFLKA) = 0
      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
* Kinetic energy of the particle (GeV)
      TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
* Particle momentum
      PMOFLK (NPFLKA) = PBEAM
*     PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
*    &                       + TWOTWO * AM (IONID) ) )

      Npercent = INT(NCASES * 0.993)

      IF (Ncase .LE. Npercent) THEN

* Particle coordinates
* Disbx and disby are the boundaries of x and y
* For beam core of Gaussian distribution, 97% is integrated by x or y from
-h*sigma to h*sigma
*      h = WHASOU(5)
*      Disbx = h*WHASOU(1)/2.355


* Obtain the x coordinate
*          CALL Cooxy(WHASOU(1),WHASOU(3),aac,bbc,coox,cooy)
*          coox = 0.0
*          cooy = 0.0


* Cosines (tx,ty,tz)
* Epx and Epy are the emittances of x and y(UNIT: Pi cm mrad)
* Epx = Pi * 3sigmax *3sigmax'
* Epy = Pi * 3sigmay *3sigmay'
* Obtain the x' and y' of beam core

       Epx = 8.1D0
       Epy = 8.1D0
* coeff is the ratio coefficiet boundary in order to make the same ratio
boundary
       coeff = WHASOU(5)/(WHASOU(1)/2.355)
       WHASOU(2) = 2.355 * (Epx / (coeff * WHASOU(5)*3.1415926d0))
       WHASOU(4) = 2.355 * (Epy / (coeff *
     &             (WHASOU(5)*WHASOU(3)/WHASOU(1))*3.1415926d0))

*      write(*,*)WHASOU(2),WHASOU(4)
*      CALL coxy(WHASOU(2),WHASOU(4),xpp,ypp,cox,coy)

*       CALL coxp(WHASOU(1),WHASOU(2),coox,xp)
*       CALL coxp(WHASOU(3),WHASOU(4),cooy,yp)
        CALL cxyp(WHASOU(1),WHASOU(2),WHASOU(3),WHASOU(4),
     & aac,bbc,coox,cooy,xpp,ypp)

      ELSE

* Obtain the x and y coordinate of beam halo
*          CALL Cooxy(WHASOU(7),WHASOU(9),aah,bbh,coox,cooy)

* Obtain the x' and y' of beam halo
*      Epx = 35.D0
*      Epy = 35.D0
*
      coeff = WHASOU(12)/(WHASOU(7)/2.355)
*      WHASOU(8) = 2.355 * (Epx / (coeff * WHASOU(12)*3.1415926d0))
*      WHASOU(10) = 2.355 * (Epy / (coeff *
*     &             (WHASOU(12)*WHASOU(9)/WHASOU(7))*3.1415926d0))
* Don't defined the emmittance of halo, make the same ratio for (sx_core /
sx_halo) and
* (xp_core/xp_halo), therefore halo ellipse and core ellipse can link
together.

      WHASOU(8) = WHASOU(7) * WHASOU(2) / WHASOU(1)
      WHASOU(10) = WHASOU(9) * WHASOU(4) / WHASOU(3)
*      write(*,*)WHASOU(2),WHASOU(4)

*       CALL coxp(WHASOU(7),WHASOU(8),coox,xp)
*       CALL coxp(WHASOU(9),WHASOU(10),cooy,yp)
        CALL cxyp(WHASOU(7),WHASOU(8),WHASOU(9),WHASOU(10),
     & aah,bbh,coox,cooy,xpp,ypp)


      END IF
*      write(*,*)coeff,WHASOU(2),WHASOU(4),WHASOU(8),WHASOU(10)
*      write(*,*) 'Current number of beam particles:', Ncase
*      write(32,FMT='(4(1X,F12.5))')  coox*10.,xpp,cooy*10.,ypp
*      write(*,*) coox, cooy,ZBEAM

*      xpp = 0.001 * xpp
*      ypp = 0.001 * ypp
*      cosx = xpp / sqrt(xpp**2+ypp**2+1.0D0)
*      cosy = ypp / sqrt(xpp**2+ypp**2+1.0D0)
      cox = 0.0d0
      coy = 0.0d0

* Cosines (tx,ty,tz)
*      TXFLK  (NPFLKA) = UBEAM
*      TYFLK  (NPFLKA) = VBEAM
*      TZFLK  (NPFLKA) = WBEAM
      TXFLK  (NPFLKA) = cox
      TYFLK  (NPFLKA) = coy
      TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
     &                       - TYFLK (NPFLKA)**2 )

* Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER

* Particle coordinates
*      XFLK   (NPFLKA) = XBEAM
*      YFLK   (NPFLKA) = YBEAM
*      ZFLK   (NPFLKA) = ZBEAM

      XFLK   (NPFLKA) = coox
      YFLK   (NPFLKA) = cooy
      ZFLK   (NPFLKA) = ZBEAM





*  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






*************************************************************************
* Give a cosines of x, y
      subroutine coxy(xdiv, ydiv, xpp, ypp, cosx, cosy)

* cosx and cosy are the cosines of x and y direction
* xp and yp are the angle divergences, xp = Px/Pz, yp =Py/Pz
* xdiv is FWHM of gaussian distribution of xp(unit: mrad)
* ydiv is FWHM of gaussian distribution of yp(unit: mrad)
c      parameter(PIPIPI = 3.141592653589793238D0)
c      parameter(TWOPIP = 2 * PIPIPI)

      double precision xdiv, ydiv, xp, yp, cosx, cosy
      double precision RGAUSx, RGAUSy,xpp,ypp

*      WRITE(*,*)xdiv,ydiv
      CALL FLNRR2(RGAUSx,RGAUSy)
* 0.001 : transform mrad to rad

       xp = 0.001 * xdiv * RGAUSx / 2.355D0
* divergent beam
       yp = 0.001 * ydiv * RGAUSy / 2.355D0
* focusing beam
*       yp = - 0.001 * ydiv * RGAUSy / 2.355D0

      cosx = xp / sqrt(xp**2+yp**2+1.0D0)
      cosy = yp / sqrt(xp**2+yp**2+1.0D0)

      xpp = xp * 1000.0
      ypp = yp * 1000.0
      return

      end





*************************************************************************
* A 2D-gaussian annular distribution with inter and outer radii of a and b
      subroutine Cooxy(fwhmx,fwhmy,a,b,coox,cooy)
* note: a<b
* a, b and inner and outer cut-off position in polar coordinate system
* fwhmx,fwhmy are FWHMs of gaussian distribution of x and y(unit: cm)
* sx,sy are the sigma of gaussian distribution
* Xmin = a * sx
* Xmax = b * sx
* Ymin = a * sy
* Ymax = b * sy


      PARAMETER(PI = 3.1415926535897932384626433832795D0)
      double precision a, b, fwhmx,sx,fwhmy,sy,r,the
      double precision XRAND1,XRAND2, coox, cooy


      sx = fwhmx / 2.355D0
      sy = fwhmy / 2.355D0

      XRAND1 = FLRNDM(XDUMMY1)

      r = DSQRT(-2.0*DLOG((1.0D0 - XRAND1)*DEXP(-a**2/2.)
     &     + XRAND1*DEXP(-b**2/2.)))


      XRAND2 = FLRNDM(XDUMMY2)
      the = 2.0 * PI * XRAND2


      coox = r * sx * DCOS(THE)
      cooy = r * sy * DSIN(THE)

      return

      end



*************************************************************************
* A 2D-gaussian annular distribution with inter and outer radii of a and b
      subroutine cxyp(fwhmx,fwhmxp,fwhmy,fwhmyp,a,b,coox,cooy,xpp,ypp)
* note: a<b
* a, b and inner and outer cut-off position in polar coordinate system
* fwhmx,fwhmy are FWHMs of gaussian distribution of x and y(unit: cm)
* sx,sy are the sigma of gaussian distribution
* Xmin = a * sx
* Xmax = b * sx
* Ymin = a * sy
* Ymax = b * sy


      PARAMETER(PI = 3.1415926535897932384626433832795D0)
      double precision a, b, fwhmx,fwhmxp,sx,fwhmy,fwhmyp,sy,r,the
      double precision XRAND1,XRAND2, coox, cooy,xp,yp,xpp,ypp


      sx = fwhmx / 2.355D0
      sy = fwhmy / 2.355D0
      xp = fwhmxp /  2.355D0
      yp = fwhmyp /  2.355D0

      XRAND1 = FLRNDM(DUMMY1)

      r = DSQRT(-2.0*DLOG((1.0D0 - XRAND1)*DEXP(-a**2/2.)
     &     + XRAND1*DEXP(-b**2/2.)))


      XRAND2 = FLRNDM(DUMMY2)
      the = 2.0 * PI * XRAND2


      coox = r * sx * DCOS(THE)
      xpp = r * xp * DSIN(THE)
      cooy = r * sy * DSIN(THE)
      ypp = r * yp * DCOS(THE)

      return

      end








