*$ 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.            *
*                                                                      *
*  K. Ott HZB/BESSY                                                    * 
*  12/2016 Neutron spectrum behind 60 cm concrete sampled isotropic    *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      INCLUDE '(EMFSTK)'
      INCLUDE '(CASLIM)'
      INCLUDE '(SCOHLP)'
*
      DIMENSION EBIN(1:268), SYCUM(1:268)  

      LOGICAL LFIRST
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*------------------------------------------------------------------
*    Fluka low energy bins (255) + 12 high energy bins
      DATA EBIN /9.209217E-13,
     &1.349028E-12, 1.976147E-12,2.894792E-12,4.240485E-12,6.211746E-12,
     &9.099382E-12, 1.332938E-11, 1.952578E-11,2.860266E-11,4.18991E-11,
     &6.13766E-11, 8.990856E-11, 1.317041E-10, 1.92929E-10,2.826153E-10,
     &4.139938E-10, 5.315785E-10,6.250621E-10,6.825603E-10,8.336811E-10,
     &8.764248E-10, 1.125352E-09, 1.44498E-09, 1.855391E-09,2.38237E-09,
     &3.059023E-09, 3.927864E-09,5.043477E-09,6.475952E-09,8.315288E-09,
     &1.067704E-08, 1.370959E-08, 1.760346E-08,2.260329E-08,2.90232E-08,
     &3.726653E-08, 4.785117E-08,6.144212E-08,7.889325E-08,1.013009E-07,
     &1.30073E-07, 1.67017E-07, 2.144541E-07, 2.753645E-07, 3.53575E-07,
     &4.539993E-07, 5.829466E-07,7.485183E-07,9.611165E-07,1.234098E-06,
     &1.363889E-06, 1.507331E-06,1.584613E-06,1.665858E-06,1.841058E-06,
     &2.034684E-06, 2.248673E-06,2.485168E-06,2.612586E-06,2.746536E-06,
     &2.863488E-06, 3.035391E-06,3.354626E-06,3.707435E-06,4.09735E-06,
     &4.307425E-06, 4.528272E-06,5.004514E-06,5.530844E-06,6.267267E-06,
     &7.101744E-06, 8.04733E-06, 9.11882E-06, 1.033298E-05, 1.17088E-05,
     &1.32678E-05, 1.503439E-05, 1.70362E-05, 1.930454E-05,2.133482E-05,
     &2.187491E-05, 2.357862E-05,2.417552E-05,2.478752E-05,2.605841E-05,
     &2.808794E-05, 3.182781E-05,3.430669E-05,3.517517E-05,3.606563E-05,
     &4.086771E-05, 4.630919E-05,5.247518E-05,5.656217E-05,5.946217E-05,
     &6.251086E-05, 6.737947E-05,7.635094E-05,8.651695E-05,9.803654E-05,
     &0.00011109, 0.0001167857, 0.0001227734, 0.0001290681,0.0001356856,
     &0.0001426423, 0.0001499558,0.0001576442,0.0001616349,0.0001657268,
     &0.0001699221, 0.0001742237, 0.0001831564,0.000192547,0.0002024191,
     &0.0002127974, 0.0002237077,0.0002351775,0.0002472353,0.0002599113,
     &0.0002732372,0.0002801543, 0.0002872464,0.0002945181,0.0003019738,
     &0.0003096183, 0.0003174564,0.0003337327,0.0003508435,0.0003688317,
     &0.0003877421, 0.000407622, 0.0004285213, 0.000450492,0.0004735892,
     &0.0004978707, 0.0005104743,0.0005233971,0.0005366469,0.0005502322,
     &0.0005784433, 0.0006081006,0.0006392786,0.0006720551,0.0007065121,
     &0.0007427358, 0.0007808167, 0.00082085, 0.0008629358,0.0009071795,
     &0.0009536916, 0.0009616402, 0.0009778343, 0.001002588,0.001053992,
     &  0.001108032, 0.001164842, 0.00119433, 0.001224564, 0.001287349,
     &  0.001353353, 0.001422741, 0.001495686, 0.00153355, 0.001572372,
     &  0.001612176, 0.001652989, 0.001737739, 0.001826835, 0.001873082,
     &  0.001920499, 0.001969117, 0.002018965, 0.00212248, 0.002231302,
     &  0.002268877, 0.002306855, 0.002345703, 0.002365253, 0.002385205,
     &  0.00242513, 0.00246597, 0.002592403, 0.002725318, 0.002865048,
     &  0.003011942, 0.00308819, 0.003166368, 0.003246525, 0.003328711,
     &  0.003499377, 0.003678794, 0.00386741, 0.004065697, 0.004274149,
     &  0.00449329, 0.004607038, 0.004723666, 0.004843246, 0.004965853,
     &  0.005220458, 0.005488116, 0.005769498, 0.005915554, 0.006065307,
     &  0.006218851, 0.006376282, 0.006537698, 0.006592384, 0.006647595,
     &  0.0067032, 0.006872893, 0.007046881, 0.007225274, 0.007408182,
     &  0.007595721, 0.007788008, 0.007985162, 0.008187308, 0.00839457,
     &  0.00860708, 0.008824969, 0.009048374, 0.009277435, 0.009512294,
     &  0.009753099, 0.01, 0.01025315, 0.01051271, 0.01077884,
     &  0.01105171, 0.01133148, 0.01161834, 0.01191246, 0.01221403,
     &  0.01252323, 0.01284025, 0.01316531, 0.01349859, 0.01384031,
     &  0.01419068, 0.01454991, 0.01491825, 0.0152959, 0.01568312,
     &  0.01608014, 0.01648721, 0.01690459, 0.01733253, 0.01777131,
     &  0.01822119, 0.01868246, 0.01915541, 0.01964033, 0.02,
     &  0.026154, 0.03420159, 0.04472542, 0.05848743, 0.07648402,
     &  0.1000182, 0.1307937, 0.171039, 0.2236677, 0.2924902,
     &  0.3824895, 0.5001815/

c     neutron spectrum  group fluence integration normalized
      DATA SYCUM /
     &  0, 0.0002054442, 0.0007492366, 0.001877797, 0.004182149,
     &  0.008896646, 0.01837331, 0.03785863, 0.07518708, 0.13959,
     &  0.2381733, 0.361028, 0.4736936, 0.5451429, 0.5778204,
     &  0.5910688, 0.5987512, 0.6033961, 0.6060977, 0.6076517,
     &  0.6109545, 0.6116745, 0.615573, 0.6194674, 0.6232362,
     &  0.6269655, 0.6307192, 0.6339805, 0.6374515, 0.6405156,
     &  0.6438449, 0.646897, 0.6497577, 0.6526895, 0.6555811,
     &  0.6582134, 0.6611298, 0.6641868, 0.6666766, 0.6694534,
     &  0.6723662, 0.6749117, 0.677416, 0.6800069, 0.6825005,
     &  0.6848962, 0.6872925, 0.6895802, 0.6918508, 0.694188,
     &  0.6966859, 0.6976558, 0.6985064, 0.6989082, 0.6993794,
     &  0.7003585, 0.7011939, 0.7021133, 0.7030997, 0.7035074,
     &  0.7038633, 0.7041231, 0.704504, 0.7053688, 0.706228,
     &  0.707042, 0.7074749, 0.7079418, 0.7089226, 0.7097114,
     &  0.7108208, 0.7119378, 0.7129914, 0.7141641, 0.7153323,
     &  0.7165159, 0.7176279, 0.7189264, 0.7200181, 0.7213412,
     &  0.7223268, 0.7225534, 0.7232957, 0.7235256, 0.7238753,
     &  0.7244125, 0.7252375, 0.7265375, 0.7272242, 0.7275097,
     &  0.7277083, 0.728978, 0.7303386, 0.7317606, 0.7324684,
     &  0.7329905, 0.7336013, 0.7344632, 0.7360867, 0.7373959,
     &  0.7389637, 0.7406399, 0.7413203, 0.7420268, 0.7427698,
     &  0.7434288, 0.7442966, 0.745124, 0.7459504, 0.7464689,
     &  0.7469049, 0.7473832, 0.7477822, 0.7484213, 0.7489113,
     &  0.7495032, 0.7500719, 0.7509061, 0.7517409, 0.7524601,
     &  0.7531787, 0.7541249, 0.7544803, 0.7549784, 0.7554789,
     &  0.7559372, 0.7564253, 0.7569811, 0.7578993, 0.7587259,
     &  0.7594662, 0.7601126, 0.7608918, 0.7614432, 0.7619361,
     &  0.7627892, 0.7640328, 0.7649488, 0.765744, 0.766584,
     &  0.7674845, 0.7693176, 0.7709838, 0.7728369, 0.7745794,
     &  0.7765322, 0.7785367, 0.7801929, 0.7819343, 0.7834524,
     &  0.7847799, 0.7857145, 0.7858581, 0.7860931, 0.7864107,
     &  0.7872248, 0.7884941, 0.7902746, 0.7914137, 0.7923765,
     &  0.7942837, 0.795822, 0.7982314, 0.8005607, 0.8016441,
     &  0.8027645, 0.803822, 0.8047753, 0.8071907, 0.8100837,
     &  0.8108844, 0.811816, 0.8128696, 0.8143191, 0.8173391,
     &  0.8211097, 0.8229166, 0.8249187, 0.8276541, 0.8293448,
     &  0.8304485, 0.8321946, 0.8333867, 0.8362051, 0.8390039,
     &  0.8414531, 0.8436366, 0.8446806, 0.845525, 0.8462254,
     &  0.8467203, 0.8477207, 0.8487785, 0.8499132, 0.8513578,
     &  0.8531026, 0.8544441, 0.8552392, 0.8561245, 0.8569959,
     &  0.8576906, 0.859157, 0.8609494, 0.8623035, 0.8629197,
     &  0.8635945, 0.8643816, 0.8650856, 0.8657741, 0.8660587,
     &  0.8663113, 0.8665292, 0.8670903, 0.8678971, 0.8684514,
     &  0.8689411, 0.8695146, 0.8700226, 0.8704776, 0.8710619,
     &  0.8715048, 0.8720357, 0.8725495, 0.8730915, 0.8736187,
     &  0.8741468, 0.8747215, 0.8752701, 0.8757725, 0.8762851,
     &  0.8768784, 0.8774869, 0.8779664, 0.8784316, 0.8789135,
     &  0.8794578, 0.8799824, 0.880528, 0.8810723, 0.881658,
     &  0.8822545, 0.8828247, 0.8834523, 0.8840625, 0.8846588,
     &  0.8853052, 0.8859331, 0.886595, 0.8872942, 0.8879929,
     &  0.8886721, 0.8893826, 0.890098, 0.8907822, 0.8914565,
     &  0.8919426, 0.9003184, 0.9105521, 0.9234138, 0.9391529,
     &  0.9564459, 0.9722523, 0.9845492, 0.9929516, 0.9974921,
     &  0.9992884, 0.9998735, 1.0/
*----------------------------------------------------------------------*
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
*======================================================================*
*                                                                      *
*                 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(lunout,*)
         write(lunout,'(a,132a)') ("*",i=1,132)
         write(lunout,*)
         write(lunout,*)
         write(lunout,'(a)') 'Neutron 60 cm source, isotropic'
         write(6,'(a)') 'Neutron 60 cm source, isotropic'
         write(lunout,*)
         write(lunout,*)
         write(lunout,'(a,132a)') ("*",i=1,132)
         write(lunout,*)
      END IF
*  |
*  +-------------------------------------------------------------------*
*     Sample the energy group
      XI = FLRNDM(DUMMY)
      DO 500 K = 1, 268
         IF(XI .LE. SYCUM(K)) THEN
            ENERGY=EBIN(K)
            GO TO 600
         END IF
 500  CONTINUE
      STOP ' Failed to sample the energy group'
 600  CONTINUE
      kount=kount+1 
*  +-------------------------------------------------------------------*
* Isotropic distribution
* Sample cosine of polar angle
      COSTHE = TWOTWO * FLRNDM(DUMMY) - ONEONE
      SINTHE = SQRT(ONEONE - COSTHE**2)
* Sample azimuthal angle 
      XI = TWOPIP * FLRNDM(DUMMY)
      COSPHI = COS(XI)
      SINPHI = SIN(XI)
*  +-------------------------------------------------------------------*
*  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)
      TKEFLK (NPFLKA) = ENERGY
*      TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IJBEAM)**2 ) - AM (IJBEAM)
* Particle momentum
*      PMOFLK (NPFLKA) = PBEAM
      PMOFLK (NPFLKA) = SQRT(ENERGY * (ENERGY + TWOTWO * AM(IJBEAM)))
* Cosines (tx,ty,tz) distribution isotropic
      TXFLK (NPFLKA) = SINTHE * COSPHI
      TYFLK (NPFLKA) = SINTHE * SINPHI
      TZFLK (NPFLKA) = COSTHE
c--------------------------------------
* Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates
      XFLK   (NPFLKA) = XBEAM
      YFLK   (NPFLKA) = YBEAM
      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

