;+ ;PROCEDURE: ; thm_part_energy_interpolate ; ;PURPOSE: ; Interpolate particle data by energy between sst & esa distributions using ; a weighted curve fitting routine. ; ;INPUTS: ; dist_sst: SST particle distribution structure in flux ; dist_esa: ESA particle distribution structure in flux ; energies: The set of target energies to interpolated the SST to. ; ;OUTPUTS: ; Replaces dist_sst with the interpolated data set ; ;KEYWORDS: ; error: Set to 1 on error, zero otherwise ; get_error: if set, interpolates scaling factor needed for error propagation ; ;NOTES: ; #1 The number of time samples and the times of those samples must ; be the same for dist_sst & dist_esa (use thm_part_time_interpolate.pro) ; #2 The number of angles and the angles of each sample must be ; the same for dist_sst & dist_esa (use thm_part_sphere_interp.pro) ; ;SEE ALSO: ; thm_part_dist_array ; thm_part_smooth ; thm_part_subtract, ; thm_part_omni_convert ; thm_part_time_interpolate.pro ; thm_part_sphere_interp.pro ; ; $LastChangedBy: pcruce $ ; $LastChangedDate: 2013-09-26 16:32:05 -0700 (Thu, 26 Sep 2013) $ ; $LastChangedRevision: 13156 $ ; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/thmsoc/trunk/idl/themis/spacecraft/particles/thm_part_energy_interpolate.pro $ ;- pro thm_part_energy_interp,dist_sst,dist_esa,energies,error=error,extrapolate_esa=extrapolate_esa,get_error=get_error;,dist_sst_counts=dist_sst_counts,dist_esa_counts=dist_esa_counts,emin=emin compile_opt idl2 error=1 dist_esa_i = 0 dist_esa_j = 0 min_flux = 1e-4 ; order of magnitude of min 1 count flux for esa/sst. Used so that we can to log/log interpolation on data with a lot of zeros ;TBD units must be FLUX ;TBD input validation ;TBD verification that number of samples in dist_sst matches the number in dist_esa... ;TBD verification that the number of angles in dist_sst match the number in dist_esa ;TBD verification of units on input distributions blankarr = (fltarr(n_elements(energies))+1) for i = 0,n_elements(dist_sst)-1 do begin extrap_num = 0 ;keep track of how many angles were extrapolated ;Note that most of the calculations below assume that variables ;are not changing along one dimension or another. If that assumption ;fails, these calculations will need to be made more complex. sst_dim = dimen((*dist_sst[i])[0].data) sst_template_out = (*dist_sst[i])[0] sst_energy_out = energies # (fltarr(sst_dim[1])+1) sst_denergy_out = blankarr # (fltarr(sst_dim[1])) sst_data_out = blankarr # (fltarr(sst_dim[1])) sst_scaling_out = blankarr # (fltarr(sst_dim[1])) sst_bins_out = blankarr # (*dist_sst[i])[0].bins[0,*] sst_phi_out = blankarr # sst_template_out.phi[0,*] sst_theta_out = blankarr # sst_template_out.theta[0,*] sst_dphi_out = blankarr # sst_template_out.dphi[0,*] sst_dtheta_out = blankarr # sst_template_out.dtheta[0,*] ;update all the supplemental variables that are required by the moment routines str_element,sst_template_out,'energy',sst_energy_out,/add_replace str_element,sst_template_out,'denergy',sst_denergy_out,/add_replace str_element,sst_template_out,'data',sst_data_out,/add_replace str_element,sst_template_out,'scaling',sst_scaling_out,/add_replace ;jmm, 2017-09-28 str_element,sst_template_out,'bins',sst_bins_out,/add_replace str_element,sst_template_out,'phi',sst_phi_out,/add_replace str_element,sst_template_out,'theta',sst_theta_out,/add_replace str_element,sst_template_out,'dphi',sst_dphi_out,/add_replace str_element,sst_template_out,'dtheta',sst_dtheta_out,/add_replace ;expand to match number of samples sst_mode_out = replicate(sst_template_out,n_elements(*dist_sst[i])) ;look over samples for this combined mode for j = 0,n_elements(*dist_sst[i])-1 do begin sample_sst = (*dist_sst[i])[j] sample_esa = (*dist_esa[dist_esa_i])[dist_esa_j] ;combined, pre-interpolated data combined_energy = [sample_esa.energy,sample_sst.energy] combined_data = [sample_esa.data,sample_sst.data] combined_scaling = [sample_esa.scaling,sample_sst.scaling] combined_bins = [sample_esa.bins,sample_sst.bins] combined_dim = dimen(combined_energy) ;Calculate energy bin widths for target energies a la moments_3d ; -uses 1/2 the separation of each 2 energy bin values ; -endpoints use the de/e from adjacent bin scaled by their energy combined_tmp = [sample_esa.energy[*,0],energies] esa_dim = dimen(sample_esa.data) tmp_dim = dimen(combined_tmp) combined_denergy = (shift(combined_tmp,-1)-shift(combined_tmp,1))/2. / combined_tmp combined_denergy[0] = combined_denergy[1] combined_denergy[tmp_dim[0]-1] = combined_denergy[tmp_dim[0]-2] combined_denergy *= combined_tmp ;combined_denergy = deriv([sample_esa.energy[*,0],energies]) sst_mode_out[j].denergy = (combined_denergy[esa_dim[0]:tmp_dim[0]-1]) # replicate(1.,sst_dim[1]) if max(sample_esa.energy,/nan) gt min(energies,/nan) then begin dprint,dlevel=1,'ERROR: ESA maximum energy(' + strtrim(max(sample_esa.energy,/nan),2) + ' eV) greater than minimum energy target(' + strtrim(min(energies,/nan),2) + ' eV)' return endif sst_mode_out[j].time=sample_sst.time ;copy time over sst_mode_out[j].end_time=sample_sst.end_time ;copy time over ;loop over look directions and interpolate for l=0,sst_dim[1]-1 do begin ;generate bins data for new bins(not needed...I think?) ;sst_mode_out[j].bins[*,l] = round(interpol(sample_sst.bins[*,l],sample_sst.energy[*,l],energies)) > 0 < 1 sst_idx = where(sample_sst.bins[*,l],c) if c eq 0 then begin ;extrapolate from ESA data if requested and no valid SST data exists ;set highest SST energy bin to zero to ensure interpolation doesn't go wild if keyword_set(extrapolate_esa) then begin extrap_num++ combined_bins[combined_dim[0]-1l,l] = 1 ;enable highest energy combined_data[combined_dim[0]-1l,l] = 0 ;set to zero ; combined_energy[-1,l] = energies[-1] sst_mode_out[j].bins[*,l] = 1 endif else begin combined_bins[combined_dim[0]-n_elements(energies):combined_dim[0]-1l,l] = 1 ;enable bins combined_data[combined_dim[0]-n_elements(energies):combined_dim[0]-1l,l] = !VALUES.D_NAN ;set to zero ;dprint,dlevel=1,'ERROR: No SST bins enabled for angle:'+strtrim(l,2) ;return endelse endif ;need to use proper bins so that disabled bins aren't included in interpolation calculations combined_idx = where(combined_bins[*,l],c) if c eq 0 then begin dprint,dlevel=1,'ERROR: No bins enabled for angle:'+strtrim(l,2) return endif ;The +min_flux -min_flux, turns alog(0) to alog(min_flux) preventing lots of -infinities in our interpolation, should work for scaling too sst_mode_out[j].data[*,l] = exp(interpol(alog(combined_data[combined_idx,l]+min_flux),alog(combined_energy[combined_idx,l]),alog(energies)))-min_flux if keyword_set(get_error) then sst_mode_out[j].scaling[*,l] = exp(interpol(alog(combined_scaling[combined_idx,l]+min_flux),alog(combined_energy[combined_idx,l]),alog(energies)))-min_flux ; sst_mode_out[j].data[*,l] = interpol(combined_data[combined_idx,l],combined_energy[combined_idx,l],energies) endfor ;dist_two should have matching time samples, but not necessarily ;matching mode transitions, the index logic below synchronizes ;iterations over the two data structures dist_esa_j++ if n_elements(*dist_esa[dist_esa_i]) eq dist_esa_j then begin dist_esa_i++ dist_esa_j=0 endif endfor ;notify user of how many angles were extrapolated if keyword_set(extrapolate_esa) then begin dprint, dlevel=2, 'Extrapolated '+strtrim(extrap_num,2)+' angles over '+strtrim(j,2)+' samples in mode' endif ;temporary routine bombs on some machines if out_dist is undefined, but not others if ~undefined(dist_out) then begin dist_out=array_concat(ptr_new(sst_mode_out,/no_copy),temporary(dist_out)) endif else begin dist_out=array_concat(ptr_new(sst_mode_out,/no_copy),dist_out) endelse endfor dist_sst=temporary(dist_out) heap_gc error=0 end