# Re: Convolution of Stick Spectra

```craigmnet@cow.physics.wisc.edu wrote:

>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
>[snip..]
>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.
>[snip..]
>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.

Thanks!
Todd

--

pro test

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.

time=systime(1)
for i=0L, n_elements( stick ) / 2 -1 do \$
convoluted[1,*] = stick[1,i]*exp(-((convoluted[0,*] - stick[0,i])^2)/ \$
(((.12*sqrt(stick[0,i]/1000))/1.6651)*10)^2) \$
+ convoluted[1,*]

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[0] ne -1 ) then \$
convoluted2[1,range] = stick[1,i]*exp(-((convoluted2[0,range] -  \$
stick[0,i])^2)/sigma^2) \$
+ convoluted2[1,range]
endfor

print, 'time old: ', time2-time
print, 'time new: ', systime(1)-time2

new = convoluted
new[1,*] = new[1,*] - convoluted2[1,*]

plot,new[0,*], new[1,*],yrange=[-10,10]

end ;; test

```