;+
; NAME:
;	FITSIMAGE
; Purpose:
;	The FITSIMAGE procedure 
;
; CALLING SEQUENCE:
;	fitsimage,image,hdr,outimage,outhdr=outhdr,rsun=rsun,solr=solr,pixr=pixr,$
;		psgrid=psgrid,range=range,normalize=normalize,cdelt=cdelt,min=min
; INPUTS:
;	  image - 2D-array 
;	  hdr - FITS header or STRUCTURE (returned by MREADFITS)
;
; OUTPUTS:
;	  outimage - output image
;
; OPTIONAL INPUT KEYWORDS:
;	psgrid - number of pixelof the image at PS devise (default=512)
;	/normalize - normalize output image to byte array
;	outhdr = FITS header with new image range
;	range - xyrange of showed region;[x0,y0,x1,y1]
;	cdelt = option to select FITS keyword. If CDELT is set to unit, then 
;		keywords CDELTn are used, otherwise RSUN is used. 
;	rsun = adjusted solar radius, default=960.arcsec
;	solr = keyword of FITS header contained the observed 
;		solar radius (in arcsec), default value is 'SOLR' or
;		it can be  floating-point value
;	pixr = keyword of FITS header contained the observed 
;		solar radius (in image pixels), default value is 'R_SUN' 
;	min  = filling of empty area of the image by minimum value of the array  
;
; FITS HEADER KEYWORDS:
;   This procedure uses the following keywords from FITS header:
;	CDELT1, CDELT2 
;	CRPIX1, CRPIX2 
;	CRVAL1, CRVAL2 
;	R_SUN - observed solar radius in pixels
;	or
;	SOLR - observed solar radius in arcsec.
;
; MODIFICATION HISTORY:
; 	Written by:	Vladimir Garaimov, May 2002
;-
pro fitsimage,img,hdr,outimg,outhdr=outhdr,rsun=rsun,solr=solr,pixr=pixr,$
    psgrid=psgrid,range=range,normalize=normalize,cdelt=cdelt,min=min

if n_params(0) lt 3 then begin
	doc_library,'fitsimage'
	return
end

if n_elements(psgrid) ne 0 then psn=psgrid(0) else psn=512
if n_elements(solr) eq 0 then _solr='SOLR' else _solr=solr(0)
if n_elements(pixr) eq 0 then _pixr='R_SUN' else _pixr=pixr(0)

position=fltarr(4)
position(0)=0
position(2)=psn
position(1)=0
position(3)=psn

x_0=position(0)
x_1=position(2)
y_0=position(1)
y_1=position(3)

xx=x_1-x_0 & yy=y_1-y_0
xdy=float(xx)/float(yy)
if n_elements(range) eq 4 then begin
  n1=float(range(2)-range(0))
  n2=float(range(3)-range(1))
  if n2 eq 0 then n1n2=1. else n1n2=n1/n2
end else begin
  nn=size(img) 
  n1n2=float(nn(1))/float(nn(2))
 end

if xdy gt n1n2 then begin
 xx=(yy*n1n2) & x_1=x_0+xx
 position(2)=position(0)+xx
endif else begin
 yy=(xx/n1n2) & y_1=y_0+yy
 position(3)=position(1)+yy
 end

swx=x_1-x_0
swy=y_1-y_0

;image calculations
thkname=strarr(30) & thkname(*)=' '
mx0=0. & my0=0. & mcdelt1=1.& mx_rad=960
i=size(hdr)
if i(i(0)+1) ge 7 then begin
 mx0=sx_par(hdr,'CRPIX1')-1.0
 my0=sx_par(hdr,'CRPIX2')-1.0
 mcdelt1=abs(sx_par(hdr,'CDELT1'))>1e-5
 mcdelt2=abs(sx_par(hdr,'CDELT2'))>1e-5
 i=sx_par(hdr,'CRVAL1')
 mx0=mx0-i/mcdelt1
 i=sx_par(hdr,'CRVAL2')
 my0=my0-i/mcdelt2
;radius
 mx_rad=sx_par(hdr,_pixr)
 if mx_rad eq 0 then begin
	i=size(_solr)
	if i(i(0)+1) eq 7 then mx_rad=sx_par(hdr,_solr) $
	else mx_rad=float(_solr)
	mx_rad=mx_rad/mcdelt1
	end
 if mx_rad eq 0 then begin
	mx_rad=960./mcdelt1
	print,'Solar Radius is not defined. use default value 960arcsec.'
	end
end

nn=size(img)-1

if n_elements(rsun) ne 0 then rs=rsun(0) else rs=960.
if not keyword_set(cdelt) then mcdelt1=float(rs)/mx_rad

xx=[0.,0.] & yy=xx
xx(0)=-mx0*mcdelt1
xx(1)=xx(0)+nn(1)*mcdelt1
yy(0)=-my0*mcdelt1
yy(1)=yy(0)+nn(2)*mcdelt1

if n_elements(range) ne 4 then begin
  img1=img & xy0=[0.,0.] 
endif else begin

 if range(0) ge range(2) or range(1) ge range(3) then begin
  print,'Range Box is wrong!' & return 
 end

 if xx(0) ge range(2) or yy(0) ge range(3) or xx(1) $
    le range(0) or yy(1) le range(1) $
 then begin&print,'Overlaing box is empty'&return&end 

img1=img 
fn=1 

 x0=xx-range([0,2])
 y0=yy-range([1,3])
 xy0=[0.,0.] & i0=0 & j0=0
 if x0(0) ge 0 then xy0(0)=x0(0) else begin
  i=-x0(0)/mcdelt1 & i0=fix(i)
  if i-i0 ne 0 and fn then begin
    xy0(0)=float(i0+1-i)*mcdelt1 & i0=i0+1 
  end else xy0(0)=0
  end
  if y0(0) ge 0 then xy0(1)=y0(0) else begin
    i=-y0(0)/mcdelt1 & j0=fix(i)
    if i-j0 ne 0 and fn then begin
      xy0(1)=float(j0+1-i)*mcdelt1 & j0=j0+1
    end else xy0(1)=0
  end

  i1=nn(1) &j1=nn(2)
  if x0(1) gt 0 then begin
    i=x0(1)/mcdelt1 & i1=i1-fix(i)
    if i-fix(i) ne 0 then i1=i1-1 
  end
  if y0(1) gt 0 then begin
     i=y0(1)/mcdelt1 & j1=j1-fix(i)
     if i-fix(i) ne 0 then j1=j1-1 
  end
  if i0 ge i1 or j0 ge j1 then return
  img1=img1(i0:i1,j0:j1)
  nn=size(img1)-1

 xx(0)=range(0) & xx(1)=range(2)
 yy(0)=range(1) & yy(1)=range(3)
 i=xx(1)-xx(0)
 xy0(0)=xy0(0)*swx/i
 i=float(nn(1))*mcdelt1/i
 swx=fix(float(swx)*i)< (x_1-x_0+1.0)
 i=yy(1)-yy(0) 
 xy0(1)=xy0(1)*swy/i
 i=float(nn(2))*mcdelt1/i
 swy=fix(float(swy)*i)< (y_1-y_0+1.0)
end

bb=congrid(img1,swx,swy,/interp,/minus_one)
nxy=size(bb)

outimg=fltarr(round(position(2)),round(position(3)))
if n_elements(min) eq 0 then outimg(*,*)=median(bb) else outimg(*,*)=min(bb)
outimg(fix(x_0+xy0(0)),fix(y_0+xy0(1)))=bb

if keyword_set(normalize) then begin
 mx=max(outimg,min=mn)
 outimg=byte(255.*(outimg-mn)/(mx-mn))
end

xdel=(xx(1)-xx(0))/float(position(2)-1)

;outhdr=hdr
sxaddpar,outhdr,'NAXIS',2
sxaddpar,outhdr,'NAXIS1',round(position(2))
sxaddpar,outhdr,'NAXIS2',round(position(3))
sxaddpar,outhdr,'CRPIX1',-xx(0)/xdel+1.
sxaddpar,outhdr,'CRPIX2',-yy(0)/xdel+1.
sxaddpar,outhdr,'CRVAL1',0.0
sxaddpar,outhdr,'CRVAL2',0.0
sxaddpar,outhdr,'CDELT1',xdel
sxaddpar,outhdr,'CDELT2',xdel
sxaddpar,outhdr,'R_SUN',rs/xdel
sxaddpar,outhdr,'SOLR',rs
sxaddpar,outhdr,'DATE_TIME',sx_par(hdr,'DATE_TIME')

range=fltarr(4)
range(0)=xx(0)
range(1)=yy(0)
range(2)=xx(1)
range(3)=yy(1)

end