;+ ;PROCEDURE: thm_esa_energy_extrapolate ;PURPOSE: Performs linear extrapolation of esa data to a new set of higher energies ; ;INPUTS: ; dist_data: ; The esa 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='peef',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=[100000,70000,50000] ; ;NOTES: ; Removes the retrace bin from ESA modes that include retrace. ; ; $LastChangedBy: pcruce $ ; $LastChangedDate: 2013-01-10 18:10:18 -0800 (Thu, 10 Jan 2013) $ ; $LastChangedRevision: 11424 $ ; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_5_0/projects/themis/spacecraft/particles/ESA/thm_esa_energy_extrapolate.pro $ ;- pro thm_esa_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) if n_elements(str_energy_dim) eq 1 then begin str_energy_dim = [str_energy_dim,1] endif ;use number of energies to detect retrace bin ;the variable below is used as the index to the first bin that will be included in the output if str_energy_dim[0] ne 32 then begin retrace=0 endif else begin retrace=1 endelse ;create source target, remove retrace energy_target = str_template.energy[retrace:str_energy_dim[0]-1,0] energy_source = str_template.energy[retrace:str_energy_dim[0]-1,*] energy_target = [reverse(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 ;if n_elements(str_energy_dim) eq 2 then begin 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[retrace:str_energy_dim[0]-1,0]]#(dblarr(str_energy_dim[1])+1) ;very quick and very dirty delta-energies ;endif else if n_elements(str_energy_dim) eq 1 then begin ; new_energy_arr = energy_target ; new_denergy_arr = [(abs(deriv(energy_target)))[0:n_elements(add_energy)-1],str_template.denergy[retrace:str_energy_dim[0]-1,0]] ;endif else begin ; dprint,dlevel=1,"ERROR: unexpected energy dimensions" ; return ;endelse 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) || (str_energy_dim[1] eq 1 && str_energy_dim[0] eq tag_dimen[0] && n_elements(tag_dimen) eq 1) 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)[retrace: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)[retrace,*] 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[retrace:str_energy_dim[0]-1,*] source_bins = (*dist_data[i])[j].bins[retrace:str_energy_dim[0]-1,*] for k = 0,str_energy_dim[1]-1 do begin source_dim = dimen(source_data) out_data[j].data[dindgen(source_dim[0])+n_elements(add_energy),k] = source_data[*,k] 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 endfor endfor *dist_data[i] = out_data endfor end