;+ ;FUNCTION: fu_fit2d(dat,ENERGY=en,ERANGE=er,EBINS=ebins,...) ; ; Select energy/angle range and perform functional fit to data ; (default Maxwellian) ;INPUTS: ; dat - data structure containing 2d data ; - e.g. "get_fa_ees, get_fa_ies, etc." ;KEYWORDS: ; NFITF integer optional number of functions for the ; fit, default = 1 ; FITF string optional function for the fit, default ; = maxwellian ; AUTO 0/1 if set, defaults are used and peak ; energy determined from the data ; ENERGY: fltarr(2), optional, min,max energy range for fit ; ERANGE: fltarr(2), optional, min,max energy bin numbers ; for fit ; EBINS: bytarr(na), optional, energy bins array for fit ; 0,1=exclude,include, ; na = dat.nenergy ; ANGLE: fltarr(2), optional, min,max pitch angle range ; for fit ; ARANGE: fltarr(2), optional, min,max angle bin numbers ; for fit ; BINS: bytarr(nb), optional, angle bins array for fit ; 0,1=exclude,include, ; nb = dat.ntheta ; BINS: bytarr(na,nb), optional, energy/angle bins array for ; fit 0,1=exclude,include ; LIMITS - A structure containing limits and display options. ; see: "options", "xlim" and "ylim", to change limits ; UNITS - convert to given data units before plotting ; F_UNITS - convert dat.data to given data units before ; passing to function "FITF" ; COLOR - array of colors to be used for each bin ; LABEL - Puts bin labels on the plot if set. ; SUMPLT - redraws the summed data and the fit after ; completion. ; MSEC - Subtitle will include milliseconds ; ;NOTES: ; This program still has bugs and may not weight the points ; properly !!!!! ; ; See "spec2d" for another means of plotting data. ; See "conv_units" to change units. ; See "time_stamp" to turn time stamping on and off. ; Future changes: fitfunct.pro, tempfit.pro, tempfit2.pro ; ; ;CREATED BY: J. McFadden 96-11-14 ;FILE: funct_fit3d.pro ;VERSION 1. ;LAST MODIFICATION: Ver 1.1 GTDelory 97-07-30 ;- function fu_fit2d,dat, $ ; NFUNCT = nfitf, $ ; FITF = fitf, $ ; AUTO = auto, $ ENERGY = en, $ ERANGE = er, $ EBINS = ebins, $ ANGLE = an, $ ARANGE = ar, $ BINS = bins, $ LIMITS = limits, $ UNITS = units, $ ; F_UNITS = f_units, $ ; COLOR = col, $ LABEL = label, $ ; SUMPLT = sumplt, $ ; EBARS = ebars, $ MSEC = msec common fit_mass,mass ; Exit if no data ; if data_type(dat) ne 8 or dat.valid eq 0 then begin print,'Invalid Data' return,[0,0,0,0] endif !y.omargin =[2,3] ; temporary fix mass=dat.mass*1.6e-22 ; dat.mass is a scaler with units of eV/(km/s)^2. Use 1.6e-12erg/eV * (km/1.e5cm)^2 to convert ; Determine angle range for fit ; default angle range is all bins ; nb=dat.nbins bins2=replicate(1b,nb) if keyword_set(an) then begin if ndimen(an) gt 1 then begin print,'Error - angle keyword must be fltarr(n)' endif else begin if dimen1(an) eq 1 then bins2=angle_to_bins(dat,[an,an]) if dimen1(an) eq 2 then bins2=angle_to_bins(dat,an) if dimen1(an) gt 2 then begin ibin=angle_to_bin(dat,an) bins2(*)=0 & bins2(ibin)=1 endif endelse endif if keyword_set(ar) then begin bins2(*)=0 if ar(0) gt ar(1) then begin bins2(ar(0):nb-1)=1 bins2(0:ar(1))=1 endif else begin bins2(ar(0):ar(1))=1 endelse endif if keyword_set(bins) then bins2=bins ; Start the plot ; if not keyword_set(units) then units='eflux' spec2d,dat, units=units, limits=limits, color=color, label=label, $ msec=msec ; Determine energy range for fit ; if keywords 'erange' or 'energy' not used then use crosshairs ; energy = dat.energy(*,0) if keyword_set(en) or keyword_set(er) then begin if keyword_set(en) then xx=en else xx=energy(er) endif else begin crosshairs,xx,yy endelse nxx = dimen1(xx) print,'nxx=',nxx ; If maxwellian include the acceleration if not keyword_set(fitf) then begin print,'Maxwellian Fit with potential drop - default function fit' epeak = xx(0) etmp = min(abs(energy - xx(0)),epeak_ind) print,'epeak=',epeak if nxx le 2 then begin er_ind = intarr(2) etmp = min(abs(energy - xx(0)),tmpindex) er_ind(0) = tmpindex if nxx eq 2 then begin if xx(1) gt xx(0) then begin etmp = min(abs(energy - xx(1)),tmpindex) er_ind(1) = tmpindex endif else er_ind(1) = 0 endif else er_ind(1) = 0 nxx = 3 print,'FLAG!' endif else begin er_ind = intarr(nxx-1) n=0 while n lt nxx-1 do begin etmp = min(abs(energy - xx(n+1)),tmpindex) er_ind(n)=tmpindex if n then begin if er_ind(n) gt er_ind(n-1) then begin tmpindex = er_ind(n) er_ind(n) = er_ind(n-1) er_ind(n-1) = tmpindex endif endif n=n+1 endwhile endelse endif else begin print,' Fitting function = ',fitf er_ind = intarr(nxx-1) n=0 while n lt nxx-1 do begin etmp = min(abs(energy - xx(n+1)),tmpindex) er_ind(n)=tmpindex if n then begin if er_ind(n) gt er_ind(n-1) then begin tmpindex = er_ind(n) er_ind(n) = er_ind(n-1) er_ind(n-1) = tmpindex endif endif n=n+1 endwhile endelse ; Determine the fitting routine ; if keyword_set(nfitf) then begin if fix(nxx/2) lt nfitf then nfitf = fix(nxx/2) endif else nfitf=fix(nxx/2) print, 'Number of functions = ',nfitf tfmaxwellian = 0 if not keyword_set(fitf) then begin fitf = 'maxwellian_'+ strcompress(string(nfitf),/remove_all) print,' Default function = ',fitf f_units='eflux' print,'FUNCT_FIT2d: Default units are '+f_units tfmaxwellian = 1 endif else begin print,'Use function = ',fitf f_units='df' tfmaxwellian = 0 endelse ; Get data for the fit if dat.nbins ne 1 then tmpdat = omni2d(dat,bins=bins2) else tmpdat=dat tmpdat = conv_units(tmpdat,f_units) xvar = tmpdat.energy yvar = tmpdat.data maxind=max(er_ind) minind=min(er_ind) xfit=reverse(xvar(minind:maxind)) yfit=reverse(yvar(minind:maxind)) ; Weighted fit according to IDL, and using ".ddata" structure returned ; from routine OMNI2D. Tends to hang CURVEFIT when .ddata is squared. wvar = (1./tmpdat.ddata)^2 wfit=reverse(wvar(minind:maxind)) ; Now fit the data ; call_procedure, fitf, xvar, afit, yvar, pder, index=er_ind, units=f_units print,'Initial guess : afit=',afit cfit = curvefit(xfit, yfit, wfit, afit, sigmaa, function_name=fitf, $ iter=iter, chi2=chi2) print,'Final fit : afit=',afit ; Print out fit information ; print,'iter=',iter print,'curvefit chi2=',chi2 chi2_2 = total((yfit-cfit)^2*wfit)/(n_elements(yfit) - $ N_elements(afit)) print,'FUNCT_FIT2d: CHI2 = '+strcompress(string(chi2_2), /remove_all) if tfmaxwellian then begin for m=0,nfitf-1 do begin print,' Population',m print,' T = ',-1./afit(2*m+1) print,' f0 = ',afit(2*m) & f0=afit(2*m) print,' source density =', $ afit(2*m)*exp(afit(2*m+1)*xx(0))*1.e-15 * $ (-1.*2.*3.1416*1.6e-12/(afit(2*m+1)*(mass)))^1.5 T0 =-1./afit(2*m+1) f0 = afit(2*m) n0 = afit(2*m)*exp(afit(2*m+1)*xx(0))*1.e-15 * $ (-1.*2.*3.1416*1.6e-12/(afit(2*m+1)*(mass)))^1.5 endfor endif if fitf eq 'maxwellian_beam' then print,'Vo = ',afit(2) ; Add the fit to the plot ; xfit=dat.energy(*,0) xfit=reverse(xfit) call_procedure, fitf, xfit, afit, yfit_max, pder, units=f_units chi2_2 = total((yfit_max-yvar)^2*wvar)/(N_elements(yvar) - $ N_elements(afit)) print,'Total data-to-fit chi2 = '+strcompress(string(chi2_2), $ /remove_all) if keyword_set(sumplt) then begin sumdat=omni2d(dat,bins=bins2) spec2d,sumdat endif else begin spec2d,dat, units=f_units, limits=limits, color=color, label=label, $ msec=msec, angle=an, arange=ar, bins=bins endelse if keyword_set(ebars) then begin error_high = tmpdat.data + tmpdat.ddata error_low = tmpdat.data - tmpdat.ddata errplot, tmpdat.energy, error_low, error_high, width=0 endif oplot,xfit,yfit_max return,[T0,n0,epeak,f0] end