;+
;PROCEDURE: thm_sst_energy_extrapolate
;PURPOSE: Performs linear extrapolation of sst data to a new set of lower energies
;
;INPUTS:
; dist_data:
; The sst data structure on which the extrapolation should be performed. (loaded by thm_part_dist_array)
;KEYWORDS:
;
; add_energy=add_energy : Adds these energies to the current set of energies for the particle data when extrapolating.
; lin_energy=lin_energy : Set this keyword to perform extrapolation on energy, not logarithmic.
; lin_counts=lin_counts: Set this keyword to perform extrapolation on counts, not logarithmic
; lsquadratic=lsquadratic: Set this keyword to perform least square quadratic extrapolation of count data(see interpol documentation in IDL help.)
; quadratic=quadratic: Set this keyword to perform quadratic extrapolation of count data.(see interpol documentation in IDL help.)
; lsquares=lsquares: Set this keyword to the number of bins that you want to use for least squares extrapolation of the count data(Uses poly_fit)
; spline=spline: Set this keyword to perform spline extrapolation of count data.(see interpol documentation in IDL help.)
; trange=trange: Set this keyword to a two element array specifying a subset of the data that the operation should be performed on.(Don't need to modify the whole thing with the same parameters)
; error=error: Returns 1 if an error occurred. Returns 0 if operation completed successfully.
;
;
; bin_select: set the bin numbers that you want to use in the extrapolation
;EXAMPLES:
; dist_data = thm_part_dist_array(probe='a',type='psef',trange=time_double(['2012-02-08/09','2012-02-08/12']))
; thm_part_conv_units,dist_data,error=e
; thm_esa_energy_extrapolate,dist_data,add_energy=[5000,7000,10000]
;
; NOTES:
;
; $LastChangedBy: pcruce $
; $LastChangedDate: 2012-12-13 16:23:14 -0800 (Thu, 13 Dec 2012) $
; $LastChangedRevision: 11355 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/tags/tdas_8_00/idl/themis/spacecraft/particles/SST/thm_sst_energy_extrapolate.pro $
;-
pro thm_sst_energy_extrapolate,dist_data,add_energy,lin_energy=lin_energy,lin_counts=lin_counts,lsquadratic=lsquadratic,quadratic=quadratic,spline=spline,lsquares=lsquares,error=error,bin_select=bin_select
error = 1
if size(dist_data,/type) ne 10 then begin
dprint,dlevel=1,"ERROR: dist_data undefined or has wrong type"
return
endif
if keyword_set(lsquares) && ~is_num(lsquares) then begin
dprint,dlevel=1,"ERROR: lsquares is set but is not of type number"
return
endif
for i = 0,n_elements(dist_data)-1 do begin
str_template=(*dist_data[i])[0]
str_energy_dim = dimen(str_template.energy)
;create source target
energy_target = str_template.energy[0:str_energy_dim[0]-1,0]
energy_source = str_template.energy[0:str_energy_dim[0]-1,*]
energy_target = [add_energy[bsort(add_energy)],energy_target]
;Just uses a trick, uses the matrix multiply with an array of ones to create a two dim array with repeated rows(or cols, depending on your notation) from the linear array
new_energy_arr = energy_target#(dblarr(str_energy_dim[1])+1)
new_denergy_arr = [(abs(deriv(energy_target)))[0:n_elements(add_energy)-1],str_template.denergy[0:str_energy_dim[0]-1,0]]#(dblarr(str_energy_dim[1])+1) ;very quick and very dirty delta-energies
new_energy_dim = dimen(new_energy_arr)
str_template.nenergy = new_energy_dim[0]
str_element,str_template,'energy',new_energy_arr,/add_replace
str_element,str_template,'denergy',new_denergy_arr,/add_replace
tag_names=tag_names(str_template)
for j=0,n_elements(tag_names)-1 do begin
;these are special cases, so we skip them.
if strlowcase(tag_names[j]) ne 'energy' and $
strlowcase(tag_names[j]) ne 'denergy'and $
strlowcase(tag_names[j]) ne 'nenergy' then begin
tag_dimen=dimen(str_template.(j))
if array_equal(tag_dimen,str_energy_dim) then begin
new_data_arr = dblarr(new_energy_dim)
str_element,str_template,tag_names[j],new_data_arr,/add_replace
endif else begin
str_element,str_template,tag_names[j],(*dist_data[i])[0].(j)
endelse
endif
endfor
out_data = replicate(str_template,n_elements(*dist_data[i]))
;take log of target & source data if requested, makes subsequent fitting code a little easier
if ~keyword_set(lin_energy) then begin
energy_target=alog(energy_target)
energy_source=alog(energy_source)
endif
;loop over distributions
for j = 0,n_elements(out_data)-1 do begin
for k=0,n_elements(tag_names)-1 do begin
;these are special cases, so we skip them.
if strlowcase(tag_names[k]) ne 'energy' and $
strlowcase(tag_names[k]) ne 'denergy'and $
strlowcase(tag_names[k]) ne 'nenergy'and $
strlowcase(tag_names[k]) ne 'data' then begin
tag_dimen=dimen(str_template.(k))
if array_equal(tag_dimen,new_energy_dim) then begin
new_data_arr = dblarr(new_energy_dim)
new_data_arr[n_elements(add_energy):new_energy_dim[0]-1,*] = (*dist_data[i])[j].(k)[0:str_energy_dim[0]-1,*]
new_data_arr[0:n_elements(add_energy)-1,*] = (dblarr(n_elements(add_energy))+1)#(*dist_data[i])[j].(k)[0,*]
out_data[j].(k) = new_data_arr
endif else begin
out_data[j].(k) = (*dist_data[i])[j].(k)
endelse
endif
endfor
source_data = (*dist_data[i])[j].data[0:str_energy_dim[0]-1,*]
source_bins = (*dist_data[i])[j].bins[0:str_energy_dim[0]-1,*]
for k = 0,str_energy_dim[1]-1 do begin
source_dim = dimen(source_data)
if ~keyword_set(lsquares) then begin
if ~keyword_set(lin_counts) then begin
;throw out disabled bins and zero bins during log extrapolation
idx = where(source_bins[*,k] and source_data[*,k] gt 0,c)
if c eq 0 then continue
out_data[j].data[*,k] = exp(interpol(alog(source_data[idx,k]),energy_source[idx,k],energy_target,lsquadratic=lsquadratic,quadratic=quadratic,spline=spline))
endif else begin
;throw out disabled bins during extrapolation
idx = where(source_bins[*,k],c)
if c eq 0 then continue
out_data[j].data[*,k] = interpol(source_data[idx,k],energy_source[idx,k],energy_target,lsquadratic=lsquadratic,quadratic=quadratic,spline=spline)
endelse
endif else begin
fit_use_bins=dindgen(lsquares < source_dim[0])
if ~keyword_set(lin_counts) then begin
idx = where(source_bins[fit_use_bins,k] and source_data[fit_use_bins,k] gt 0,c)
if c eq 0 then continue
;note that status keyword needs to be set, even if unused, so that singular matrixes will not cause program halt
fit = poly_fit(energy_source[fit_use_bins[idx],k],alog(source_data[fit_use_bins[idx],k]),1,/double,status=s)
if s ne 0 then continue
out_data[j].data[dindgen(n_elements(add_energy)),k] = exp(fit[0]+fit[1]*energy_target[dindgen(n_elements(add_energy))])
endif else begin
idx = where(source_bins[fit_use_bins,k],c)
if c eq 0 then continue
;note that status keyword needs to be set, even if unused, so that singular matrixes will not cause program halt
fit = poly_fit(energy_source[fit_use_bins[idx],k],source_data[fit_use_bins[idx],k],1,/double,status=s)
if s ne 0 then continue
out_data[j].data[dindgen(n_elements(add_energy)),k] = fit[0]+fit[1]*energy_target[dindgen(n_elements(add_energy))]
endelse
endelse
out_data[j].data[dindgen(source_dim[0])+n_elements(add_energy),k] = source_data[*,k]
endfor
endfor
*dist_data[i] = out_data
endfor
end