Viewing contents of file '../idllib/contrib/tappin/graffer/gr_fit_funct.pro'
function Gr_fit_funct, pdefs, ftype, npt, slice, fset, pr, resid

;+
; GR_FIT_FUNCT
;	Make a least squares fit to a dataset
;
; Usage:
;	gr_fit_funct, pdefs
;
; Argument:
;	pdefs	struct	in/out	The GRAFFER plot structure.
;
; History:
;	Original: 28/8/96; SJT
;	Add support for median fits and simplify representation of fit
;	types: 23/9/96; SJT
;	Change to make this the inner routine & allow higher degree
;	fits: 22/11/96; SJT
;	Shorten name: 25/11/96; SJT
;	Add CHECK error handling for fitting calls thus eliminating
;	the need for local SVDFIT version and making safer: 29/1/97; SJT
;	Add options to do "piecewise linear" fits: 4/2/97; SJT
;-

handle_value, pdefs.data, data, /no_copy
handle_value, data(fset).xydata, xy, /no_copy

if (ftype(2) eq 0) then case (data(fset).type) of
    0:                          ; No errs
    1: wy = transpose(xy(2, *)) ; Y
    2: wy = total(xy(2:3, *), 1)/2. ; +-Y
    3: wx = transpose(xy(2, *)) ; X
    4: wx = total(xy(2:3, *), 1)/2. ; +-X
    5: begin
        wx = transpose(xy(2, *))
        wy = transpose(xy(3, *)) ; XY
    end
    6: begin
        wx = transpose(xy(2, *))
        wy = total(xy(3:4, *), 1)/2. ; X+-Y
    end
    7: begin
        wx = total(xy(2:3, *), 1)/2.
        wy = transpose(xy(4, *)) ; +-XY
    end
    8: begin
        wx = total(xy(2:3, *), 1)/2.
        wy = total(xy(4:5, *), 1)/2. ; +-X+-Y
    end
endcase

x = transpose(xy(0, *))
y = transpose(xy(1, *))

handle_value, data(fset).xydata, xy, /no_copy, /set

if (ftype(1)) ne 0 then begin
    temp = temporary(x)
    x = temporary(y)
    y = temporary(temp)
    if (N_elements(wx) ne 0) then w = wx
    fftype = -2
    var = 'y'
endif else begin
    fftype = -1
    if (N_elements(wy) ne 0) then w = wy
    var = 'x'
endelse

if (ftype(0) ne 4) then begin
    if (ftype(0) and 2) ne 0 then begin
        x = alog(x)
        var = 'alog('+var+')'
    endif
    
    if (ftype(0) and 1) ne 0 then begin
        if (n_elements(w) ne 0) then w = w/y
        y = alog(y)
    endif
endif

if (slice ne '') then begin
    junk1 = execute('x = x'+slice)
    junk2 = execute('y = y'+slice)
    if (n_elements(w) ne 0) then junk3 = execute('w = w'+slice) $
    else junk3 = 1
    if (junk1+junk2+junk3 ne 3) then begin
        graff_msg, resid, 'Bad slice specification -- fit ' + $
          'not performed'
        pr = 0.
        handle_value, pdefs.data, data, /no_copy, /set
        return, ''
    endif
endif

if (n_elements(w) ne 0) then w = 1./w

sing = 0
c2 = 0.

                                ; It is quite possible to make SVDFIT
                                ; do a belly up so we need an error
                                ; handler to catch the error otherwise
                                ; GRAFFER will crash ignominiously
                                ; with the probable loss of pdefs.

catch, ecode
if (ecode ne 0) then begin
    widget_control, resid, set_value = 'Fitting error - degree ' + $
      'too high?'
    handle_value, pdefs.data, data, /no_copy, /set
    pr = 0.
    return, ''
endif

if (ftype(0) eq 4) then fit = gr_pwfit(x, y, w, ftype(3), c2) $
else if (ftype(2)) then fit = ladfit(x, y) $
else fit = svdfit(x, y, ftype(3)+1, weight = w, sing = sing, chisq = c2)

catch, /cancel                  ; Should only need the Catcher for SVDFIT

func = gr_str_fun(fit, var, pieces = ftype(0) eq 4)
if (ftype(0) and 1) then func = 'exp('+func+')'


xydata = {graff_funct,  $
          Range: [0., 0.],  $
          Funct: func}

data(pdefs.cset).ndata = npt
data(pdefs.cset).type = fftype

handle_value, data(pdefs.cset).xydata, xydata, /no_copy, /set
handle_value, pdefs.data, data, /no_copy, /set

free = n_elements(x) - n_elements(fit)

catch, ecode                    ; CHISQR_PDF can also fall down
if (ecode ne 0) then pr = 0. $
else pr = 1. - chisqr_pdf(c2, free)
catch, /cancel                  ; Should only need the Catcher for CHISQR_PDF

return, func

end