' fsctlib.bas - jackord@kw.igs.net - revised 27 Oct 02 - Liberty Basic v3.01
' Rutherford scattering from a fixed force center for impact parameters
'   less than five times distance of closest approach at head-on impact
' SEQ - 26 particles at regular impact parameter increments
' RND - 26 particles randomly distributed across the area of the beam
' CNT - angular distribution for 500 randomly distributed particles

' Initialize Window
  nomainwin
  button#1, "SEQ", [seq], UL, 145, 1, 30, 20
  button#1, "RND", [ran], UL, 185, 1, 30, 20
  button#1, "CNT", [cnt], UL, 225, 1, 30, 20
  WindowWidth=410     ' pixel scale 0-400
  WindowHeight=329    ' pixel scale 0-300
  UpperLeftX=100: UpperLeftY=100
  open "Rutherford Scattering" for graphics_nsb as #1   ' no slide bars
  #1 "trapclose [quit]"
  dim cc(20): dim yy(316)
  gosub [grid]

[waitHere]
  wait

[grid]
  #1 "cls ; down; backcolor green ; color green"
  #1 "place 200 150 ; circlefilled 10"
  #1 "color blue ; backcolor white"
return

[hit]                                                  ' Feynman algorithm
  x=0-200: dt=1: vx=.4416: vy=0: x1=0: y1=0            '   with trajectory
  r3=((x-200)*(x-200)+(y-150)*(y-150))^1.5
  vx=vx+(x-200)/r3*dt/2: vy=vy+(y-150)/r3*dt/2
  phase=1
  while phase<3 and r3<66000000
    x=x+vx*dt: y=y+vy*dt
    r3=((x-200)*(x-200)+(y-150)*(y-150))^1.5
    vx=vx+(x-200)/r3*dt: vy=vy+(y-150)/r3*dt
    if phase=2 then
      #1 "goto "; int(x+.5); " "; int(y+.5)
      if x>400 or y>300 then phase=3
    end if
    if phase=1 and x>=0 then #1 "place "; int(x+.5); " "; int(y+.5): phase=2
  wend
return

[ctr]                                                  ' Feynman Algorithm
  x=0-200: dt=1: vx=.4416: vy=0                        '   no trajectory
  r3=((x-200)*(x-200)+(y-150)*(y-150))^1.5
  vx=vx+(x-200)/r3*dt/2: vy=vy+(y-150)/r3*dt/2
  while r3<66000000
    x=x+vx*dt: y=y+vy*dt
    r3=((x-200)*(x-200)+(y-150)*(y-150))^1.5
    vx=vx+(x-200)/r3*dt: vy=vy+(y-150)/r3*dt
  wend
  tt=atn(vy/vx)/atn(1)*45
  if vx<0 and vy>=0 then tt=tt+180
  if vx<0 and vy<0 then tt=tt-180
  tt=abs(tt)
return


[gety]                                                 ' Get random y
  test=0
  while test=0
    x=rnd(1): y=rnd(1)
    r=((x-.5)*(x-.5)+(y-.5)*(y-.5))^.5
    if r<=.5 then test=1
  wend
  r=50*r: if y<.5 then r=0-r
  y=150+r
return

[seq]                                                  ' Stepped impact
  gosub [grid]                                         '   parameter
  for i=125 to 175 step 2
    y=i: gosub [hit]
  next i
goto [waitHere]

[ran]                                                  ' Random impact
  gosub [grid]                                         '   parameter
  for i=0 to 25
    gosub [gety]: gosub [hit]
  next i
goto [waitHere]

[cnt]                                                  ' 500 particle counter
  redim cc(20)
  xorg=30: yorg=280
  redim yy(316)
  pi=4*atn(1): t=2*atn(.2)/pi*180: u=tan(pi/180*t/2)
  yy(0)=yorg-int(10*pi/9*(1+u*u)/u/u/u+.5)
  for i=46 to 360
    u=tan(pi/180*i/4)
    yy(i-45)=yorg-int(10*pi/9*(1+u*u)/u/u/u+.5)
  next i
  #1 "cls ; color black ; down"
  for i=30 to 390 step 20
    #1 "line "; i; " 20 "; i; " "; yorg
  next i
  for i=20 to 280 step 20
    #1 "line "; xorg; " "; i; " 390 "; i
  next i
  #1 "place "; xorg-4; " "; yorg+16: #1 "\0"
  #1 "place "; 150-8; " "; yorg+16: #1 "\60"
  #1 "place "; 270-12; " "; yorg+16: #1 "\120"
  #1 "place "; 390-12; " "; yorg+16: #1 "\180"
  #1 "place 180 "; yorg+17: #1 "\Theta  deg"
  #1 "place 0 136": #1 "\Cnt"
  #1 "place "; xorg-12; " "; yorg+6: #1 "\0"
  #1 "place "; xorg-26; " "; yorg+6-100: #1 "\100"
  #1 "place "; xorg-26; " "; yorg+6-200: #1 "\200"
  #1 "color blue"                                ' Analytic function
  for i=0 to 315
    #1 "line "; xorg+45+i; " "; yy(i); " "; xorg+45+i; " "; yorg
  next i
  #1 "place "; xorg+20; " "; yorg+16: #1 "\N(<60) = 440"
  #1 "place "; xorg+265; " "; yorg+16: #1 "\N(>60) = 60"
  timer 2000, [cont]
  wait
[cont]
  timer 0
  #1 "color red ; backcolor red"
  for i=1 to 500
    gosub [gety]: gosub [ctr]: x1=int(tt/10): cc(x1)=cc(x1)+1: x2=xorg+20*x1
    y2=yorg-cc(x1): #1 "line "; x2; " "; y2; " "; x2+20; " "; y2
    scan
  next i
  x1=0: x2=0
  for i=2 to 5: x1=x1+cc(i): next i
  for i=6 to 17: x2=x2+cc(i): next i
  #1 "backcolor white"
  #1 "place "; xorg+20; " "; yorg+16: #1 "\"; space$(12)
  #1 "place "; xorg+20; " "; yorg+16: #1 "\N(<60) = "; x1
  #1 "place "; xorg+265; " "; yorg+16: #1 "\"; space$(12)
  #1 "place "; xorg+265; " "; yorg+16: #1 "\N(>60) = "; x2
goto [waitHere]

[quit]
  close #1
end