diff --git a/doc/manual.pdf b/doc/manual.pdf
index f2591ec8..3a47bd15 100644
Binary files a/doc/manual.pdf and b/doc/manual.pdf differ
diff --git a/doc/manual.tex b/doc/manual.tex
index 31a22416..50ce6627 100644
--- a/doc/manual.tex
+++ b/doc/manual.tex
@@ -1303,7 +1303,73 @@ \subsection{Radiation Seed Input File}
 
 The seed file can be in 2 forms - one describing a homogeneous seed of Gaussian or flat top profiles in each of the 3 spatial directions, or a list of hdf5 files containing the values of the $x$ and $y$ fields at the nodes of the radiation mesh.
 
-As with the beam file, the seed file contains 2 namelists - one specifying which of the 2 input methods is being used and the number of seeds, and the other containing the specification for the selected method.
+As with the beam file, the seed file contains 2 namelists - one specifying which of the 2 input methods is being used and the number of seeds, and the other containing the specification for the selected method. The namelist {\bf nslist} tells Puffin which of these options to expect.
+
+Please note that the mesh is only imported as is. It is unfortunately not processed or interpolated onto the Puffin mesh. It is read as-is.
+
+
+\subsubsection{Namelist: nslist}
+
+{\bf dtype}
+
+String, specifying type of radiation input. Can be either {\it simple} or {\it h5} (for importing a full mesh)
+
+{\bf nseeds}
+
+Integer. The number of electron beams to set up. For a simple beam, multiple seeds can be setup. For the mesh/h5 input, this parameter is ignored.
+
+If {\bf dtype} $=$ {\it simple}, then a namelist called {\bf slist} must be used, and if {\bf dtype} $=$ {\it h5} the namelist {\bf sh5list} must be entered.
+
+\subsubsection{Namelist: slist}
+
+For the simple seed file, one can set up a seed with a simple homogeneous distribution in every dimension. The radiation seed is assumed to have either a Gaussian or flat-top distribution profile in each dimension. The beam is therefore specified by giving the widths and so forth of the beam. Multiple radiation seeds can be defined, so each parameter is an array, with the $i^{th}$ element describing the $i^{th}$ radiation seed.
+
+
+{\bf freqf}
+
+The frequency of each seed, as a fraction of the reference resonant frequency calculated in the scaled frame parameters in the main input file.
+
+{\bf ph\_sh}
+
+The phase of each seed ($0-2\pi$).
+
+{\bf sA0\_X}
+
+The peak intensity of the $x$-polarized field vector. In dimensionless scaled units of $|A|^2$ ($\sim 1$ for FEL steady-state saturation for an electron beam with properties determined by the scaled frame), and in $Wm^{-2}$ for unscaled (SI) units.
+
+{\bf sA0\_Y}
+
+The peak intensity of the $y$-polarized field vector. In dimensionless scaled units of $|A|^2$ ($\sim 1$ for FEL steady-state saturation for an electron beam with properties determined by the scaled frame), and in $Wm^{-2}$ for unscaled (SI) units.
+
+{\bf sSigmaF}
+
+The standard deviation $\sigma$ of the intensity profile of the seed in each dimension, in the order $x, y, z2$. For multiple seeds, specify all of the $x$ values first, then all the $y$ values, then all the $\bar{z}_2$ values. If a flat top distribution is specified (see {\bf qFlatTop} below), then this is the radius \textit{i.e.} half the total length of the radiation seed in $\bar{z}_2$. The profile in the transverse $x$ and $y$ directions is always Gaussian. In scaled $\bar{x}, \bar{y}, \bar{z}_2$ (dimensionless) units for {\bf qScaled $=$ .true.} in the main input file, if {\bf .false.} then must be specified in meters for $x$ and $y$, and in seconds for the third $t$ dimension.
+
+
+{\bf qFlatTop}
+
+Array of logicals of size {\bf nseeds}, describing whether the intensity/power distribution in the longditudinal/temporal direction is flat-top ($=${\bf .true.}) or Gaussian ($=${\bf .false.}).
+
+
+{\bf meanZ2}
+
+The mean of the seed power profile in the temporal dimension. If outside the left edge of the mesh ($\bar{z}_2 = 0$), then the front edge of the seed is shifted to the front edge of the mesh. So if {\bf meanZ2} is undefined or $=0$, then the seed mean will be recalculated so that the front of the seed is aligned with the front of the mesh.
+
+
+
+{\bf qRndFj\_G}
+
+If {\bf qFlatTop} $=$ {\bf .true.}, then this controls whether the flat top profile is rounded off smoothly at the edges with a Gaussian. See {\bf sSigFj\_G} below to control the width of the gaussian tails. If {\bf qFlatTop} $=$ {\bf .false.}, then this is ignored.
+
+
+{\bf sSigFj\_G}
+
+RMS width of the Gaussian tails, in the case of a flat top radiation seed pulse.
+
+{\bf qMatchS\_G}
+
+If {\bf .true.}, then the value for the \textit{rms} $\sigma$ of the transverse field distributions above will be ignored, and instead the field distribution in the transverse plane will be matched to the electron beam transverse charge distribution.
+
 
 
 \newpage
diff --git a/inputs/simple/1D/2colour/beam_file.in b/inputs/simple/1D/2colour/beam_file.in
index 72bf5fdf..d39c45d2 100644
--- a/inputs/simple/1D/2colour/beam_file.in
+++ b/inputs/simple/1D/2colour/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -35,7 +35,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =0.1, 0.1, 1E8, 1.0, 1.0, 1E8
 sLenE = 1E-6, 1E-6, 50.0, 1E-6, 1E-6, 1e-4
-iNumElectrons = 1, 1, 73000, 1, 1, 1
+iNumMPs = 1, 1, 73000, 1, 1, 1
 emitx = 6.75599500154485e-01
 emity = 6.75599500154485e-01
 sQe = 100E-12
diff --git a/inputs/simple/1D/2colour/wchicanes/beam_file.in b/inputs/simple/1D/2colour/wchicanes/beam_file.in
index 3c67ab4f..3256db52 100644
--- a/inputs/simple/1D/2colour/wchicanes/beam_file.in
+++ b/inputs/simple/1D/2colour/wchicanes/beam_file.in
@@ -14,7 +14,7 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
 ! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
@@ -35,7 +35,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 1E8
 sLenE = 1E-6, 1E-6, 50.0, 1E-6, 1E-6, 1e-4
-iNumElectrons = 1, 1, 73000, 1, 1, 1
+iNumMPs = 1, 1, 73000, 1, 1, 1
 emitx = 6.75599500154485e-01
 emity = 6.75599500154485e-01
 sQe = 100E-12
diff --git a/inputs/simple/1D/2colour/wdrifts/beam_file.in b/inputs/simple/1D/2colour/wdrifts/beam_file.in
index ce3d65ca..3256db52 100644
--- a/inputs/simple/1D/2colour/wdrifts/beam_file.in
+++ b/inputs/simple/1D/2colour/wdrifts/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -35,8 +35,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 1E8
 sLenE = 1E-6, 1E-6, 50.0, 1E-6, 1E-6, 1e-4
-iNumElectrons = 1, 1, 73000, 1, 1, 1
-sEmit_n = 6.75599500154485e-01
+iNumMPs = 1, 1, 73000, 1, 1, 1
+emitx = 6.75599500154485e-01
+emity = 6.75599500154485e-01
 sQe = 100E-12
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/1D/CLARA/beam_file.in b/inputs/simple/1D/CLARA/beam_file.in
index 1eb6edf9..d529d39a 100644
--- a/inputs/simple/1D/CLARA/beam_file.in
+++ b/inputs/simple/1D/CLARA/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =0.084, 0.084, 1E8, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 50.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 1, 1, 10000, 1, 1, 19
+iNumMPs = 1, 1, 10000, 1, 1, 19
 emitx = 1.0
 emity = 1.0
 sQe = 250E-12
diff --git a/inputs/simple/1D/OptCommV165pp65-70/fig1/beam_file.in b/inputs/simple/1D/OptCommV165pp65-70/fig1/beam_file.in
index feda6492..26a0b385 100644
--- a/inputs/simple/1D/OptCommV165pp65-70/fig1/beam_file.in
+++ b/inputs/simple/1D/OptCommV165pp65-70/fig1/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 1.0
 sLenE = 1E-6, 1E-6, 2.0, 1E-6, 1E-6, 1E-6
-iNumElectrons = 1, 1, 800, 1, 1, 1
+iNumMPs = 1, 1, 800, 1, 1, 1
 emitx = 1.0
 emity = 1.0
 sQe = 1.0E-9
diff --git a/inputs/simple/1D/OptCommV165pp65-70/fig1/unscaled/beam_file.in b/inputs/simple/1D/OptCommV165pp65-70/fig1/unscaled/beam_file.in
index a9737b49..b525ab51 100644
--- a/inputs/simple/1D/OptCommV165pp65-70/fig1/unscaled/beam_file.in
+++ b/inputs/simple/1D/OptCommV165pp65-70/fig1/unscaled/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 1.0
 sLenE = 1E-6, 1E-6, 2.1224E-13, 1E-6, 1E-6, 1E-6
-iNumElectrons = 1, 1, 800, 1, 1, 1
+iNumMPs = 1, 1, 800, 1, 1, 1
 emitx = 1.0
 emity = 1.0
 sQe = 1.0E-9
diff --git a/inputs/simple/1D/OptCommV165pp65-70/fig2/beam_file.in b/inputs/simple/1D/OptCommV165pp65-70/fig2/beam_file.in
index 34aeac93..de3efc60 100644
--- a/inputs/simple/1D/OptCommV165pp65-70/fig2/beam_file.in
+++ b/inputs/simple/1D/OptCommV165pp65-70/fig2/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,8 +34,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE = 0.1, 0.1, 1E8, 1.0, 1.0, 1.0
 sLenE = 1E-6, 1E-6, 10.0, 1E-6, 1E-6, 1E-6
-iNumElectrons = 1, 1, 4000, 1, 1, 1
-sEmit_n = 1.0
+iNumMPs = 1, 1, 4000, 1, 1, 1
+emitx = 1.0
+emity = 1.0
 sQe = 1.0E-9
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/1D/PREv64pp046503/fig4/beam_file.in b/inputs/simple/1D/PREv64pp046503/fig4/beam_file.in
index 3a0ddedf..a76e5c86 100644
--- a/inputs/simple/1D/PREv64pp046503/fig4/beam_file.in
+++ b/inputs/simple/1D/PREv64pp046503/fig4/beam_file.in
@@ -14,7 +14,7 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam
 ! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =0.1, 0.1, 1E8, 1.0, 1.0, 1.0
 sLenE = 1E-6, 1E-6, 40.0, 1E-6, 1E-6, 1E-6
-iNumElectrons = 1, 1, 16000, 1, 1, 1
+iNumMPs = 1, 1, 16000, 1, 1, 1
 emitx = 1.0
 emity = 1.0
 sQe = 2.0E-9
diff --git a/inputs/simple/1D/PREv64pp046503/fig4ROUNDH/beam_file.in b/inputs/simple/1D/PREv64pp046503/fig4ROUNDH/beam_file.in
index c965fa05..3f9c2c83 100644
--- a/inputs/simple/1D/PREv64pp046503/fig4ROUNDH/beam_file.in
+++ b/inputs/simple/1D/PREv64pp046503/fig4ROUNDH/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =0.1, 0.1, 1E8, 1.0, 1.0, 1E-4
 sLenE = 1E-6, 1E-6, 40.0, 1E-6, 1E-6, 6E-4
-iNumElectrons = 1, 1, 32000, 1, 1, 19
+iNumMPs = 1, 1, 32000, 1, 1, 19
 emitx = 1.0
 emity = 1.0
 sQe = 2.0E-9
diff --git a/inputs/simple/1D/PREv64pp046503/fig4ROUNDP/beam_file.in b/inputs/simple/1D/PREv64pp046503/fig4ROUNDP/beam_file.in
index 8b96335f..de66a925 100644
--- a/inputs/simple/1D/PREv64pp046503/fig4ROUNDP/beam_file.in
+++ b/inputs/simple/1D/PREv64pp046503/fig4ROUNDP/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 1.0
 sLenE = 1E-6, 1E-6, 40.0, 1E-6, 1E-6, 1E-6
-iNumElectrons = 1, 1, 32000, 1, 1, 1
+iNumMPs = 1, 1, 32000, 1, 1, 1
 emitx = 1.0
 emity = 1.0
 sQe = 2.0E-9
diff --git a/inputs/simple/1D/PRSTABVv6pp070701/fig4/beam_file.in b/inputs/simple/1D/PRSTABVv6pp070701/fig4/beam_file.in
index 923001a9..04b7038f 100644
--- a/inputs/simple/1D/PRSTABVv6pp070701/fig4/beam_file.in
+++ b/inputs/simple/1D/PRSTABVv6pp070701/fig4/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -34,7 +34,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =0.1, 0.1, 1E8, 1.0, 1.0, 0.004
 sLenE = 1E-6, 1E-6, 6.0, 1E-6, 1E-6, 0.024
-iNumElectrons = 1, 1, 2400, 1, 1, 9
+iNumMPs = 1, 1, 2400, 1, 1, 9
 emitx = 1.0
 emity = 1.0
 sQe = 1.0E-9
diff --git a/inputs/simple/1D/multi-beam/beam_file.in b/inputs/simple/1D/multi-beam/beam_file.in
index deac97ff..29064f37 100644
--- a/inputs/simple/1D/multi-beam/beam_file.in
+++ b/inputs/simple/1D/multi-beam/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -36,8 +36,8 @@ sSigmaE(1,:) = 0.1, 0.1, 1E8, 1.0, 1.0, 1.0
 sSigmaE(2,:) = 0.1, 0.1, 1E8, 1.0, 1.0, 1.0
 sLenE(1,:) = 1E-6, 1E-6, 3.0, 1E-6, 1E-6, 1E-6
 sLenE(2,:) = 1E-6, 1E-6, 3.0, 1E-6, 1E-6, 1E-6
-iNumElectrons(1,:) = 1, 1, 2500, 1, 1, 1
-iNumElectrons(2,:) = 1, 1, 2500, 1, 1, 1
+iNumMPs(1,:) = 1, 1, 2500, 1, 1, 1
+iNumMPs(2,:) = 1, 1, 2500, 1, 1, 1
 emitx = 1.0, 1.0
 emity = 1.0, 1.0
 sQe = 1.0E-9, 4.0E-9
diff --git a/inputs/simple/3D/CLARA/beam_file.in b/inputs/simple/3D/CLARA/beam_file.in
index 1ed8e211..ea4486f6 100644
--- a/inputs/simple/3D/CLARA/beam_file.in
+++ b/inputs/simple/3D/CLARA/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -37,7 +37,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =7.213715e-5, 2.4323e-5, 250E-15, 1.0, 1.0, 0.0001
 sLenE = 1E-6, 1E-6, 1.75E-12, 1E-6, 1E-6, 0.0006
-iNumElectrons = 20, 20, 105000, 7, 7, 7
+iNumMPs = 20, 20, 105000, 7, 7, 7
 sQe = 2.5E-10
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/CLARA/single-slice/beam_file.in b/inputs/simple/3D/CLARA/single-slice/beam_file.in
index 6a8edf31..0cd6923d 100644
--- a/inputs/simple/3D/CLARA/single-slice/beam_file.in
+++ b/inputs/simple/3D/CLARA/single-slice/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -37,7 +37,7 @@ dtype = 'simple'
 &BLIST
 sSigmaE =7.213715e-5, 2.4323e-5, 1E8, 1.0, 1.0, 0.0001
 sLenE = 1E-6, 1E-6, 3.329781820520314e-16, 1E-6, 1E-6, 0.0006
-iNumElectrons = 20, 20, 4, 7, 7, 7
+iNumMPs = 20, 20, 4, 7, 7, 7
 sQe = 1.3319127282081256e-13 ! 1.0655301825665005e-12 !2.5E-10
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7/beam_file.in b/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7/beam_file.in
index 2b803383..6eadce40 100755
--- a/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7/beam_file.in
+++ b/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -38,8 +38,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 6.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 9, 9, 900, 9, 9, 5
-sEmit_n = 1.0
+iNumMPs = 9, 9, 900, 9, 9, 5
+emitx = 1.0
+emity = 1.0
 sQe = 1.0E-9
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7a/beam_file.in b/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7a/beam_file.in
index b2b57ef1..e817cc53 100755
--- a/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7a/beam_file.in
+++ b/inputs/simple/3D/PhyOfPlasmasV19pp093119/fig7a/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -38,8 +38,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 6.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 9, 9, 900, 9, 9, 5
-sEmit_n = 1.0
+iNumMPs = 9, 9, 900, 9, 9, 5
+emitx = 1.0
+emity = 1.0
 sQe = 1.0E-9
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/diffraction/beam_file.in b/inputs/simple/3D/diffraction/beam_file.in
index 815e25f4..d928359e 100644
--- a/inputs/simple/3D/diffraction/beam_file.in
+++ b/inputs/simple/3D/diffraction/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -37,8 +37,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 10.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 7, 7, 2050, 7, 7, 7
-sEmit_n = 1.0
+iNumMPs = 7, 7, 2050, 7, 7, 7
+emitx = 1.0
+emity = 1.0
 sQe = 1.774e-10
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/shortCSE/beam_file.in b/inputs/simple/3D/shortCSE/beam_file.in
index 56f1a4d6..0c2a0b3a 100755
--- a/inputs/simple/3D/shortCSE/beam_file.in
+++ b/inputs/simple/3D/shortCSE/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -37,8 +37,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 0.222, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 2.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 7, 7, 400, 7, 7, 7
-sEmit_n = 1.0
+iNumMPs = 7, 7, 400, 7, 7, 7
+emitx = 1.0
+emity = 1.0
 sQe = 3.5908E-11
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/single_spike_SASE/beam_file.in b/inputs/simple/3D/single_spike_SASE/beam_file.in
index 9c9031c6..71abe313 100755
--- a/inputs/simple/3D/single_spike_SASE/beam_file.in
+++ b/inputs/simple/3D/single_spike_SASE/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -37,8 +37,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 10.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 7, 7, 2050, 7, 7, 7
-sEmit_n = 1.0
+iNumMPs = 7, 7, 2050, 7, 7, 7
+emitx = 1.0
+emity = 1.0
 sQe = 1.774e-10
 bcenter = 0.0
 gammaf = 1.0
diff --git a/inputs/simple/3D/test1/beam_file.in b/inputs/simple/3D/test1/beam_file.in
index ea52d492..8157274d 100644
--- a/inputs/simple/3D/test1/beam_file.in
+++ b/inputs/simple/3D/test1/beam_file.in
@@ -14,8 +14,8 @@
 ! sSigmaE - gaussian std dev in each dimension - x, y, z2, px, py, gamma, then for additional beams
 ! sLenE   - Total length of beam modelled in each dimension - x, y, z2, px, py, gamma...
 ! bcenter - Center of beam in z2
-! iNumElectrons  -  Number of macroparticles in each dimension used to model the beam 
-! sEmit_n - Scaled transverse beam emittance
+! iNumMPs  -  Number of macroparticles in each dimension used to model the beam 
+! emitx, emity - Scaled transverse beam emittance
 ! sQe - Beam charge 
 ! gammaf - Ratio of average beam energy to reference beam energy gamma / gamma_r 
 ! chirp  - Energy chirp in z2 i.e. dgamma/dz2
@@ -37,8 +37,9 @@ dtype = 'simple'
 &BLIST
 sSigmaE =1.0, 1.0, 1E8, 1.0, 1.0, 0.001
 sLenE = 1E-6, 1E-6, 13.0, 1E-6, 1E-6, 0.006
-iNumElectrons = 20, 20, 1100, 7, 7, 7
-sEmit_n = 1.0
+iNumMPs = 20, 20, 1100, 7, 7, 7
+emitx = 1.0
+emity = 1.0
 sQe = 4.5E-9
 bcenter = 0.0
 gammaf = 1.0
diff --git a/source/FreadData.f90 b/source/FreadData.f90
index e2a56c59..49629058 100644
--- a/source/FreadData.f90
+++ b/source/FreadData.f90
@@ -280,7 +280,7 @@ subroutine read_in(zfilename, &
                  beam_file, sElectronThreshold, &
                  iNumNodesY, iNumNodesX, &
                  nodesPerLambdar, sFModelLengthX, &
-               sFModelLengthY, sFModelLengthZ2, &
+                 sFModelLengthY, sFModelLengthZ2, &
                  iRedNodesX, iRedNodesY, sFiltFrac, &
                  sDiffFrac, sBeta, seed_file, srho, &
                  sux, suy, saw, sgamma_r, &
@@ -528,7 +528,7 @@ subroutine read_in(zfilename, &
                      qMatched_A,qOKL)
 
   CALL read_seedfile(seed_file,nseeds,sSigmaF,sA0_Re,sA0_Im,freqf,&
-                     ph_sh, qFlatTopS,SmeanZ2,field_file,qOKL)
+                     ph_sh, qFlatTopS,SmeanZ2,field_file,qscaled,qOKL)
 
   call FileNameNoExtension(beam_file, zBFile_G, qOKL)
   
@@ -919,7 +919,7 @@ END SUBROUTINE read_beamfile
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 SUBROUTINE read_seedfile(se_f, nseeds,sSigmaF,sA0_X,sA0_Y,freqf,ph_sh,&
-                         qFlatTop, meanZ2,field_file,qOK)
+                         qFlatTop, meanZ2,field_file,qsc, qOK)
 
   IMPLICIT NONE
 
@@ -934,6 +934,7 @@ SUBROUTINE read_seedfile(se_f, nseeds,sSigmaF,sA0_X,sA0_Y,freqf,ph_sh,&
   INTEGER(KIND=IP), INTENT(INOUT) :: nseeds
 
   LOGICAL, ALLOCATABLE, INTENT(OUT) :: qFlatTop(:)
+  logical, intent(in) :: qsc
   LOGICAL, INTENT(OUT) :: qOK
 
 !                     LOCAL ARGS
@@ -972,7 +973,15 @@ SUBROUTINE read_seedfile(se_f, nseeds,sSigmaF,sA0_X,sA0_Y,freqf,ph_sh,&
 
 !  Default value
 
-  sSigmaF = 1.0_wp
+  if (qsc) then
+    sSigmaF(:,1) = 0.1
+    sSigmaF(:,2) = 0.1
+    sSigmaF(:,3) = 1.0
+  else
+    sSigmaF(:,1) = 50.0e-6
+    sSigmaF(:,1) = 50.0e-6
+    sSigmaF(:,1) = 1.0e-12
+  end if
   freqf = 1.0_wp
   ph_sh = 0.0_wp
   sA0_X = 0.0_wp
@@ -980,7 +989,7 @@ SUBROUTINE read_seedfile(se_f, nseeds,sSigmaF,sA0_X,sA0_Y,freqf,ph_sh,&
   qFlatTop = .false.
   meanZ2 = 0.0_wp
   qRndFj_G = .false.
-  sSigFj_G = 0.01_wp
+  sSigFj_G = sSigmaF(:,3) / 100.0_wp
   qMatchS_G = .true.
   if (dtype == 'simple') then 
     if (se_f .ne. '') then
diff --git a/source/Jsetupcalcs.f90 b/source/Jsetupcalcs.f90
index 5fea6954..3965f3c1 100644
--- a/source/Jsetupcalcs.f90
+++ b/source/Jsetupcalcs.f90
@@ -628,7 +628,7 @@ END SUBROUTINE SetUpInitialValues
 subroutine scaleParams(sEleSig, sLenEPulse, sSigEdge, &
                        beamCenZ2, chirp, sEmit, emitx, emity, gamFrac, &
                        sFieldModelLength, sLengthofElm, &
-                       sSeedSigma)
+                       sSeedSigma, sAx, sAy, scr)
 
     real(kind=wp), intent(inout) :: sEleSig(:,:), sLenEPulse(:,:), &
                                     sSigEdge(:), beamCenZ2(:), &
@@ -636,7 +636,8 @@ subroutine scaleParams(sEleSig, sLenEPulse, sSigEdge, &
                                     sFieldModelLength(:), &
                                     sLengthofElm(:), &
                                     sSeedSigma(:,:), &
-                                    emitx(:), emity(:)
+                                    emitx(:), emity(:), &
+                                    sAx(:), sAy(:), scr(:)
 
     real(kind=wp), intent(in) :: gamFrac(:)
 
@@ -699,6 +700,10 @@ subroutine scaleParams(sEleSig, sLenEPulse, sSigEdge, &
       call scaleX(sSeedSigma(is,iX_CG), lg_G, lc_G)
       call scaleX(sSeedSigma(is,iY_CG), lg_G, lc_G)
       call scaleT(sSeedSigma(is,iZ2_CG), lc_G)
+      call scaleIntensity(sAx(is), lg_G, lc_G, sGammaR_G, sKappa_G)
+      call scaleIntensity(sAy(is), lg_G, lc_G, sGammaR_G, sKappa_G)
+      call scaleT(scr(is), lc_G)
+      call scaleT(sSigFj_G(is), lc_G)
 
     end do
 
@@ -1510,8 +1515,8 @@ SUBROUTINE getSeed(NN,sig,cen,magx,magy,qFT,qRnd, &
 
 !     x and y polarized fields in z2
 
-  oscx = -z2env * cos(fr * z2nds / (2.0_WP * rho) - ph_sh)
-  oscy = z2env * sin(fr * z2nds / (2.0_WP * rho) - ph_sh)
+  oscx = z2env * sin(fr * z2nds / (2.0_WP * rho) - ph_sh)! + 4.0_wp*(cos(10_wp * z2nds)))
+  oscy = z2env * cos(fr * z2nds / (2.0_WP * rho) - ph_sh)!+ 4.0_wp*(cos(10_wp * z2nds)))
 
 !     Full 3D field
 
diff --git a/source/Msetup.f90 b/source/Msetup.f90
index 608dbb67..b3d8bb82 100644
--- a/source/Msetup.f90
+++ b/source/Msetup.f90
@@ -174,10 +174,12 @@ subroutine init(sZ, qOK)
     call scaleParams(sEleSig, sLenEPulse, sSigEj_G, &
                      beamCenZ2, chirp, sEmit_n, emitx, emity, gamma_d, &
                      sFieldModelLength, sLengthofElm, &
-                     sSeedSigma)
+                     sSeedSigma, sA0_Re, sA0_Im, SmeanZ2)
   end if
 
 
+  sA0_Re = sqrt(2.0_wp*sA0_Re)
+  sA0_Im = sqrt(2.0_wp*sA0_Im)  ! Convert intensity to peak field magnitude
 
 
   call calcSamples(sFieldModelLength, iNodes, sLengthofElm, &
diff --git a/source/scale.f90 b/source/scale.f90
index 07fd0fdb..63b828bf 100644
--- a/source/scale.f90
+++ b/source/scale.f90
@@ -12,45 +12,52 @@ module scale
 implicit none
 
 
-INTERFACE scaleG
-    MODULE PROCEDURE scaleG_single, scaleG_array
-END INTERFACE
+interface scaleG
+    module procedure scaleG_single, scaleG_array
+end interface
 
 
-INTERFACE unscaleG
-    MODULE PROCEDURE unscaleG_single, unscaleG_array
-END INTERFACE
+interface unscaleG
+    module procedure unscaleG_single, unscaleG_array
+end interface
 
 
-INTERFACE scaleX
-    MODULE PROCEDURE scaleX_single, scaleX_array
-END INTERFACE
+interface scaleX
+    module procedure scaleX_single, scaleX_array
+end interface
 
 
-INTERFACE unscaleX
-    MODULE PROCEDURE unscaleX_single, unscaleX_array
-END INTERFACE
+interface unscaleX
+    module procedure unscaleX_single, unscaleX_array
+end interface
 
 
-INTERFACE scalePx
-    MODULE PROCEDURE scalePx_single, scalePx_array
-END INTERFACE
+interface scalePx
+    module procedure scalePx_single, scalePx_array
+end interface
 
 
-INTERFACE unscalePx
-    MODULE PROCEDURE unscalePx_single, unscalePx_array
-END INTERFACE
+interface unscalePx
+    module procedure unscalePx_single, unscalePx_array
+end interface
 
 
-INTERFACE scaleT
-    MODULE PROCEDURE scaleT_single, scaleT_array
-END INTERFACE
+interface scaleT
+    module procedure scaleT_single, scaleT_array
+end interface
 
 
-INTERFACE unscaleT
-    MODULE PROCEDURE unscaleT_single, unscaleT_array
-END INTERFACE
+interface unscaleT
+    module procedure unscaleT_single, unscaleT_array
+end interface
 
+interface scaleIntensity
+    module procedure scaleIntensity_single, scaleIntensity_array
+end interface
+
+interface unscaleIntensity
+    module procedure unscaleIntensity_single, unscaleIntensity_array
+end interface
 
 contains
 
@@ -367,7 +374,7 @@ end subroutine unscaleT_array
 
 subroutine scaleZ(sZ, Lg)
 
-	real(kind=wp), intent(inout) :: sZ
+	  real(kind=wp), intent(inout) :: sZ
     real(kind=wp), intent(in) :: Lg
 
     sZ = sZ / Lg
@@ -388,4 +395,68 @@ end subroutine unscaleZ
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Scaling of I -> |A|^2
+
+subroutine scaleIntensity_single(intensity, Lg, Lc, gammar, kappa)
+
+    real(kind=wp), intent(out) :: intensity
+    real(kind=wp), intent(in) :: Lg, Lc, gammar, kappa
+
+    real(kind=wp) :: intScale
+
+    intScale = c * e_0 * ((gammar * m_e * c**2.0_wp ) / &
+                 (q_e * kappa * Lg ))**2.0_wp
+
+	  intensity = intensity / intScale
+
+end subroutine scaleIntensity_single
+
+subroutine scaleIntensity_array(intensity, Lg, Lc, gammar, kappa)
+
+    real(kind=wp), intent(out) :: intensity(:)
+    real(kind=wp), intent(in) :: Lg, Lc, gammar, kappa
+
+    real(kind=wp) :: intScale
+
+    intScale = c * e_0 * ((gammar * m_e * c**2.0_wp ) / &
+                 (q_e * kappa * Lg ))**2.0_wp
+
+	  intensity = intensity / intScale
+
+end subroutine scaleIntensity_array
+
+subroutine unscaleIntensity_single(intensity, Lg, Lc, gammar, kappa)
+
+    real(kind=wp), intent(out) :: intensity
+    real(kind=wp), intent(in) :: Lg, Lc, gammar, kappa
+
+    real(kind=wp) :: intScale
+
+    intScale = c * e_0 * ((gammar * m_e * c**2.0_wp ) / &
+                 (q_e * kappa * Lg ))**2.0_wp
+
+	  intensity = intensity * intScale
+
+end subroutine unscaleIntensity_single
+
+subroutine unscaleIntensity_array(intensity, Lg, Lc, gammar, kappa)
+
+    real(kind=wp), intent(out) :: intensity(:)
+    real(kind=wp), intent(in) :: Lg, Lc, gammar, kappa
+
+    real(kind=wp) :: intScale
+
+    intScale = c * e_0 * ((gammar * m_e * c**2.0_wp ) / &
+                 (q_e * kappa * Lg ))**2.0_wp
+
+	  intensity = intensity * intScale
+
+end subroutine unscaleIntensity_array
+
+!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+
 end module scale