PRO GRS_IntInt, cube, outputimage, vel_low, vel_high
;;
;; RYS attempt to remake Simon's moment maker, but now on FITS cubes.
;; goal is to run FASTER!
;;
; Changed to gaussian smoothed 11/08 by LDA
; Updated with convol for speed
;
; Updated JMR : 24 March 2006
; changed input values to be 'filename of cube', 'filename of output image', and 'velocity in km/s'
; added readfits for input cube and writefits for final image created
; added calculation of correct image min and max values to be added to the image header
; fixed calculation of Integrated intensity and final coordinates.

IF n_params() LT 4 THEN BEGIN
  print,'----------------------------------------------------------------------'
  print,'PROCEDURE GRS_IntInt, input_cube, output_map, vel_low, vel_high '
  print, ''
  print,' Make an Integrated Intensity image from a GRS data cube'
  print, ''
  print,' Input velocity should be in km/s.'
  print,'----------------------------------------------------------------------'
  return
endif

cubefilename=strcompress(cube,/remove_all)
output=strcompress(outputimage,/remove_all)
vel_low_init=vel_low*1e3
vel_high_init=vel_high*1e3

; Read in the fits cube
image=readfits(string(cubefilename),header,/SILENT)

; Read in header values and find image size
dim = size(image, /Dimensions)
n_x = dim[0]
n_y = dim[1]
n_z = dim[2]

crpix_x = SXPAR(header, 'CRPIX1')
cdelt_x = SXPAR(header, 'CDELT1')
crval_x = SXPAR(header, 'CRVAL1')
naxis_x = SXPAR(header, 'NAXIS1')

crpix_y = SXPAR(header, 'CRPIX2')
cdelt_y = SXPAR(header, 'CDELT2')
crval_y = SXPAR(header, 'CRVAL2')
naxis_y = SXPAR(header, 'NAXIS2')

crpix_z = SXPAR(header, 'CRPIX3')
cdelt_z = SXPAR(header, 'CDELT3')
crval_z = SXPAR(header, 'CRVAL3')
naxis_z = SXPAR(header, 'NAXIS3')

; If longitude < 39 then velocity ranges to 135000. Otherwise to 85000
x0 = (0.-crpix_x)*cdelt_x+crval_x < (n_x-crpix_x)*cdelt_x+crval_x
IF x0 lt 39 then low=1 else low=0

; If no values are input, use all values
if N_ELEMENTS(vel_low_init) EQ 0 then vel_low = 20000
if N_ELEMENTS(vel_high_init) EQ 0 then begin
    if low EQ 1 then vel_high_init = 100000 else vel_high_init = 60000
endif

; If values input incorrectly, switch
vel_low = vel_low_init < vel_high_init
vel_high = vel_low_init > vel_high_init

; Compute the values to be used in the arrays
vel_low_arr = Round((vel_low-crval_z)/cdelt_z+crpix_z)
vel_high_arr = Round((vel_high-crval_z)/cdelt_z+crpix_z)
    
; distance_arr is just the distance from the center of a 3x3 matrix
distance_arr=[[2^.5, 1, 2^.5], [1,0,1], [2^.5, 1, 2^.5]]

; The 46 is in arcseconds. cdelt_x is in degrees. FWHM = 2.354 * sigma. variance = sigma^2
beam = 46/(ABS(cdelt_x)*3600)
variance = (beam/2.354)^2
gauss_arr = (2*!pi*variance)^(-.5)*exp(-.5*(distance_arr^2)/variance)

; Compute moment with convol
scale_factor = total(gauss_arr)
subcube = total_1d(image[*,*,vel_low_arr:vel_high_arr], 3)*(cdelt_z/1000)
moment = convol(subcube, gauss_arr, scale_factor, /nan)

;--------------------------------------------------------------------------------
; Update the header
;--------------------------------------------------------------------------------
    moment_hdr = header

    n_x_out = N_ELEMENTS(moment[*,0])
    n_y_out = N_ELEMENTS(moment[0,*])

; Update the header
    SXADDPAR, moment_hdr, 'CRPIX1', crpix_x+2
    SXADDPAR, moment_hdr, 'CRPIX2', crpix_y+2
    SXADDPAR, moment_hdr, 'NAXIS', 2
    SXADDPAR, moment_hdr, 'NAXIS1', n_x_out
    SXADDPAR, moment_hdr, 'NAXIS2', n_y_out

    SXADDPAR, moment_hdr, 'BUNIT', 'K km/s'    

    SXDELPAR, moment_hdr, 'CTYPE3'
    SXDELPAR, moment_hdr, 'NAXIS3'
    SXDELPAR, moment_hdr, 'CRPIX3'
    SXDELPAR, moment_hdr, 'CRVAL3'
    SXDELPAR, moment_hdr, 'CDELT3'
    SXDELPAR, moment_hdr, 'CROTA3'

   sxaddpar,moment_hdr,'DATAMIN',min(moment)
   sxaddpar,moment_hdr,'DATAMAX',max(moment)

writefits,output,moment,moment_hdr,/CHECKSUM

return
end