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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2020      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  03-Jan-20    by             Paola Sala           *
*                                                Infn - Milan          *
*                                                                      *
*  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)'
*
      LOGICAL LFIRST, LISNUT
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*
* ===========================     my definitions  .....
*

**      PARAMETER (NMAX=70000)           ! - size of my input file...
**      PARAMETER (NMAX=3407000)           ! - for C1  or  C2 file
      PARAMETER (NMAX=82000)           ! - for C1  or  C2 file
*      PARAMETER (NMAX=6813000)           ! - for BIG-input, total

*      COMMON /MYCOM/ NNN, IDET(NMAX), MYIJ(NMAX), ERG(NMAX), WGT(NMAX), 
*     &                    TT(NMAX), XXX(NMAX), YYY(NMAX), ZZZ(NMAX),
*     &                    COSX(NMAX), COSY(NMAX), COSZ(NMAX)
      COMMON /MYCOM/ ERG(NMAX), WGT(NMAX), TT(NMAX), 
     &               XXX(NMAX), YYY(NMAX), ZZZ(NMAX),
     &               VCOSX(NMAX), VCOSY(NMAX), VCOSZ(NMAX),
     &               IDET(NMAX), MYIJ(NMAX), NNN

**
**  NNN - line number ...
**
**  IDET - detector : 5 (C1), 6 (C2)
**  MYIJ   - part. type                 - replace  IJBEAM here
**  ERG - energy, GeV 
**  WGT - ev. weight
**  TT  - ev. time, sec
**  XXX, YYY, ZZZ - coordinates,  cm
*** COSX, COSY, COSZ - normalized angles (cos)
**

*      SAVE XXX, YYY, ZZZ
*      SAVE COSX, COSY, COSZ
*      SAVE MYIJ, ERG, WGT, TT

      SAVE /MYCOM/

      INTEGER NENTRY                   ! counter
      SAVE    NENTRY

      INTEGER MYLUN           ! channel to print
 
***      INTEGER DUM, I

**      double precision SFLOODR  ! radius of the sphere (sflood)
**      double precision ZMIN  ! m
**      integer N0       ! number of trials to call sflood for a given event
**      integer IETYPE, FLXGEN


* ============================     my  definitions - end .....

*
*  Statement function:
      LISNUT (IJ) = INDEX ( PRNAME (IJ), 'NEUTRI' ) .GT. 0
*======================================================================*
*                                                                      *
*                 BASIC VERSION            ...and my  inouts           *
*                                                                      *
*======================================================================*
*
* it is the code  to read in my input data for neutrons and photons
* my lines are commented with -     ! my...
*
* NMAX - it corresponds to 2e+06 beam partciles (P)
*         OR               4e+06
*
*       see ==>       SUBROUTINE READ_FILE 
*
*  Det - Type -  Nev - Ekin - Ltot- L.Orig - Time - w - X - Y - Z  - 
*                                                      cosX - cosY - cosZ 
*  NEWREG, JTRACK, NCASE, Ekin,  CTRACK, CMTRCK, ATRACK, - 
*                       WTRACK, XSCO, YSCO, ZSCO   CXTRCK, CYTRCK, CZTRCK 
*
* det = 5  for C1
*     = 6  for C2
* typ = particle type
* nev  = event number
* ekin = particle kinetic energy , GeV
* ltot, lorig - 2 type of distancies,  for C1 
*               more natural is Lorig
* time - time in sec
* ww   = particle weight
* xx-yy-zz - coordinates
*  cosx-cosy-cosz - angles ...

*======================================================================
*
*  use new format (corrected in Analiz-C1/
*            idet, itype, nev, ekin, ltot, lorig, told, time, ww,
*                                  xx, yy, zz, xcos, ycos, zcos);
*
*======================================================================


      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:

      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
*
*                    .... my start - read file ...
* --------------------------------------------------------------

      MYLUN= 90


*      OPEN ( UNIT = 90, FILE = 'my-90infor', STATUS = 'NEW', FORM =
       OPEN ( UNIT = MYLUN, FILE = 'my-infor-45', STATUS = 'NEW', FORM =
     &          'FORMATTED' )

       write (MYLUN,*) "...my file start.. : LUNOUT= , IODRAW= :", 
     &                 LUNOUT,IODRAW
       WRITE(MYLUN,'(A32)') RUNTIM
       write(MYLUN,*) 
*       ------------------

         CALL READ_FILE          ! --- read file and fill common block

         NENTRY = 0

*                        ---------------------------
*                    .... my start - end of read file ...
      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

*  +-------------------------------------------------------------------*
*
*            .... my part - select a random ev-line in my arrays ...
*                           to organyze FLUKA cycles...
*                          write all selected numbers 
*                      -------------------------
      NENTRY = NENTRY + 1

      RNDSIG = FLRNDM (RNDSIG)              ! FLUKA rndm

      N = INT(NNN * RNDSIG) + 1             ! select line from array

      IF ( NENTRY .LT. 100) THEN            ! type some lines
*       WRITE (LUNOUT,*) '*** N=',N, " ", IJ(N)," ", ERG(N)," ", WGT(N)
*        WRITE (LUNOUT,*) '*** N=',N,  MYIJ(N), ERG(N), XXX(N), COSX(N)
*        WRITE (45,*) '*** N=',N,  MYIJ(N), ERG(N), XXX(N), COSX(N)
        WRITE (MYLUN,*) '*** N=',N,  MYIJ(N), ERG(N), XXX(N), VCOSX(N)
      END IF 

* ************************
*        convert my variables into local ones
*           --------------

      IJBEAM = MYIJ(N)           ! part. type

      UBEAM = VCOSX(N)            !  vector of the particles Cosines (tx,ty,tz)
      VBEAM = VCOSY(N)
      WBEAM = VCOSZ(N)

      XBEAM = XXX(N)
      YBEAM = YYY(N)
      ZBEAM = ZZZ(N)

* *****************************

*
*            .... my part - end of the part ...
*
*  +-------------------------------------------------------------------*


      NPFLKA = NPFLKA + 1
*  Wt is the weight of the particle
***      WTFLK  (NPFLKA) = ONEONE
      WTFLK  (NPFLKA) = WGT(N)                ! my - put my weight,, not 1

      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 = IPROM  * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA
         IJHION = IJHION * 100    + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         EEXSTK (NPFLKA) = EXENRG (IONID)
         TMNSTK (NPFLKA) = TMNLF  (IONID)
         ILVSTK (NPFLKA) = IEXLVL (IONID)
         LFRPHN (NPFLKA) = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:

      ELSE IF ( IJBEAM .EQ. -2 ) THEN

         IJHION = IPROM  * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA
         IJHION = IJHION * 100    + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         EEXSTK (NPFLKA) = EXENRG (IONID)
         TMNSTK (NPFLKA) = TMNLF  (IONID)
         ILVSTK (NPFLKA) = IEXLVL (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
*  |  Parent radioactive isotope:
         IRDAZM (NPFLKA) = 0
*  |  Particle age (s)
         AGESTK (NPFLKA) = +ZERZER
*  |  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) ) )
         LFRPHN (NPFLKA) = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:        -  ! it is my case : N  and G ....

      ELSE
         IONID = IJBEAM

         EEXSTK (NPFLKA) = EXENRG (IJBEAM)
         TMNSTK (NPFLKA) = TMNLF  (IJBEAM)
         ILVSTK (NPFLKA) = IEXLVL (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
*  |  Parent radioactive isotope:
         IRDAZM (NPFLKA) = 0

*  |  Particle age (s)
***         AGESTK (NPFLKA) = +ZERZER
         AGESTK (NPFLKA) = TT(N)             ! my : time-history > 0  !

*  |  Kinetic energy of the particle (GeV)
***         TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 )
***     &                   - AM (IONID)
      TKEFLK (NPFLKA) = ERG(N)              ! my :  particle energy...

*  |  Particle momentum
***         PMOFLK (NPFLKA) = PBEAM
         PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (IONID) ) )

*        PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
*    &                          + TWOTWO * AM (IONID) ) )

*  |  +----------------------------------------------------------------*
*  |  |  Check if it is a neutrino, if so force the interaction
*  |  |  (unless the relevant flag has been disabled)
         IF ( LISNUT (IJBEAM) .AND. LNUFIN ) THEN
            LFRPHN (NPFLKA) = .TRUE.
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |  Not a neutrino
         ELSE
            LFRPHN (NPFLKA) = .FALSE.
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
      END IF

*  |
*  +-------------------------------------------------------------------*
*  From this point .....
*  Particle generation (1 for primaries)
      LOFLK  (NPFLKA) = 1
*  User dependent flag:
      LOUSE  (NPFLKA) = 0
*  No channeling:
      KCHFLK (NPFLKA) = 0
      ECRFLK (NPFLKA) = ZERZER
*  Extra infos:
      INFSTK (NPFLKA) = 0
      LNFSTK (NPFLKA) = 0
      ANFSTK (NPFLKA) = ZERZER
*  Parent variables:
      IPRSTK (NPFLKA) = 0
      EKPSTK (NPFLKA) = ZERZER
*  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
*  +-------------------------------------------------------------------*

      AKNSHR (NPFLKA) = -TWOTWO

*  Cosines (tx,ty,tz)
      TXFLK  (NPFLKA) = UBEAM
      TYFLK  (NPFLKA) = VBEAM
      TZFLK  (NPFLKA) = WBEAM
*     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

*  Calculate the total kinetic energy of the primaries: don't change
*  +-------------------------------------------------------------------*
*  |  (Radioactive) isotope:
      IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      ELSE IF ( ILOFLK (NPFLKA) .EQ. -2 .OR.
     &          ILOFLK (NPFLKA) .GT. 100000 ) THEN
         TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
*  |
*  +-------------------------------------------------------------------*
*  |  Standard particle:
      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 =========================================*
*                  ... my : end of job .....
*** 9999 CONTINUE
***      CALL FLABRT('SOURCE','Error reading source file')
*                  ..........................
      END




**************************************************
*=======================================================


      SUBROUTINE READ_FILE 
*               ---------
*       use my common and fill it with input data
*

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*

**      PARAMETER (NMAX=70000)           ! - size of my input file...
**      PARAMETER (NMAX=3407000)           ! - for C1  or  C2 file
      PARAMETER (NMAX=82000)           ! - for C1  or  C2 file
*      PARAMETER (NMAX=6813000)           ! - for BIG-input, total

*      COMMON /MYCOM/ NNN, IDET(NMAX), IJ(NMAX), ERG(NMAX), WGT(NMAX), 
*     &                    TT(NMAX), XXX(NMAX), YYY(NMAX), ZZZ(NMAX),
*     &                    COSX(NMAX), COSY(NMAX), COSZ(NMAX)
      COMMON /MYCOM/ ERG(NMAX), WGT(NMAX), TT(NMAX), 
     &               XXX(NMAX), YYY(NMAX), ZZZ(NMAX),
     &               VCOSX(NMAX), VCOSY(NMAX), VCOSZ(NMAX),
     &               IDET(NMAX), MYIJ(NMAX), NNN

**
**  NNN - line number ...
**
**  IDET - detector : 5 (C1), 6 (C2)
**  MYIJ   - part. type                - replace  IJBEAM here
**  ERG - energy, GeV 
**  WGT - ev. weight
**  TT  - ev. time, sec
**  XXX, YYY, ZZZ - coordinates,  cm
*** COSX, COSY, COSZ - normalized angles (cos)
**

*      SAVE XXX, YYY, ZZZ
*      SAVE COSX, COSY, COSZ
*      SAVE MYIJ, ERG, WGT, TT

      SAVE /MYCOM/

      CHARACTER*250 LINE

*      LOGICAL FIRST
**
*      SAVE FIRST
*      DATA FIRST / .TRUE. /

       INTEGER MYLUN

*======================================================================
*   information about input file:
**************************************
* NMAX - it corresponds to 2e+06 beam partciles (P)
*            or 4e+06 for BIG-input  file
*                         ----------
*  Det - Type -  Nev - Ekin - Ltot- L.Orig - Time - w - X - Y - Z  - 
*                                                      cosX - cosY - cosZ 
*  NEWREG, JTRACK, NCASE, Ekin,  CTRACK, CMTRCK, ATRACK, - 
*                       WTRACK, XSCO, YSCO, ZSCO   CXTRCK, CYTRCK, CZTRCK 
*
* det = 5  for C1
*     = 6  for C2
* typ = particle type
* nev  = event number
* ekin = particle kinetic energy , GeV
* ltot, lorig - 2 type of distancies,  for C1 
*               more natural is Lorig
* time - time in sec
* ww   = particle weight
* xx-yy-zz - coordinates
*  cosx-cosy-cosz - angles ...
*
*======================================================================
*
*  use new format (corrected in Analiz-C1/
*            idet, itype, nev, ekin, ltot, lorig, told, time, ww,
*                                  xx, yy, zz, xcos, ycos, zcos);
*
*
*======================================================================

*
*                    .... my start - read file ...
*                      -------------------------
****  |  WHASOU(1) = UNIT number (File name is linked with OPEN card)
***         LUNRD = NINT(WHASOU(1))
****         POSMIN = WHASOU(2)
****         POSMAX = WHASOU(3)

*-----------------------------
**          
**            WRITE(LUNOUT,*) "-my- ",LUNRD
**
***            LUNOUT= 91    ! ...see parameter
**
***            OPEN(UNIT = 81, FILE = "../sr.spectrum", 
**            OPEN(LUNRD, FILE = "test.dat",   IOSTAT= IOS,
**     &         STATUS = "OLD")
**
**            WRITE(LUNOUT,*) "-my- open: iostat = ",IOS
**
***            WRITE(LUNOUT,'(A,A,A)') "    sr white spectrum ",
***     &      "(file sr.spectrum)" 
**
**
**
********************************************

       MYLUN = 90

*       OPEN ( UNIT = MYLUN, FILE = 'my-infor-45',STATUS = 'UNKNOWN',FORM =
       OPEN ( UNIT = MYLUN, FILE = 'my-infor-45',STATUS = 'NEW',FORM =
     &          'FORMATTED' )

         WRITE (MYLUN,*) "...my file start.. : LUNOUT= , IODRAW= :", 
     &                 LUNOUT,IODRAW
         WRITE(MYLUN,'(A32)') RUNTIM
         write(MYLUN,*) 


*----------------------------
*    another method to open file 
*

        WRITE(MYLUN,*)         !  put empty line in *out  file

        CALL OAUXFI('../../Analiz-C1/BIG-10-out', LUNRD, 'OLD', IERR)
*        CALL OAUXFI('../../Analiz-C1/BIG-30-out', LUNRD, 'OLD', IERR)

        WRITE(MYLUN,*)
        WRITE(MYLUN,*) " -my- : ierr= ", IERR
        WRITE(MYLUN,'(A,132A)') ("*",I=1,132)   ! line with **/...
        WRITE(MYLUN,*)

*----------------------------

*         WTOT = ZERZER               ! internal  - total weight
         NNN = 0                     ! internal ... counter

 10      CONTINUE                   ! start loop
************************** 

*            READ( LUNRD, '(A)', ERR=9999, END=20 ) LINE
            READ( LUNRD, '(A)', ERR=9999, END=20, IOSTAT= IOS) LINE

            IF (NNN .LT. 10) THEN
              WRITE(MYLUN,*) "-my- read : ios =",IOS
              WRITE(MYLUN,'(A)') LINE
***            WRITE(LUNOUT,*) "***", LINE
            ENDIF


***  now read line ...
***  ****************************
*            READ (LINE,*,ERR=10) 
*            READ (LINE,*, IOSTAT= IOS,ERR=10) 
*     &          IDD, IP, NN, EK, DUM, DUM, T W, X, Y, Z, XC, YC, ZC
*     &         I, IDD, DUM, DUM, IP, EK, AL, T, W, X, Y, Z, XC, YC, ZC
*            WRITE(LUNOUT,*) "***", IOS
*   **********************************8

***  new format !
***             idet, itype, nev, ekin, ltot, lorig, told, time, ww,
***                   xx, yy, zz, xcos, ycos, zcos
***  Told - FULKA time,  Tnew = + Trndm for bunch
***  X = is corrected == 6.6 (not 6.7)
***  ****************************
*            READ (LINE,*,ERR=10) 
            READ (LINE,*, IOSTAT= IOS,ERR=10) 
     &    IDD, IP, NN, EK, DUM, DUM, Told,Tnew, W, X, Y, Z, XC, YC, ZC

*            WRITE(LUNOUT,*) "***", IOS
*   **********************************8

*           ... now save good event ....
*
            NNN = NNN + 1

            IF (NNN.GT.NMAX) CALL FLABRT('SOURCE','Increase NMAX')

            IDET(NNN) = IDD              ! 5 (C1), 6 (C2)
            MYIJ(NNN)  = IP                ! part. type
            ERG(NNN) = EK                !  kinetic energy, GeV  
            WGT(NNN) = W                 ! weight

*            TT(NNN)  = T                 ! time, sec 
            TT(NNN)  = Told                 ! time, sec - FLUKA time 
*            TT(NNN)  = Tnew                 ! time, sec - with Trndm

            XXX(NNN) = X                 ! coordinates
            YYY(NNN) = Y
            ZZZ(NNN) = Z

*           ... make normalization of forward direction vector:
*               input angles are :   YC  and ZC,,,, XC is calculated !
*
            VCOSY(NNN) = YC
            VCOSZ(NNN) = ZC
            VCOSX(NNN) = SQRT ( ONEONE - YC**2 - ZC**2)


            IF (NNN .LT. 10) THEN
          WRITE(MYLUN,* ) "-N- :", NNN, MYIJ(NNN), ERG(NNN), TT(NNN), 
     &                    XXX(NNN), YYY(NNN), ZZZ(NNN), VCOSX(NNN)
            ENDIF


         GOTO 10                    ! read next ev.
*********************************

 20      CONTINUE                   ! end of ev. loop

         IF (NNN.EQ.0) CALL FLABRT('SOURCE','Error reading file')
         WRITE (MYLUN,*)
         WRITE (MYLUN,*) '*** input file:  ',NNN,' particles loaded',
     &' - total weight: ',WTOT
         WRITE (MYLUN,*)

*                        ---------------------------
*                    .... my start - end of read file ...
*
      RETURN

*=== End of subroutine READ_FILE =========================================*
*                  ... my : end of job .....
 9999 CONTINUE
      CALL FLABRT('SOURCE-READ_FILE','Error reading source file')
*                  ..........................


      END
*===================================================


