From 4785bd61cb63b100f2b4a555037fc6fe2cba50ea Mon Sep 17 00:00:00 2001
From: knodlseder <>
Date: Tue, 28 Jun 2011 03:26:44 +0000
Subject: [PATCH] initial version

---
 test/wcs_test.pro | 64 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 64 insertions(+), 0 deletions(-)
 create mode 100644 test/wcs_test.pro

diff --git a/test/wcs_test.pro b/test/wcs_test.pro
new file mode 100644
index 0000000..88a4e95
--- /dev/null
+++ b/test/wcs_test.pro
@@ -0,0 +1,64 @@
+;==============================================================================
+; Script: wcs_test
+; ----------------------------------------------------------------------------
+; Create WCS test map
+; ----------------------------------------------------------------------------
+; Parameters:
+;  ra       - Right Ascension of map centre (deg)
+;  dec      - Declination of map centre (deg)
+;  pixsize  - Pixelsize for FITS map
+;  proj     - Map projection
+;  filename - Filename
+;==============================================================================
+pro wcs_test, ra=ra, dec=dec, pixsize=pixsize, proj=proj, filename=filename
+
+
+;====================
+; Default parameters
+;=====================
+if n_elements(ra)       eq 0 then ra       = 83.63
+if n_elements(dec)      eq 0 then dec      = 22.0
+if n_elements(pixsize)  eq 0 then pixsize  = 0.02
+if n_elements(size)     eq 0 then size     = 200
+if n_elements(proj)     eq 0 then proj     = 'CAR'
+if n_elements(filename) eq 0 then filename = 'map.fits'
+
+
+;==============================
+; Build astronometry structure
+;==============================
+galactic = 0
+equinox  = 2000.0
+cdelt    = [pixsize, pixsize]
+crpix    = [double(size)/2.0+0.5, double(size)/2.0+0.5]
+crval    = [ra, dec]
+ctype    = ['RA---'+proj,'DEC--'+proj]
+make_astr, astr, delt=cdelt, crpix=crpix, crval=crval, ctype=ctype
+
+
+;=========
+; Set map
+;=========
+map = dblarr(size, size)
+for ix=0,size-1 do begin
+  x = replicate(double(ix),size)
+  y = dindgen(size)
+  xy2ad, x, y, astr, a, d
+  gcirc, 0, 83.63*!dtor, 22.0*!dtor, a*!dtor, d*!dtor, distance
+  distance = distance/!dtor
+  inx = where(distance le 0.1, nvalues)
+  if (nvalues gt 0) then map[ix,inx] = 1.0
+endfor
+
+
+;==========
+; Save map
+;==========
+mkhdr, hdr, map
+putast, hdr, astr ;, cd_type=1
+;sxdelpar, hdr, 'LONPOLE'
+;sxdelpar, hdr, 'LATPOLE'
+;sxdelpar, hdr, 'PV2_1'
+mwrfits, map, filename, hdr, /create
+
+end
--
libgit2 0.21.4