;Procedure to construct an image cube of depth nslices=1 from an input
;2D fits image, trim as necessary, and alter its header to include WCS
;coordinate system information.
;
; inputs: 
; inname - name of input FITS file from model (string)
; xmin,xmax,ymin,ymax (integers) minimum and maximum x and y values
; for the output image (starting counting from 0)
;
; output: twodmodel.fits - output FITS file for reading into CASA
; 
; Example: IDL>make_2dimage,'einstein.fits',0,299,27,326
;
;
; the image is pointed at the CDFS, at 245GHz by default, with
; 0.043" pixels
;

pro make_2dimage,inname,xmin,xmax,ymin,ymax

	nslices=1

	im=readfits(inname, header)

        imnew = float(im[xmin:xmax,ymin:ymax,0])

        crp1 = float(xmax-xmin)/2.0
        crp2 = float(ymax-ymin)/2.0

	size=size(imnew)	

	im2d=fltarr(size[1],size[2],nslices)

	im2d[*,*,0] = imnew

	sxaddpar, header, 'CTYPE1','RA---SIN','First axis is RA'

	sxaddpar, header, 'CTYPE2','DEC--SIN','Second is DEC'

	sxaddpar, header, 'CTYPE3','FREQ','Third is frequency'

	sxaddpar, header, 'CRVAL1',52.5,'RA ref pos in Chandra deep field south'
	sxaddpar,header,'NAXIS3',1

	sxaddpar, header, 'CRVAL2',-28.5,'DEC ref position'

	sxaddpar, header, 'CRVAL3',245e9,'Reference frequency; obs at 245GHz'

	sxaddpar, header, 'CDELT1',-1.183e-5,'RA step = .043arcsec'

	sxaddpar, header, 'CDELT2',1.185e-5,'DEC step = .043arcsec'

;	sxaddpar, header, 'CDELT3',.006e9,'FREQ step = .006GHz slices'

	sxaddpar, header, 'CRPIX1',crp1,'RA reference point pixel'

	sxaddpar, header, 'CRPIX2',crp2,'DEC reference point pixel'

	sxaddpar, header, 'CRPIX3',1.0,'FREQ ref image is the first one'

	sxaddpar,header,'EPOCH',2000.0

	print,'header appended'

	writefits, 'twodmodel.fits',im2d,header

	print,'fits file written again'

end