pro fit_cum2,xhist,yhist,xh,yh,yfit,ind_fit,x0,x2,a,da,chi_cum,noise,dev,nbin_cum,opt
;+
; Project     : all data
;
; Name        : POWERLAWFIT
;
; Category    : Event statistics
;
; Explanation : Fitting of powerlaw function with threshold effects
;
; Syntax      : IDL>powerlawfit,x,
;
; Inputs      : x=fltarr(n)=array of values of n events  
;
; Outputs     ; a	=powerlaw slope 
;
; History     : 20-Apr-2015, Version 1 written by Markus J. Aschwanden
;
; Contact     : aschwanden@lmsal.com
;-

common powerlaw_var,xh_,yh_,ind_fit_,yfit_,x0_,noise_,dev_
x0_	=x0

;_______________________CUMULATIVE HISTOGRAM___________________________
x1      =min(xhist)
x2      =max(xhist)
x1log   =alog10(x1)
dhlog   =1./float(nbin_cum)
nh      =long((alog10(x2)-alog10(x1))/dhlog+1)
xh      =10.^(x1log+dhlog*findgen(nh))
yh      =fltarr(nh)
for i=0,nh-1 do begin
 ind    =where((xhist ge xh(i)),nhist)		;cumulative number
 nbin	=total(yhist(ind))
 yh(i)  =nbin
endfor
ynext	=shift(yh,+1)
ind	=where((yh ne ynext) and (yh ne 0),nh)
xh	=xh(ind)
yh	=yh(ind)

;_______________________POWERLAW SLOPE GLOBAL SEARCH__________________
xlog	=alog10(xh)
nlog	=alog10(yh)
x1log	=xlog(0)
x2log	=xlog(nh-1)
n0	=max(yh)
nvar	=100
devmin	=1.e10
for i=1,nvar-1 do begin
 a	=1.0+2.0*float(i)/float(nvar-1)
 yfit_  =n0*((x2+x0)^(1.-a)-(xh+x0)^(1.-a))/((x2+x0)^(1.-a)-(x1+x0)^(1.-a))
 ind	 =where((yh ge 1) and (yh lt 0.9*n0),nind)
 dev    =sqrt(total((alog10(yfit_(ind))-nlog(ind))^2)/float(nind))
 if (dev lt devmin) then begin
  devmin=dev
  a_best =a
  n0_best=n0
  yfit_best=yfit_
 endif
endfor

;________________________POWELL OPTIMIZATION___________________________
xh_     =xh
yh_     =yh
yfit_	=yfit_best
a	=a_best
n0	=n0_best
nmax	=max(yh)
ind_fit =where((yh ge 1.) and (yh lt 0.9*nmax),nf)
ind_fit_=ind_fit
nev     =max(yh)
coeff   =[alog10(n0),a_best]
ncoeff  =n_elements(coeff)
xi      =fltarr(ncoeff,ncoeff)
for i=0,ncoeff-1 do xi(i,i)=0.1
ftol    =1.e-4
if (opt eq 1) then powell,coeff,xi,ftol,chi_cum,'fit_cum_powell'
chi_cum =fit_cum_powell(coeff)         ;-->updates yfit_,noise_,dev_,chi_cum
yfit	=yfit_
ind_fit	=ind_fit_
n0	=10.^coeff(0)
a       =coeff(1)
da	=a/sqrt(nev)
noise   =noise_
dev	=dev_
end