;------------------------------------------
;Model: gaussian with a given sigma
function ks_model,x
   common ks_model_com,sigma_com

;the only reason to use auxillary variable SIGMA_COM
;is to avoid accidental conflict between DIFFERENT routines,  
;caused because SIGMA is so often used variable name...


   sigma=sigma_com
   res=gauss_pdf(x/sigma)
   return,res
end

;------------------------------------------
;main program KS_EXAMPLE
;test the hypothesis that radial velocities in 
;a planetary ring simulation follow a Gaussian
;------------------------------------------
   common ks_model_com,sigma_com

   ps=0.1
   program='ks_example'
   psdirect,program,ps,/color

;---------------------------
;1000 particles
   restore,'tb_iimodel2_a100tau010.pos'
;   restore,'tb_iimodel1_a100tau400.pos'
   data=vx/omega



   sigma=sqrt(mean(data^2))
   SIGMA_COM=SIGMA


   !p.multi=[0,2,2]
   nwin
   ff='(f6.4)'

   !x.range=[-40,40]
   !p.title='1000 radial velocity points from a simulation'
   histo_f,data,-30,30,1,/plot,/gauss,psym=10
   xyouts,10,0.07,'sigma='+string(sigma,ff)
   !p.title='Output of ksone.pro'

   ksone,data,'ks_model',d,prob,/plot,xtitle='vx',ytitle='CFD'
   xyouts,-35,0.8,'D='+string(d,ff)
   xyouts,-35,0.7,'Prob='+string(prob,ff)
   
;----------------------------------------- 
;take just 30 points from the distribution
   data=data(0:29)
   
   nwin
   !p.title='sample 30 out of 1000 points'
   histo_f,data,-30,30,1,/plot,/gauss,psym=10
   !p.title='Output of ksone.pro'

   ksone,data,'ks_model',d,prob,/plot,xtitle='vx',ytitle='CFD'
   xyouts,-35,0.8,'D='+string(d,ff)
   xyouts,-35,0.7,'Prob='+string(prob,ff)
   
   !p.multi=0
   defplot
   psdirect,program,ps,/stop
   
;-------------------------------------------------
;THE DATA FITS THE MODEL TOO WELL
;what happens if somebody forges the data to be too good?

;take the same original datapoints with their sigma
   restore,'tb_iimodel2_a100tau010.pos'
   data=vx/omega
   sigma=sqrt(mean(data^2))
   SIGMA_COM=SIGMA
;replace true data by values from the distribution and add some noise
;samples cumulative density distribution evenly:
   for i=0,999 do begin 
       data(i)=sigma*gauss_cvf(1.-0.001*(i+0.5)) 
   endfor
;noise
   data=data+randomn(seed,1000)*.5
 


   psdirect,program+'_too_good',ps
   !p.multi=[0,2,1]
   !x.range=[-40,40]
   nwin
   !p.title='TOO GOOD TO BE TRUE?'
   histo_f,data,-30,30,1,/plot,/gauss,psym=10
   !p.title='Output of ksone.pro'

   ksone,data,'ks_model',d,prob,/plot,xtitle='vx',ytitle='CFD'
   xyouts,-35,0.8,'D='+string(d,ff)
   xyouts,-35,0.7,'Prob='+string(prob,ff)
   
   !p.multi=0
   psdirect,program+'_too_good',ps,/stop
end