Title :DCJAEMMC Keywords :Monte Carlo Simulations, XEDS Resolution Computer :IBM PC and compatibles Operating System :MS DOS v 2.0 or greater Programming Language :TurboPascal Hardware Requirements :Math Chip Optional but recommended Author(s) :David C. Joy Correspondence Address :University of Tenn :Em Facility, Dept. of Zoology, Knoxville, TN :(615)-974-3638 :BITNET:JOY@UTKVX1 AEMMC - uses a single scattering Monte Carlo model to investigate the spatial resolution of X ray production in thin films. After computing and displaying the electron trajectories the program plots two pieces of data. The first is the cumulative profile of Xray generation in the solid, calibrated in angstroms and the 0% and 90% values are drawn in to enable the value to be measured easily. The second display calculates the Xray signal that would be observed from a thin unit width plane-parallel layer such as a monolayer of some material sandwiched between two grains. ------------------------------------------------------------------------ Title :DCJAEMMC Keywords :Monte Carlo Simulations, XEDS Resolution Computer :IBM PC and compatibles Operating System :MS DOS v 2.0 or greater Programming Language :TurboPascal Hardware Requirements :Math Chip Optional but recommended Author(s) :David C. Joy Correspondence Address :University of Tenn :Em Facility, Dept. of Zoology, Knoxville, TN :(615)-974-3638 :BITNET:JOY@UTKVX1 Monte Carlo Simulations of Electron Beam - Solid Interactions The disc supplied with these notes in designed for use with an IBM PC or most compatible clones. For some programs you will need a graphics adapter. These programs will automatically determine which type of graphics card is in use and set themselves up appropriately. It is also desirable that the computer has a maths co-processor chip installed (i.e an 8087 or 80x87) as this will greatly speed up the operation of the program. The program will sense whether or not a maths chip is installed and use it if it is available. Any version of MS-DOS above 2.0 is supported. The basic programs are supplied in two forms - as source code files in Version 5.0 Turbo Pascal, and identified by the .PAS extension; and as executable files identified by the .EXE extension. The disc also contains the graphics drivers, .BGI extension. If you copy this disc to a hard disc all of these files should be in the same directory. You will need a minimum 256k memory to run these demonstrations. To modify any of the code you will also require Version 5.0 or higher Turbo Pascal. The programs on the disc are: AEMMC - uses a single scattering Monte Carlo model to investigate the spatial resolution of Xray production in thin films Note that when typing in these names it does not matter whether you use upper case, lower case or a mixture of both. These programs are copyrighted (C) by David C Joy 1990. They may be freely copied and distributed but only for non-commercial purposes. If you use them in a publication I would appreciate an acknowledgement. Any suggestions for improvements are welcomed, and any problems or bugs that you identify will also be gratefully acknowledged. I can be reached at (615)-974-3638 or on BITNET as JOY@UTKVX1 or on FAX at (615)-974-3642 Using the programs The attached note describes the basic physics that are built into these programs. Reading this section before using the programs will probably be helpful in giving you some idea of what to expect. To run a program: (1) Start your computer in the usual way. If you wish to be able to generate hard copy screen dumps of a CGA screen on to your printer then when the ">" prompt appears type "graphics" and hit return. If you are using a VGA or EGA graphics card then screen dumps require a special utility program. (2) If you booted your computer from a floppy disc then when the ">" prompt reappears put the Monte Carlo disc in the floppy disc drive from which you booted-up (usually drive A). If you start your computer from a hard disc then put the Monte Carlo disc into floppy drive A as before but at the "C>" prompt type "A: ". The prompt will now read "A>". In either case now choose the program that you wish to run, e.g AEMMC and then type AEMMC and follow this with the return key. Alternatively you may wish to copy the Monte Carlo disc to a sub-directory of your hard disc. Be careful to put the runnable program (.EXE files) and the graphics driver files (.BGI) in the same directory. If you choose your sub-directory name to be (say) MC then to run a program, at the "C>" prompt type "cd\MC" then when the "C>" returns type "aemmc " . (3) If the program loads but fails to run you may be using a non-standard video card. In that case try a program that does not require graphics. (4) The programs may prompt you for various pieces of input data such as the atomic number, atomic weight and density of your sample. The list below contains data for a few selected common elements and materials. If you make a mistake when entering data hit the escape key then reenter the correct value. DATA TABLE FOR MC PROGRAMS Material Z At.Wt Density Carbon 6 12 2.2 Aluminum 13 27 2.7 Silicon 14 28 2.34 Chromium 26 52 7.1 Copper 29 63.55 8.96 Silver 47 107.9 10.5 Gold 79 193 19.6 GaAs 32 72.2 5.3 Examples of use for each program AEMMC This program uses the single scattering code to investigate the spatial resolution of Xray microanalysis in thin films. The program initially looks like SS_MC but asks, in addition, for the diameter of the electron probe that is to be used. Typically this will be in the range from 5 angstroms for a dedicated STEM to about 75A for a conventional TEM/STEM. If in doubt put the probe size equal to zero. After computing and displaying the electron trajectories as before the program freezes until the return key is hit when it performs some numerical integrations (with sound effects) and plots two pieces of data. The first is the cumulative profile of Xray generation in the solid - for example if the electron beam is scanned towards a plane parallel grain boundary then the graph shows the Xray signal from the grain on the far side of the boundary would vary with the position of the beam. By analogy with "resolution" in the SEM we can define the Xray resolution as the distance the beam must travel to take the signal from 10% to 90% of its maximum value. The scale is calibrated in angstroms and the 10% and 90% values are drawn in to enable the value to be measured easily. Hitting the return key again produces a second display. This calculates the Xray signal that would be observed from a thin unit width plane-parallel layer such as a monolayer of some material sandwiched between two grains. The maximum level of the signal relative to the matrix value is shown - note that this does not account for effects such as relative ionization cross-sections or fluorescent yields. The dotted line on the display plots the corresponding signal from the surrounding matrix. In this case the resolution might be taken as the full width half maximum of the distribution. Suggested trial experiments: (a) Using copper as the material, compute the X-ray resolution for an incident energy of 100keV, for a beam diameter of 30A and for sample thicknesses of 250, 500, 750, 1000, 1500 and 2000A using a minimum of 500 trajectories per run. How does the resolution vary with the foil thickness? Is the probe size ever a limiting factor? What happens at 200keV? (b) Repeat experiment (a) for foils of carbon, silicon, silver and gold. What is the effect of making the foil of a lighter or heavier element. Compare the computed Xray resolution with the 'beam broadening' predicted by the simple formula derived by Goldstein et al (SEM 1977 Proceedings, ed O Johari, 1, p315). (c) For cases (a) and (b) how does the spatial resolution for an interfacial layer vary ? How is it related to the probe size and the resolution determined previously?  Title :DCJAEMMC Keywords :Monte Carlo Simulations, XEDS Resolution Computer :IBM PC and compatibles Operating System :MS DOS v 2.0 or greater Programming Language :TurboPascal Hardware Requirements :Math Chip Optional but recommended Author(s) :David C. Joy Correspondence Address :University of Tenn :Em Facility, Dept. of Zoology, Knoxville, TN :(615)-974-3638 :BITNET:JOY@UTKVX1 Program AEMMC; { Thin Film Monte Carlo Simulation with graphical plot of trajectories and computation of Xray generation (C) DCJ 04/16/90 } {$N+} {turn on numeric coprocessor} {$E+} {turn on emulator package} { This program is in the public domain and may be copied and distributed. Commercial use is prohibited without written permission. If you use this code for any publication an acknowledgement would be appreciated. } uses CRT,DOS,GRAPH; {resources required} label exit; const two_pi=6.28318; e_min=0.5; {cutoff energy in keV in bulk case} var inc_energy,at_num,at_wht,density,mn_ion_pot:extended; sp,cp,ga,bk_sct,al_a,sg_a,lam_a,er,fe_val:extended; x,y,z,xn,yn,zn,ca,cb,cc,cx,cy,cz:extended; step,m_step,m_t_step,del_E,s_en:extended; plot_scale,hplot_scale,m_f_p,thick,ll,radial:extended; gliset,glgset,probe_size:extended; rstep,zstep:extended; num,traj_num,k,count,dum,bottom,answer:integer; center,top:integer; GraphDriver:Integer; GRAPHMODE:Integer; ErrorCode:Integer; Color,Xasp,Yasp:word; Xray:array[0..51,0..20] of extended; Strip:array[0..50] of extended; Function power(mantissa,exponent:real):real; {this is necessary because PASCAL does not have an exponentiation function} begin if mantissa<=0 then power:=0 else power:=exp(ln(mantissa)*exponent); end; Function stop_pwr(energy:real):real; {this computes the stopping power in keV/gm/cm2 using the modified Bethe expression } var temp:extended; begin if energy<0.05 then energy:=0.05; temp:=ln(1.166*(energy+0.85*mn_ion_pot)/mn_ion_pot); stop_pwr:=temp*78500*at_num/(at_wht*energy); end; Function gasdev:real; {this generates a gaussian deviate of zero mean and unit variance to simulate the finite incidence probe diameter. The code is adapted from that in Numerical Recipes, by Press et al., Cambridge University Press, 1986, page 203} var fac,r,v1,v2:real; begin if gliset=0 then begin repeat v1:=2.0*random-1.0; {pick a random number} v2:=2.0*random-1.0; {and another} r:=(v1*v1)+(v2*v2); {are they in unit circle?} until r<1.0; {if not try again} fac:=sqrt(-2.0*ln(r)/r);{Box-Muller transform} glgset:=v1*fac; {get two normal deviates -save one} gasdev:=v2*fac; {and return the other} gliset:=1; end else {we already have a spare value} begin gasdev:=glgset; {so return it} gliset:=0; {and reset the flag to zero} end; end; Procedure get_constants; {computes some constants needed by the program} begin al_a:=power(at_num,0.67)*3.43E-3; {relativistically correct the beam energy for use up to 500keV} er:=(inc_energy+511.0)/(inc_energy+1022.0); er:=er*er; lam_a:=at_wht/(density*6.0E23); {lambda in cm} lam_a:=lam_a*1.0E8; {put into angstroms} sg_a:=at_num*at_num*12.56*5.21E-21*er; end; Function lambda(energy:real):extended; {computes elastic MFP for single scattering model} var al,ak,sg:real; begin al:=al_a/energy; ak:=al*(1.+al); {giving sg cross-section in cm2 as} sg:=sg_a/(energy*energy*ak); {and lambda in angstroms is} lambda:=lam_a/sg; end; Function broadening:extended; {computes Goldstein et al beam broadening approximation} var var1,var2,var3:extended; begin var1:=sqrt(density/at_wht); var2:=at_num/(inc_energy*1000); var3:=thick*1E-8; var3:=power(var3,1.5); broadening:=6.25E13*var2*var1*var3; {in angstroms} end; Procedure set_up_screen; {gets the input data to run this program} begin ClrScr; {erases all previous data from screen} GoToXY(22,1); Writeln('Single Scattering Monte Carlo Simulation' ); {Having set up the screen now get input data} GoToXY(1,5); Write('Input beam energy in keV '); Readln(Inc_Energy); GoToXY(1,7); Write('Target Atomic Number is '); Readln(at_num); { Calculate J the mean ionization potential mn_ion_pot using the Berger-Selzer analytical fit } mn_ion_pot:=(9.76*at_num + (58.5/power(at_num,0.19)))*0.001; GoToXY(1,9); Write('Target Atomic Weight is '); Readln(At_wht); GoToXY(1,11); Write('Target density in gm/cc is '); Readln(Density); GoToXY(40,5); write('Probe diameter (angstroms) ?'); readln(probe_size); if probe_size<1 then probe_size:=1; {convert this to unit variance for GasDev routine} probe_size:=probe_size/2.24; GoToXY(40,7); write('Foil thickness (A) '); readln(thick); if thick<20 then thick:=20; {set up scales for the radial and depth counters} zstep:=thick/20; {use Goldstein et al approximation for scattering diameter} rstep:=sqrt(broadening*broadening+probe_size*probe_size)/50; {to get a sensible width to the screen display round-off this number} rstep:=trunc(rstep)*1.0; if rstep<1.0 then rstep:=1.0; { get the number of trajectories to be run in this simulation } GoToXY(40,9); write('Number of trajectories required '); readln(traj_num); end; Procedure initialize; {sets up the graphics drivers for version 5.0 Turbo Pascal} var InGraphicsMode:Boolean; PathToDriver:String; begin DirectVideo:=False; PathToDriver:=''; GraphDriver:=detect; InitGraph(GraphDriver,GRAPHMODE,PathToDriver); SetViewPort(0,0,GetMaxX,GetMaxY,True); Color:=GetMaxColor; center:=trunc(GetMaxX/2); {beam entry coordinate} top:=trunc(GetMaxY*0.1); {plot position of entry surface} end; Procedure xyplot(a,b,c,d:real); {this displays the trajectories on the pixel screen} var iy,iz,iyn,izn:integer; begin iy:=center + trunc(a*hplot_scale); {plotting coordinate #1} iz:=top + trunc(b*plot_scale); {plotting coordinate #2} iyn:=center + trunc(c*hplot_scale); {plotting coordinate #3} izn:=top + trunc(d*plot_scale); {plotting coordinate #4} if d=99 then izn:=top-2; {BS plotting limit} if d=999 then izn:=bottom+2;{transmitted plotting limit} {plot this vector on the screen} line(iy,iz,iyn,izn); end; Procedure set_up_graphics; {draws in the surfaces, beam location, and action thermometers} var a,b:integer; begin a:=GetMaxX-20; {adjust to suit your screen} b:=20; {ditto} Line(b,top,a,top); {plot in top surface} {find position of bottom surface} bottom:=top + trunc(thick*GetMaxY/1000); if bottom>GetMaxY-40 then bottom:=GetmaxY-40; plot_scale:=(bottom-top)/thick; {hence get pixels/angstrom} GetAspectRatio(Xasp,Yasp); hplot_scale:=(Yasp/Xasp)*plot_scale; {compensate for screen aspect ratio} {also plot in exit surface} Line(b,bottom,a,bottom); {plot bottom surface} Line(center,1,center,top); {plot incident beam} OutTextXY(trunc(center-80),bottom+15,'0%.......50%......100%'); OutTextXY(trunc(center-78),bottom+28,'Trajectories completed'); end; Procedure reset_coordinates; {resets coordinates at start of each trajectory} begin s_en:=inc_energy; x:=probe_size*GasDev; {set up a Gaussian probe} y:=probe_size*GasDev; {ditto} z:=0; cx:=0; cy:=0; cz:=1; end; Procedure zero_counters; {since PASCAL does not zero variables on start-up we must do this} var j,k:integer; begin bk_sct:=0; num:=0; gliset:=0; for j:=0 to 51 do begin for k:=0 to 20 do begin xray[j,k]:=0; end; end; end; Procedure s_scatter(energy:real); {calculates scattering angle using screened Rutherford cross-section} var R1,al:real; begin al:=al_a/energy; R1:=random; cp:=1.0-((2.0*al*R1)/(1.0+al-R1)); sp:=sqrt(1.0-cp*cp); {and get the azimuthal scattering angle} ga:=two_pi*random; end; Procedure new_coord(step:real); {gets xn,yn,zn from x,y,z and scattering angles} var an_n,an_m,v1,v2,v3,v4:real; begin { find the transformation angles } if cz=0 then cz:=0.000001; an_m:=(-cx/cz); an_n:=1.0/sqrt(1+(an_m*an_m)); { save computation time by getting all the transcendentals first } v1:=an_n*sp; v2:=an_m*an_n*sp; v3:=cos(ga); v4:=sin(ga); { find the new direction cosines} ca:=(cx*cp) + (v1*v3) + (cy*v2*v4); cb:=(cy*cp) + (v4*(cz*v1 - cx*v2)); cc:=(cz*cp) + (v2*v3) - (cy*v1*v4); { and get the new coordinates } xn:=x + step*ca; yn:=y + step*cb; zn:=z + step*cc; end; Procedure straight_through; {handles case where initial entry exceeds thickness} var zpos:integer; begin num:=num+1; xyplot(y,0,y,999); zpos:=trunc(thick*random/zstep); {Xray produced at random depth} xray[0,zpos]:=xray[0,zpos]+1; end; Procedure back_scatter; {handles case of backscattered electrons} var rval:extended; rpos:integer; begin num:=num+1; bk_sct:=bk_sct+1; xyplot(y,z,yn,99); rval:=sqrt(x*x+y*y); rpos:=trunc(rval/rstep); if rpos>50 then rpos:=50; xray[rpos,0]:=xray[rpos,0]+1; end; Procedure transmit_electron; {handles case of transmitted electron} var ll:extended; begin num:=num+1; ll:=(thick-z)/cc; yn:=y+ll*cb; {exit y-coordinate} xyplot(y,z,yn,999); end; Procedure add_xray; var xm,ym,zm,rpos:extended; rval,zval:integer; {general procedure to add Xray yield} begin {find mid-point of trajectory} xm:=(x+xn)/2; ym:=(y+yn)/2; zm:=(z+zn)/2; {find radial position} rpos:=sqrt(xm*xm+ym*ym); {and express as integer} rval:=trunc(rpos/rstep); {check for out of bounds} if rval<0 then rval:=0; if rval>50 then rval:=51; {put it out of bounds} {corresponding depth variable} zval:=trunc(zm/zstep); {check out of bounds} if zval<0 then zval:=0; if zval>20 then zval:=20; {now add the yield term} xray[rval,zval]:=xray[rval,zval]+1; end; Procedure reset_next_step; {resets variables for next trajectory step} begin xyplot(y,z,yn,zn); cx:=ca; cy:=cb; cz:=cc; x:=xn; y:=yn; z:=zn; {find the energy lost on this step} del_E:=step*stop_pwr(s_en)*density*1E-8; {so the current energy is} s_en:=s_en - del_E; end; Procedure show_traj_num; {updates thermometer display for % of trajectories done} var a,b:integer; begin a:=trunc(center-80); {adjust to suit your screen} b:=bottom+23; {ditto} line(a,b,a+trunc(165*(num/traj_num)),b); end; Procedure get_strips; {integrates xray yield to get intensity in strip perpendicular to y axis and running through thickness of foil} var j,k,l,irval:integer; S,T:string; freq:real; begin freq:=100; S:='Please be patient...I am integrating data'; T:='-'; OutTextXY(trunc(GetMaxX/4),trunc(GetMaxY/2),S); {initialize the array we will use} for j:=0 to 50 do begin strip[j]:=0; end; {now do the integration} for j:=0 to 50 do begin for k:=-50 to 50 do begin irval:=trunc(sqrt(j*j+k*k)); if irval>50 then irval:=50; {boundary check} for l:=0 to 20 do begin {remember to normalize by intercept of annulus} strip[j]:=strip[j]+(xray[irval,l]/(irval+1)); end; end; T:=concat('-',T); OutTextXY(trunc(GetMaxX/4)-30,trunc(GetMaxY/2)+5,T); Sound(trunc(freq)); Delay(150); NoSound; freq:=freq*1.03; end; end; Procedure show_BS_coeff; {displays BS coefficient on thermometer scale} label hang; var a,b:integer; begin a:=GetMaxX-180; {adjust to suit your screen} b:=bottom+23; {ditto} OutTextXY(a,b-8,'0....0.25...0.5'); OutTextXY(a+5,b+5,'BS coefficient'); OutTextXY(15,b+5,' to continue'); Line(a,b,trunc(a+(bk_sct/traj_num)*220),b); hang: {this loop freeze the display on the screen} If (not keypressed) then goto hang; end; Procedure plot_results; label freeze; var a,b,c,d,J,K:integer; xpos,ypos,width:integer; s:string[11]; profile:array[0..100] of extended; vscale,hscale,yval,dummy:extended; begin {by plotting the box} a:=trunc(GetMaxX*0.1); b:=trunc(GetMaxX*0.9); c:=trunc(GetMaxY*0.1); d:=trunc(GetMaxY*0.9); line(a,c,b,c); {top} line(a,d,b,d); {bottom} line(a,c,a,d); {LHS} line(b,c,b,d); {RHS} {scale the display for plotting} hscale:=(b-a)/100; vscale:=(d-c); SetlineStyle(DottedLn,0,NormWidth); for J:= 1 to 9 do {ticks on bottom axis} begin line(a+trunc(J*10*hscale),d,a+trunc(J*10*hscale),c); end; {put in 10% and 90% markers} line(a,d-trunc(0.1*vscale),b,d-trunc(0.1*vscale)); line(a,d-trunc(0.9*vscale),b,d-trunc(0.9*vscale)); line(trunc((a+b)/2),d,trunc((a+b)/2),c); SetLineStyle(SolidLn,0,NormWidth); OutTextXY(a-25,d-trunc(0.1*vscale)-5,'10%'); OutTextXY(a-25,d-trunc(0.9*vscale)-5,'90%'); OutTextXY(trunc(GetMaxX/3)-20,c-10,'Integrated Xray signal vs position'); {calibrate display width} width:=trunc(50*rstep); Str(width,s); s:=concat(s,' A'); OutTextXY(a-3,d+5,'-'+s); OutTextXY(b-7,d+5,'+'+s); OutTextXY(trunc((a+b)/2-4),d+5,'0'); {now do the profile integration} dummy:=0; for J:=0 to 100 do begin profile[J]:=0; K:=abs(50-J); dummy:=dummy+strip[k]; profile[J]:=dummy; end; {normalize this to maximum value and plot} MoveTo(a,c); {move pointer to corner} for J:=0 to 100 do begin yval:=profile[J]/profile[100]; ypos:=d-trunc(yval*vscale); xpos:=a+trunc(J*hscale); LineTo(xpos,ypos); end; readln; freeze: if (not keypressed) then goto freeze; end; Procedure plot_results2; label freeze; var a,b,c,d,J,K:integer; xpos,ypos,width:integer; s:string[11]; t:string[11]; profile:array[0..100] of extended; vscale,hscale,rescale,yval,dummy:extended; begin {by again plotting the box} ClearViewPort; a:=trunc(GetMaxX*0.1); b:=trunc(GetMaxX*0.9); c:=trunc(GetMaxY*0.1); d:=trunc(GetMaxY*0.9); line(a,c,b,c); {top} line(a,d,b,d); {bottom} line(a,c,a,d); {LHS} line(b,c,b,d); {RHS} {scale the display for plotting} hscale:=(b-a)/100; vscale:=(d-c-2); SetlineStyle(DottedLn,0,NormWidth); for J:= 1 to 9 do {ticks on bottom axis} begin line(a+trunc(J*10*hscale),d,a+trunc(J*10*hscale),c); end; {put in 20%, 50% and 80% markers} line(a,d-trunc(0.2*vscale),b,d-trunc(0.2*vscale)); line(a,trunc((d+c)/2),b,trunc((d+c)/2)); line(a,d-trunc(0.8*vscale),b,d-trunc(0.8*vscale)); line(trunc((a+b)/2),d,trunc((a+b)/2),c); SetLineStyle(SolidLn,0,NormWidth); OutTextXY(a-25,d-trunc(0.2*vscale)-5,'20%'); OutTextXY(a-25,trunc((d+c)/2)-5,'50%'); OutTextXY(a-25,d-trunc(0.8*vscale)-5,'80%'); OutTextXY(trunc(GetMaxX/3)-20,c-10,'Interface and Matrix Profiles'); {calibrate display width} width:=trunc(50*rstep); Str(width,s); s:=concat(s,' A'); OutTextXY(a-3,d+5,'-'+s); OutTextXY(b-7,d+5,'+'+s); OutTextXY(trunc((a+b)/2-4),d+5,'0'); {now compute and plot the two displays} dummy:=0; for J:=0 to 100 do {total Xray signal} begin K:=abs(50-J); dummy:=dummy+strip[k]; end; {plot the interface signal component} MoveTo(a,c); for J:=1 to 49 do begin {by smoothing the data} strip[J]:=(strip[J-1]+strip[J]+strip[J+1])/3; end; rescale:=strip[0]/(dummy*rstep); for J:=0 to 100 do begin k:=abs(50-J); yval:=strip[k]/strip[0]; {normalized value} ypos:=d-1-trunc(yval*vscale*rescale); xpos:=a+trunc(J*hscale); LineTo(xpos,ypos); end; Str((100*rescale):4:2,t); t:=concat('Max% = ',t); OutTextXY(trunc((a+b)/2+20),d-trunc(vscale*0.5),t); {now plot the matrix signal} for J:=0 to 100 do begin K:=abs(50-J); yval:=(dummy-strip[K])/dummy; ypos:=d-1-trunc(yval*vscale); xpos:=a+trunc(J*hscale); PutPixel(xpos,ypos,color); end; freeze: if (not keypressed) then goto freeze; end; { ******************************************************** * this is the start of the main program * ********************************************************} begin set_up_screen; {get input data and find J value} get_constants; {other parameters needed later on} { reset the random number generator} randomize; {set up the graphics display for plotting } Initialize; set_up_graphics; { ******************************************************************* * the Monte Carlo Loop * *******************************************************************} zero_counters; while num < traj_num do begin reset_coordinates; {allow initial entrance of electron} step:=-lambda(s_en)*ln(random); zn:=step; if zn>thick then {this one is transmitted} begin straight_through; goto exit; end else {plot this position and reset coordinates} begin xyplot(y,0,y,zn); xray[0,1]:=xray[0,1]+1; {add Xray yield} z:=zn; end; {now start the single scattering loop} repeat step:=-lambda(s_en)*ln(random); s_scatter(s_en); new_coord(step); {decide what happens next} if zn<=0 then {this one is backscattered} begin back_scatter; goto exit; end; if zn>=thick then {this one is transmitted} begin transmit_electron; goto exit; end; {otherwise we go round again} add_xray; reset_next_step; until s_en<=e_min; {end of repeat loop - the energy drops below e_min} num:=num+1; {increment counter} exit: {end of goto jumps} { *********************************************************** * end of the Monte Carlo Loop * *********************************************************** } show_traj_num; if num mod 100=0 then randomize; end; show_BS_coeff; ClearViewPort; ClearDevice; get_strips; ClearViewPort; plot_results; {Xray generation profile} plot_results2; {interface and matrix data} readln; readln; CloseGraph; end.