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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*----------------------------------------------------------------------*
*
c      DIMENSION CFCUM(0:72), ENEDGE(73), ANG(0:40), CUMANG(0:40)   
      DIMENSION CFCUM(0:53), ENEDGE(53), ANG(0:40), CUMANG(0:40)   
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(CASLIM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
      LOGICAL LFIRST
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*----------------------------------------------------------------------*
*   Neutron energy PTB
      DATA ENEDGE /
     & 4.210000000E-05,1.100000000E-04,3.300000000E-04,5.400000000E-04,
     & 7.500000000E-04,9.700000000E-04,1.180000000E-03,1.400000000E-03,
     & 1.610000000E-03,1.820000000E-03,2.040000000E-03,2.250000000E-03,
     & 2.470000000E-03,2.680000000E-03,2.900000000E-03,3.110000000E-03,
     & 3.320000000E-03,3.540000000E-03,3.750000000E-03,3.970000000E-03,
     & 4.180000000E-03,4.390000000E-03,4.610000000E-03,4.820000000E-03,
     & 5.040000000E-03,5.250000000E-03,5.470000000E-03,5.680000000E-03,
     & 5.890000000E-03,6.110000000E-03,6.320000000E-03,6.540000000E-03,
     & 6.750000000E-03,6.960000000E-03,7.180000000E-03,7.390000000E-03,
     & 7.610000000E-03,7.820000000E-03,8.030000000E-03,8.250000000E-03,
     & 8.460000000E-03,8.680000000E-03,8.890000000E-03,9.110000000E-03,
     & 9.320000000E-03,9.530000000E-03,9.750000000E-03,9.960000000E-03,
     & 1.018000000E-02,1.039000000E-02,1.060000000E-02,1.082000000E-02,
     & 1.100000000E-02/
*
*   Cumulative Am-Be  spectrum PTB
      DATA CFCUM /0.0D0,
     & 0.000000000E+00,1.439707739E-02,4.779029857E-02,7.908394596E-02,
     & 1.071782428E-01,1.321731688E-01,1.535688255E-01,1.733648069E-01,
     & 1.908612552E-01,2.100573584E-01,2.323528324E-01,2.538484688E-01,
     & 2.763439022E-01,2.991392747E-01,3.286332874E-01,3.642260621E-01,
     & 4.011185729E-01,4.357115506E-01,4.664053197E-01,4.963992310E-01,
     & 5.232937714E-01,5.518879667E-01,5.836815127E-01,6.143752818E-01,
     & 6.476685233E-01,6.780623533E-01,7.054567923E-01,7.287520633E-01,
     & 7.493478824E-01,7.675441885E-01,7.852405962E-01,8.056364558E-01,
     & 8.239327417E-01,8.402294334E-01,8.569260440E-01,8.737226343E-01,
     & 8.925188187E-01,9.109150842E-01,9.278116542E-01,9.422087316E-01,
     & 9.518867670E-01,9.584054437E-01,9.626645791E-01,9.663338342E-01,
     & 9.701430610E-01,9.752020340E-01,9.814507655E-01,9.869696452E-01,
     & 9.916486953E-01,9.953479444E-01,9.981273801E-01,9.996370737E-01,
     & 1.000000000E+00/

*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
*   Angle cosines
c      DATA ANG/ -1.0D0, -0.95D0, -0.9D0, -0.85D0, -0.8D0, -0.75D0,
c     &          -0.7D0, -0.65D0, -0.6D0, -0.55D0, -0.5D0, -0.45D0,
c     &          -0.4D0, -0.35D0, -0.3D0, -0.25D0, -0.2D0, -0.15D0,
c     &          -0.1D0, -5.D-02,  0.0D0,  5.D-02,  0.1D0,  0.15D0,
c     &           0.2D0,  0.25D0,  0.3D0,  0.35D0,  0.4D0,  0.45D0,
c     &           0.5D0,  0.55D0,  0.6D0,  0.65D0,  0.7D0,  0.75D0,
c     &           0.8D0,  0.85D0,  0.9D0,  0.95D0,  1.0D0 /
*   Cumulative angular distribution
c      DATA CUMANG / 0.0D0, 4.3787904871D-4, 9.309108266D-4,
c     &    1.5080318037D-3, 2.1835838416D-3, 2.9743548826D-3,
c     &    3.8909576476D-3, 4.8923191479D-3, 5.9862770681D-3,
c     &    7.1813938468D-3, 8.4870236975D-3, 9.9133858209D-3,
c     &    1.1471644393D-2, 1.3097707545D-2, 1.4776914040D-2,
c     &    1.6511000723D-2, 1.8301761200D-2, 2.0151047694D-2,
c     &    2.2060772967D-2, 2.4032912292D-2, 2.6069505497D-2,
c     &    2.9084203728D-2, 3.3546756876D-2, 4.0152519402D-2,
c     &    4.9930800378D-2, 6.4405250301D-2, 8.5831276094D-2,
c     &    0.1175474784D0,  0.1569148680D0,  0.1981130258D0,
c     &    0.2412270912D0,  0.2863461629D0,  0.3335634830D0,
c     &    0.3829766299D0,  0.4346877199D0,  0.4919529577D0,
c     &    0.5559660349D0,  0.6275220820D0,  0.7075099237D0,
c     &    0.8189694928D0,  1.0D0 /
*----------------------------------------------------------------------*
*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
*----------------------------------------------------------------------*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
	 open(99,file='sourceout.txt',status='unknown')
      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)') "Am-Be Neutron Source. Version PTB" 
        write(lunout,*)
        write(lunout,*)
        write(lunout,'(a,132a)') ("*",i=1,132)
        write(lunout,*)
      END IF
*  |
*  +-------------------------------------------------------------------*
*  Sample the energy group
      XI = FLRNDM(DUMMY)
      write(99,*)XI
        DO 500 K = 1,53
        IF(XI .LE.CFCUM(K)) THEN
*          NEUGRP = 0
          ENERGY = ENEDGE(K-1)+
     &    (XI-CFCUM(K-1))*(ENEDGE(K)-ENEDGE(K-1))/(CFCUM(K)-CFCUM(K-1))
c      write(99,*)XI,K,ENEDGE(k),CFCUM(k)
c       write(99,*) energy
        GO TO 600          
        END IF
 500  CONTINUE
c      write(99,*)XI,K,ENEDGE(k),CFCUM(k)
      STOP ' Failed to sample the energy group'
 600  CONTINUE
      kount=kount+1
*
*  Sample the cosine of the polar angle
c      XI = FLRNDM(XI)
c      DO 300 K = 1, 40
c        IF(XI .LE. CUMANG(K)) THEN
c          COSTHE = ANG(K-1) + 
c     &    (XI - CUMANG(K-1))*(ANG(K)-ANG(K-1))/(CUMANG(K)-CUMANG(K-1))
c          GO TO 400          
c        END IF
c 300  CONTINUE
c      STOP ' Failed to sample the polar angle cosine'
c 400  CONTINUE
c      SINTHE = SQRT(ONEONE - COSTHE**2)
c      IF(FLRNDM(COSTHE) .LE. HLFHLF) SINTHE = -SINTHE
*  Samples the azimuthal angle
c     PHI = TWOPIP * FLRNDM(COSTHE)
c      COSPHI = COS(PHI)
*  Npflka is the stack counter: of course any time source is called it
*  must be =0      
      COSTHE = ONEONE - TWOTWO*FLRNDM(YDUMMY)
      SINTHE = SQRT(ONEONE - COSTHE**2)
      IF(FLRNDM(YDUMMY).LE. HLFHLF) SINTHE = -SINTHE
* Samples the azimuthal angle
      PHI = TWOPIP * FLRNDM(ZDUMMY)
      COSPHI = COS(PHI)
      SINPHI = SQRT(ONEONE - COSPHI**2)
c      if (PHI.gt.hlfhlf) sinphi=-sinphi
      UBEAM = SINTHE * COSPHI
      VBEAM = SINTHE * SINPHI
      WBEAM = COSTHE
      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
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IPROZ  * 1000 + IPROA
         IJHION = IJHION * 100 + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         ILOFLK (NPFLKA) = IJHION
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID  = IJBEAM
         ILOFLK (NPFLKA) = IJBEAM
      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
c      IGROUP (NPFLKA) = NEUGRP
* Kinetic energy of the particle (GeV)
      TKEFLK (NPFLKA) = ENERGY
* Particle momentum
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                  + TWOTWO * AM (ILOFLK(NPFLKA)) ) )
* Cosines (tx,ty,tz) (make sure they are properly normalized)
c      TZFLK (NPFLKA) = COSTHE
c      TYFLK (NPFLKA) = SINTHE * COSPHI
c      TXFLK (NPFLKA) = SQRT(ONEONE-TZFLK(NPFLKA)**2-TYFLK(NPFLKA)**2)
c      IF(FLRNDM(COSPHI) .LE. HLFHLF) TXFLK (NPFLKA) = -TXFLK (NPFLKA)
*Cosines (tx,ty,tz) (make sure they are properly normalized)
      TZFLK (NPFLKA) = WBEAM
      TYFLK (NPFLKA) = VBEAM
      TXFLK (NPFLKA) = UBEAM 

* 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
*  Flag this is prompt radiation
      LRADDC (NPFLKA) = .FALSE.
      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

