Source code for shapelets_real_lens.pro:

FUNCTION shapelets_real_lens, decomp, A, D, PIXSCALE=pixscale

;
;  Given a decomp (and the information it contains regarding beta,
;  n_pixels, centroid, coefficients etc.), and the A & D shear and flexion
;  transformation matrices, outputs an image of size decomp.n_pixels
;  as specified by the decomp containing a lensed image.
;
;  The lensed image is made by sampling the source plane (shapelet
;  model) on an irregular grid as specified by the lens equation,
;  corresonding to a regular grid in the image plane.
;
;  Oversampling (to do more accurate sampling of the image for more
;  accurate cross-pixel integration) should be managed by
;  appropriately setting decomp.beta and decomp.n_pixels.
;
;  The PIXSCALE keyword must be set, as it puts the dimensionful D
;  into pixels
;
;  A & D should be set as follows:
;k  = "kappa"
;s1 = shear[0]
;s2 = shear[1]
;A = [[1. - k - s1, -s2], [-s2, 1. - k + s1]]

;f1 = fflex[0]
;f2 = fflex[1]
;g1 = gflex[0]
;g2 = gflex[1]
;D = fltarr(2, 2, 2)
;D[0, 0, 0] = -.5 * (3. * f1 + g1)
;D[0, 0, 1] = -.5 * (f2 + g2)
;D[0, 1, 0] = -.5 * (f2 + g2)
;D[0, 1, 1] = -.5 * (f1 - g1)
;D[1, 0, 0] = -.5 * (f2 + g2)
;D[1, 0, 1] = -.5 * (f1 - g1)
;D[1, 1, 0] = -.5 * (f1 - g1)
;D[1, 1, 1] = -.5 * (3. * f2 - g2)

if not keyword_set(pixscale) then message, "PIXSCALE keyword must be set, flexion is a dimensional quantity"


Ddim = 0.5 * D * pixscale ; put the D_ijk matrix into the correct 
                          ; dimensions (i.e. pixel^-1)

; Make an array of x and y positions at the location of each pixel in
; the output image plane...
shapelets_make_xarr, [decomp.n_pixels[0], decomp.n_pixels[1]], $
                     xlens, ylens, X0=decomp.x

; Map these into a non-regular grid in the source plane, using the weak
; lensing terms in the lens equation
;
xsource = A[0, 0]*xlens + A[0, 1]*ylens + Ddim[0, 0, 0]*xlens*xlens $
                                        + Ddim[0, 0, 1]*xlens*ylens $
                                        + Ddim[0, 1, 0]*ylens*xlens $
                                        + Ddim[0, 1, 1]*ylens*ylens
ysource = A[1, 0]*xlens + A[1, 1]*ylens + Ddim[1, 0, 0]*xlens*xlens $
                                        + Ddim[1, 0, 1]*xlens*ylens $
                                        + Ddim[1, 1, 0]*ylens*xlens $
                                        + Ddim[1, 1, 1]*ylens*ylens

; Then get the values of the surface brightness at these locations
image = simage_shapelets_surface_brightness(xsource, ysource, decomp)
return, image

END

Return to the shapelets web page or the code help menu.