[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Convolution of Stick Spectra
>After looking at your problem it looks like the widths of the lines
>are as wide as your energy range. I guess this would be appropriate
>so many exponentiations. You could have gotten a pretty big savings
>if the lines were narrow, and you could restrict the computation to a
>narrow region around the line center (say +/- 10 sigma with WHERE). I
>recommend that anyway to get rid of underflow errors.
>If for example your sigma term were
>(((.12*sqrt(stick[0,i]/1000))/1.6651)*10) [ note the last factor is
>smaller ] then things start to look interesting.
That's because there's always a danger in taking code x and modifying it
into simpler example y! As you suggested, I actually DO have a sigma term
with a *10 instead of a *1000.
I did what you suggested with the +/- 10*sigma with where (code below),
and it drastically improved running time. For the 5000 element array, the
running time went from 62.3 seconds to 4.9 seconds! I plot the error
associated with the method, and I was actually able to go down to +/-
4*sigma before there was any noticable error.
convoluted = fltarr( 2, 2048 )
convoluted[0,*] = findgen( 2048 ) / 2047. * 4.
convoluted2 = convoluted
;; Let's fake a stick spectrum, we usually have at least this many elements
stick = abs(randomn(10, 2, 5000 ))
stick[0,*] = stick[0,*] * 4.
stick[1,*] = stick[1,*] * 1000.
for i=0L, n_elements( stick ) / 2 -1 do $
convoluted[1,*] = stick[1,i]*exp(-((convoluted[0,*] - stick[0,i])^2)/ $
time2 = systime(1)
for i=0L, n_elements( stick ) / 2 -1 do begin
;; Find range
sigma = (((.12*sqrt(stick[0,i]/1000))/1.6651)*10)
range = where( convoluted2[0,*] ge (stick[0,i]-4*sigma) and $
convoluted2[0,*] le (stick[0,i]+4*sigma) )
if( range ne -1 ) then $
convoluted2[1,range] = stick[1,i]*exp(-((convoluted2[0,range] - $
print, 'time old: ', time2-time
print, 'time new: ', systime(1)-time2
new = convoluted
new[1,*] = new[1,*] - convoluted2[1,*]
end ;; test