;------------------------------------------
;like KS_EXAMPLE but using XSQ
;---------------------------

;read radial velocity of 1000 simulation particles
   restore,'tb_iimodel2_a100tau010.pos'
   data=vx/omega
   sigma=sqrt(mean(data^2))
;zero mean


;bin the data to get observed frequancies and
;calculate the expected values for the same bins
;do not worry about too small or zero frequancies

   xtable=-40.+(findgen(80)+0.5) ;roughly +/- 5 sigma, certainly enough
   dx=1.
   obfreq=xtable*0.
   exfreq=xtable*0.
   
   for i=0,n_elements(xtable) -1 do begin       
       x1=xtable(i)-dx/2.
       x2=xtable(i)+dx/2.       
       ind=where(data gt x1 and data le x2,count)
       if(count ge 1) then obfreq(i)=count       
;model: 
;difference of the cumulative frequency distribution
;at the end of the interval
       exfreq(i)=(gauss_pdf(x2/sigma)-gauss_pdf(x1/sigma))*1000.
   endfor

;xsq_test makes rebinning if bins contain too few data
;rebinned data returned in obcell, excell
   Res = XSQ_TEST(Obfreq, Exfreq , EXCELL=excell, OBCELL=obcell)

   ff='(f6.2)'
   print,'N=1000,   chi2='+string(res(0),ff)+$
     '   Prob of chi2 or greater'+string(res(1),ff)

   print,'obcell'
   print,obcell
   print,'excell'
   print,excell


;----------------------------------------- 
;take just 30 points from the distribution

   data=data(0:29)
   obfreq=xtable*0.
   exfreq=xtable*0.
   
   for i=0,n_elements(xtable) -1 do begin       
       x1=xtable(i)-dx/2.
       x2=xtable(i)+dx/2.       
       ind=where(data gt x1 and data le x2,count)
       if(count ge 1) then obfreq(i)=count       
       exfreq(i)=(gauss_pdf(x2/sigma)-gauss_pdf(x1/sigma))*30.
   endfor
   Res = XSQ_TEST(Obfreq, Exfreq , EXCELL=excell, OBCELL=obcell)

   print,'N=30,     chi2='+string(res(0),ff)+$
     '   Prob of chi2 or greater'+string(res(1),ff)

   print,'obcell'
   print,obcell
   print,'excell'
   print,excell
end