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