; $Id: gen_fit.pro,v 1.4 2003/04/22 19:10:16 jmspark Exp $
;*******************************************************************************
;23456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************
;
; Gaussian = p[0] exp( - 4 ln(2) ( x - p[1] )^2 / p[2]^2 )
;
; p[0] = Amplitude (Peak Value)
; p[1] = Center of the Peak
; p[2] = Full Width Half Maximum
;
; Area = p[0] * p[2] * sqrt( !pi / alog(2.) ) / 2.
;
PRO gau ,x ,par ,y ,pdpar ,i

	z=(x-par(1))/par(2)

	con=4.*ALOG(2.)
	ez=FLTARR(n_elements(x))
	index = WHERE(abs(z) LE 7.,count)

	IF (count GT 0) THEN ez(index)=exp(-con*z(index)^2)
	y=par(0)*ez

	IF n_params(0) LE 3 THEN RETURN

	pdpar(*,3*i)=ez
	pdpar(*,3*i+1)=8.*alog(2.)*par(0)*ez*z/par(2)
	pdpar(*,3*i+2)=pdpar(*,3*i+1)*z

END
;*******************************************************************************
;23456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
;       10        20        30        40        50        60        70        80
;*******************************************************************************

PRO lor,x,par,y,pdpar,i
z=(x-par(1))/par(2)
lz=1./(1+z^2.)
y=par(0)*lz
IF n_params(0) LE 3 THEN RETURN
pdpar(*,3*i)=lz
pdpar(*,3*i+1)=2.*lz^2.*z/par(2)
pdpar(*,3*i+2)=pdpar(*,3*i+1)*z
END

PRO lorgauconv,x,par,y,pdpar,i
COMMON gaussval, sg	;need to get the gaussian FWHM
;The convolution of a lorentzian with a gaussian using complex error functions.
;This is more accurate than Gauss-Hermite quadrature in any case and should not fail
;for any values of the width parameters.
;See documentation on cerrorf.pro and references therein about how this is used.
;Paul Schleger, Riso national Laboratory, Sept. 1993

;some constants needed
PI = 3.141592654
n1=SQRT(ALOG(2.))
n2=SQRT(PI*ALOG(2.))
rootpi=SQRT(PI)

;Map the fit parameters:
x0 = par(1)
sl = 2.*par(2)
ymax = par(0)


;See Nuclear Instruments and Methods Vol. 137, 179-181 (1976) for
;what's going on here....
eh=fltarr(n_elements(x))
eh(*)=n1*sl/sg
v=2.*n1*(x-x0)/sg
w=cerrorf(COMPLEX(v,eh))
y=ymax*n2*sl*FLOAT(w)/sg
IF n_params(0) LE 3 THEN RETURN

;Calculate partial derivatives:
dwdv=2.*(-v*FLOAT(w)+eh*IMAGINARY(w))
dwda=2.*(-1./rootpi+v*IMAGINARY(w)+eh*FLOAT(w))
pdpar(*,3*i) = sl*n2*FLOAT(w)/sg
pdpar(*,3*i+1) = -2.*ymax*sl*rootpi*ALOG(2.)*dwdv/(sg*sg)
pdpar(*,3*i+2) =  ymax*n2*(FLOAT(w)+eh*dwda)/sg
END

pro general,x,b,f,pder
common fitp,ngau,nlor,a,vary,pnames,gw,ufunct,nufun
common gaussval, local_sg

;First generate the gauss-hermite convolution parameters for the
;user defined function (if it's defined)
IF (gw NE 0.0) AND (ufunct NE '') THEN BEGIN
  xgh=[0.3811869902,1.1571937124,1.9816567566,2.9306374202]
  wgh=[6.6114701255e-1,2.0780232581e-1,1.7077983007e-2,1.9960407221e-4]
  xgh=xgh*gw/(2.*alog(2.)^0.5)
  wgh=wgh/!pi^0.5
END

;Map the variable parameters to the actual parameters passed to the functions
c=where(vary,count)
IF count NE 0 THEN for i=0,total(vary)-1 do a(c(i))=b(i)

;I don't know why there is a check here, but, OK perhaps I'll figure it out soon
IF ((ngau GT 0) OR (nlor GT 0) OR (nufun GT 0)) THEN BEGIN
	tpd=fltarr(n_elements(x),3*ngau+3*nlor+nufun)
	ghpd=fltarr(n_elements(x),3*ngau+3*nlor+nufun)
END

;Start with the background
f=a(0)+a(1)*x+a(2)*x^2

;Add the gausians. Ignore resolution function value, since we can seperate it out later
FOR i=0,ngau-1 DO BEGIN
  IF (gw EQ 0.0) THEN BEGIN
  	IF (n_params(0) LE 3) THEN gau,x,a(3+3*i:5+3*i),y ELSE $
     	gau,x,a(3+3*i:5+3*i),y,tpd,i
  END ELSE IF (n_params(0) LE 3) THEN $
  		gau,x,[a(3+3*i),a(4+3*i),SQRT(gw^2+(a(5+3*i))^2)],y ELSE $
   		gau,x,[a(3+3*i),a(4+3*i),SQRT(gw^2+(a(5+3*i))^2)],y,tpd,i
  f=f+y
END

;Add the lorentzians
FOR i=ngau,ngau+nlor-1 DO BEGIN
  IF (gw eq 0.0) THEN BEGIN
    IF (n_params(0) le 3) THEN lor,x,a(3+3*i:5+3*i),y ELSE $
     lor,x,a(3+3*i:5+3*i),y,tpd,i
    f=f+y
  END ELSE BEGIN
   local_sg = gw
   IF (n_params(0) le 3) THEN lorgauconv,x,a(3+3*i:5+3*i),y ELSE $
     lorgauconv,x,a(3+3*i:5+3*i),y,tpd,i
    f=f+y
  END
END

;Finally, add the userfunction
IF (ufunct NE '') THEN BEGIN
  IF (gw EQ 0.0) THEN BEGIN 	;no resolution function to convolve
    IF (n_params() LE 3) THEN $
      fstr=ufunct(0)+',x,a(6+3*(ngau+nlor-1):5+3*(ngau+nlor-1)+nufun),y' ELSE $
      fstr=ufunct(0)+',x,a(6+3*(ngau+nlor-1):5+3*(ngau+nlor-1)+nufun),y,tpd,(ngau+nlor)'
    run=execute(fstr)
    f=f+y
  END ELSE BEGIN		;convolve the resolution function (Gauss-Hermite??)
    FOR gh=0,3 DO BEGIN		;We seem to be doing a 4th order G-H
      IF (n_params() LE 3) THEN $	;First we do the x-xgh
        fstr=ufunct(0)+',x-xgh(gh),a(6+3*(ngau+nlor-1):5+3*(ngau+nlor-1)+nufun),y' ELSE $
        fstr=ufunct(0)+',x-xgh(gh),a(6+3*(ngau+nlor-1):5+3*(ngau+nlor-1)+nufun),y,ghpd,(ngau+nlor)'
      run=execute(fstr)
      f=f+y*wgh(gh)
      IF (n_params() GT 3) THEN $
        tpd(*,3*(ngau+nlor):3*(ngau+nlor)+nufun-1) = $
        	tpd(*,3*(ngau+nlor):3*(ngau+nlor)+nufun-1) + $
        	wgh(gh)*ghpd(*,3*(ngau+nlor):3*(ngau+nlor)+nufun-1)
      IF (n_params() LE 3) THEN $	;Then we do the x+xgh
        fstr=ufunct(0)+',x+xgh(gh),a(6+3*(ngau+nlor-1):5+3*(ngau+nlor-1)+nufun),y' ELSE $
        fstr=ufunct(0)+',x+xgh(gh),a(6+3*(ngau+nlor-1):5+3*(ngau+nlor-1)+nufun),y,ghpd,(ngau+nlor)'
      run=execute(fstr)
      f=f+y*wgh(gh)
      IF (n_params() gt 3) THEN $
        tpd(*,3*(ngau+nlor):3*(ngau+nlor)+nufun-1) = $
        	tpd(*,3*(ngau+nlor):3*(ngau+nlor)+nufun-1) + $
        	wgh(gh)*ghpd(*,3*(ngau+nlor):3*(ngau+nlor)+nufun-1)
    END
  END
END
IF (n_params(0) LE 3) THEN RETURN ;Don't want the partial derivatives

;Map the partial derivatives to their appropriate places
pder = fltarr(n_elements(x),total(vary))
pind=0
if (vary(0) eq 1) then begin
  pder(*,pind) = 1.
  pind=pind+1
  endif
if (vary(1) eq 1) then begin
  pder(*,pind) = x
  pind=pind+1
  endif
if (vary(2) eq 1) then begin
  pder(*,pind) = x^2
  pind=pind+1
  endif
for i=0,ngau+nlor-1 do begin
  if (vary(3+3*i) eq 1) then begin
    pder(*,pind)=tpd(*,3*i)
    pind=pind+1
    endif
  if (vary(4+3*i) eq 1) then begin
    pder(*,pind)=tpd(*,3*i+1)
    pind=pind+1
    endif
  if (vary(5+3*i) eq 1) then begin
    pder(*,pind)=tpd(*,3*i+2)
    pind=pind+1
    endif
  endfor
for i=3*(ngau+nlor),3*(ngau+nlor)+nufun-1 do $
  if (vary(3+i) eq 1) then begin
    pder(*,pind)=tpd(*,i)
    pind=pind+1
    endif
return
end


function gen_fit,x,y,w,sigma=sigma,minuit=minuit
;+
; NAME:
;       GEN_FIT
; PURPOSE:
;       fit y=f(x) where:
;       f(x) = a0+a1*x+a2*x^2+arbitrary number of gaussians+lorentzians
; CATEGORY:
;       ?? - fitting
; CALLING SEQUENCE:
;       yfit = gen_fit(x,y,w)
; INPUTS:
;       x = independent variable, must be a vector.
;       y = dependent variable, must have the same number of points
;               as x.
;       quiet = set to inhibit printing curfit iterations.
; OUTPUTS:
;       chisquared.
; OPTIONAL OUTPUT PARAMETERS:
;       a = coefficients. a three element vector as described above.
;
; COMMON BLOCKS:
;       None.
; SIDE EFFECTS:
;       None.
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;       Thom Mason, Riso National Laboratory, 1993
;-
;
common fitp,ngau,nlor,a,vary,pnames,gw,ufunct,nufun
COMMON com_minchoice, minchoice, minpath

IF NOT keyword_set(sigma) THEN sigma=0 ELSE sigma = fltarr(n_elements(a))
PRINT,MINCHOICE
IF (keyword_set(minuit) OR minchoice EQ 'MFit' OR minchoice EQ 'Simp' $
   OR minchoice EQ 'MFit+') THEN $
;IF 1 eq 1 then $
	chisquared=minuitfit(x,y,w,sigma=sigma,minimizer=minchoice)  $
ELSE BEGIN
	n=n_elements(y)
	w=replicate(0.,n)+w
	index=where(vary,count)
	IF count NE 0 THEN BEGIN
		b=a(index)
		nerd=gcurfit(x,y,w,b,sigmab,/quiet) ;call curfit
		chisquared=TOTAL(w*(y-nerd)^2)/(n_elements(x)-TOTAL(vary))
		IF keyword_set(sigma) THEN sigma(index)=sigmab
	END ELSE BEGIN
		nerd=0.0
		chisquared=0
	END
END
return, chisquared
END





