 PROGRAM Beam
 IMPLICIT NONE

 DOUBLE PRECISION, PARAMETER ::     &
   C_Light         = 299792458.0D0, &
   Electron_Charge = 1.6021892D-19, &
   Electron_Mass   = 0.9109534D-30

 REAL :: Radius, Amps
 REAL, DIMENSION(3) :: Position, Velocity
 REAL, PARAMETER :: Frequency= 2.85e9

 INTEGER :: Number, Nr_Of_Clouds_So_Far, iSeed, Nr_of_clouds

 integer :: jx,jy,nx,ny,iColour,jjColour,nTime,jTime,iTime
 integer :: iPass
 real :: x,y,charge_over_mass,charge,timestep,cloudsize,time
 real :: dx, dy
 real :: deltaxposition,deltayposition,deltazposition
 real :: z0, p90, FirstTime, LastTime


 !
 ! Read the parameters from the standard input unit:
 !
    READ (*,*) iTime               ! The first timestep
    READ (*,*) TimeStep            ! The width of the timesteps
    READ (*,*) CloudSize           ! The size of the macroparticles
    READ (*,*) Nr_Of_Clouds_So_Far ! The number of macroparticles
                                   ! which were already emitted
                                   ! in earlier timesteps.
!
! Extract 1 nC of electrons from the cathode, located in the center
! of the cavity wall (z=1 mm in the notation of the RF-gun.geo file)
! and having a 1 mm spot-size.
! The photoelectrons are extracted for 10 ps with a uniform time
! distribution (the transverse distribution is also
! uniform with cylindrical symmetry). The RF peak field on the cathode
! is 140 MV/m and the launching phase (beginning of the electron
! extraction) is 30 RF deg (90 deg is the max of the electric field).
!

    iSeed= -1     ! Initialise the random number generator

    Z0= 1.1e-3    ! The z-position where the macroparticles
                  ! come into existance.
                  ! The lowest z-coordinate in the cavity is 1mm
    Radius= 10e-3    ! The radius of the beam
!
! If the field at t=0 is the electric field at its negative maximum,
! 'p90' must be 90, to inject the particles at 30 degrees.
!
!!    p90= 0
    p90= 90
    FirstTime= ((p90+30.0)/360.0) / Frequency
    LastTime= FirstTime + 10 * 1e-12 * (Radius / 0.5e-3)
!  LastTime= FirstTime + 50 * TimeStep

    Amps= 1e-9 / (LastTime-FirstTime)   ! 1nC Gesamtladung

 !
 ! The output of this program are the properties of the emitted
 ! particles. For each timestep, there is a dataset which describes
 ! the macroparticles, which shall come into existance at that timestep.
 !
 ! The number of timesteps for which this program writes datasets
 ! can be larger than 1.
 !

    NTime= 100      ! Nr of timesteps for which we write the dataset

    Number= Nr_Of_Clouds_So_Far

    WRITE (*,*) NTime,' # Nr of timesteps..'

    For_NTimes: &
    DO jTime= iTime, iTime+NTime-1, 1
       Time= jTime*TimeStep

       !
       ! Particles are only emitted in a small time.
       !
       IF ((Time .LT. FirstTime) &
      .OR. (Time .GT. LastTime)) THEN
          Nr_of_Clouds = 0
          WRITE (*,*) Nr_of_Clouds,' # Nr of Clouds to follow.'
          CYCLE
       ENDIF

       !
       ! The first record in a dataset must contain the number
       ! of particles which are described by that dataset.
       ! We compute that number by generating the data of the particles
       ! twice. In the first pass, we count that particles, and
       ! in the second run, we write the data of the particles.
       !
       Nr_of_Clouds= 0

       For_two_passes : &
       DO iPass= 1, 2, 1

          IF (iPass .EQ. 2) THEN
             WRITE (*,*) Nr_of_Clouds,' # Nr of clouds to follow.'
          ENDIF

          !
          ! We want to have more than one particle per gridcell
          !
          dx= Cloudsize / 2
          dy= dx

          jjColour= 1
          nx= MAX(2, NINT( 2*Radius / dx ) )
          ny= nx

          For_jy_Column : &
          DO jy= -(ny/2), (ny/2), 1

             For_jx_Row : &
             DO jx= -(nx/2), (nx/2), 1
                x= jx*dx
                y= jy*dy

                !
                ! If the position (x,y) really is inside of the circle
                ! with radius 'radius'
                !
                Inside_Circle : &
                IF (x**2 + y**2 .LT. Radius**2) THEN
                   jjColour= jjColour+1

                   Velocity(1)= 0.0
                   Velocity(2)= 0.0
                   Velocity(3)= 0.0

                   !
                   ! We eject the particles at random positions.
                   ! The amplitude of the fluctuations shall be such
                   ! that the average current density is uniform.
                   !  => the x-amplitudes shall be 'dx'
                   !  => the y-amplitudes shall be 'dy'
                   !  => the z-amplitudes shall be such that the z-width
                   !     of the sheet of particle created in the current
                   !     timestep is the same size as the pathlength
                   !     of the particles in one timestep.
                   !
                   ! We cannot take 'Velocity(3)' here, since that is zero.
                   ! We guess that the particles are accelerated ver fast.
                   !
                   DeltaXPosition= dx
                   DeltaYPosition= dy
                   DeltaZPosition= 0.5 * C_Light * TimeStep
                   Position(1)= x + DeltaXPosition * (0.5-wbran1(iSeed))
                   Position(2)= y + DeltaYPosition * (0.5-wbran1(iSeed))
                   Position(3)= Z0 + DeltaZPosition * wbran1(iSeed)

                   IF (iPass .EQ. 1) THEN
                      !
                      ! In the first pass, we only count the number of
                      ! of macroparticles which will be created.
                      !
                      Nr_of_Clouds= Nr_of_Clouds + 2
                   ELSE
                      !
                      ! In the second pass, we know how many particles
                      ! will be created: 'Nr_of_Clouds'.
                      ! We adjust the charge of the macroparticles such that
                      ! the total current equals 'Amps'
                      !
                      Charge= (1.0/Nr_of_Clouds) * Amps * TimeStep
                      Charge_Over_Mass= Electron_Charge / Electron_Mass

                      !
                      ! In addition to the physical properties that each
                      ! macroparticles has, they also have to nonphysical
                      ! properties, which are: 'iColour' 'Number'.
                      ! Particles with the same value 'iColour' are plotted
                      ! by the postprocessor in the same colour.
                      ! The 'Number' property is not yet used by GdfidL.
                      !
                      iColour= MOD(jjColour, 30)

                      !
                      ! 'Number' was initialised to the total number
                      ! of clouds so far.
                      ! We count up the 'Number', so each macroparticle
                      ! has a unique value for 'Number'.
                      !
                      Number= Number + 1

                      !
                      ! We write the properties of the current macroparticle.
                      !
                      WRITE (*,*) ' # Next Cloud..'    ! Could be empty line
                      WRITE (*,*) Charge_Over_Mass     ! A real number
                      WRITE (*,*) Charge               ! A real number
                      WRITE (*,*) Position(:)          ! Three real numbers
                      WRITE (*,*) Velocity(:)          ! Three real numbers
                      WRITE (*,*) iColour              ! Any positive integer
                      WRITE (*,*) Number               ! Any positive integer

                      !
                      ! We also create an anti particle at exactly
                      ! the same position, but with an initial velocity
                      ! downwards.
                      !
                      Velocity(1)= 0
                      Velocity(2)= 0
                      Velocity(3)= -0.8*C_Light

                      Number= Number + 1
                      WRITE (*,*) ' # Next Cloud..'
                      WRITE (*,*) -Charge_Over_Mass
                      WRITE (*,*) -Charge
                      WRITE (*,*) Position(:)
                      WRITE (*,*) Velocity(:)
                      WRITE (*,*) iColour
                      WRITE (*,*) Number

                   ENDIF
                ENDIF Inside_Circle
             ENDDO For_jx_Row
          ENDDO For_jy_Column
       ENDDO For_two_Passes
    ENDDO For_NTimes

    !
    ! This is a main program.
    ! The datasets for NTIME timesteps have been written.
    ! The program stops now, and will be executed again by GdfidL
    ! after NTIME timesteps have been performed.
    !
    STOP

 CONTAINS
    !
    ! This is a random number generator.
    !
    FUNCTION wbran1(IDUM)
    IMPLICIT REAL (a-h, o-z)
    IMPLICIT INTEGER (i-n)

    INTEGER, INTENT(INOUT) :: IDUM

    PARAMETER (IA= 16807, IM= 2147483647, AM= 1./IM, IQ= 127773,IR= 2836, &
               NTAB= 32, NDIV=1+(IM-1)/NTAB, EPS=1.2E-7, RNMX= 1.-EPS)

    ! Minimal ranmdom generator of Park and Miller with
    ! Bays-Durham shuffle and added safeguards.
    ! Returns a uniform random deviate between 0.0 and 1.0
    ! (exclusive of the endpoint values).
    ! Call with idum a negative integer to initialise;
    ! Thereafter, do not alter idum between successive deviatives
    ! in a sequence. RNMX should approximate the largest floating value
    ! that is less than 1.

    INTEGER, DIMENSION(NTAB), SAVE :: iv= 0
    INTEGER, SAVE :: iy= 0

      IF ((idum .LE. 0) .OR. (iy .EQ. 0)) THEN
         !
         ! Initialise
         !
         ! Be sure to prevent idum=0
         !
         idum= MAX(-idum,1)
         DO j= NTAB+8, 1, -1
            !
            ! Load the shuffle table (after 8 warm-ups)
            !
            k= idum/iq
            idum= ia*(idum-k*iq)-ir*k
            IF (idum .LT. 0) idum= idum+im
            IF (j .LE. NTAB) THEN
               iv(j)= idum
            ENDIF
         ENDDO
         iy= iv(1)
      ENDIF

      k= idum/iq

      !
      !     Compute idum=mod(IA*idum,IM) without overflows by Schrage's method
      !
      idum= ia*(idum-k*iq)-ir*k
      IF (idum .LT. 0) idum= idum+im

      !
      ! j will be in the range 1:NTAB
      !
      j= 1+iy/ndiv
      IF ((j .LT. LBOUND(iv, DIM= 1)) &
     .OR. (j .GT. UBOUND(iv, DIM= 1))) THEN
          WRITE (*,*) 'wbran1: j out of range: j, iy, ndiv:', j, iy, ndiv
         STOP
      ENDIF

      !
      ! Output previously stored value and refill the shuffle table
      !
      iy= iv(j)
      iv(j)= idum

      !
      ! Because users don't expect endpoint values
      !
      wbran1= MIN(am*iy,rnmx)
    END FUNCTION wbran1

 END PROGRAM Beam
