dimensional headache [message #43322] |
Fri, 01 April 2005 08:49 |
Margrethe
Messages: 4 Registered: March 2005
|
Junior Member |
|
|
I'm developing a two-dimensional model and to save computational time
I've tried to avoid for-loops whereever I could and used matrices
instead. I thought I had it all figured out, but I'm doing something
wrong...I'd be grateful if someone could tell me what's wrong with the
following:
The model lies on a grid xgrid=findgen(160), ygrid=findgen(160). I want
to calculate the numerical value of the model at the coordinate
xmod,ymod
Here's my old code using for loops
for i=0,159 do begin
for j=0,159 do begin
dx_sqr = ( xmod - xgrid (i) )^2.
dy_sqr = ( ymod - ygrid (j) )^2.
d = sqrt ( dx_sqr + dy_sqr )
; don't take the wings of the psf into account
indx = where (d lt (3. * (fwhm/2.35) ) ,ct)
if (ct gt 0) then begin
vel_p = ( vel_sys + v_rad(i,j) )
psf = gaussfunc(st_dev,xgrid(i),ygrid(j),xmod,ymod )
lineprof(i,j) = lineprof(i,j) + $
total ( psf * image(i,j) * $
exp( (vel - vel_p(i,j))^2./(-2.*v2)) )
endif
endfor
endfor
The above example works fine, but of course takes a long time to
compute.
Here's my attempt with matrices, which doesn't give the correct output.
dx_sqr = (xmod - findgen(160))^2.
dy_sqr = (ymod - findgen(160))^2.
; expand one dimension
dx_sqr_expand = rebin(dx_sqr,160,160)
dy_sqr_expand = transpose(rebin(dy_sqr,160,160))
d = sqrt(dx_sqr_expand + dy_sqr_expand)
; don't take the wings of the psf into account
WhereToMulti,d,where(d lt (3.*(fwhm/2.35)),ct),colindx,rowindx
vel_p = (vel_sys + v_rad(colindx,rowindx))
psf = gaussfunc(st_dev,xgrid(colindx),ygrid(rowindx),xmod,ymod)
lineprof(k,l) = lineprof(k,l) + $
total( psf * image(colindx,rowindx) * $
exp( (vel - vel_p(colindx,rowindx))^2./(-2.*v2)) )
|
|
|