;MSW 2/11/07
;routine to convert IFU cube into RSS format and output corresponding position file

PRO cube2rss, infile, outfile, pf
if (n_params() ne 3) then $
  message,'Usage: cube2rss, "infile.fits", "outfile.fits", position_file.fits"'

; read in the fits file selected and extract some dimension info
fits_read,infile,z,hdr
dim_size = size(z,/dimensions)
dim_size0 = dim_size[0]
dim_size1 = dim_size[1]
dim_size2 = dim_size[2]

print,'dimensions   no of x values'
print,strtrim(string(dim_size0),2),'x',strtrim(string(dim_size1),2), dim_size2

; read in lambda_start and delta_lambda from FITS header 
crval1=sxpar(hdr,'crval1')
cdelt1=sxpar(hdr,'cdelt1')
crval2=sxpar(hdr,'crval2')
cdelt2=sxpar(hdr,'cdelt2')
crval3=sxpar(hdr,'crval3')
cdelt3=sxpar(hdr,'cdelt3')
; check whether crval1 and cdelt1 are valid header params
if (n_elements(crval1) eq 0) then $
        message,'crval1 not defined in FITS header'
if (n_elements(cdelt1) eq 0) then $
        message,'cdelt1 not defined in FITS header'
if (n_elements(crval2) eq 0) then $
        message,'crval2 not defined in FITS header'
if (n_elements(cdelt2) eq 0) then $
        message,'cdelt2 not defined in FITS header'
if (n_elements(crval3) eq 0) then $
        message,'crval3 not defined in FITS header'
if (n_elements(cdelt3) eq 0) then $
        message,'cdelt3 not defined in FITS header'
print,'      crval3       cdelt3'
print,crval3,cdelt3

; generate wavelength grid based on header keywords
wave=findgen(dim_size2)*cdelt3+crval3

;generate x and y values and write position table
delt=0.2
y=findgen(dim_size0)*delt;+crval1
x=findgen(dim_size1)*delt;+crval2
xrev=reverse(x)   ;reverse the x axis
k=1
openw,1,pf
printf,1,'#C bla R 0.2 -0.2 0 test'
printf,1,'#I 0.2 0.2'
printf,1,'#P S'
for i=0,n_elements(x)-1 do begin
  for j=0,n_elements(y)-1 do begin
    printf,1,strtrim(string(k),2),' ',strtrim(string(xrev[i]),2),' ',strtrim(string(y[j]),2),' 1'
    k=k+1
  endfor
endfor
close,1

;reverse x axis of array as well
zrev=reverse(z,1)
;reform array into 2D
zrev=reform(zrev,dim_size0*dim_size1,dim_size2,/overwrite)
ztra=transpose(zrev)
zdims = size(ztra)
dim_size = size(ztra,/dimensions)
dim_size0 = dim_size[0]
dim_size1 = dim_size[1]
print,'reformed array'
print,'        x dimension  no of spax'
print,dim_size0,dim_size1


;do some editing of the header info
sxaddpar,hdr,'naxis',zdims[0]
sxaddpar,hdr,'naxis1',dim_size1
sxaddpar,hdr,'naxis2',dim_size0
sxdelpar,hdr,'naxis3'
sxaddpar,hdr,'crval2','0.'
sxaddpar,hdr,'cdelt2','1.'
crval1=sxpar(hdr,'crval3')
cdelt1=sxpar(hdr,'cdelt3')
sxdelpar,hdr,'crval3'
sxdelpar,hdr,'cdelt3'
sxdelpar,hdr,'cd3_3'
sxaddpar,hdr,'crval1',crval1
sxaddpar,hdr,'cdelt1',cdelt1
sxaddpar,hdr,'cd1_1',cdelt1
sxaddpar,hdr,'dispaxis','1'

;write out to fits file
fits_write,outfile,ztra,hdr


END