PROGRAM scatter
! compute differential and total cross section
DIM dN(0 to 360)
LIBRARY "csgraphics"
CALL initial(x0,vx0,R,delta_b,dt,nbin,dtheta)
LET N = 0                         ! number of incident particles
LET b = 0
! make sure that b = R actually done, but b might exceed R slightly
DO while b <= R + delta_b/2
   CALL trajectory(vx,vy,vx0,x0,b,dt)
   CALL detector(dN(),vx,vy,b,dtheta,N)
   LET b = b + delta_b
LOOP
SET COLOR "black/white"
SET CURSOR 22,1
PRINT "Hit any key for the differential cross section."
DO
LOOP until key input
CALL output(dN(),N,nbin,dtheta,R)
END

SUB initial(x0,vx0,R,delta_b,dt,nbin,dtheta)
    INPUT prompt "kinetic energy of incident particles = ": E
    INPUT prompt "incremental change in impact parameter = ": delta_b
    LET R = 1.0                   ! beam radius
    LET vx0 = sqr(2*E)            ! mass = 1
    LET dt = 0.01                 ! time step
    LET nbin = 9                  ! number of detectors
    LET dtheta = 180/nbin         ! degrees
    ! initial position of beam relative to target
    LET x0 = -3
    CALL compute_aspect_ratio(abs(x0),xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
    SET COLOR "blue"
    LET Rs = 1.0                  ! scattering region radius
    BOX CIRCLE -Rs,Rs,-Rs,Rs
    PLOT 0,0
    SET COLOR "red"
END SUB

SUB trajectory(vx,vy,vx0,x0,b,dt)
    LET y = b
    LET x = x0
    LET vx = vx0
    LET vy = 0
    DO
       CALL EulerRichardson(x,y,vx,vy,dt)   ! one time step
       PLOT POINTS: x,y
    LOOP until x*x + y*y > x0*x0 + b*b
END SUB

SUB detector(dN(),vx,vy,b,dtheta,N)
    DECLARE DEF arctan
    LET theta = int(arctan(vx,vy)/dtheta)   ! scattering angle
    LET dN(theta) = dN(theta) + b
    LET N = N + b
END SUB

SUB EulerRichardson(x,y,vx,vy,dt)
    ! could use arrays to reduce # of statements
    CALL force(x,y,fx,fy)
    ! mass is unity in program units
    LET vxmid = vx + 0.5*fx*dt
    LET vymid = vy + 0.5*fy*dt
    LET xmid = x + 0.5*vx*dt
    LET ymid = y + 0.5*vy*dt
    CALL force(xmid,ymid,fxmid,fymid)
    LET vx = vx + fxmid*dt
    LET vy = vy + fymid*dt
    LET x = x + vxmid*dt
    LET y = y + vymid*dt
END SUB

SUB force(x,y,fx,fy)
    LET r2 = x*x + y*y
    LET r = sqr(r2)
    ! force for simple model of Hydrogen atom
    ! range of force set equal to unity
    IF r2 <= 1 then
       LET f = 1/r2 - r
    ELSE
       LET f = 0
    END IF
    LET fx = f*x/r
    LET fy = f*y/r
END SUB

SUB output(dN(),N,nbin,dtheta,R)
    CLEAR
    LET dtheta_rad = dtheta*pi/180     ! radians
    LET target_density = 1/(pi*R*R)
    PRINT "theta","dN/N","sigma(theta)"
    PRINT
    FOR i = 0 to nbin
        LET theta = (i + 0.5)*dtheta_rad
        IF dN(i) > 0 then
           LET d_Omega = 2*pi*sin(theta)*dtheta_rad
           LET diffcross = dN(i)/(d_Omega*N*target_density)
           LET totalcross = totalcross + diffcross*d_Omega
           PRINT (i+1)*dtheta,dN(i)/N,diffcross
        END IF
    NEXT i
    PRINT "total cross section ="; totalcross
END SUB

FUNCTION arctan(x,y)
    ! convert -pi/2 to pi/2 radians to 0 to 180 degrees
    IF x > 0  then
       LET theta = atn(abs(y/x))
    ELSE IF  y > 0 then
       LET theta = atn(y/x) + pi
    ELSE
       LET theta = pi - atn(y/x)    ! corrected from text
    END IF
    LET arctan = theta*180/pi
END DEF