;Collision model using a Maxwell-Boltzman distribution for the ;collision gas energy and accomodating collisions from behind ;the ion ;Created by Dave Dahl defa _Mean_Free_Path_mm 4.0 ; mean free path (MFP) in mm defa _Steps_per_MFP 20.0 ; avg number of int steps in MFP defa _Collision_Gas_Mass_amu 4.0 ; assume helium defa _Temperature_K 273.0 ; gas temperature defa _Mark_Collisions 0 ; mark collisions if not equal zero adefa v0_dist 1200 ;Maxwell-Boltzmann v distribution array adefa totals 500 adefa counts 500 seg initialize rcl ion_number 1 xy goto skip ;skip recording if total below next record rcl recordstep + ;inc nextrecord by recordstep sto nextrecord ;save next recording objective ; rcl total ;verifying printout ; rcl v ; message ; v = #, total = # rcl v ;rcl v at current step rcl index 1 + sto index ;inc array index asto v0_dist ;save v into array at index lbl skip ;skipped recording entry rcl v ;get current velocity rcl stepsize + ;increment it by a small step sto v ;store v's new value rcl total ;get current total rcl maxrecord ;get max total to stay in loop x>=y goto loop ;keep looping if total is less message ;Maxwell-Boltzmann dist array Done! lbl skip_init ;continue any normal initialize functions exit seg tstep_adjust ; insure time steps are short enough rcl ion_vz_mm ; get ion's speed rcl ion_vy_mm ; load velocity components rcl ion_vx_mm >p3d ; convert velocity to polar coords sto v ; save speed in temporary variable rcl _Mean_Free_Path_mm ;get MFP rcl v / rcl _Steps_per_MFP / ;time step for MFP / _Steps_per_MFP sto tmax ;remember rcl ion_time_step ;get proposed time step xp3d ; convert velocity to polar coords sto v ; save in temporary variables rlup sto az rlup sto el ; key? ; x=0 goto skipkey ; rcl ion_mass ; rcl v ; >ke ; message ;ion ke = # eV ; lbl skipkey ; test for collision rcl v rcl ion_time_step * ; compute distance from tstep * v rcl effective_free_path / chs e^x 1 x<>y - ;(1-e(-d/fp)) rand ;get random number from 0 - 1 x>y goto skip1 ; skip if no collision ;*********** collision model ************ ; collision -- assume variable position hit 0.999999999 rand * ; insure collision (r < 1.0) abs x>0 sqrt ; equal area weighted r value asin ; collision angle in radians sto impact_angle_rad ; save in local variable rcl v rcl impact_angle_rad cos * sto vr ; ion's impact velocity rcl v rcl impact_angle_rad sin * sto vt ; ion's normal velocity rcl _temperature_k 3 * 1.380658e-23 * ;get thermal v0 of collion gas molecule rcl _Collision_Gas_Mass_amu abs ;get abs mass of collision gas 1.6605402e-27 * / ;vrms = sqrt(3kT/m) sqrt 1000 / ;root mean squared thermo velocity in mm/usec sto v0thermal ;store as thermal v0 rand 998 * 1 + sto n ;get array pointer fp index (1-999) int sto low_index ;int and save as low array index 1 + ;inc for high array index arcl v0_dist ;get high v array factor rcl low_index ;get low array index arcl v0_dist sto low_factor ;get low v array factor - ;delv = (vhigh - vlow) rcl n frac * rcl low_factor + ;get interpolated v value rcl v0thermal * ;scale with v0thermal entr * sto v_squared ;save velocity squared of ion ;vx will hold impact ke vsquared rand 0.5 - sto vx ;rand vx (-0.5 to +0.5) rand 0.5 - sto vy ;rand vy (-0.5 to +0.5) rand 0.5 - sto vz ;rand vy (-0.5 to +0.5) rcl vx abs rcl vy abs rcl vz abs + + ;sum of abs vals of vx, vy, and vz sto vsum ;sto sum of abs vals as vt rcl vx rcl vsum / rcl v_squared * sto vx ;signed vx component of ke x=0 1 sto temp ;save 1 if vx =0 else save vx rcl vx abs x=0 1 ;save 1 if abs vx = 0 else save abs vx rcl temp / sto vx_sign ;get sign of vx component (-1 or +1) rcl vx abs sqrt rcl vx_sign * sto vx ;compute impact velocity of gas molecule ;Note: if vx is greater than vr -- a hit from behind is assumed rcl ion_mass rcl _collision_gas_mass_amu - rcl ion_mass rcl _collision_gas_mass_amu + / ; (mi - mc)/(mi + mc) rcl vr * sto vr ; attenuated impact velocity rcl _collision_gas_mass_amu entr + rcl ion_mass rcl _collision_gas_mass_amu + / ; (mc + mc)/(mi + mc) rcl vx * rcl vr + sto vr ; add to vr component rcl vr rcl vr * rcl vt rcl vt * + sqrt sto v ; compute resulting velocity 1.0e-10 * sto dr ; protect against divide by zero rcl vr x<0 goto neg_vr ; jump if vr is negative ; D_impact_angle = impact_angle - new_angle rcl impact_angle_rad rcl vt rcl vr rcl dr + / atan - >deg ; convert to degrees 90 + ; swing additional 90 deg ccw goto el_done ; continue lbl neg_vr ; when vr is negative ; resulting angle = impact_angle + new_angle rcl impact_angle_rad rcl vt rcl vr abs rcl dr + / atan + >deg ; convert to degrees 90 - ; swing additional 90 deg cw lbl el_done ; continue 360 rand * ; azimuth angle rcl v ; velocity >r3d ; compute rect assuming vertical is on original line -90 >elr ; el rotate back from 90 vertical rcl el >elr ;el rotate to initial el rcl az >azr ;az rotate to initial az sto ion_vx_mm ;store back into user variables rlup sto ion_vy_mm rlup sto ion_vz_mm rldn rldn >p3d rcl ion_mass x<>y >ke rcl ion_number arcl totals + rcl ion_number asto totals rcl ion_number arcl counts 1 + rcl ion_number asto counts rcl _Mark_Collisions abs x!=0 mark ;mark ion's collision if requested lbl skip1 seg terminate rcl ion_number arcl totals rcl ion_number arcl counts sto count 0.000001 + / rcl count rcl ion_number message ; ion = # collisions = #, avg ke = # eV