pro prob8_3,ss,mbar,vbar,a,b

;simulates sheet 8/problem 3
; if a=0, the analytic (first order, linear) relation between e and p
; is overplotted (problem 3c)

;example: prob8_3,1000,4.,5.,0.,0.3
;example: prob8_3,1000,4.,5.,0.,0.05
;example: prob8_3,1000,4.,-5.,0.,0.3
;example: prob8_3,1000,4.,-5.,0.1,0.1
;example: prob8_3,1000,4.,5.,0.3,0.3
;example: prob8_3,1000,100.,-5.,0.,0.8
  
;=>as long as a and b small, everything OK.
;  Otherwise, non-linear effects take over.  
  
sig=mbar*a
m=randomn(seed,ss)*sig+mbar

sig=vbar*b
v=randomn(seed,ss)*sig+vbar

p=m*v
e=m*v^2/2.

plot,p,e,psym=1,xtitle='p',ytitle='E'

pbar=mbar*vbar
ebar=mbar*vbar^2/2.

pmin=min(p)
pmax=max(p)
emin=min(e)
emax=max(e)

plots,[pbar,pbar],[0.,emax],lines=1
plots,[pmin,pmax],[ebar,ebar],lines=1

print,'correlation coefficient'
print,'mc:  ',correlate(p,e)
print,'calc:',(a^2+2.*b^2)/sqrt((a^2+b^2)*(a^2+4.*b^2))*pbar/abs(pbar)

if(a eq 0.) then begin
emin=-ebar+2.*ebar/pbar*pmin
emax=-ebar+2.*ebar/pbar*pmax
plots,[pmin,pmax],[emin,emax],lines=2
endif

return
end