pro mk_img, _extra=_extra

if ~exist(xsiz) then xsiz = 768d ; in pixels
if ~exist(ysiz) then ysiz = 512d ; in pixels
aspect_ratio = ysiz/xsiz
if ~exist(bigdot_xsiz) then bigdot_xsiz = 5d
if ~exist(bigdot_ysiz) then bigdot_ysiz = 5d
if ~exist(lildot_xsiz) then lildot_xsiz = 3d
if ~exist(lildot_ysiz) then lildot_ysiz = 3d

if ~exist(n_xdots) then n_xdots = 8d+1d
if ~exist(n_ydots) then n_ydots = n_xdots*aspect_ratio
xdelt_grid = xsiz/(n_xdots-1)
ydelt_grid = ysiz/(n_ydots-1)

xgrid = ((lindgen(n_xdots)-((n_xdots-1)/2))*xdelt_grid) + 0.5d
ygrid = ((lindgen(n_ydots)-((n_ydots-1)/2))*ydelt_grid) + 0.5d

if ~exist(interp) then interp = 0
if ~exist(pivot) then pivot = 0
if ~exist(mag) then mag = 1.25
if ~exist(missing) then missing = 0

if ~exist(wavelnth) then wavelnth = 171
if ~exist(crota2) then crota2 = 30d ; deg
if ~exist(cdelt1) then cdelt1 = 0.5d ; arcsec/pixel
if ~exist(cdelt2) then cdelt2 = 0.5d ; arcsec/pixel
naxis1 = xsiz
naxis2 = ysiz

img = dblarr(xsiz, ysiz)

xcen_pix = xsiz/2d + 0.5d
ycen_pix = ysiz/2d + 0.5d

; Create test pattern in image:
for i=0, n_xdots-1 do $
   for j=0, n_ydots-1 do $
      img[(xgrid[i]-(lildot_xsiz-1)/2):(xgrid[i]+(lildot_xsiz-1)/2), $
          (ygrid[j]-(lildot_ysiz-1)/2):(ygrid[j]+(lildot_ysiz-1)/2)] = 1d

img[(xcen_pix-(bigdot_xsiz-1)/2):(xcen_pix+(bigdot_xsiz-1)/2), $
    (ycen_pix-(bigdot_ysiz-1)/2):(ycen_pix+(bigdot_ysiz-1)/2)] = 2d

; Define remaining necessary keywords:

;naxis1 - defined
;naxis2 - defined

;cdelt1 - defined
;cdelt2 - defined

if ~exist(xcen) then xcen = (1024d)*0.6d ; AIA pixel coord of fov center
                                         ; converted to arcsec
if ~exist(ycen) then ycen = (1280d)*0.6d ; AIA pixel coord of fov center
                                         ; converted to arcsec

crpix1 = xcen/cdelt1 + 0.5d
crpix2 = ycen/cdelt2 + 0.5d
crval1 = 0.0d
crval2 = 0.0d

; Optionally embed sub image in final image array:

if keyword_set(do_embed) then begin

   if ~exist(crota2_embed) then crota2_embed = 10d ; degrees
   delta_angle = crota2_embed - crota2            ; TODO: Check this for sign error

   if not exist(cdelt11_embed) then cdelt1_embed = .7d ; arcsec/pixel
   if not exist(cdelt12_embed) then cdelt2_embed = .7d ; arcsec/pixel

   mag = cdelt1/cdelt1_embed

; Transform CRPIX1/2 from subimg frame to embedding frame:

; First rotate the coordinates to align the axes:
   sin0 = sin(delta_angle/!radeg)
   cos0 = cos(delta_angle/!radeg)
   crpix1_embed = crpix1*cos0 - crpix2*sin0
   crpix2_embed = crpix1*sin0 + crpix2*cos0
; Now translate coordinates to common origin:
   crpix1_embed = crpix1_embed + (crpix1_embed - 

   
; Remap sub image into embedding array:
rotimg = rot(img, delta_angle, mag, x0, y0, $
             interp=interp, cubic=cubic, missing=missing, pivot=0, $
             _extra=_extra)















subimg = rotimg[((xsiz_img/2+1)-xsiz/2):((xsiz_img/2+1)+xsiz/2-1), $
                ((ysiz_img/2+1)-ysiz/2):((ysiz_img/2+1)+ysiz/2-1)]

wdef, 0, 1024, 1024 & tvscl, rebin(img, 1024, 1024, /samp)
STOP, 'STOPPING after insertion of test pattern into original image.'

; Transform and extract desired sub-image:
rotimg = rot(img, angle, mag, x0, y0, $
             interp=interp, cubic=cubic, missing=missing, pivot=0, $
             _extra=_extra)
subimg = rotimg[((xsiz_img/2+1)-xsiz/2):((xsiz_img/2+1)+xsiz/2-1), $
                ((ysiz_img/2+1)-ysiz/2):((ysiz_img/2+1)+ysiz/2-1)]
; NB: [xsiz, ysiz] are actually [xsiz_subimg, ysiz_subimg]

wdef, 0, xsiz, ysiz & tvscl, subimg
STOP, 'STOPPING after extraction and display of sub-image.'
; Now update WCS keywords for sub-image:
subhdr = hdr
subhdr.crval1 = 0d
subhdr.crval2 = 0d
subhdr.crpix1 = x0
subhdr.crpix2 = y0
subhdr.crota2 = angle

subhdr.cdelt1 = hdr.cdelt1/mag
subhdr.cdelt2 = hdr.cdelt1/mag

subhdr.naxis = 2l
subhdr.naxis1 = xsiz
subhdr.naxis2 = ysiz

; Note - I should really be using Bill Thompson's WCS routine suite here

index2map, subhdr, subimg, submap

wdef, 0, 1024
erase & wait, 0.5
plot_map, submap, /limb, grid=15 ; , center=[0,0], fov=[40, 40]

help, index0, data0
STOP, 'STOPPING after creation and display of extracted raw sub-map.'

aia_prep, index0, data0, oindex0, odata0, /use_ref, /cutout, index_ref=subhdr, /no_pnt_update

index2map, oindex0, odata0, omap0

erase & wait, 0.5
plot_map, omap0, /limb, grid=15, ref_map=submap

STOP, "STOPPING after creation and display of extracted aia_prep'd sub-map."

STOP, 'Stopping before return.'

end