[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Convolution of Stick Spectra

mole6e23@hotmail.com (Todd Clements) writes:

> 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

You're welcome.  But call me a pedant too.  You can save yet more
computations by precomputing the gaussian argument.


xx = convoluted(0,*)
yy = convoluted(1,*)
for i = 0L, n_elements(stick)/2 -1 do begin
  z1 = (xx-stick(0,i))/(((.12*sqrt(stick[0,i]/1000))/1.6651)*10)^2
  wh = where(z1 LT (2.*4^2), ct)  ;; Four sigma (remember the 1/2!)
  if ct GT 0 then yy(wh) = yy(wh) + stick(1,i)*exp(-z1)
convoluted(1,*) = yy

I hope I got it all right, but you get the idea.  By the way, as I
note above, the definition of the gaussian is exp(-x^2/(2*sigma^2)).
Is the 1/2 buried somewhere in your definition?

Craig B. Markwardt, Ph.D.         EMAIL:    craigmnet@cow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response