;+
; NAME:
;	CREATE_WCS_MAP
; PURPOSE:
;	Create WCS map from GALEX data format of Murthy (2014) ApJS 213, 32
; EXPLANATION:
;	This program is designed to take the data from the output of
;	get_galex_hlsp_data and to make a map based on the WCS
;	coordinates. The data variable must be a structure with at least 
;	the following elements:
;	GLON:		Galactic longitude in degrees
;	GLAT:		Galactic latitude in degrees
;	FUV_FINAL:	FUV flux in photon units
;	FUV_EXP_TIME:	FUV exposure time in seconds
;	NUV_FINAL:	NUV flux in photon units
;	NUV_EXP_TIME:	NUV exposure time in seconds
;
;	Notes:  The data are added together into bins defined by the WCS projection.
;		Each data point is weighted by its exposure time. Only non-zero
;		points are added. Other than the binning, no further smoothing
;		is done.
;	
; CALLING SEQUENCE
;	create_wcs_map, hdr, data, fuv_grid, fuv_var, nuv_grid, nuv_var, glon, glat, $
;			fov, cdelt, proj = proj, wcs = wcs
;
; Inputs:
;	DATA -	Structure as defined above. The input may be taken directly 
;		from get_gslp_data.pro.
; Optional Inputs:
;	GLON	  -	Galactic longitude of center of field in degrees.
;	GLAT	  -	Galactic latitude of center of field in degrees.
;	FOV	  -	Field of view in degrees (radius).
;	CDELT	  -	Degrees per pixel.
; Outputs:
;	HDR  -		FITS header. Note the variable is rewritten.
;	FUV_GRID  -	2-dim array with binned FUV data. 
;	FUV_VAR	  -	2-dim array with total FUV exposure time for each pixel.
;	NUV_GRID  -	2-dim array with binned NUV data. 
;	NUV_VAR	  -	2-dim array with total NUV exposure time for each pixel.
; Optional Keywords
;	PROJ	  -	WCS projection. Any projection recognized by adxy. Unless
;			this keyword is specifically defined, a TAN projection
;			is assumed.
;	WCS	  -	Structure with WCS parameters. If empty, WCS parameters
;			used are passed out. If not empty, the WCS structure is 
;			filled.
;
; Notes: Documentation and links at 
;	 http://www.iiap.res.in/personnel/murthy/Jayant_Murthy/GALEX_Data.html
;	 Points more than 0.5 degree from center are unreliable so best to use
;	 those points for which (x - 24)^2 + (y - 24)^2 < 225
;
; Procedures Used:
;	ADXY
;Revision History
; REVISION HISTORY
;	Written by Jayant Murthy (jmurthy@yahoo.com): 	May 21, 2014
;	Minor corrections and documentation:		Jul 24, 2014
; COPYRIGHT
;	Simplified BSD license copied from Wikipedia 
;	(https://en.wikipedia.org/wiki/BSD_licenses)
;	Copyright (c) 2014, Jayant Murthy
;	All rights reserved.
;	Redistribution and use in source and binary forms, with or without
;	modification, are permitted provided that the following conditions are met:
;
;	1. Redistributions of source code must retain the above copyright notice, this
;	   list of conditions and the following disclaimer. 
;	2. Redistributions in binary form must reproduce the above copyright notice,
;	   this list of conditions and the following disclaimer in the documentation
;	   and/or other materials provided with the distribution.

;	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
;	ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
;	WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
;	DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
;	ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
;	(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
;	LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
;	ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
;	(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
;	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
;
;	The views and conclusions contained in the software and documentation are those
;	of the authors and should not be interpreted as representing official policies, 
;	either expressed or implied, of the FreeBSD Project.
;-


function check_wcs, wcs, glon, glat, fov, cdelt, proj
;
;Function to check whether the WCS parameters are defined and, if not, to define
;them. The relevant parameters are held in the wcs structure defined below.

if (n_elements(wcs) eq 0)then begin
	wcs = {naxis1: 0, naxis2: 0, crval1: 0., crval2: 0., crpix1: 0. , crpix2:0., $
	cdelt1: 0., cdelt2: 0., ctype1: 'GLON', ctype2: 'GLAT'}

;All input parameters are in degrees
	if (n_elements(glon) eq 0)then $
		read, "Enter Galactic Longitude: ", 		glon
	if (n_elements(glat) eq 0)then $
		read, "Enter Galactic Latitude: ",  		glat
	if (n_elements(fov) eq 0)then $
		read, "Enter Field of View in degrees: ",	fov
	if (n_elements(cdelt) eq 0)then $
		read, "Enter spatial resolution in degrees: ",  cdelt
	if (n_elements(proj) eq 0)then proj = '-TAN'
		
	wcs.crval1 = glon
	wcs.crval2 = glat
	wcs.cdelt1 = -cdelt
	wcs.cdelt2 = cdelt
	nx = fix(fov/cdelt) + 1
	ny = fix(fov/cdelt) + 1
	wcs.naxis1 = nx
	wcs.naxis2 = ny
	wcs.crpix1 = nx/2
	wcs.crpix2 = ny/2
	wcs.ctype1 = wcs.ctype1 + proj
	wcs.ctype2 = wcs.ctype2 + proj
endif
return,1
end

pro create_wcs_map, hdr, data, fuv_grid, fuv_var, $
		    nuv_grid, nuv_var, glon, glat, fov, cdelt,  $
		    proj = proj, wcs = wcs

;Check to see if wcs structure is defined. If not, we read the relevant
;parameters from the terminal and populate the structure.    
	tst = check_wcs(wcs, glon, glat, fov, cdelt, proj)

;Define the arrays
	fuv_grid = fltarr(wcs.naxis1, wcs.naxis2)
	nuv_grid = fuv_grid
	fuv_var  = fuv_grid
	nuv_var  = fuv_grid

;Setup the FITS header. 	
	mkhdr, hdr, fuv_grid
	sxaddpar, hdr, "CRVAL1", wcs.crval1
	sxaddpar, hdr, "CRVAL2", wcs.crval2
	sxaddpar, hdr, "CDELT1", wcs.cdelt1
	sxaddpar, hdr, "CDELT2", wcs.cdelt2
	sxaddpar, hdr, "CTYPE1", wcs.ctype1
	sxaddpar, hdr, "CTYPE2", wcs.ctype2
	sxaddpar, hdr, "CRPIX1", wcs.crpix1
	sxaddpar, hdr, "CRPIX2", wcs.crpix2
	hist = "Written on "+ systime(0)
	sxaddhist, hist, hdr
	
;Convert the galactic coordinates to x and y indices in the grid
	nlines = n_elements(data)
	adxy, hdr, data.glon, data.glat, x, y
	x = fix(x)
	y = fix(y)
	
;Select those points which fit inside the grid and where the data points are
;greater than 0.
	q = where((x ge 0) and (x lt wcs.naxis1) and $
		  (y ge 0) and (y lt wcs.naxis2) and $
		  (data.fuv_final gt 0) and $
		  (data.nuv_final gt 0), nq)
;Cut down the arrays to include only the valid points.
	x = x(q)
	y = y(q)
	fuv = data(q).fuv_final
	nuv = data(q).nuv_final
	ftime = data(q).fuv_exp_time
	ntime = data(q).nuv_exp_time
	
;Start the binning
	for i = 0L, nq - 1 do begin
	
;Convert the integers into float because of the damn IDL integer issues.
		fuv_val  = float(fuv(i))
		fuv_time = float(ftime(i))
		
;If the data point is > 0, add into the array weighting by the exposure time
		if (fuv_val gt 0)then begin
			fuv_grid(x(i), y(i)) = fuv_grid(x(i), y(i)) + $
				fuv_val*fuv_time
			fuv_var(x(i), y(i)) = fuv_var(x(i), y(i)) + $
				fuv_time
		endif
		
;And now do the same for the NUV
		nuv_val  = float(nuv(i))
		nuv_time = float(ntime(i))
		if (nuv_val gt 0)then begin
			nuv_grid(x(i), y(i)) = nuv_grid(x(i), y(i)) + $
				nuv_val*nuv_time
			nuv_var(x(i), y(i)) = nuv_var(x(i), y(i)) + $
				nuv_time
		endif
endfor

;Divide by the exposure time to convert back into photon units.
	q=where(fuv_var gt 0)
	fuv_grid(q)=fuv_grid(q)/fuv_var(q)
	q=where(nuv_var gt 0)
	nuv_grid(q)=nuv_grid(q)/nuv_var(q)
	
			
end