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 :: DCcurrent= 0.1
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, Particles_per_Area
real :: deltaxposition,deltayposition,deltazposition
real :: z0
! This program gets its parameters via the standard input unit,
! and writes its output to the standard output unit.
!
! 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.
iSeed= -1 ! Initialise the random number generator
Z0= -0.05 ! The z-position where the macroparticles
! come into existance.
Amps= DCcurrent
Radius= 10e-2 ! The radius of the beam
!
! 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
!
! 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
Particles_per_Area= 100
dx= SQRT( (2*Radius)**2 / Particles_per_Area )
dy= dx
jjColour= 1
nx= NINT(SQRT( Particles_per_Area ))
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 + y2 .LT. Radius*2) THEN
jjColour= jjColour+1
Velocity(1)= 0.0
Velocity(2)= 0.0
Velocity(3)= 0.5 * C_Light
!
! 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.
!
DeltaXPosition= dx
DeltaYPosition= dy
DeltaZPosition= Velocity(3) * 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 + 1
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
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