--[[
This is a lua workbench file of a quadrupole cylindrical ion trap.
This lua exists out of three parts, namely: The changing electric field, the drag due to the gas, and thermal fluctuations.
This lua imports three saperate lua's for each of the effects as given above.
Designed by Salem Najja, bachelor student at the radboud university, student number: s4331648
--]]


--simion.workbench_program()

--[[
===================================================================================
The electric field.
===================================================================================
==]]

local a = 0.1
local seed = 80
math.randomseed(seed)
local       R0  =  4.172E-3           --Radius of the effective area in which the ions can move freely mm
--local      mc  =  ion_mass            --mass ions
adjustable   PE_U_E_S  =   00.05    --updating each micro second.
local    Last_PE_  =   00.00        --Starting time
adjustable T_MAX = 500.000000  -- maximum aantal microseconden dat we willen meten.adjustable 
adjustable T_MAX_RIT = 50.000000  -- maximum aantal microseconden dat we willen meten.




function segment.other_actions()
	if
		abs(ion_time_of_flight - Last_PE_) >= PE_U_E_S
	then
		Last_PE_    =  ion_time_of_flight
		sim_update_pe_surface =  1
	end
end


--[[
===================================================================================
Drag solid sphere
===================================================================================
==]]
local      SSC         	=             {}
SSC.segment            	=             {}
local       damping    	=         1         	--Drag ON=1/OFF=0
local       T          	=        293    	     	--temperatuur in gass chamber
local       kb         	=        1.381E-23      --bolzmanconstant Joules/Kelvin
local       R_gas      	=       188             --radius gas in picometer (Argon)
local       M_gas      	=       39.95           --mass in amu.
local       kg_amu     	=       1.6605402e-27	-- (kg/amu) conversion factor
local       eV_J       	=       6.2415095e+18   -- (eV/J) conversion factor
local    	R  			=  (481.36/2+R_gas)   --1018/2
local    	l  			=    -1

local    	sigma  		=  math.pi*(R*1E-12)^2

adjustable   SpMFP  =     20.0

adjustable   G_vx  =   0.0      --(mm/microseconden)
adjustable   G_vy  =   0.0      --(mm/microseconden)
adjustable   G_vz  =   0.0      --(mm/microseconden)


--[[
-----------------------------------------
Desired registration and some other stuff
-----------------------------------------
--]]

adjustable _mark_collisions	=    0      -- collision markers: 0= off / 1 = on
-- How much trace data (average KE) to output.
-- (0=none, 1=at each splat, 2=at each collision)
adjustable _trace_level   	=    0
--  adjustable _trace_skip		=    100      -- This reduces the verbosity of the trace.


local  function erf(z)
 local z_2     =   abs(z)
 local   t      =   (1+0.32759109962*z_2)^(-1)
 local res     =   -1.061405429*t
   res     =    (res+1.453152027 ) * t
   res     =    (res + 0.2844966736) * t
   res     =   ((res - 0.254829592 ) * t) * exp(-(z_2)^2)
   res     =   res + 1
    if z < 0 then res = -res end
    return res
end

local function Gaus_rand()
 local   s     =   1
 local   v_1, v_2

 while   s       >=   1
 do
     v_1    =   2*math.random() -1
     v_2    =   2*math.random() -1
     s     =    (v_1)^2+(v_2)^2
 end
 local rand_1    =   v_1*sqrt(-2*ln(s)/s)
 return  rand_1
end

local  KE_AV    =   {}
local  L_col_time   =   {}

function segment.initialize()

end
SSC.segment.initialize = segment.initialize

function segment.tstep_adjust()
  if
   ion_time_of_flight<T_MAX
  then
   ion_time_step   = min(ion_time_step, T_MAX - ion_time_of_flight)
  end
end
SSC.segment.tstep_adjust  =  segment.tstep_adjust
local is_splat_next = false
function segment.other_actions()

 if
  is_splat_next then ion_splat = -1

 end
 if ion_time_of_flight >= T_MAX then
  is_splat_next = true
 end

 if
  ion_splat ~= 0 then
  is_splat_next = false end

 if ion_instance == 3 then
 -- Temporarily translate ion velocity (mm/us) frame of
    -- reference such that mean background gas velocity is zero.
    -- This simplifies the subsequent analysis.
	if ion_time_of_flight >= T_MAX_RIT then
		is_splat_next = true
	end
	local vx = ion_vx_mm - G_vx
    local vy = ion_vy_mm - G_vy
    local vz = ion_vz_mm - G_vz

    -- Obtain ion speed (relative to mean background gas velocity).
    local speed_ion = sqrt(vx^2 + vy^2 + vz^2)
    if speed_ion < 1E-7 then
         speed_ion = 1E-7  -- prevent divide by zero and such effects later on
    end

    -- Compute mean-free-path.
    -- > See notes.pdf for discussion on the math.
    if l > 0 then -- explicitly specified
        effective_mean_free_path_mm = l
    else  -- calculate from current ion velocity
        -- Note: in some cases, mean-free-path will not change significantly, so
        -- we don't need to recompute it on every time step.  But it is simpler
        -- and less error prone to do so and doesn't affect run times much.
        if last_ion_number ~= ion_number or
                abs(speed_ion / last_speed_ion - 1) > 0.05  -- changed
        then
            -- Compute mean gas speed (mm/us)
            local c_bar_gas = sqrt(8*kb*T/math.pi/(M_gas * kg_amu)) / 1000

            -- Compute median gas speed (mm/us)
            local c_star_gas = sqrt(2*kb*T/(M_gas * kg_amu)) / 1000

            -- Compute mean relative speed (mm/us) between gas and ion.
            local s = speed_ion / c_star_gas
            local c_bar_rel = c_bar_gas * ((s + 1/(2*s)) * 0.5 * sqrt(math.pi) * erf(s) + 0.5 * exp(-s*s))
            
			-- Compute mean-free-path (mm)
            -- print(p)
			
			
			effective_mean_free_path_mm = 1000 * kb * T * speed_ion  / (c_bar_rel * p * sigma)
            --effective_mean_free_path_mm = (speed_ion / c_bar_rel) / (p * sigma)
			

            -- Store data about this calculation.
            last_speed_ion = speed_ion
            last_ion_number = ion_number

            --print("DEBUG:ion[c],gas[c_bar],c_bar_rel,MFP=",
            --      speed_ion, c_bar_gas, c_bar_rel, effective_mean_free_path_mm)

            -- Note: The following is a simpler and almost as suitable
            -- approximation for c_bar_rel, which you may used instead:
            -- c_bar_rel = sqrt(speed_ion^2 + c_bar_gas^2)
        end

     end
 
  -- Limit time-step size to a fraction of the MFP.
    max_timestep = effective_mean_free_path_mm / speed_ion / SpMFP

    -- Compute probability of collision in current time-step.
    -- > For an infinitesimal distance (dx) traveled, the increase in the
    --   fraction (f) of collided particles relative to the number
    --   of yet uncollided particles (1-f) is equal to the distance
    --   traveled (dx) over the mean-free-path (lambda):
    --     df/(1-f) = dx / lambda
    --   Solving this differential equation gives
    --     f = 1 - exp(- dx / lambda) = 1 - exp(- v dt / lambda)
    --   This f can be interpreted as the probability that a single
    --   particle collides in the distance traveled.
    local collision_prob = 1 -
        exp(- speed_ion * ion_time_step / effective_mean_free_path_mm)

    -- Test for collision.
    if math.random() > collision_prob then
        return -- no collision
    end

    ----- Handle collision.

    -- Compute standard deviation of background gas velocity in
    -- one dimension (mm/us).
    -- > From kinetic gas theory (Maxwell-Boltzmann), velocity in
    --   one dimension is normally distributed with standard
    --   deviation sqrt(kT/m).
    local vr_stdev_gas =
        sqrt(kb * T / (M_gas * kg_amu)) / 1000

    -- Compute velocity of colliding background gas particle.
    -- > For the population of background gas particles that collide with the
    --   ion, their velocities are not entirely Maxwell (Gaussian) but
    --   are also proportional to the relative velocities the ion and
    --   background gas particles:
    --     p(v_gas) = |v_gas - v_ion| f(v_gas)
    --   See notes.pdf for discussion.
    -- > To generate random velocities in this distribution, we may
    --   use a rejection method (http://en.wikipedia.org/wiki/Rejection_sampling)
    --   approach:
    --   > Pick a gas velocity from the Maxwell distribution.
    --   > Accept with probability proportional to its
    --     speed relative to the ion.
    local vx_gas, vy_gas, vz_gas -- computed velocities
    -- > scale is an approximate upper-bound for "len" calculated below.
    --   We'll use three standard deviations of the three dimensional gas velocity.
    local scale = speed_ion + vr_stdev_gas * 1.732 * 3  --sqrt(3)=~1.732
    repeat
        vx_gas = Gaus_rand() * vr_stdev_gas
        vy_gas = Gaus_rand() * vr_stdev_gas
        vz_gas = Gaus_rand() * vr_stdev_gas
        local len = sqrt((vx_gas - vx)^2 + (vy_gas - vy)^2 + (vz_gas - vz)^2)
        --assert(len <= scale) -- true at least ~99% of the time.
    until math.random() < len / scale

    -- Alernately, for greater performance and as an approximation, you might
    -- replace the above with a simple Maxwell distribution:
    --  vx_gas = Gaus_rand() * vr_stdev_gas
    --  vy_gas = Gaus_rand() * vr_stdev_gas
    --  vz_gas = Gaus_rand() * vr_stdev_gas

    -- Translate velocity reference frame so that colliding
    -- background gas particle is stationary.
    -- > This simplifies the subsequent analysis.
    vx = vx - vx_gas
    vy = vy - vy_gas
    vz = vz - vz_gas

    -- > Notes on collision orientation
    --   A collision of the ion in 3D can now be reasoned in 2D since
    --   the ion remains in some 2D plane before and after collision.
    --   The ion collides with an gas particle initially at rest (in the
    --   current velocity reference frame).
    --   For convenience, we define a coordinate system (r, t) on the
    --   collision plane.  r is the radial axis through the centers of
    --   the colliding particles, with the positive direction indicating
    --   approaching particles.  t is the tangential axis perpendicular to r.
    --   An additional coordinate theta defines the the rotation of the
    --   collision plane around the ion velocity axis.

    -- Compute randomized impact offset [0, 1) as a fraction
    -- of collisional cross-section diameter.
    -- 0 is a head-on collision; 1 would be a near miss.
    -- > You can imaging this as the gas particle being a stationary
    --   dart board of radius 1 unit (representing twice the actual radius
    --   of the gas particle) and the ion center is a dart
    --   with velocity perpendicular to the dart board.
    --   The dart has equal probability of hitting anywhere on the
    --   dart board.  Since a radius "d" from the center represents
    --   a ring with circumference proportional to "d", this implies
    --   that the probability of hitting at a distance "d" from the
    --   center is proportional to "d".
    -- > Formally, the normalized probability density function is
    --   f(d) = 2*d for d in [0,1].  From the fundamental transformation
    --   law of probabilities, we have
    --   integral[0..impact_offset] f(d) dd = impact_offset^2 = U,
    --   where U is a uniform random variable.  That is,
    --   impact_offset = sqrt(U).  Decrease it it slightly
    --   to prevent overflow in asin later.
    local impact_offset = sqrt(0.999999999 * math.random())

    -- Convert impact offset to impact angle [0, +pi/2) (radians).
    -- Do this since the target is really a sphere (not flat dartboard).
    -- This is the angle between the relative velocity
    -- between the two colliding particles (i.e. the velocity of the dart
    -- imagined perpendicular to the dart board) and the r axis
    -- (i.e. a vector from the center of the gas particle to the location
    -- on its surface where the ion hits).
    -- 0 is a head-on collision; +pi/2 would be a near miss.
    local impact_angle = asin(impact_offset)

    -- In other words, the effect of the above is that impact_angle has
    -- a distribution of p(impact_angle) = sin(2 * impact_angle).

    -- Compute randomized angle [0, 2*pi] for rotation of collision
    -- plane around radial axis.  The is the angle around the
    -- center of the dart board.
    -- Note: all angles are equally likely to hit.
    -- The effect is that impact_theta has a distribution
    -- of p(impact_theta) = 1/(2*pi).
    local impact_theta = 2*math.pi*math.random()

    -- Compute polar coordinates in current velocity reference frame.
    local speed_ion_r, az_ion_r, el_ion_r = rect3d_to_polar3d(vx, vy, vz)

    -- Compute ion velocity components (mm/us).
    local vr_ion = speed_ion_r * cos(impact_angle)    -- radial velocity
    local vt_ion = speed_ion_r * sin(impact_angle)    -- normal velocity

    -- Attenuate ion velocity due to elastic collision.
    -- This is the standard equation for a one-dimensional
    -- elastic collision, assuming the other particle is initially at rest
    -- (in the current reference frame).
    -- Note that the force acts only in the radial direction, which is
    -- normal to the surfaces at the point of contact.
    local vr_ion2 = (vr_ion * (ion_mass - M_gas))
                  / (ion_mass + M_gas)

    -- Rotate velocity reference frame so that original ion velocity
    -- vector is on the +y axis.
    -- Note: The angle of the new velocity vector with respect to the
    -- +y axis then represents the deflection angle.
    vx, vy, vz = elevation_rotate(90 - deg(impact_angle), vr_ion2, vt_ion, 0)

    -- Rotate velocity reference frame around +y axis.
    -- This rotates the deflection angle and in effect selects the
    -- randomized impact_theta.
    vx, vy, vz = azimuth_rotate(deg(impact_theta), vx, vy, vz)

    -- Rotate velocity reference frame back to the original.
    -- For the incident ion velocity, this would have the effect
    -- of restoring it.
    vx, vy, vz = elevation_rotate(-90 + el_ion_r, vx, vy, vz)
    vx, vy, vz = azimuth_rotate(az_ion_r, vx, vy, vz)

    -- Translate velocity reference frame back to original.
    -- This undoes the prior two translations that make velocity
    -- relative to the colliding gas particle.
    vx = vx + vx_gas + G_vx
    vy = vy + vy_gas + G_vy
    vz = vz + vz_gas + G_vz

    -- Set new velocity vector of deflected ion.
    ion_vx_mm, ion_vy_mm, ion_vz_mm = vx, vy, vz

    -- Now lets compute some statistics...

    -- Compute running average of KE.  This is for statistical reporting only.
    -- At thermal equilibrium, KE of the ion and KE of the gas would
    -- be approximately equal according to theory.
  --[[  if _trace_level  ~= 0 then record_ke_other_actions() end

    if _mark_collisions ~= 0 then
        mark() -- draw dot at collision point
    end
 --]]
     if _trace_level >= 1 then
        -- Compute new ion speed and KE.
        local speed_ion2 = sqrt(ion_vx_mm^2 + ion_vy_mm^2 + ion_vz_mm^2)
        local ke2_ion = speed_to_ke(speed_ion2, ion_mass)

        -- To average ion KE somewhat reliably, we do a running (exponential decay)
        -- average of ion KE over time.  The reset time of the exponential decay
        -- is set to some fraction of the total time-of-flight, so the average
        -- will become more steady as the run progresses (we assume this is a
        -- system that approaches equilibrium).
        -- Note: exp(-x) can be approximated with 1-x for small x.

        -- time between most recent collisions
        local dt = ion_time_of_flight - (last_collision_times[ion_number] or 0)
        -- average over some fraction of TOF
        reset_time = ion_time_of_flight * 0.5
        -- weight for averaging.
        local w = 1 - (dt / reset_time)  -- ~= exp(-dt / reset_time)
        -- update average ion KE
        KE_AV[ion_number] = w * (KE_AV[ion_number] or ke2_ion)
                                + (1-w) * ke2_ion
        if _trace_level >= 2 then -- more detail
            local T_ion = KE_AV[ion_number] / eV_J / (1.5 * kb)
            if trace_count % _trace_skip == 0 then
                print(string.format(
                    "n=,%d,TOF=,%0.3g,ion KE (eV)=,%0.3e,ion mean KE (eV)=," ..
                    "%0.3e,ion mean temp (K)=,%0.3e",
                    ion_number, ion_time_of_flight, ke2_ion,
                    KE_AV[ion_number], T_ion))
            end
            trace_count = (trace_count + 1) % _trace_skip
        end
        last_collision_times[ion_number] = ion_time_of_flight
    end

    if _mark_collisions ~= 0 then
        mark() -- draw dot at collision point
    end
end
end
SSC.segment.other_actions = segment.other_actions

function segment.accel_adjust()


    ion_ax_mm = (ion_ax_mm - damping * (ion_vx_mm - G_vx)) -- ion_v[xyz]_mm is particle velocity in mm/usec; ion_a[xyz]_mm is particle acceleration in mm/usec^2.
    ion_ay_mm = (ion_ay_mm - damping * (ion_vy_mm - G_vy))
    ion_az_mm = (ion_az_mm - damping * (ion_vz_mm - G_vz))
end