Description of Radiative Dense Plasma Focus Computation Package RADPFV5.16

 and Downloads -  Lee model code   


·       Numerical Experimental Facility

·       Simulate any Mathers-type plasma focus, computes dynamics

·       Design new plasma focus machines

·       Thermodynamics included; H2, D2, Ne, Ar, Xe, He, N2, Kr and D-T

·       Model parameters to fit experimental axial, radial phase times

·       Radiative phase computes line radiation, recombination and total yield. Computes neutron yield for deuterium operation; based on an improved beam-target model and calibrated at an experimental point of 0.5MA Yn=5.6x10^9 (extended to D-T). Thermonuclear fusion neutrons also computed for D and D-T operation. Termination of computation of slow phase: transit time for slow disturbance speed.  Plasma Self-absorption based on revised equations presented in File 3; appendix by N A D Khattak.

Also includes:

Time guard feature

Choice of Tapered electrode [note: if this feature is not wanted ensure taper? y=1, n=0 cell is filled with a ‘0’]

Quick choice of specified machines; one click loading of chosen machine; at present 3 machines may be click-loaded: the UNU/ICTP PFF, the NX2 and the PF1000 [note: if this feature is not wanted ensure the relevant cell in top RH corner of sheet is filled with a ‘0’]

Acknowledgement: Scroll down to see acknowledgement.

There are altogether 4 files in this package.

File1: PDF File "Description of Radiative Dense Plasma  Focus Computation Package": This file

File2: PDF File "Theory of Radiative Plasma Focus Model"                 

File3: PDF file "Appendix by N A D Khattak".   


"Radiative PF Standard Code" RADPFV5.16 and  RADPFV5.16FIB (added Fast Ion Beam section)

(If difficulty with these two above links, try right click, select open link in new window; then click on the reload this page (circular arrow) sign left of the URL address window ; the EXCEL file will be downloaded, appearing in the bottom left of your screen, also in your download folder.)


In addition, there are files for the computation of thermodynamic data needed for this code.



Introductory description


A simple 2 phase (axial and radial) model was developed by S.Lee in 1985 as a component of a 3kJ plasma focus experimental package which became known as the UNU/ICTP PFF. This network of basically identical 3kJ PF machines, with different experimental and application emphases, is now operated by groups in countries including Singapore, Malaysia, Thailand, Indonesia, India, Pakistan,, Egypt and Zimbabwe.


The model was written as a 3 phase (non-radiative) model (in GWBASIC) for an experimental program at the 1991 Spring College in Plasma Physics at the ICTP.


The present 5-phase package (axial, radial inward shock, radial reflected shock, slow compression radiative and expanded large column phase) is re-written in Microsoft EXCEL VISUAL BASIC in order to make it available for wider usage.


The model may be configured to any conventional Mather-type plasma focus by inputting machine parameters: inductance, capacitance, electrode radii and length.  And operating parameters: charging voltage and fill gas pressure. The thermodynamics (specific heat ratio and charge number as functions of temperature) are included for gases including hydrogen, deuterium, neon, argon, helium, nitrogen, Kr and Xe also D-T(50:50).  The gases may be selected by simply inputting atomic number, molecular weight and dissociation number (2 for deuterium and hydrogen, 1 for the others).


The results are the following: waveforms for the total discharge current and tube voltage, axial phase trajectory and speed, radial trajectories for the shock front, current sheath and column length and the corresponding speeds, plasma temperature and radiation yields (bremsstrahlung, line and recombination) and power; and thermodynamic quantities such as specific heat ratios and charge numbers. These are output in graphical as well as tabular forms. Also computed are plasma pinch current and neutron yield, and any distributions, if required.


Are the results any good?

But are there any indications that our computed results are anywhere near the actual results that may be measured on the device in actual operation?


NOT if we just guess the model parameters fm, fc, fmr, fcr; as suggested above. Then the results are just hypothetical; although with experience we may assign some reasonable values of the model parameters for the particular machine in its particular operating conditions. And the results may be valid for planning or designing purposes.


How do we make the results realistic?


The standard practice that we pioneered from the early 90’s is to fit the computed total current waveform to an experimentally measured total current waveform by variations of 2 sets of model parameters (for mass swept up and effective current fractions) one for the axial phase and the other for the radial phase. This technique has turned out to be effective in thus accounting for all physical and machine mechanisms at play during a particular discharge and simulating the effects into the code. The code than runs realistically like that machine for that shot. Some machine effects or even physical mechanisms may change from shot- to- shot, then the fitted model parameters change correspondingly. Thus, even shot-to-shot global variations of the same machine are captured by the code. With realization that the power of this technique is all-encompassing, this technique is being emulated by several other codes. A detailed explanation follows.


          From experience it is known that the current trace of the focus is one of the best indicators of gross performance. The axial and radial phase dynamics and the crucial energy transfer into the focus pinch are among the important information that is quickly apparent from the current trace.


The exact time profile of the total current trace is governed by the bank parameters namely capacitance Co, external, or static inductance Lo and circuit resistance ro, by the focus tube geometry namely electrode radii, outer ‘b’ and inner anode ‘a’, and the anode length ‘zo’; and on the operational parameters which are the charging voltage Vo and the fill pressure Po and the fill gas. It also depends on the fraction of mass swept-up and the fraction of sheath current and the variation of these fractions through the axial and radial phases. These parameters determine the axial and radial dynamics, specifically the axial and radial speeds which in turn affect the profile and magnitudes of the discharge current.  The detailed profile of the discharge current during the pinch phase will also reflect the joule heating and radiative yields. At the end of the pinch phase the total current profile will also reflect the sudden transition of the current flow from a constricted pinch to a large column flow. Thus the discharge current powers all the dynamic, electrodynamic, thermodynamic and radiation processes in the various phases of the plasma focus. Conversely all the dynamic, electrodynamic, thermodynamic and radiation processes in the various phases of the plasma focus affect the discharge current.

It is then no exaggeration to say that the discharge current waveform contains information on all the dynamic, electrodynamic, thermodynamic and radiation processes that occurs in the various phases of the plasma focus.


The discharge current of a plasma focus is one of the simplest measurement that is routinely carried out. Yet this easily measured waveform carries all the information of the physical processes in the plasma focus.


Our standard practice for any existing plasma focus is to obtain a measured current trace. Then we fit the computed current trace to the measured current trace. The fitting process involves adjusting the model parameters fm, fc, fmr, fcr one by one, or in combination until the computed current waveform fits the measured current waveform.


Once this fitting is done our experience is that the other computed properties including dynamics, energy distributions and radiation are all realistic


Fitting computed current trace to experimental current trace of existing machine:


The main model parameters are the tube current flow factor CURRF (known to be 0.7 for most machines) and the mass swept-up factor (MASSF or fm, for axial & MASSFR or fmr, for radial).  These have been pre-selected in the model, but may be adjusted so that the time of focus, and the radial inward shock transit time, fit the experimentally observed times for each machine. The computed current trace may be compared with the experimental current trace.


Features for comparison include current risetime and rising shape, peak current, current 'roll off' and dip, both shape and amplitude. Absolute values should be compared. Our experience with a number of machines shows that the fit is usually very good, occasionally almost exact.

The machine parameters and operating conditions should already have been determined and inputted into the active sheet. The model parameters are then adjusted, one by one, or in combination until best fit is obtained between the computed current trace and the experimental current trace.

Besides the model parameters, sometimes (when all else fails in the fitting process) the inductance (as published or given by the experimenters) needs to be adjusted. The reason is that the inductance Lo may be given as the short circuit bank inductance whereas it should be the ‘static’ inductance of the plasma focus; ie the inductance of the PF before the current sheet moves. Usually also the value of stray resistance r0 needs to be guessed at as few experimenters determine this carefully. We usually start with the value of r0 as 0.1 of (Lo/Co)0.5; and make small adjustment as necessary. Sometimes, especially for PF’s using very low values of Co, it may also be necessary (when all else fails) to adjust the value of Co (for sub-uF capacitor banks, the closely spaced connecting parallel plates and parallel connecting cables may actually significantly change the value of Co. In one or two cases where there is very good fit in current profiles but the absolute values of currents don’t match, it has been reasonable to suspect that the calibration constant for the current profile has been given wrongly by the experimenter.


Designing a new plasma focus


If a machine has not been built the model may be used to aid design. First use the following rule of thumb procedure [use SI units].


What capacitance ( C ) are you planning?

How low is the inductance (L) you expect to attain?

What maximum voltage (V) do you expect to operate?

Enter these values into the appropriate spaces for these machine parameters.


For the stray (circuit) resistance, take 1/4 the value of (L/C)1/2 .

Estimate the undamped peak current using the formula I=V/(L/C)1/2. 

Use (I/a)=250kA max undamped current per cm to assign the value of centre electrode radius 'a'.

Put in double this value for outer electrode radius 'b'.

The length of the electrode may be assigned as 5 times the value of 1.6(LC)1/2 . This length is in cm when the value of (LC)1/2 is expressed in microsecond. (This gives a length which will provide an average axial speed of 5 cm/us which usually gives a peak speed at end of axial phase of 8cm/us. For operation in H2 or D2 20% longer may be better; for neon operation to get suitable line radiation (12-13.5 A) for SXR microlithography purposes, 20% shorter may be better as we require a lower speed to get to the correct level ionization stages. The focus is normally operated so that the start of current dip (signifying the end of the axial phase) occurs at or just after peak current. For argon to generate SXR a high speed is required so use the same length as for H2. For xenon, if the aim is for EUV (around 13 nm) for experiments for NGL the model has predicted a requirement for very low speeds, around 1.3cm/us! So it appears one needs very short anodes, at least 5 times shorter than that needed for D2. However there is not much experimental experience accumulated so far.


 For pressure values assign as follows: D2: 4 torr;     neon: 1.5 torr;     Ar: 0.7 torr.

For xenon, runs on the code suggest several torr to go with the very short anode length.


From the above rule of thumb design parameters, is your PF fat or thin? (according to the ratio length of the centre electrode divided by diameter; for NX2 this ratio is 1.2, fat; for UNU/ICTP PFF this ratio is 17, thin )


If it is fat use the model parameters suggested for the NX2 (These suggested values are tabulated at the top right of the active sheet which appears when you open RADPFV5.008).

If it is thin assign the model parameters closer to the UNU/ICTP PFF.


Run the computation and from results make adjustment to 'a', 'b', length (V may also easily be varied, especially reduced since we have started with max V; C also, use more or less capacitors; careful with L, normally make L as small as possible, but be realistic). Adjust parameters for best results over a range of pressures and gases. Best results could mean strong current dip or biggest line emission in the case of neon, which is useful for developing microlithography SXR sources.


Model may be adapted to suit requirements


The axial phase (trajectory) going into the radial phase (trajectories of shock front, current sheet and length of the focussing column) is portrayed reasonably well. As the radial inward shock goes on-axis, a reflected shock phase follows. The reflected shock moves outward until it hits the incoming radial piston which was moving behind the radial inward shock. Now follows the slow compression, radiative phase.


The radiative phase is the interesting phase, which is presented in the package in a form which gives reasonable results.  The slow piston motion is coupled to the rate of change of current, the elongation and power gain/loss due to joule heating, Bremsstrahlung and line emission and recombination losses.  Thus radiation collapse (critical current of 1.6 MA for deuterium, but much reduced to possibly below 100kA for neon and argon under certain conditions) is included into the modelling.  Reasonable line radiation yields (all lines) are computed for the UNU/ICTP PFF. One may wish to include emission in specific lines. Plasma self-absorption effects (adapted from equations discussed by Khattak in File 3) have been implemented in the code.


There is room for further interesting development.  For example, allowing the radiation collapse to couple to a ‘piston’ motion will lead to a huge voltage spike in a ‘high pressure’ regime which is not observed experimentally. In this model, this effect is ‘artificially’ restricted by ‘house keeping’ procedures in the package.  It should be looked into further.  Instabilities could be introduced into the package by insertion of a suitable ‘anomalous’ resistance time function. This should be coupled into the voltage and current equation; but not into the ‘piston’ equation. The current would rapidly diffuse as this ‘anomalous’ resistance kicks in, causing the necessary abandonment of the concept of a ‘piston’.


The radiative phase is followed by an expanded large column phase, in which the current flows in a large column with the radius of the centre electrode.


The theoretical basis (with all equations) is given in the next (separate) PDF file; File 2., with an appendix by N A D Khattak in PDF file 3.


Most recently a reasonable beam-target mechanism has been incorporated into the code.




This model and code is the result of almost 40 years of experimental work on shock waves and the Plasma Focus. It incorporates the best experience (as understood by me) shared with each and every one of the students and Fellows I have trained and each and every one of the colleagues I have worked with. Some of them are named here: Prof (Dr) Mohamad Akel, Prof (Dr) Tsai Shih Tung, Prof (Dr) V A Gribkov,  Dr Chow Sai Pew,  (the late) Dr Chen You Hor, Yong Yeow Chin, Prof (Dr) Tou Teck Yong, Dr Liu Mahe, Dr Bing Shan, Assoc Prof (Dr) Paul C K Lee, Prof (Dr) Rajdeep Singh Rawat, Dr Moo Siew Pheng, Assoc Prof (Dr) Kwek Kuan Hiang, Prof (Dr) Wong Chiow San, Dr Chew Ah Chuan, Dr Adrian Serban, Dr Susetyo Mulyodrono, Assoc Prof Zhang Guixin, Prof (Dr) Jalil Ali, Assoc Prof (Dr) Saw Sor Heoh, Dr Alin Patran, Tan Thian Hoo, Dr Suresh Kumar, Chin You Hon, Prof (Dr.) C.K.Chakrabarty, Prof (Dr) M Zakaullah, Prof (Dr) A V Gholap, Dr. Ashutosh Srivastava, Prof (Dr) N A D Khattak, Dr. Vahid Damideh, Prakash Gautam, Assoc Prof Arwinder Singh, Dr. Federico Roy, Dr Teh Thiam Oun, Assoc Prof (Dr) Yap Seong Ling, Dr Lim Ling Hong.


Description of Programme Package


The package consists of an ACTIVE SHEET in EXCEL linked to a MACRO (where the basic programme is written).  The package may be operated from the sheet, without going into the MACRO.  The machine parameters may be inputted directly onto the sheet, as may the operating conditions and gas.  The model parameters (CURRF, MASSF and MASSFR), if required to be adjusted can also be directly inputted onto the sheet.


After downloading the programme, the EXCEL Sheet appears.  The first section, first 19 rows, contain essential information and inputted quantities.  Parameters which may be changed directly on the sheet are in bold underlined.


The programme (as downloaded) contains all machine parameters for the UNU/ICTP PFF. Thermodynamic data (for the 3 gases D, Ne and Ar; and xenon) are also preloaded.  The pre-selected operational conditions and gas are shown.


 The programme may be operated directly from the sheet. Place the cursor in any empty cell; then press “Ctrl + A”


The computation should take less than 1 minute (depending on speed of machine and input parameters).


Results are outputted on the sheet in columns as follows (starting row 20):




A                     time in microsecond

B                     current in kA

C                     tube voltage in kV

D                     axial position in cm

E                      axial speed in cm/ms

F                       time in microsecond (starting on row 20, radial phase)

G                      time in nanosecond, referenced to start of radial phase

H                      current in kA (radial phase data only)

I                        tube voltage in kV (radial phase data only)

J                      radial shock position, referenced to axis, in mm

K                      radial piston position, referenced to axis, in mm

L                      axial position of focus column, referenced to anode end, in mm

M                     radial shock speed in cm/ms

N                      radial piston speed in cm/ms

O                      elongation speed of column in cm/ms

P                       reflected shock radial position in mm

Q                      temperature in oK

R                      Joule heating power in watts

S                       Bremsstrahlung emission power in watts   Ab corrected)

T                     recombination emission power in watts (Ab corrected)

U                      line emission power in watts (Ab corrected)

V                      sum of S,T & U (Absorption Corrected)

W                    sum of power ie of R & V (Absorption Corrected)

X                     time-integrated Joule heating in Joules

Y                    time-integrated Bremsstrahlung in Joules (Absorption Corrected)

Z                     time-integrated recombination emission in Joules (Absorption Corrected)

AA                  time-integrated line emission in Joules (Absorption Corrected)

AB                   Total radiation in Joules (Absorption Corrected)

AC                   Total energy gain/loss in J, ie X+AB (Absorption Corrected)

AD                  plasma self-absorption corrected coefficient

AE                   plasma radiation power (if black body) in W

AF                   specific heat ratio

AG                  effective plasma charge number

AH                  Thermal neutron yield

AI                    Beam-target produced neutron yield

AJ                   Total neutron yield

AK                  plasma ion density per m^3

AL                   Volume radiation power in W

AM                 Surface radiation power in W

AN                  plasma self-absorption corrected coefficient

AO                  radial phase piston work in J


The number of data rows may go up to 7000.  The data is also presented near the top of the sheet in graphical forms.  The Figures may be reconfigured from time to time.


The horizontal axis (for Figs 1 & 2) is time in ms.  The other Figures display computed data of the radial phases. Radial phase time scale is in ns, referenced to the start of the radial phase’



series 1

series 2

series 3

series 4






Fig 1 (top left)

Circuit Current




Fig 2 (top right)

axial position

axial speed



Fig 3

Radial shock position

radial piston position

axial focus length


Fig 4





Fig 5

Radial shock speed

radial piston speed

Elongation speed

line radiation energy

Fig 6

Plasma temperature




Fig 7

Joule heat energy

Bremsstrahlung energy

recombination energy

Line radiation

Fig 8

Joule power

Bremstrahlung power

recombination Power

Line radiation power

Inset                            Sp Heat Ratio        charge number


To get into the code (containing the programme lines)

(from top toolbar) click on 'Tools'  or ‘View’

select 'Macro' and select 'View Macros'

select 'radpf005' by clicking on it

 On the right hand panel click on  'step into'

That gets you into the code. You may then modify the code as required.


When finished with the code

click 'red cross' on top right hand corner

click 'OK' when message "This command will stop the Debugger" appears

This gets you back to the active sheet.




Summary Steps


The complete computer package (active sheet and macro code) is provided as an EXCEL File RADPFV5.16 (latest) or with added FIB section RADPFV5.16FIB. (Either may be downloaded from

Open the file.

When the EXCEL Sheet appears, press ‘Ctrl + A”. The programme will run and present data in dataline of row 17 and in Columns (as described) and in the 8 Figures and inset.



It is recommended that you keep a reference copy of the current standard version, which you can refer to.  Modified programmes should be kept under another title such as RadPFTrialVersion. (so that you can get back to the original latest version should your trial code not work)


The model has been tested for the UNU/ICTP PFF 3kJ plasma focus, in the following ranges, operating at 14kV.


            Deuterium      0.5 to 19 torr

            Neon               0.1 to 5.5 torr

            Argon             0.1 to 2.5 torr


Any plasma focus will only be able to operate properly within a range of parameters. For example if the parameters are such that the current sheet moves too slowly in the axial phase, by the time the radial phase starts the drive current may have dropped to too low a value and the radial phase cannot complete.


In other words there needs to be a matching between the sum of the characteristic axial & radial times and the characteristic capacitor discharge time.


The code has recently incorporated a Time Match Safeguard. Before axial phase computation is started the code checks that there is suitable matching within a suitable range. This is done by checking that the ratio ALT  in the code is LESS than a certain value. We have set the lower limit of ALT at 0.68 for D2 and 0.65 for neon and argon. If too large a pressure is set for selected capacitor voltage and the value of ALT falls below this set value an error message will appear recommending to set the pressure lower or the voltage higher. (For some purpose you may wish to disable this time check)


If this happens click the 'red cross' on the upper right hand corner, click OK on the 'debugger' message, getting you back to the active sheet and adjust the pressure and or voltage accordingly.


(The value of ALT for each computation is shown in the active sheet.)


Furthermore if you have changed drastically the values of the parameters such as the capacitance C e.g. changing the value of C from 30 uF to 3 uF, you would have reduced the capacitor discharge time by 3 times (square root of 10), hence it may be difficult to match the axial transit time (by simply increasing voltage or reducing pressure) unless you also reduce the anode length by a similar factor. Likewise if you drastically alter the value of the circuit inductance L or the anode radius a you would have to adjust other parameters accordingly. Time matching is crucial for proper computation just as it is in actual operation in the laboratory.


If machines parameters of another machine are entered, keep in the middle (or lower end) of the pressure range.  If parameters are inappropriately chosen the Time Match Safeguard will stop the programme execution and give appropriate instructions to remedy. 


Another safeguard for inappropriately low operating pressures is also incorporated into the code which stops the execution and gives you appropriate warning.


You may also need to adjust the source data of the figures. Source data has been set for a maximum of 7000 points for Figs 1 & 2, and 6000 points for the other Figures.


The Active sheet comes pre-loaded with parameters of the UNU/ICTP PFF

 Operating parameters set as 13kV in neon at 3 torr



If you have difficulties you cannot solve, please e-mail me parameters of your machine.