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 + y**2 .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