Commit 4785bd61cb63b100f2b4a555037fc6fe2cba50ea

Authored by knodlseder
1 parent 40a2c333

initial version

Showing 1 changed file with 64 additions and 0 deletions
test/wcs_test.pro 0 → 100644
  1 +;==============================================================================
  2 +; Script: wcs_test
  3 +; ----------------------------------------------------------------------------
  4 +; Create WCS test map
  5 +; ----------------------------------------------------------------------------
  6 +; Parameters:
  7 +; ra - Right Ascension of map centre (deg)
  8 +; dec - Declination of map centre (deg)
  9 +; pixsize - Pixelsize for FITS map
  10 +; proj - Map projection
  11 +; filename - Filename
  12 +;==============================================================================
  13 +pro wcs_test, ra=ra, dec=dec, pixsize=pixsize, proj=proj, filename=filename
  14 +
  15 +
  16 +;====================
  17 +; Default parameters
  18 +;=====================
  19 +if n_elements(ra) eq 0 then ra = 83.63
  20 +if n_elements(dec) eq 0 then dec = 22.0
  21 +if n_elements(pixsize) eq 0 then pixsize = 0.02
  22 +if n_elements(size) eq 0 then size = 200
  23 +if n_elements(proj) eq 0 then proj = 'CAR'
  24 +if n_elements(filename) eq 0 then filename = 'map.fits'
  25 +
  26 +
  27 +;==============================
  28 +; Build astronometry structure
  29 +;==============================
  30 +galactic = 0
  31 +equinox = 2000.0
  32 +cdelt = [pixsize, pixsize]
  33 +crpix = [double(size)/2.0+0.5, double(size)/2.0+0.5]
  34 +crval = [ra, dec]
  35 +ctype = ['RA---'+proj,'DEC--'+proj]
  36 +make_astr, astr, delt=cdelt, crpix=crpix, crval=crval, ctype=ctype
  37 +
  38 +
  39 +;=========
  40 +; Set map
  41 +;=========
  42 +map = dblarr(size, size)
  43 +for ix=0,size-1 do begin
  44 + x = replicate(double(ix),size)
  45 + y = dindgen(size)
  46 + xy2ad, x, y, astr, a, d
  47 + gcirc, 0, 83.63*!dtor, 22.0*!dtor, a*!dtor, d*!dtor, distance
  48 + distance = distance/!dtor
  49 + inx = where(distance le 0.1, nvalues)
  50 + if (nvalues gt 0) then map[ix,inx] = 1.0
  51 +endfor
  52 +
  53 +
  54 +;==========
  55 +; Save map
  56 +;==========
  57 +mkhdr, hdr, map
  58 +putast, hdr, astr ;, cd_type=1
  59 +;sxdelpar, hdr, 'LONPOLE'
  60 +;sxdelpar, hdr, 'LATPOLE'
  61 +;sxdelpar, hdr, 'PV2_1'
  62 +mwrfits, map, filename, hdr, /create
  63 +
  64 +end
... ...