;+
;
; NAME       : newkirk (function)
; PURPOSE    :
;	perform Newkirk filter operation on 10-cm coronagraph data
; CATEGORY :
;	idl/lib/nkr
; CALLING SEQUENCE :
;  img = newkirk(img0,norder=norder,dh=dh,nrep=nrep,verbose=ver, 	$
;	flat=flat,od=od,round=round,skyint=skyint,sclbyt=sclbyt,	$
;	trend=trend,com=com,sbt=sbt,sample=sample,skycoeff=skycoeff,	$
;	fltcoeff=fltcoeff,i0=i0,enlarge=enlarge)
; INPUTS :
; 	img0     -- image data assumed to be FIX type
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
; 	norder  -- order of polinomial (default=5)
; 	dh      -- skin height above the limb (default=2)
;	nrep    -- repeat no. of fitting procedure for lower edge of sky
;			 (default=0)
; 	/ver    -- display fitting curve
; 	flat    -- flat image if make correction before newkirk
; 	/od     -- use occulting disk edge to detect the limb
; 	skyint  -- return sky intensity
;	sclbyt(2) -- scale 0-255 with range of sky*[(1+sclbyt(0))~(1+sclbyt(1))]
;	trend   -- remove linear trend of sky and set height from the limb
;	com     -- return ox,oy,r information
;	/sbt    -- make subtraction other than division
;	sample  -- sampling step of sky (default=10)
;	skycoeff -- return coefficients of sky fitting polinomial
;	fltcoeff -- coefficients of sky fitting polinomial for flat data
; 	      if sbt is set, re-correct flat inhomogeneity
;	i0      -- return central intensity for calibration
;	/enlarge -- operation for enlarged image ('left', 'right')
; OUTPUTS :
;	none
; COMMON BLOCKS : 	none
; SIDE EFFECTS :
;	if /ver is set, fitting curve is displayed
; RESTRICTIONS :	none
; PROCEDURE :
; MODIFICATION HISTORY :
;	K.I. 1991/12/31	
;	K.I. 1992/04/22
;-

;*************************************************************************
;NAME       : center (procedure)
;FUNCTION   : get circle center from limb points x,y
;=========================================================================
pro center,x,y,ox,oy,r

nn=n_elements(x)

; ------------  get ox & oy from 3-points ----------
x1=x(0)	& x2=x(nn/2)	& x3=x(nn-1)
y1=y(0)	& y2=y(nn/2)	& y3=y(nn-1)
a12=x1-x2
b12=(y1-y2)
c12=(y1^2-y2^2+x1^2-x2^2)/2.
a23=x2-x3
b23=(y2-y3)
c23=(y2^2-y3^2+x2^2-x3^2)/2.
ox=+(c12*b23-c23*b12)/(a12*b23-a23*b12)
oy=-(c12*a23-c23*a12)/(a12*b23-a23*b12)

rsam=total(sqrt((x-ox)^2+(y-oy)^2))
r=rsam/nn

end

;*************************************************************************
;NAME       : odcenter2 (procedure)
;FUNCTION   : obtain center of occulting disk & its radius 
;			for enlarged 10cm image
;=========================================================================
pro odcenter2,img,ox,oy,r,part=part
; part  -- which limb? 'left' or 'right'

nn=5	; no. of edge detection cross line
s=size(img)	& nx=s(1)	& ny=s(2)
dy=10
ycross=(ny-2*dy)/(nn-1)*(indgen(nn))+dy
x1=fltarr(nn)
x2=fltarr(nn)
for n=0,nn-1 do begin
	if n eq nn/2 then get2edges,img(*,ycross(n)),edge1,edge2,region=0.3 $
	else  get2edges,img(*,ycross(n)),edge1,edge2
	x1(n)=edge1
	x2(n)=edge2
endfor
if keyword_set(part) then begin
	if part eq 'left' then center,x1,ycross,ox,oy,r	$
	else if part eq 'right' then center,x2,ycross,ox,oy,r
endif else begin
	ox=fltarr(2)
	oy=fltarr(2)
	r=fltarr(2)
	center,x1,ycross,ox0,oy0,r0
	ox(0)=ox0	& oy(0)=oy0	& r(0)=r0
	center,x2,ycross,ox0,oy0,r0
	ox(1)=ox0	& oy(1)=oy0	& r(1)=r0
endelse

end

;*************************************************************************
;NAME       : limbcenter (procedure)
;FUNCTION   : obtain center of the limb where intensity is icrit
;=========================================================================
pro limbcenter,img,ox,oy,r,rxy=rxy
; rxy   --  ration between x & y-radius

nn=3	; no. of edge detection cross line = nn*2+1
x1=intarr(2*nn+1)
x2=intarr(2*nn+1)
yc=intarr(2*nn+1)
y1=intarr(2*nn+1)
y2=intarr(2*nn+1)
xc=intarr(2*nn+1)

s=size(img)
nx=s(1)
ny=s(2)
nx2=nx/2
ny2=ny/2

;----- O.D. edge on central cross section -----
get2edges,img(*,ny2:ny2),xmin,xmax,region=0.33
get2edges,transpose(img(nx2:nx2,*)),ymin,ymax,region=0.33
dx=(xmax-xmin)/(nn+1)/2
dy=(ymax-ymin)/(nn+1)/2
;----- get the intensity of limb determination -----
imax=intarr(4*(nn*2+1))
for n=0,nn*2 do begin
	yc(n)=ny2+dy*(n-nn)
	imax(4*n)   = max(img(0:nx2,yc(n):yc(n)),m)
	imax(4*n+1) = max(img(nx2:nx-1,yc(n):yc(n)),m)
	xc(n)=nx2+dx*(n-nn)
	imax(4*n+2) = max(transpose(img(xc(n):xc(n),0:ny2)),m)
	imax(4*n+3) = max(transpose(img(xc(n):xc(n),ny2:ny-1)),m)
endfor
limb=min(imax)
sky =(total(img(*,0:0))+total(img(*,ny-1:ny-1)))/nx/2
icrit=sky+(limb-sky)*0.9
;----- obtain limb positions where intensity is icrit -----
for n=0,nn*2 do begin
	wp=where(img(*,yc(n):yc(n)) gt icrit)
	x1(n)=wp(0)
	x2(n)=wp(n_elements(wp)-1)
	wp=where(transpose(img(xc(n):xc(n),*)) gt icrit)
	y1(n)=wp(0)
	y2(n)=wp(n_elements(wp)-1)
endfor

;----- obtain center pos. & radius -----
ox = float(total(x1)+total(x2))/float(2*(nn*2+1))
oy = float(total(y1)+total(y2))/float(2*(nn*2+1))
r = ( total(sqrt((x1-ox)^2+(yc-oy)^2))+total(sqrt((x2-ox)^2+(yc-oy)^2))	  $ 
    + total(sqrt((xc-ox)^2+(y1-oy)^2))+total(sqrt((xc-ox)^2+(y2-oy)^2)) ) $ 
			/float(4*(nn*2+1))
;----- obtain ratio between x & y radius -----
if keyword_set(rxy) then begin
    for n=0,nn*2 do begin
	wp=where(img(*,oy-nn+n:oy-nn+n) gt icrit)
	x1(n)=wp(0)
	x2(n)=wp(n_elements(wp)-1)
	wp=where(transpose(img(ox-nn+n:ox-nn+n,*)) gt icrit)
	y1(n)=wp(0)
	y2(n)=wp(n_elements(wp)-1)
    endfor
    rx=float(total(x2)-total(x1))/float(2*(nn*2+1))
    ry=float(total(y2)-total(y1))/float(2*(nn*2+1))
    rxy=rx/ry
endif

end

;*************************************************************************
;NAME       : mkround (function)
;FUNCTION   : make coronagraph data to be exact round
;-------------------------------------------------------------------------
function mkround,img

s=size(img)
nx=s(1)
ny=s(2)
nn=3	; no of cross line to determine edges = 2*nn+1

odcenter,img,ox,oy,r  	; determin the center of O.D.
if (ox le 0) or (ox ge nx) or (oy le 0) or (oy ge ny) then begin
	print,'cannot make makeround operation!!'
	return,img
endif

xmins=0
xmaxs=0
ymins=0
ymaxs=0
for i=0,nn*2 do begin
	get2edges,img(*,oy-nn+i:oy-nn+i),xmin,xmax,region=0.33
	xmins=xmins+xmin
	xmaxs=xmaxs+xmax
	get2edges,transpose(img(ox-nn+i:ox-nn+i,*)),ymin,ymax,region=0.33
	ymins=ymins+ymin
	ymaxs=ymaxs+ymax
endfor
wx=float(xmaxs-xmins)/(nn*2+1)
wy=float(ymaxs-ymins)/(nn*2+1)
return,congrid(img,nx,ny*wx/wy,/interp)

end

;*************************************************************************
;NAME       : skypoint (procedure)
;FUNCTION   : get relation between height from the limb & sky intensity
;=========================================================================
pro skypoint,img,ox,oy,r,height,sky,sample=sample
; ox,oy	  - center of sun
; r       - radius of the sun
; height  - height from the limb
; sky     - sky intensity
; sample  - sampling step of sky

s=size(img)
nx=s(1)
ny=s(2)
di=10	; sampling step in x
dj=10	; sampling step in y
if keyword_set(sample) then begin
	di=sample
	dj=sample
endif
ni=(nx-1)/di
r2=r^2

i = di*where( (findgen(ni)*di-ox)^2+(0.-oy)^2 gt r2 )
height=sqrt((float(i)-ox)^2+(0.-oy)^2)-r
sky=img(i,0)
for j=dj,ny-1,dj do begin
	i = di*where( (findgen(ni)*di-ox)^2+(j-oy)^2 gt r2 , count)
	if(count ne 0l) then begin
		height=[height,sqrt((float(i)-ox)^2+(float(j)-oy)^2)-r]
		sky=[sky,img(i,j)]
	endif
endfor

end

;*************************************************************************
;NAME       : trend1 (procedure)
;FUNCTION   : get sky distribution around the limb
;=========================================================================
pro trend1,img,ox,oy,r,height=height
;INPUT
; img     - 10-cm image
; ox,oy	  - center of sun
; r       - radius of the sun
;OUTPUT
; img     - position

if keyword_set(height) then rh=height	$
else rh=10	; height from limb at which linear trend is determined
nn=360	; no. of points for plane fit
s=size(img)
nx=s(1)
ny=s(2)
th=2.*!pi/nn*indgen(nn)
x = ox+(r+rh)*cos(th)
y = oy+(r+rh)*sin(th)
i = where((x gt 0) and (x lt nx-2) and (y gt 0) and (y lt nx-2))
x=x(i)
y=y(i)
z = img(fix(x(i)),fix(y(i)))

aa=dblarr(3,3)
bb=dblarr(3)
aa(0,0)=total(x^2) &	aa(1,0)=total(x*y) &	aa(2,0)=total(x)
aa(0,1)=total(x*y) &	aa(1,1)=total(y^2) &	aa(2,1)=total(y)
aa(0,2)=total(x) &	aa(1,2)=total(y) &	aa(2,2)=n_elements(x)
bb(0)=total(x*z) &	bb(1)=total(y*z) &	bb(2)=total(z)
cc=invert(aa)
c = cc # bb
fit=intarr(nx,ny)
for j=0,ny-1 do begin
	fit(*,j) = c(0)*indgen(nx) + c(1)*j + c(2)
endfor
img=fix(float(img)/float(fit)*c(2))

end

;*************************************************************************
;NAME       : newkirk (function)
;FUNCTION   : perform Newkirk filter operation on 10-cm coronagraph data
;DATE 	    : 1991/12/31 k.i. 1992/04/22
;=========================================================================
function newkirk,img0,norder=norder,dh=dh,nrep=nrep,verbose=ver,	$
	flat=flat,od=od,round=round,skyint=skyint,sclbyt=sclbyt,	$
	trend=trend,com=com,sbt=sbt,sample=sample,skycoeff=skycoeff,	$
	fltcoeff=fltcoeff,i0=i0,enlarge=enlarge
; img0     -- image data assumed to be FIX type
; norder  -- order of polinomial (default=5)
; dh      -- skin height above the limb (default=2)
; nrep    -- repeat no. of fitting procedure for lower edge of sky (default=0)
; /ver    -- display fitting curve
; flat    -- flat image if make correction before newkirk
; /od     -- use occulting disk edge to detect the limb
; /round  -- make exact round correction (ny->ny*rx/ry)
; skyint  -- return sky intensity
; sclbyt(2) -- scale 0-255 with range of sky*[(1+sclbyt(0))~(1+sclbyt(1))]
; trend   -- remove linear trend of sky and set height from the limb
; com     -- return ox,oy,r information
; /sbt    -- make subtraction other than division
; sample  -- sampling step of sky (default=10)
; skycoeff -- return coefficients of sky fitting polinomial
; fltcoeff -- coefficients of sky fitting polinomial for flat data
; 	      if sbt is set, re-correct flat inhomogeneity
; i0      -- return central intensity for calibration
; /enlarge -- operation for enlarged image ('left', 'right')
;-------------------------------------------------------------------------
if not keyword_set(img0) then begin
	print,'NEWKIRK: make newkirk operation for coronagraph data'
	return,0
endif

img=img0
s=size(img)
nx=s(1)
ny=s(2)
ny0=ny
if not keyword_set(norder) then norder=5 
if not keyword_set(dh) then dh=2.	 
if not keyword_set(nrep) then nrep=0
if keyword_set(od) then od=1 else od=0
if keyword_set(round) then round=1 else round=0
if keyword_set(ver) then ver=1 else ver=0

;=========================================================================
if keyword_set(enlarge) then begin	; +++++++++ enlarged image  ++++++

; ------------  get center position and radius of O.D. ---------- 
odcenter2,img,oxa,oya,ra

wcal=nx/20	; width for searching calib. spot
nwcal=3		; sampling width for calib.
calimg=img(nx/2-wcal:nx/2+wcal,ny/2-wcal:ny/2+wcal)
mm=max(calimg,msb)
y_cal=msb/(2*wcal+1)
x_cal=msb-y_cal*(2*wcal+1)
if (x_cal-nwcal lt 0 or x_cal+nwcal ge 2*wcal+1) $
    or (y_cal-nwcal lt 0 or y_cal+nwcal ge 2*wcal+1) then $
	i0sum=total(img(nx/2-nwcal:nx/2+nwcal,ny/2-nwcal:ny/2+nwcal)) $
else 	i0sum=total(calimg(x_cal-nwcal:x_cal+nwcal,y_cal-nwcal:y_cal+nwcal))
i0=i0sum/(2*nwcal+1)^2
if keyword_set(ver) then begin
	print,format='("calibration spot:  x=",i3,"    y=",i3,"    i0=",i5)', $
		x_cal,y_cal,i0
endif


if string(enlarge) eq 'right' or string(enlarge) eq 'left' then  both=0 $
else both=1
if string(enlarge) eq 'right' then  side='right'	$
else side='left'

loop1:

if side eq 'right' then begin
	ox=oxa(1)	& oy=oya(1)	& r=ra(1)+dh
endif else begin
	ox=oxa(0)	& oy=oya(0)	& r=ra(0)+dh
endelse
com='nk: ox='+string(ox,format="(f6.2)")	$
	+' oy='+string(oy,format="(f6.2)") $
	+' r='+string(r,format="(f6.2)")	$
	+' ry='+string(float(ny)/ny0,format="(f5.3)")

if keyword_set(ver) then begin
   print,'side: ',side
   print,format='("O.D.fitting:  ox=",f7.2,"   oy=",f7.2,"   r=",f7.2,"   i0=",f7.2)',ox,oy,r,i0
endif

; -----------------  give flat & trend corrections  --------------------
if(n_elements(flat) gt 0) then begin	; flat correction
	if keyword_set(ver) then print,'flat correction:'
	fact=total(flat)/n_elements(flat)
	if side eq 'left' then begin
		half=img0(0:nx/2-1,*)
		half_flat=flat(0:nx/2-1,*)
		x00=0
	endif else if side eq 'right' then begin
		half=img0(nx/2:nx-1,*)
		half_flat=flat(nx/2:nx-1,*)
		x00=nx/2
	endif
	half=fix((float(half)/float(half_flat>1)*fact))
	if keyword_set(sbt) then begin	; re-correct flat slope
	    if not keyword_set(fltcoeff) then begin
		if keyword_set(ver) then print,'getting flat slope:'
		fltcoeff=1
		dumy=newkirk(flat,od=od,dh=dh,norder=norder	$
			,skycoeff=fltcoeff,ver=ver,enlarge=enlarge)
	    endif
	    if keyword_set(ver) then print,'re-correct flat slope:'
    	    for j=0,ny-1 do begin
		height=sqrt((findgen(nx/2)+x00-ox)^2+(float(j)-oy)^2) - r
		ws=where(height gt 0, count)
		if count ne 0l then  half(ws,j:j) = fix( float(half(ws,j:j)) $
				*poly(alog10(height(ws)+1),fltcoeff)/fact )
    	    endfor
	endif
	if side eq 'left' then img(0:nx/2-1,*)=half	$
	else if side eq 'right' then img(nx/2:nx-1,*)=half
endif


; ------------------  get sky distribution  -------------------
skypoint,img,ox,oy,r,height,sky,sample=sample ;	get sky distribution
coeff =poly_fit( alog10(height+1), sky, norder ) ; fitting of sky
for i = 1,nrep do begin
	wu=where( poly(alog10(height+1),coeff) gt sky )
	height=height(wu)
	sky=sky(wu)
	coeff =poly_fit( alog10(height+1), sky, norder )
	if keyword_set(ver) then $
		print,'repeat_fitting:',i,'   no of points=',n_elements(height)
endfor
if keyword_set(ver) then begin
	print,format='("skyfit:  coeff=",10f9.1)',transpose(coeff)
	plot,height,sky,psym=4,symsize=0.1
	skyfit=poly( alog10(lindgen(200)+1), coeff )
	oplot,indgen(200),skyfit
endif
if keyword_set(skycoeff) then begin
	skycoeff=coeff
	if keyword_set(ver) then print,'return coeff. of sky-fitting:'
	return,0
endif
fact=total(sky)/n_elements(sky)
skyint=fix(fact)

; --------------------  newkirk operation  --------------------------
if side eq 'left' then begin
	half=img(0:nx/2-1,*)
	x00=0
endif else if side eq 'right' then begin
	half=img(nx/2:nx-1,*)
	x00=nx/2
endif
if keyword_set(sbt) then begin
    if keyword_set(ver) then print,'newkirk operation with subtraction:'
    for j=0,ny-1 do begin
	height=sqrt((findgen(nx/2)+x00-ox)^2+(float(j)-oy)^2) - r
	ws=where(height gt 0, count)
	if count ne 0l then  half(ws,j:j) = fix( half(ws,j:j)	$
			- fix(poly(alog10(height(ws)+1),coeff)) )
	ws=where(height le 0, count)
	if count ne 0l then half(ws,j:j) = 0
    endfor
endif else begin
    if keyword_set(ver) then print,'newkirk operation with division:'
    for j=0,ny-1 do begin
	height=sqrt((findgen(nx/2)+x00-ox)^2+(float(j)-oy)^2) - r
	ws=where(height gt 0, count)
	if count ne 0l then  half(ws,j:j) = fix( float(half(ws,j:j))	$
			/poly(alog10(height(ws)+1),coeff)*fact )
	ws=where(height le 0, count)
	if count ne 0l then half(ws,j:j) = fix(fact)
    endfor
endelse

if both eq 1 then begin
    if side eq 'left' then begin
	img(0:nx/2-1,*)=half
	side='right'
	goto,loop1
    endif
    img(nx/2:nx-1,*)=half
endif else begin
    img=half
endelse

if n_elements(sclbyt) eq 2 then 	$
	img=bytscl(img,min=(skyint*(1.+sclbyt(0))),max=(skyint*(1.+sclbyt(1))))
return,img

;=========================================================================
endif else begin	; ++++++++++++++     round image     +++++++++++++

; ---------------  round correction  ----------------
ny0=ny
if keyword_set(round)  then begin	; modify y-scale to make round O.D.
	img=mkround(img)
	nx0=nx & ny0=ny
	s=size(img)	& nx=s(1)	& ny=s(2)
	if keyword_set(ver) then print,'ny modification : ',ny0,' --->',ny
endif

; ------------  get center position and radius od O.D. ----------
if keyword_set(od)  then begin
	odcenter,img,ox,oy,r  	; determin the center of O.D.
	r=r+dh
endif else begin
	limbcenter,img,ox,oy,r	; determin the center of limb
	r=r+dh-2.
endelse
i0=total(img(fix(ox)-1:fix(ox)+1,fix(oy)-1:fix(oy)+1))/3/3

com='nk: ox='+string(ox,format="(f6.2)")+' oy='+string(oy,format="(f6.2)") $
	+' r='+string(r,format="(f6.2)")+	$
	' ry='+string(float(ny)/ny0,format="(f5.3)")

if keyword_set(ver) then 	$
   print,format='("O.D.fitting:  ox=",f7.2,"   oy=",f7.2,"   r=",f7.2,"   i0=",f7.2)',ox,oy,r,i0
;circle,ox,oy,r
if (ox le 0) or (ox ge nx) or (oy le 0) or (oy ge ny) or (r gt ny/2) then begin
	print,'cannot make NEWKIRK operation!!'
	return,img0
endif

; -----------------  give flat & trend corrections  --------------------
if(n_elements(flat) gt 0) then begin	; flat correction
	if keyword_set(ver) then print,'flat correction:'
	fact=total(flat)/n_elements(flat)
	img=fix((float(img0)/float(flat>1)*fact))
	if keyword_set(round)  then img=congrid(img,nx,ny,/interp)
	if keyword_set(sbt) then begin	; re-correct flat slope
	    if not keyword_set(fltcoeff) then begin
		if keyword_set(ver) then print,'getting flat slope:'
		fltcoeff=1
		dumy=newkirk(flat,od=od,round=round,trend=trend,dh=dh,	$
			norder=norder,skycoeff=fltcoeff,ver=ver)
	    endif
	    if keyword_set(ver) then print,'re-correct flat slope:'
    	    for j=0,ny-1 do begin
		height=sqrt((findgen(nx)-ox)^2+(float(j)-oy)^2) - r
		ws=where(height gt 0, count)
		if count ne 0l then  img(ws,j:j) = fix( float(img(ws,j:j)) $
				*poly(alog10(height(ws)+1),fltcoeff)/fact )
    	    endfor
	endif
endif

if keyword_set(trend) then begin	; trend correction
	if keyword_set(ver) then print,'removing linear trend:'
	trend1,img,ox,oy,r,height=trend
endif

; ------------------  get sky distribution  -------------------
skypoint,img,ox,oy,r,height,sky,sample=sample	 ; get sky distribution
coeff =poly_fit( alog10(height+1), sky, norder ) ; fitting of sky
for i = 1,nrep do begin
	wu=where( poly(alog10(height+1),coeff) gt sky )
	height=height(wu)
	sky=sky(wu)
	coeff =poly_fit( alog10(height+1), sky, norder )
	if keyword_set(ver) then $
		print,'repeat_fitting:',i,'   no of points=',n_elements(height)
endfor
if keyword_set(ver) then begin
	print,format='("skyfit:  coeff=",10f7.1)',transpose(coeff)
	plot,height,sky,psym=4,symsize=0.1
	skyfit=poly( alog10(lindgen(200)+1), coeff )
	oplot,indgen(200),skyfit
endif
if keyword_set(skycoeff) then begin
	skycoeff=coeff
	if keyword_set(ver) then print,'return coeff. of sky-fitting:'
	return,0
endif
fact=total(sky)/n_elements(sky)
skyint=fix(fact)

; --------------------  newkirk operation  --------------------------
if keyword_set(sbt) then begin
    if keyword_set(ver) then print,'newkirk operation with subtraction:'
    for j=0,ny-1 do begin
	height=sqrt((findgen(nx)-ox)^2+(float(j)-oy)^2) - r
	ws=where(height gt 0, count)
	if count ne 0l then  img(ws,j:j) = fix( img(ws,j:j)	$
			- fix(poly(alog10(height(ws)+1),coeff)) )
	ws=where(height le 0, count)
	if count ne 0l then img(ws,j:j) = 0
    endfor
endif else begin
    if keyword_set(ver) then print,'newkirk operation with division:'
    for j=0,ny-1 do begin
	height=sqrt((findgen(nx)-ox)^2+(float(j)-oy)^2) - r
	ws=where(height gt 0, count)
	if count ne 0l then  img(ws,j:j) = fix( float(img(ws,j:j))	$
			/poly(alog10(height(ws)+1),coeff)*fact )
	ws=where(height le 0, count)
	if count ne 0l then img(ws,j:j) = fix(fact)
    endfor
endelse

if keyword_set(round)  then 	$	; recover original y-scale 
	img=congrid(img,nx0,ny0,/interp) 

if n_elements(sclbyt) eq 2 then 	$
	img=bytscl(img,min=(skyint*(1.+sclbyt(0))),max=(skyint*(1.+sclbyt(1))))
return,img

endelse

end