# collision_hs1.sl # A hard-sphere, elastic, ion-neutral collision model for SIMION 7 (or 6). # # Note to SIMION 8 Users: A slightly more updated and documented # Lua version of this code is included in SIMION 8 # (collision_hs1 example). It's recommended you use that version instead. # # The code implements a rather complete hard-sphere collision model. # Collision models are useful for simulating non-vacuum conditions, in # which case ions collide against a background gas and are deflected # randomly. # # Features and assumptions of the model: # - Ion collisions follow the hard-sphere collision model. # - Ion collisions are elastic. # - Background gas is assumed neutral in charge. # - Background gas velocity follows the Maxwell-Boltzmann distribution. # - Background gas mean velocity may be non-zero. # - Kinetically cooling and heating collisions are simulated. # - Background gas as a whole is unaffected by ion collisions. # # Note on time-steps: each individual ion-gas collision is modeled, # which requires the time-step to be some fraction of mean-free-path. # Therefore, simulations with frequent collisions (i.e. higher # pressure) can be computationally intensive. # # This code does not account for absorptions (e.g. when using electrons # rather than ions). That can be easily supported by setting ion_splat, # likely as a function of impact_offset. # # The code has been influenced by a variety of prior SIMION hard-sphere # collision models: # # [Dahl] _Trap/INJECT.PRG in SIMION 7.0 # [Dahl2] http://www.simion.com/examples/dahl_drag.prg # [Appelhans2001] http://dx.doi.org/10.1016/S1387-3806(02)00627-9 # [Ding2002] http://dx.doi.org/10.1016/S1387-3806(02)00921-1 # [Ling1997] http://dx.doi.org/10.1002/(SICI)1097-0231(19970830)11:13<1467::AID-RCM54>3.0.CO;2-X # # HISTORY: # REV-4 made this correction: # Issue I362 - HS1 collision model does not accurately thermalize. # http://www.simion.com/issue/362 # # David Manura. REV-4 (200702) # (c) Scientific Instrument Services, Inc. 2006-2007. # Mean free path (MFP) (mm) # Set to -1 to calculate this automatically from pressure and temperature. adjustable _mean_free_path_mm = -1 # Mean number of time steps per MFP. adjustable _steps_per_MFP = 20.0 # Mass of background gas particle (amu) adjustable _gas_mass_amu = 4.0 # Background gas temperature (K) adjustable _temperature_k = 273.0 # Collision marker flag. # If non-zero, markers will be placed at the collisions. adjustable _mark_collisions = 1 # Background gas pressure (Pa) # Note: (Pa/mtorr) = 0.13328. adjustable _pressure_pa = 0.53 # Collision-cross section (m^2) # (The collision diameter is roughly the sum of the diameters # of the colliding ion and buffer gas particles.) # (2.1E-19 is roughly for two Helium atoms--Atkins1998-Table 1.3) adjustable _sigma_m2 = 2.27E-18 # Mean background gas velocity (mm/usec) in x,y,z directions. # Normally, these are zero. adjustable _vx_bar_gas_mmusec = 0 adjustable _vy_bar_gas_mmusec = 0 adjustable _vz_bar_gas_mmusec = 0 #-- Internal variables # Statistics adjustable[100] ion_ke_totals adjustable[100] ion_distance_totals adjustable[100] ion_collision_totals # Currently used mean-free path (-1 = undefined). static effective_mean_free_path_mm = -1 # Last known ion speed (-1 = undefined). static last_speed_ion = -1 # Last known ion number (-1 = undefined). static last_ion_number = -1 # Error function (erf). # erf(z) = (2/sqrt(PI)) * integral[0..z] exp(-t^2) dt # This algorithm is quite accurate. It is based on # "Q-Function Handout" by Keith Chugg: # http://tesla.csl.uiuc.edu/~koetter/ece361/Q-function.pdf # See also http://www.theorie.physik.uni-muenchen.de/~serge/erf-approx.pdf # I also find that the following makes a reasonable approximation: # 1 - exp(-(2/sqrt(PI))x - (2/PI)x^2) sub erf(z) returns(res) z2 = abs(z) t = 1 / (1 + 0.32759109962 * z2) res = ( - 1.061405429 ) * t res = (res + 1.453152027 ) * t res = (res - 1.421413741 ) * t res = (res + 0.2844966736) * t res =((res - 0.254829592 ) * t) * exp(-z2*z2) res = res + 1 if z < 0 res = -res endif endsub sub initialize endsub sub tstep_adjust # Ensure time steps are sufficiently small. They should be some # fraction of mean-free-path so that collisions are not missed. if effective_mean_free_path_mm > 0 # Get speed (speed, unused, unused) = rect3d_to_polar3d(ion_vx_mm, ion_vy_mm, ion_vz_mm) # Limit time-step tmax = effective_mean_free_path_mm / speed / _steps_per_MFP if ion_time_step > tmax ion_time_step = tmax endif endif unused = unused endsub sub other_actions if _pressure_pa == 0 # if collisions disabled effective_mean_free_path_mm = -1 exit endif # Define constants k = 1.3806505e-23 # Boltzmann constant (J/K) # R = 8.3145 # Ideal gas constant (J/(mol*K)) kg_amu = 1.6605402e-27 # (kg/amu) conversion factor PI = 3.1415926535 # PI constant #eV_J = 6.2415095e+18 # (eV/J) conversion factor # Translate ion velocity (mm/us) frame of reference such # that mean background gas velocity is zero. # This simplifies the subsequent analysis. vx = ion_vx_mm - _vx_bar_gas_mmusec vy = ion_vy_mm - _vy_bar_gas_mmusec vz = ion_vz_mm - _vz_bar_gas_mmusec # Convert ion velocity vector to polar coordinates (mm/us). speed_ion = sqrt(vx*vx + vy*vy + vz*vz) if speed_ion < 1E-7 speed_ion = 1E-7 # prevent divide by zero and such effects later on endif #=== Notes on Mean-Free-Path === # The mean-free-path (lambda) is generally a function of the mean # relative speed (c_bar_rel) between the ion and moving background # gas: # lambda = (c_ion/c_bar_rel) / (n * sigma) # where # c_ion is ion speed # c_bar_rel is the mean relative speed between ion and moving # background gas # n is the number of gas particles per unit volume # sigma is the collision cross section (roughly, the area with # a diameter equal to the sum of the diameters of the # colliding ion and gas particles) # # Mean relative speed (c_bar_rel) is formally determined by # c_bar_rel = tripple_integral |v_ion - v_gas| * f(v_gas) dv_gas # where # f(v) is the product of the one-dimensional Maxwell # in all three dimensions: # f(v) = (m_gas / (2*pi*k*T))^(3/2) * exp(-m*v^2 / (2*k*T)) # which evaluates to # c_bar_rel = c_bar_gas * (s + 1/(2*s)) * 0.5 * sqrt(PI) * erf(s) + # (1/2)*exp(-s^2)) (s > 0) # where # s = c_ion/c_star_gas # c_bar_gas = sqrt(8*k*T_b/pi/m_b) is the mean gas speed # c_star_gas = sqrt(2*k*T_b/m_b) is the median gas speed # # This approach is recommended by Ding2002. # # Ling1997 uses a simpler (and almost as suitable) # approximation for c_bar_rel: # c_bar_rel ~= sqrt(c_ion^2 + c_bar_gas^2) #=== # Compute effective mean-free-path (or use a specific value). if _mean_free_path_mm < 0 #...using current ion velocity # Note: only recompute mean-free-path if speed_ion has # changed significantly. This is intended to speed up the # calculation a bit. Handles flying ions by groups too. if last_ion_number != ion_number or abs(speed_ion / last_speed_ion - 1) > 0.05 # Compute mean gas speed (mm/us) c_bar_gas = sqrt(8*k*_temperature_k/PI/(_gas_mass_amu * kg_amu)) / 1000 # Compute median gas speed (mm/us)z c_star_gas = sqrt(2*k*_temperature_k/(_gas_mass_amu * kg_amu)) / 1000 # Compute mean relative speed (mm/us) s = speed_ion / c_star_gas c_bar_rel = c_bar_gas * ( (s + 1/(2*s)) * 0.5 * sqrt(pi) * erf(s) + 0.5 * exp(-s*s)) # Compute mean-free-path (mm) effective_mean_free_path_mm = 1000 * k * _temperature_k * (speed_ion / c_bar_rel) / (_pressure_pa * _sigma_m2) last_speed_ion = speed_ion last_ion_number = ion_number endif else #...using specified value effective_mean_free_path_mm = _mean_free_path_mm endif #print("DEBUG:ion[c=#],gas[c_bar=#],c_bar_rel=#, MFP=#", # speed_ion, c_bar_gas, c_bar_rel, effective_mean_free_path_mm) # Compute collision probability per distance traveled collision_prob = 1 - exp(- speed_ion * ion_time_step / effective_mean_free_path_mm) # Was there a collision? if rand() > collision_prob exit endif #--- collisions # Compute standard deviation of gas velocity in one dimension (mm/us). # The following is derived from kinetic gas theory. vr_stdev_gas = sqrt(k * _temperature_k / (_gas_mass_amu * kg_amu)) / 1000 # Compute a normalized Gaussian random variable (-inf, +inf). # This uses the Box-Muller algorithm. # rand1-3 are Gaussian random variables. s = 1 while s >= 1 v1 = 2*rand() - 1 v2 = 2*rand() - 1 s = v1*v1 + v2*v2 endwhile # (assume divide by zero improbable?) rand1 = v1*sqrt(-2*ln(s) / s) rand2 = v2*sqrt(-2*ln(s) / s) s = 1 while s >= 1 v1 = 2*rand() - 1 v2 = 2*rand() - 1 s = v1*v1 + v2*v2 endwhile rand3 = v1*sqrt(-2*ln(s) / s) vx_gas = rand1 * vr_stdev_gas vy_gas = rand2 * vr_stdev_gas vz_gas = rand3 * vr_stdev_gas # Or a slightly more correct thing might be to make probability # of (vx_gas,vy_gas,vz_gas) proportional to # |v_gas - v_ion| as well (see Lua version) # 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. # The probability of a given impact_offset is made # proportional to impact_offset^2. # Note: 0 is a head-on collision; 1 would be a near miss. impact_offset = sqrt(0.999999999 * rand()) # Compute randomized impact angle [0, +PI/2) (radians) # between ion velocity vector and radial axis. # Note: 0 is a head-on collision; +PI/2 would be a near miss. impact_angle = asin(impact_offset) # Compute randomized angle [0, 2*PI] for rotation of collision # plane around radial axis. # Note: all angles are equally likely. impact_theta = 2*PI*rand() # Compute polar coordinates in current velocity reference frame. (speed_ion_r, az_ion_r, el_ion_r) = rect3d_to_polar3d(vx, vy, vz) # Compute ion velocity components (mm/us). # Note: this choice of coordinates ensures that the vector is # always in the first (+/+) quadrant. vr_ion = speed_ion_r * cos(impact_angle) #.. radial velocity 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. vr_ion2 = (vr_ion * (ion_mass - _gas_mass_amu)) / (ion_mass + _gas_mass_amu) # Rotate velocity frame of reference 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-degrees(impact_angle), vr_ion2, vt_ion, 0) # Rotate velocity frame of reference around +y axis. # This rotates the deflection angle and in effect chooses the # collision plane (impact_theta), which was left unchosen before. (vx, vy, vz) = azimuth_rotate(degrees(impact_theta), vx, vy, vz) # Rotate velocity frame of reference back to the original. (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 frame of reference back to the original. # This undoes the prior two translations that make velocity # relative to the colliding gas particle. vx = vx + vx_gas + _vx_bar_gas_mmusec vy = vy + vy_gas + _vy_bar_gas_mmusec vz = vz + vz_gas + _vz_bar_gas_mmusec # Set new velocity vector. (ion_vx_mm, ion_vy_mm, ion_vz_mm) = (vx, vy, vz) # Now lets compute some statistics... # Calculate new ion speed and KE. (speed_ion2, unused, unused) = rect3d_to_polar3d(ion_vx_mm, ion_vy_mm, ion_vz_mm) ke2_ion = speed_to_ke(speed_ion2, ion_mass) # Compute mean gas KE #ke_bar_gas = ( # (3/2) * k * _temperature_k + # (1/2) * (_gas_mass_amu * kg_amu) * ( # _vx_bar_gas_mmusec*_vx_bar_gas_mmusec + # _vy_bar_gas_mmusec*_vy_bar_gas_mmusec + # _vz_bar_gas_mmusec*_vz_bar_gas_mmusec # ) * 1e+6 #) * eV_j #print("DEBUG:ion[ke=#],gas[ke=#]", ke2_ion, ke_bar_gas) # Record KE after collisions. This is later used to compute average KE. if ion_number <= 100 ion_ke_totals[ion_number] = ion_ke_totals[ion_number] + ke2_ion ion_collision_totals[ion_number] = ion_collision_totals[ion_number] + 1 endif if _mark_collisions != 0 mark() endif unused = unused endsub sub terminate # Display some statistics. # Note: At equilibrium, the ion and gas KE become roughly equal. if ion_number <= 100 k = 1.3806505e-23 # Boltzmann constant (J/K) eV_J = 6.2415095e+18 # (eV/J) conversion factor ke_bar = ion_ke_totals[ion_number] / (ion_collision_totals[ion_number] + 1E-10) T_bar = ke_bar / eV_J / (1.5 * k) print("ion=#, collisions=#, mean KE=# eV, mean T=# K", ion_number, ion_collision_totals[ion_number], ke_bar, T_bar) endif endsub #sub efield_adjust # # For testing, apply a quadratic potential well # # to trap ions in. The kinetic cooling of the buffer # # gas causes ions to collect near the center of the well. # # V(x,y,z) = x*x + y*y* + z*z = r*r # # E(x,y,z) = -(2*x, 2*y, 2*z) # r_max = 100 # radius # V_max = 10 # voltage at r_max # a = 2 * V_max / (r_max * r_max) # ion_dvoltsx_gu = ion_px_gu * a # ion_dvoltsy_gu = ion_py_gu * a # ion_dvoltsz_gu = ion_pz_gu * a #endsub