; $Id: ft_data__define.pro,v 1.9 2005/02/08 18:49:48 dimeo Exp $
;###############################################################################
;
; NAME:
;  FT_DATA__DEFINE
;
; PURPOSE:
;  Class definition for the data used in the Fourier
;  Transform toolkit.
;
; CATEGORY:
;  DAVE, Data Analysis, FT Toolkit.
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
; NAME:
;
;  FT_DATA__DEFINE
;
; PURPOSE:
;
;  Object class for data in the Fourier Transform application.
;  Class enables user to perform forward and reverse FFT of
;  the DAVE data.  Also various filters can be used on the
;  the data such as: boxcar averaging (using the IDL SMOOTH
;  function), Gaussian smoothing, and the Savitzky-Golay
;  filter.
;
; AUTHOR:
;
;  Robert Dimeo
;  National Institute of Standards and Technology
;  Center for Neutron Research
;  100 Bureau Drive, Mail Stop 8562
;  Gaithersburg, MD 20899
;  Tel: (301) 975-8135
;  Email: robert.dimeo@nist.gov
;
; CATEGORY:
;
;  DAVE application
;
; INPUT PARAMETERS
;
;  None
;
; INPUT KEYWORDS
;
;  FILENAME -- Name of a DAVE file saved to disk (optional)
;
; METHODS
;
;  INIT (f), CLEANUP (p) -- Standard lifecycle methods
;  GET_PROPERTY (f) -- Accessor method using keywords
;  SET_PROPERTY (f) -- Accessor method using keywords
;  READ_DAVE (f) -- Restores a DAVE pointer previously saved on disk
;  SAVE_DAVE (f) -- Four files are saved to disk:  the four files
;                   correspond to the imaginary, real, magnitude, and
;                   phase of the time-correlation function.
;
; REQUIREMENTS
;
;  IDL 6.0 or higher
;
; ADDITIONAL REQUIRED PROGRAMS
;
;  GET_DAVE_PTR_CONTENTS.PRO
;  SHOWPROGRESS__DEFINE.PRO
;
; MODIFICATION HISTORY:
;
;  Written -- 3/3/04 (RMD)
; -
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ft_data::cleanup
if ptr_valid(self.dave_ptr) then heap_free,self.dave_ptr
ptr_free,self.qty_time_ptr,self.qty_energy_ptr
ptr_free,self.time_ptr,self.energy_ptr,self.plot_ptr
ptr_free,self.dqty_energy_ptr,self.dqty_time_ptr
ptr_free,self.yvals_ptr
ptr_free,self.dqty_time_orig_ptr,self.qty_time_orig_ptr
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::clone
oclone = obj_new('ft_data')

return,oclone
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::erase_energy
ptr_free,self.dqty_energy_ptr,self.qty_energy_ptr
ptr_free,self.energy_ptr
self.dqty_energy_ptr = ptr_new(/allocate_heap)
self.qty_energy_ptr = ptr_new(/allocate_heap)
self.energy_ptr = ptr_new(/allocate_heap)
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::save_dave,filename
if n_params() eq 0 then return,0
if ~ptr_valid(self.dave_ptr) then return,0
filename = filename + '_time'
;  Here add on nonsensical error bars
qty = *self.qty_time_ptr
qty_re = float(qty)
qty_im = imaginary(qty)

error_bars = n_elements(*self.dqty_time_ptr) ne 0

if error_bars then begin
   dqty = *self.dqty_time_ptr
   dqty_re = float(dqty)
   dqty_im = imaginary(dqty)
endif else begin
   strout = strarr(3)
   strout[0] = 'You must use the Monte-Carlo error estimate'
   strout[1] = 'for both the sample and resolution prior to'
   strout[2] = 'saving as DAVE files.  Returning...'
   void = dialog_message(dialog_parent = self.wid, strout)
   return,0
endelse

xvals = *self.time_ptr
n = n_elements(xvals)
n21 = n/2+1
nqty = size(qty,/dimensions)

if n_elements(nqty) eq 2 then begin
   for i = 0,nqty[1]-1 do begin
      qty[*,i] = shift(qty[*,i],-n21)
      qty_re[*,i] = shift(qty_re[*,i],-n21)
      qty_im[*,i] = shift(qty_im[*,i],-n21)
      if error_bars then begin
         dqty[*,i] = shift(dqty[*,i],-n21)
         dqty_re[*,i] = shift(dqty_re[*,i],-n21)
         dqty_im[*,i] = shift(dqty_im[*,i],-n21)
      endif
   endfor
endif else begin
   qty = shift(qty,-n21)
   qty_re = shift(qty_re,-n21)
   qty_im = shift(qty_im,-n21)
   if error_bars then begin
      dqty = shift(dqty,-n21)
      dqty_re = shift(dqty_re,-n21)
      dqty_im = shift(dqty_im,-n21)
   endif
endelse
qty_mag = sqrt(qty_re^2+qty_im^2)
qty_phase = (!radeg)*atan(qty_im/qty_re)
if error_bars then begin
   dqty_mag = abs(dqty)
   dqty_phase = qty_phase*sqrt((dqty_re/qty_re)^2+(dqty_im/qty_im)^2)
endif

pos = where(xvals ge 0.0,count_pos)
if count_pos gt 0 then begin
   if n_elements(nqty) eq 2 then begin
      qty_re = qty_re[pos,*]
      qty_im = qty_im[pos,*]
      qty_mag = qty_mag[pos,*]
      qty_phase = qty_phase[pos,*]
      if error_bars then begin
         dqty_re = dqty_re[pos,*]
         dqty_im = dqty_im[pos,*]
         dqty_mag = dqty_mag[pos,*]
         dqty_phase = dqty_phase[pos,*]
      endif
   endif else begin
      qty_re = qty_re[pos]
      qty_im = qty_im[pos]
      qty_mag = qty_mag[pos]
      qty_phase = qty_phase[pos]
      if error_bars then begin
         dqty_re = dqty_re[pos]
         dqty_im = dqty_im[pos]
         dqty_mag = dqty_mag[pos]
         dqty_phase = dqty_phase[pos]
      endif
   endelse
endif
err_re = dqty_re
err_im = dqty_im
err_mag = dqty_mag
err_phase = dqty_phase
xvals = xvals[pos]

yvals = (*(*(*self.dave_ptr).datastrptr).commonstr.histptr).y
hist_str = {qty:qty_mag,err:err_mag,x:xvals,y:yvals}
*(*(*self.dave_ptr).datastrptr).commonstr.histptr = hist_str
(*(*self.dave_ptr).datastrptr).commonstr.xlabel = 'time (ps)'
(*(*self.dave_ptr).datastrptr).commonstr.histlabel = 'ABS[I(t)]'
(*(*self.dave_ptr).datastrptr).commonstr.xunits = 'time:ps'
mag_filename = filename+'_mag.dave'
daveptr = self.dave_ptr
save,filename = mag_filename,daveptr

hist_str = {qty:qty_re,err:err_re,x:xvals,y:yvals}
*(*(*self.dave_ptr).datastrptr).commonstr.histptr = hist_str
(*(*self.dave_ptr).datastrptr).commonstr.xlabel = 'time (ps)'
(*(*self.dave_ptr).datastrptr).commonstr.histlabel = 'Re[I(t)]'
(*(*self.dave_ptr).datastrptr).commonstr.xunits = 'time:ps'
real_filename = filename+'_real.dave'
daveptr = self.dave_ptr
save,filename = real_filename,daveptr

hist_str = {qty:qty_im,err:err_im,x:xvals,y:yvals}
*(*(*self.dave_ptr).datastrptr).commonstr.histptr = hist_str
(*(*self.dave_ptr).datastrptr).commonstr.xlabel = 'time (ps)'
(*(*self.dave_ptr).datastrptr).commonstr.histlabel = 'Im[I(t)]'
(*(*self.dave_ptr).datastrptr).commonstr.xunits = 'time:ps'
im_filename = filename+'_imaginary.dave'
daveptr = self.dave_ptr
save,filename = im_filename,daveptr

hist_str = {qty:qty_phase,err:err_phase,x:xvals,y:yvals}
*(*(*self.dave_ptr).datastrptr).commonstr.histptr = hist_str
(*(*self.dave_ptr).datastrptr).commonstr.xlabel = 'time (ps)'
(*(*self.dave_ptr).datastrptr).commonstr.histlabel = 'Phase[I(t)]'
(*(*self.dave_ptr).datastrptr).commonstr.xunits = 'time:ps'
phase_filename = filename+'_phase.dave'
daveptr = self.dave_ptr
save,filename = phase_filename,daveptr

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::read_dave
; Loads in the contents of a DAVE data file saved on disk.
; This will only deal with
if ~file_test(self.filename) then return,0
restore,filename = self.filename
; The next statement is absolutely critical so that we do not
; leak any memory accidentally!!!
if ptr_valid(self.dave_ptr) then heap_free,self.dave_ptr
self.dave_ptr = daveptr ; pointer name on disk is daveptr

ret_val =   get_dave_ptr_contents(     $
            self.dave_ptr,             $
            QTY = qty,                 $
            ERR = err,                 $
            XVALS = xvals,             $
            XUNITS = xunits,           $
            YVALS = yvals,             $
            TREATMENT = treatment,     $
            ERMSG = errmsg             )

; If the units of x are not energy then free up the local
; DAVE_PTR and get out.
units = strupcase(xunits)
if (units ne 'ENERGY:UEV') and $
   (units ne 'ENERGY:MEV') then begin
   if ptr_valid(self.dave_ptr) then heap_free,self.dave_ptr
   return,0
endif

; If the bin widths are not equal then get out
nx = n_elements(xvals)
dx = xvals[1:nx-1] - xvals[0:nx-2]
dx_sig = sqrt((moment(dx))[1])
if dx_sig gt 10.0*(machar()).eps then begin
   strout = 'The spacing between energy points must be equal'
   void = dialog_message(dialog_parent = self.wid,strout)
   if ptr_valid(self.dave_ptr) then heap_free,self.dave_ptr
   return,0
endif

qty_size = size(qty,/dimensions)
if n_elements(qty_size) eq 2 then begin
   n = qty_size[0]
;   if ((n mod 2) eq 1) then begin
;      qty = qty[0:n-2,*]
;      err = err[0:n-2,*]
;      xvals = xvals[0:n-2]
;   endif
   n = n_elements(xvals)
   n21 = n/2 + 1
   qty = shift(qty,-n21,0)
   err = shift(err,-n21,0)
endif
if n_elements(qty_size) eq 1 then begin
   qty = reform(qty)
   n = n_elements(qty)
;   if ((n mod 2) eq 1) then begin
;      qty = qty[0:n-2]
;      err = err[0:n-2]
;      xvals = xvals[0:n-2]
;   endif
   n = n_elements(xvals)
   n21 = n/2 + 1
   qty = shift(qty,-n21)
   err = shift(err,-n21)
endif

*self.yvals_ptr = yvals
ny = n_elements(yvals)
*self.sub_ptr = fltarr(2,ny)
*self.qty_energy_ptr = qty
*self.dqty_energy_ptr = err
if strupcase(xunits) eq 'ENERGY:UEV' then xvals = 1.e-3*xvals
*self.energy_ptr = xvals
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::get_property,  filename = filename,          $
                                 display_name = display_name,  $
                                 nx = nx,                      $
                                 ny = ny,                      $
                                 energy = energy,              $
                                 qty_energy = qty_energy,      $
                                 qty_time = qty_time,          $
                                 time = time,                  $
                                 qvals = qvals,                $
                                 dqty_energy = dqty_energy,    $
                                 dqty_time = dqty_time,        $
                                 nmc = nmc,                    $
                                 wid = wid,                    $
                                 sub_array = sub_array,        $
                                 range = range
filename = self.filename
nmc = self.nmc
range = self.range
wid = self.wid
if arg_present(sub_array) and n_elements(*self.sub_ptr) ne 0 then $
   sub_array = *self.sub_ptr

if arg_present(dqty_time) then begin
   if n_elements(*self.dqty_time_ptr) eq 0 then $
      dqty_time = 0.0 else $
      dqty_time = *self.dqty_time_ptr
endif
if arg_present(qvals) then begin
   if n_elements(*self.yvals_ptr) eq 0 then $
      qvals = 0.0 else $
      qvals = *self.yvals_ptr
endif
if arg_present(dqty_energy) then begin
   if n_elements(*self.dqty_energy_ptr) eq 0 then $
      dqty_energy = 0.0 else $
      dqty_energy = *self.dqty_energy_ptr
endif
if arg_present(time) then begin
   if n_elements(*self.time_ptr) eq 0 then $
      time = 0.0 else $
      time = *self.time_ptr
endif
if arg_present(qty_energy) then begin
   if n_elements(*self.qty_energy_ptr) eq 0 then $
      qty_energy = 0.0 else $
      qty_energy = *self.qty_energy_ptr
endif
if arg_present(qty_time) then begin
   if n_elements(*self.qty_time_ptr) eq 0 then $
      qty_time = 0.0 else $
      qty_time = *self.qty_time_ptr
endif
if arg_present(display_name) then begin
   if self.filename ne '' then begin
      display_name = file_basename(self.filename,'.dave')
   endif else begin
      display_name = ''
   end
endif
if arg_present(nx) or arg_present(ny) then begin
   if n_elements(*self.qty_energy_ptr) eq 0 then return,0
   qty_size = size(*self.qty_energy_ptr,/dimensions)
   nx = qty_size[0]
   if n_elements(qty_size) eq 2 then ny = qty_size[1] else $
      ny = 1
endif
if arg_present(energy) then begin
   if n_elements(*self.energy_ptr) eq 0 then $
      energy = 0.0 else $
      energy = *self.energy_ptr
endif
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_calc_time,etemp
energy = etemp/4.136
nenergy = n_elements(energy)
de = energy[1]-energy[0]
n21 = nenergy/2 + 1
; The array of subscripts:
t = indgen(nenergy)
; Insert negative frequencies in elements F(N/2 +1), ..., F(N-1):
if (nenergy mod 2) eq 0 then val = 2 else val = 1
t[n21] = n21 - nenergy + findgen(n21-val)
; Compute T0 frequency:
t = shift((t/(nenergy*de)),-n21)
return,t
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_calc_energy,ttemp
return,ft_calc_time(ttemp)
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::draw,_Extra = extra
if n_elements(*self.plot_ptr) eq 0 then return,0
x = (*self.plot_ptr).x
y = (*self.plot_ptr).y
names = tag_names(*self.plot_ptr)
dy_index = where(strupcase(names) eq 'DY',count_dy)
if count_dy ne 0 then dy = (*self.plot_ptr).dy

xtitle = (*self.plot_ptr).xtitle
ytitle = (*self.plot_ptr).ytitle
if self.autoscale then begin
   yhi = max(y,min = ylo)
   xhi = max(x,min = xlo)
   dely = (yhi-ylo)
   delx = (xhi-xlo)
   fx = 0.1 & fy = 0.1
   self.xrange = [xlo-fx*delx,xhi+fx*delx]
   self.yrange = [ylo-fy*dely,yhi+fy*dely]
   if keyword_set(extra.ylog) then $
      self.yrange[0] = ylo > 1.e-9
   if keyword_set(extra.xlog) then begin
      xpos = where(x gt 0.0)
      self.xrange[0] = 0.1*x[xpos[0]]
   endif
endif
color_exist = where(strupcase(tag_names(extra)) eq 'COLOR',count_colors)
pos1 = [0.1,0.15,0.9,0.85]
if ~keyword_set(extra.overplot) then begin
   tags = tag_names(extra)
   xlog_there = where(tags eq 'XLOG',count_xlog)
   if count_xlog ne 0 then begin
      if extra.xlog then begin
         pos = where(x gt 0.0)
         x = x[pos] & y = y[pos]
         if n_elements(dy) ne 0 then dy = dy[pos]
      endif
   endif
   plot,x,y,xtitle = xtitle,ytitle = ytitle,xrange = self.xrange, $
      yrange = self.yrange,/xsty,/ysty,title = self.title, _Extra = extra, $
      position = pos1
   if n_elements(dy) ne 0 and count_colors gt 0 then $
      color_errplot,x,y-dy,y+dy,width = 0.0,color = extra.color
endif else begin
   tags = tag_names(extra)
   xlog_there = where(tags eq 'XLOG',count_xlog)
   if count_xlog ne 0 then begin
      if extra.xlog then begin
         pos = where(x gt 0.0)
         x = x[pos] & y = y[pos]
         if n_elements(dy) ne 0 then dy = dy[pos]
      endif
   endif
   oplot,x,y,_Extra = extra
   if n_elements(dy) ne 0 and count_colors gt 0 then $
      color_errplot,x,y-dy,y+dy,width = 0.0,color = extra.color
endelse
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::get_plot_info, energy = energy,        $
                                 time = time,            $
                                 index = index,          $
                                 phase = phase,          $
                                 magnitude = magnitude,  $
                                 real = real,            $
                                 imag = imag,            $
                                 range = range

range = self.range
energy = self.energy
time = self.time
index = self.index
phase = self.phase
magnitude = self.magnitude
real = self.real
imag = self.imag
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::set_plot_info, energy = energy,        $
                                 time = time,            $
                                 index = index,          $
                                 phase = phase,          $
                                 magnitude = magnitude,  $
                                 real = real,            $
                                 imag = imag,            $
                                 xrange = xrange,        $
                                 yrange = yrange,        $
                                 autoscale = autoscale

if n_elements(xrange) ne 0 then self.xrange = xrange
if n_elements(yrange) ne 0 then self.yrange = yrange
if n_elements(autoscale) ne 0 then self.autoscale = autoscale

if n_elements(index) eq 0 then index = 0
if n_elements(index) ne 0 then begin
   self.index = index
   format = '(f8.2)'
   qvalue = strtrim(string((*self.yvals_ptr)[index],  $
      format = format),2)
   sym = '!3!sA!r!u!9 %!3!n!E-1!N'
   title = 'Q = '+qvalue+' '+sym
   self.title = title
endif
if keyword_set(energy) then begin
   self.energy = 1B & self.time = 0B
   if n_elements(*self.qty_energy_ptr) eq 0 then return,0
   qty = *self.qty_energy_ptr
   dqty = *self.dqty_energy_ptr
   x = *self.energy_ptr
   qty_size = size(qty,/dimensions)

   if n_elements(qty_size) eq 1 then begin
      y = qty
      dy = dqty
   endif else begin
      y = reform(qty[*,index])
      dy = reform(dqty[*,index])
   endelse

   n = n_elements(x)
   n21 = n/2 + 1
   if n_elements(dy) ne 0 then dy = shift(dy,-n21)
   y_im = shift(imaginary(y),-n21)
   y_re = shift(float(y),-n21)
   xtitle = 'E (meV)'

   if keyword_set(real) then begin
      self.real = 1B & self.phase = (self.magnitude = (self.imag = 0B))
      yp = y_re
      dyp = float(dy)
      ytitle = 'Re(I(E))'
   endif
   if keyword_set(imag) then begin
      self.imag = 1B & self.phase = (self.magnitude = (self.real = 0B))
      yp = y_im
      dyp = imaginary(dy)
      ytitle = 'Im(I(E))'
   endif
   if keyword_set(phase) then begin
      self.phase = 1B & self.real = (self.magnitude = (self.imag = 0B))
      yp = (!radeg)*(atan(y_im/y_re))
      dy_re = float(dy) & dy_im = imaginary(dy)
      dyp = (!radeg)*yp*sqrt((dy_im/y_im)^2+(dy_re/y_re)^2)
      ytitle = '<I(E)'
   endif
   if keyword_set(magnitude) then begin
      self.magnitude = 1B & self.phase = (self.real = (self.imag = 0B))
      yp = sqrt(y_im^2+y_re^2)
      dyp = abs(dy)
      ytitle = 'ABS(I(E))'
   endif
   *self.plot_ptr = {x:x,y:yp,dy:dyp,xtitle:xtitle,ytitle:ytitle}
endif
if keyword_set(time) then begin
   self.time = 1B & self.energy = 0B
   if n_elements(*self.qty_time_ptr) eq 0 then return,0
   if n_elements(*self.time_ptr) eq 0 then return,0
   qty = *self.qty_time_ptr
   x = *self.time_ptr
   qty_size = size(qty,/dimensions)
   n = n_elements(x)
   n21 = n/2 + 1
   if n_elements(qty_size) eq 1 then begin
      y = qty
      if n_elements(*self.dqty_time_ptr) ne 0 then begin
         err_t = *self.dqty_time_ptr
         dy_re = shift(float(err_t),-n21)
         dy_im = shift(imaginary(err_t),-n21)
      endif
   endif else begin
      y = qty[*,index]
      if n_elements(*self.dqty_time_ptr) ne 0 then begin
         err_t = reform((*self.dqty_time_ptr)[*,index])
         dy_re = shift(float(err_t),-n21)
         dy_im = shift(imaginary(err_t),-n21)
      endif
   endelse
   xtitle = 't (ps)'


   y_im = shift(imaginary(y),-n21)
   y_re = shift(float(y),-n21)

   if keyword_set(phase) then begin
      self.phase = 1B & self.real = (self.magnitude = (self.imag = 0B))
      yp = (!radeg)*(atan(y_im/y_re))
      if n_elements(dy_re) ne 0 then dyp = yp*sqrt((dy_re/y_re)^2+(dy_im/y_im)^2)
      ytitle = '<I(t)'
   endif
   if keyword_set(magnitude) then begin
      self.magnitude = 1B & self.phase = (self.real = (self.imag = 0B))
      yp = sqrt(y_im^2+y_re^2)
      if n_elements(dy_re) ne 0 then dyp = sqrt(dy_re^2+dy_im^2)
      ytitle = 'ABS(I(t))'
   endif
   if keyword_set(real) then begin
      self.real = 1B & self.phase = (self.magnitude = (self.imag = 0B))
      yp = y_re
      if n_elements(dy_re) ne 0 then dyp = dy_re
      ytitle = 'Re(I(t))'
   endif
   if keyword_set(imag) then begin
      self.imag = 1B & self.phase = (self.magnitude = (self.real = 0B))
      yp = y_im
      if n_elements(dy_im) ne 0 then dyp = dy_im
      ytitle = 'Im(I(t))'
   endif
   pos = where(x ge 0.0,count_pos)
   if count_pos gt 0 then begin
      x = x[pos]
      yp = yp[pos]
      if n_elements(dyp) ne 0 then dyp = dyp[pos]
   endif
   if n_elements(dyp) gt 0 then begin
      *self.plot_ptr = {x:x,y:yp,dy:dyp,xtitle:xtitle,ytitle:ytitle}
   endif else begin
      *self.plot_ptr = {x:x,y:yp,xtitle:xtitle,ytitle:ytitle}
   endelse
endif

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::reverse_fft
; This method computes I(t) from I(E)
if ~ptr_valid(self.dave_ptr) then return,0
if n_elements(*self.qty_energy_ptr) eq 0 then return,0

qty = *self.qty_energy_ptr
dqty = *self.dqty_energy_ptr
energy = *self.energy_ptr
range = self.range
if n_elements(range) eq 0 then begin
   count_energy = n_elements(energy)
   e_indices = indgen(count_energy)
endif else begin
   elo = range[0] & ehi = range[1]
   e_indices = where(((energy) ge elo) and ((energy) le ehi), $
               count_energy)
endelse
; Now check how many elements have been selected.  After all we
; do not want to calculate the inverse Fourier Transform for a
; function that is composed of less than 3 points (arbitrarily
; chosen of course).
if count_energy le 2 then return,0
; ...and check if the number of elements in the calculation of
; the FFT are even.  If not then add an element to either end
; depending on which side has been truncated.
;if count_energy mod 2 ne 0 then begin
;   if e_indices[0] ne energy[0] then begin
;      e_indices = [e_indices[0]-1,e_indices]
;   endif else begin
;      e_indices = [e_indices,e_indices[count_energy-1]+1]
;   endelse
;endif

; We need to flip qty and dqty prior to performing the
; range selection (added next 5 lines by RMD 7/12/04)
n = n_elements(energy)
n21 = n/2
new_qty = shift(qty,-n21)
new_dqty = shift(dqty,-n21)
qty = new_qty[e_indices,*]
dqty = new_dqty[e_indices,*]


qty = qty[e_indices,*]
dqty = dqty[e_indices,*]

; Now flip them back....(added next two lines by RMD 7/12/04)
qty = shift(qty,n21)
dqty = shift(dqty,n21)

energy = energy[e_indices]
direction = 1
yfft = fft(qty,direction,dimension = 1,double = 1)

t = ft_calc_time(energy)
tsort = sort(t) & t = t[tsort]
yfft = yfft[tsort,*]
*self.qty_time_ptr = yfft
*self.time_ptr = t
*self.qty_time_orig_ptr = yfft
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::forward_fft
; This method computes I(E) from I(t)
if ~ptr_valid(self.dave_ptr) then return,0
if n_elements(*self.qty_time_ptr) eq 0 then return,0

qty = *self.qty_time_ptr
if n_elements(*self.dqty_time_ptr) ne 0 then $
   dqty = *self.dqty_time_ptr
time = *self.time_ptr

direction = -1
yfft = fft(qty,direction,dimension = 1,double = 1)

e = ft_calc_energy(time)
esort = sort(e) & e = e[esort]
yfft = yfft[esort,*]
*self.qty_energy_ptr = yfft
*self.energy_ptr = e

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::mc_estimate_errors
; Use Fanning's progress bar utility
msg = 'Estimating errors using the "bootstrap" Monte-Carlo method'
title = 'Progress'
progressBar = Obj_New("SHOWPROGRESS",message = msg,title = title, $
   self.wid)
progressBar->Start

nmc = self.nmc
yfft = *self.qty_time_ptr
qty = *self.qty_energy_ptr
err = *self.dqty_energy_ptr
energy = *self.energy_ptr
range = self.range

if n_elements(range) eq 0 then begin
   count_energy = n_elements(energy)
   e_indices = indgen(count_energy)
endif else begin
   elo = range[0] & ehi = range[1]
   e_indices = where(((energy) ge elo) and ((energy) le ehi), $
               count_energy)
endelse
energy = energy[e_indices]
sz = size(qty)
if sz[0] eq 2 then begin
   qty = qty[e_indices,*]
   err = err[e_indices,*]
endif else begin
   qty = qty[e_indices]
   err = err[e_indices]
endelse
sz = size(yfft,/dimensions)
if n_elements(sz) gt 1 then begin
   qty_re_save = fltarr(nmc,sz[0],sz[1])
   qty_im_save = fltarr(nmc,sz[0],sz[1])
endif else begin
   qty_re_save = fltarr(nmc,sz[0])
   qty_im_save = fltarr(nmc,sz[0])
endelse

for i = 0L,nmc-1 do begin
   if n_elements(sz) gt 1 then begin
      r_normal = randomn(s,sz[0],sz[1])
   endif else begin
      r_normal = randomn(s,sz[0])
   endelse
   derr = err*r_normal
   this_qty = qty+derr
   ret = self->set_property(qty_energy = this_qty)
   ret = self->reverse_fft()
   ret = self->get_property(qty_time = qty_time)
   qt_re = float(qty_time)
   qt_im = imaginary(qty_time)
   if n_elements(sz) eq 2 then begin
      qty_re_save[i,*,*] = qt_re
      qty_im_save[i,*,*] = qt_im
   endif else begin
      qty_re_save[i,*] = qt_re
      qty_im_save[i,*] = qt_im
   endelse
   progressBar->Update, (i+1)*100/nmc
endfor

dy_re = 0.*err
dy_im = 0.*err
if n_elements(sz) eq 2 then begin
   for j = 0L,sz[1]-1 do begin
      for i = 0L,sz[0]-1 do begin
         re_result = moment(qty_re_save[*,i,j])
         im_result = moment(qty_im_save[*,i,j])
         dy_re[i,j] = sqrt(re_result[1])
         dy_im[i,j] = sqrt(im_result[1])
      endfor
   endfor
endif else begin
   for i = 0L,sz[0]-1 do begin
      re_result = moment(qty_re_save[*,i])
      im_result = moment(qty_im_save[*,i])
      dy_re[i] = sqrt(re_result[1])
      dy_im[i] = sqrt(im_result[1])
   endfor
endelse
ret = temporary(qty_re_save) & ret = temporary(qty_im_save)
*self.dqty_time_ptr = complex(dy_re,dy_im)
ret = self->set_property(qty_energy = qty,qty_time = yfft)

progressBar->Destroy
Obj_Destroy, progressBar

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::remove_error_bars
ptr_free,self.dqty_time_ptr
self.dqty_time_ptr = ptr_new(/allocate_heap)
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::set_property,  filename = filename,       $
                                 time = time,               $
                                 energy = energy,           $
                                 qty_energy = qty_energy,   $
                                 qty_time = qty_time,       $
                                 dqty_time = dqty_time,     $
                                 nmc = nmc,                 $
                                 range = range,             $
                                 sub_array = sub_array,     $
                                 wid = wid

if n_elements(sub_array) ne 0 then *self.sub_ptr = sub_array
if n_elements(wid) ne 0 then self.wid = wid
if n_elements(dqty_time) ne 0 then *self.dqty_time_ptr = dqty_time
if n_elements(range) ne 0 then self.range = range
if n_elements(nmc) ne 0 then self.nmc = nmc
if n_elements(time) ne 0 then *self.time_ptr = time
if n_elements(energy) ne 0 then *self.energy_ptr = energy
if n_elements(qty_energy) ne 0 then $
   *self.qty_energy_ptr = qty_energy
if n_elements(qty_time) ne 0 then $
   *self.qty_time_ptr = qty_time
if n_elements(filename) ne 0 then self.filename = filename
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::boxcar_average,width,   time = time, $
                                          energy = energy
; Error checking
if ~ptr_valid(self.dave_ptr) then return,0
if n_params() eq 0 then return,0
if width lt 2 then return,0
if keyword_set(time) and keyword_set(energy) then return,0
if ~keyword_set(time) and ~keyword_set(energy) then return,0
if keyword_set(time) then begin
   if n_elements(*self.qty_time_ptr) eq 0 then return,0
   qty = *self.qty_time_ptr
   x = *self.time_ptr
endif
if keyword_set(energy) then begin
   qty = *self.qty_energy_ptr
   x = *self.energy_ptr
endif

qty_size = size(qty,/dimensions)
if n_elements(qty_size) eq 2 then begin
   for i = 0,qty_size[1]-1 do begin
      qty[*,i] = smooth(qty[*,i],width)
   endfor
endif else begin
   qty[*] = smooth(qty[*],width)
endelse

if keyword_set(time) then *self.qty_time_ptr = qty
if keyword_set(energy) then *self.qty_energy_ptr = qty
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::gaussian_smooth,width,  time = time,   $
                                          energy = energy
; Width is in pixels
; Error checking
if ~ptr_valid(self.dave_ptr) then return,0
if n_params() eq 0 then return,0
if width lt 3 then return,0
if keyword_set(time) and keyword_set(energy) then return,0
if ~keyword_set(time) and ~keyword_set(energy) then return,0
if keyword_set(time) then begin
   if n_elements(*self.qty_time_ptr) eq 0 then return,0
   qty = *self.qty_time_ptr
   x = *self.time_ptr
endif
if keyword_set(energy) then begin
   qty = *self.qty_energy_ptr
   x = *self.energy_ptr
endif

qty_size = size(qty,/dimensions)
sigma = width/2.354
n = n_elements(x)
nn = findgen(3*width)
num_n = n_elements(nn)
filter = (1.0/sqrt(2.0*!dpi*sigma^2))*exp(-0.5*((nn-(0.5*num_n))/sigma)^2)

qty_size = size(qty,/dimensions)
if n_elements(qty_size) eq 2 then begin
   for i = 0,qty_size[1]-1 do begin
      qty[*,i] = convol(qty[*,i], filter, /edge_truncate)
   endfor
endif else begin
   qty[*] = convol(qty[*], filter, /edge_truncate)
endelse

if keyword_set(time) then *self.qty_time_ptr = qty
if keyword_set(energy) then *self.qty_energy_ptr = qty
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::hanning_smooth,width,   time = time,   $
                                          energy = energy
; Width is in pixels
; Error checking
if ~ptr_valid(self.dave_ptr) then return,0
if n_params() eq 0 then return,0
if width lt 3 then return,0
if keyword_set(time) and keyword_set(energy) then return,0
if ~keyword_set(time) and ~keyword_set(energy) then return,0
if keyword_set(time) then begin
   if n_elements(*self.qty_time_ptr) eq 0 then return,0
   qty = *self.qty_time_ptr
   x = *self.time_ptr
endif
if keyword_set(energy) then begin
   qty = *self.qty_energy_ptr
   x = *self.energy_ptr
endif

qty_size = size(qty,/dimensions)
filter = hanning(width)
qty_size = size(qty,/dimensions)
if n_elements(qty_size) eq 2 then begin
   for i = 0,qty_size[1]-1 do begin
      qty[*,i] = convol(qty[*,i], filter, /edge_truncate)
   endfor
endif else begin
   qty[*] = convol(qty[*], filter, /edge_truncate)
endelse

if keyword_set(time) then *self.qty_time_ptr = qty
if keyword_set(energy) then *self.qty_energy_ptr = qty
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::savitzky_golay,nleft,nright,degree, $
                                 time = time,         $
                                 energy = energy
if ~ptr_valid(self.dave_ptr) then return,0
if n_params() lt 3 then return,0
if keyword_set(time) and keyword_set(energy) then return,0
if ~keyword_set(time) and ~keyword_set(energy) then return,0

if keyword_set(time) then begin
   if n_elements(*self.qty_time_ptr) eq 0 then return,0
   qty = *self.qty_time_ptr
   x = *self.time_ptr
endif
if keyword_set(energy) then begin
   qty = *self.qty_energy_ptr
   x = *self.energy_ptr
endif

qty_size = size(qty,/dimensions)
filter = savgol( nleft, nright, 0, degree)

if n_elements(qty_size) eq 2 then begin
   for i = 0,qty_size[1]-1 do begin
      qty[*,i] = convol(qty[*,i], filter, /edge_truncate)
   endfor
endif else begin
   qty[*] = convol(qty[*], filter, /edge_truncate)
endelse

if keyword_set(time) then *self.qty_time_ptr = qty
if keyword_set(energy) then *self.qty_energy_ptr = qty
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::original_scaling
error_bars = n_elements(*self.dqty_time_ptr) ne 0
if error_bars then begin
   ret = self->set_property(qty_time = *self.qty_time_orig_ptr, $
      dqty_time = *self.dqty_time_orig_ptr)
endif else begin
   ret = self->set_property(qty_time = *self.qty_time_orig_ptr)
endelse
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::scale_to_one
; This method scales the data to be displayed in the time
; domain to 1 at the point t=0.
error_bars = n_elements(*self.dqty_time_ptr) ne 0
if error_bars then begin
   ret = self->get_property(qty_time = qty_time,time = time, $
               dqty_time = dqty_time)
endif else begin
   ret = self->get_property(qty_time = qty_time,time = time)
endelse

n = n_elements(time)
n21 = n/2 + 1
tmin = min(abs(time),t_min_index)
nqty = size(qty_time,/dimension)
if n_elements(nqty) eq 2 then begin
   qty_time = shift(qty_time,-n21,0)
   if error_bars then dqty_time = shift(dqty_time,-n21,0)
   ux = 1+bytarr(nqty[0])
   uy = 1+bytarr(nqty[1])
   norm_factors = abs(reform(qty_time[t_min_index,*]))
   bad_ones = where(norm_factors eq 0.0,count_bad)
   if count_bad gt 0 then norm_factors[bad_ones] = 1.0
   normalization = ux#norm_factors
   qty_time = shift(qty_time/normalization,n21,0)
   if error_bars then $
      dqty_time = shift(dqty_time/normalization,n21,0)
endif else begin
   qty_time = shift(qty_time,-n21)
   if error_bars then dqty_time = shift(dqty_time,-n21)
   normalization = abs(qty_time[t_min_index])
   bad_ones = where(normalization eq 0.0,count_bad)
   if count_bad gt 0 then normalization[count_bad] = 1.0
   if error_bars then $
      dqty_time = shift(dqty_time/normalization,n21)
   qty_time = shift(qty_time/normalization,n21)
endelse
if ~error_bars then $
   ret = self->set_property(qty_time = qty_time) $
else $
   ret = self->set_property(qty_time = qty_time,dqty_time = dqty_time)
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::map_values,OBJ_REF
; This function maps the "energy" data of SELF onto the energy
; grid of OBJ_REF
if n_params() eq 0 then return,0
if obj_class(obj_ref) ne 'FT_DATA' then return,0
ret = obj_ref->get_property(energy = e_2,qty_energy = qty_2)
err2 = 1.0*qty_2
ret = self->get_property(energy = e_1,qty_energy = qty_1)
err1 = 1.0*qty_1

; Determine the domain of the common energy values
elo = max([min(e_1),min(e_2)])
ehi = min([max(e_1),max(e_2)])

e2_indices = where((e_2 ge elo) and (e_2 le ehi),count_energy_2)
e1_indices = where((e_1 ge elo) and (e_1 le ehi),count_energy_1)

if count_energy_2 lt 2 then return,0
if count_energy_1 lt 2 then return,0

; Now rebin

x_in = e_1[e1_indices]
x_out = e_2[e2_indices]
z_in = qty_1[e1_indices,*]
dz_in = err1[e1_indices,*]

drebin,x_in,z_in,dz_in,x_out,z_out,dz_out,err=err,emsg=emsg,$
   /points,/to_points

if err ne 0 then begin
   print,emsg
   return,0
endif

ret = self->set_property(energy = x_out,qty_energy = z_out)
ret = obj_ref->set_property(energy = x_out,qty_energy = qty_2[e2_indices,*])

return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function ft_data::init,filename = filename,wid = wid

if n_elements(filename) eq 0 then self.filename = '' else $
   self.filename = filename
; If the filename is valid then read it in.
self.xrange = [0.0,1.0] & self.yrange = [0.0,1.0]
self.range = [0.0,1.0]
self.autoscale = 1B
self.yvals_ptr = ptr_new(/allocate_heap)
self.plot_ptr = ptr_new(/allocate_heap)
self.qty_time_ptr = ptr_new(/allocate_heap)
self.qty_energy_ptr = ptr_new(/allocate_heap)
self.dqty_energy_ptr = ptr_new(/allocate_heap)
self.dqty_time_ptr = ptr_new(/allocate_heap)
self.qty_time_orig_ptr = ptr_new(/allocate_heap)
self.dqty_time_orig_ptr = ptr_new(/allocate_heap)
self.energy_ptr = ptr_new(/allocate_heap)
self.time_ptr = ptr_new(/allocate_heap)
self.sub_ptr = ptr_new(/allocate_heap)
self.energy = 1B
self.time = 0B
self.index = 0
self.phase = 0B
self.magnitude = 1B
self.real = 0B
self.imag = 0B
self.nmc = 50L
if n_elements(wid) ne 0 then self.wid = wid else $
   self.wid = 0L
self.title = ''
if self.filename ne '' then ret = self->read_dave()
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ft_data__define
void =   {ft_data,   dave_ptr:ptr_new(),           $
                     filename:'',                  $
                     xrange:fltarr(2),             $
                     yrange:fltarr(2),             $
                     range:fltarr(2),              $
                     autoscale:1B,                 $
                     energy:1B,                    $
                     time:0B,                      $
                     index:0,                      $
                     phase:0B,                     $
                     magnitude:1B,                 $
                     real:0B,                      $
                     imag:0B,                      $
                     nmc:0L,                       $
                     title:'',                     $
                     yvals_ptr:ptr_new(),          $
                     sub_ptr:ptr_new(),            $
                     plot_ptr:ptr_new(),           $
                     qty_time_ptr:ptr_new(),       $
                     qty_energy_ptr:ptr_new(),     $
                     dqty_energy_ptr:ptr_new(),    $
                     dqty_time_ptr:ptr_new(),      $
                     qty_time_orig_ptr:ptr_new(),  $
                     dqty_time_orig_ptr:ptr_new(), $
                     time_ptr:ptr_new(),           $
                     wid:0L,                       $
                     energy_ptr:ptr_new()          }

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;