function getscan,			$	; Input
 	image,				$	; *
	index=index,			$	; * (Yohkoh images)
 	in_xc,				$	; *
 	in_yc,				$	; *
 	in_radius,			$	; *
 	threepoints	= threepoints,	$
 	offset		= offset,	$
 	rsun		= rsun,		$	; * (non-Yohkoh images)
 	resolution	= resolution,	$	; * (non-Yohkoh images)
 	dr		= dr_in,	$	; *
 	center		= center,	$	; *
 	evans		= evans,	$	; *
 	dpa		= dpa_in,	$	; *
 	pa		= pa,		$	; *
 	clockwise	= clockwise,	$	; *

 	nscan		= nscan,	$	; Output
 	minmax		= minmax,	$

 	plotscan	= plotscan,	$ 	; Plot control
 	new		= new,		$
 	preview		= preview,	$
 	flip		= flip,		$
 	oplot		= oplot,	$
  	hp		= hp

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	GETSCAN
; PURPOSE:
;	Extract a spherical scan from an image
; CALLING SEQUENCE:
;	yscan = getscan(image, xc, [yc, radius,	$
;		dr=dr, low=low, dpa=dpa, pa=pa])
;	yscan = getscan(data,index=index,radius,...)
; INPUTS:
;	The units of xc and yc is the pixel size. The units of rc depend
;	on the setting of the keywords `rsun': if `rsun' is set then rc is
;	in units of `rsun'; if not, then rc is given in pixels.
;
;	For non-Yohkoh images:
;
;	image	2D image array, or 3D image cube
;	xc, [yc,radius]
;		defines the circle on the image where the scan is extracted:
;		(xc,yc) is the center of the circle.
;	If only 1 argument is specified it should be a 3-element array with
;		center coordinates and radius, or a scalar containing the
;		radius (the center of the image will be taken as the center for
;		the scan).
;	If 2 arguments are specified, the first should be a 2-element array
;		with the center coordinates; the 2nd should be a scalar for the
;		radius
;	If 3 arguments are specified, they should be all scalars for center
;		coordinates and radius
;
;	For Yohkoh images:
;
;	index=index	standard Yohkoh index structure
;	data		Yohkoh image corresponding to the index structure
;	radius		scan radius (default=1.15) (in units of the solar radius)
; OPTIONAL INPUTS:
;	threepoints=[x1,y1,x2,y2,x3,y3]
;		6-element vector specifying three points on the scan circle
;		in the form [x1,y1,x2,y2,x3,y3]. This keyword overrides the
;		arguments xc,yc,radius.
;	offset=[xoff,yoff]
;		2-element vector defining an offset to be applied to the
;		center (xc,yc). Useful if a portion of a larger image is
;		manipulated; if the portion (x1:x2,y1:y2) is extracted the
;		offset would be (x1,y1), while the (xc,yc) values would still
;		be the center coordinates in the full image (obviously the
;		offset can also be applied directly to (xc,yc))
;	resolution=resolution	(non-Yohkoh images only)
;		image resolution (arcsec/pixels)
;		!!! for Yohkoh images this keyword is internally set
;	rsun=rsun		(non-Yohkoh images only)
;		solar radius in
;		pixels (if `resolution' not set), or
;		arcsec (if `resolution' is  set, or input is Yohkoh image)
;		!!! for Yohkoh images this keyword is internally set
;	dr=dr	radial width of the scan (default=0.1), specified in units of
;		the scan radius (if `rsun' not set), or
;		the solar radius (if `rsun' is  set, or input is Yohkoh image)
;		Unless the keyword /center is set, pixels in [radius,radius+dr]
;		are extracted.
;	/center	if set, pixels in [radius-0.5*dr,radius+0.5*dr] are extracted
;	pa=pa	array of position angles in [0,360.], measured counterclockwise
;		from the top of the scan circle. If not specified then
;		[0,1,2,...,360] is used (see PROCEDURE).
;	/clockwise
;		if set, the pa angle is measured clockwise from the top
;		('north') of the image
;	dpa=dpa	angular width (in degrees) of position angle bins. If not
;		specified the value pa(1)-pa(0) is used (see PROCEDURE).
;	/evans	If set, a scan with the same characteristics as the Evans scans
;		is made:
;		- position angles pa = 3*indgen(120) deg
;		- position angle width dpa = 3 deg
;		- radial width dr = 1.1' (centered on the specified scan radius)
;		This keyword overrides the keywords pa,dpa,dr and /center
;		!!! for non-Yohkoh images the `resolution' keyword MUST be
;			specified too
;	/preview
;		displays the image with edges of scan circle overplotted
;		(the intensities along the scan are not calculated)
;	/oplot	suppresses display of image; in combination with /preview
;		this can be used to overplot multiple scan circles on the
;		the same image.
;	/plotscan
;		displays image; calculates scan intensities and displays
;		results in separate window
;	/new	opens new windows to plot results (by default the windows from
;		a previous call are used again).
;	/flip	will put the image on the screen after rotating 90 degrees
;		clockwise
; OUTPUTS:
;	yscan	1-dim array of intensity values along the scan.
; OPTIONAL OUTPUTS:
;	nscan=nscan
;		# pixels used in the averaging for each bin in the scan
;	minmax=minmax
;		minimum and maximum value in each bin
; RESTRICTIONS:
;	Image cubes are processed internally by calling GetScan recursively
;	for each of the images in the cube. Only a limited number of arguments
;	will be passed to the recursive GetScan call (marked by * in the
;	declaration below). In particular none of the plot control keywords
;	will be passed. The only output is the intensity array.
; PROCEDURE:
;	For each of the position angles in `pa' all pixels inside the
;	image area defined by radius-center*0.5*dr+[0,dr] and
;	pa+0.5*[dpa,-dpa] are averaged to obtain the intensity.
; MODIFICATION HISTORY:
;	FEB-1995, Paul Hick (UCSD)
;	MAY-1995, Paul Hick (UCSD), added processing of image cube (NOT tested)
;-

plotscan= keyword_set(plotscan)
preview	= keyword_set(preview)
oplot	= keyword_set(oplot)
hp	= keyword_set(hp)

center	= keyword_set(center)
plot	= keyword_set(plot)
flip	= keyword_set(flip)
evans	= keyword_set(evans)

clockwise = keyword_set(clockwise)
pmclock  = 1B-2B*clockwise	; plus/minus 1

rsun_set = n_elements(rsun) ne 0
reso_set = n_elements(resolution) ne 0

yohkoh = n_elements(index) ne 0

a = size(image)
case a(0) of
2: begin			; Single 2D image
    nX = a(1)  &  nY = a(2)	; Image dimensions
end
3: begin			; Image cube processing
    message, /info, 'Image cube processing has not been tested properly

; this can be made more general by syncing some of the arguments before
; making the recursive call to getscan, I think

    for img=0,a(3)-1 do begin	; Loop over all images
	case yohkoh of
	0: y=getscan(				$
		image(*,*,img),			$
		in_xc,				$
		in_yc,				$
		in_radius,			$
		rsun		= rsun,		$
		resolution	= resolution,	$
		dr		= dr_in,	$
		center		= center,	$
		evans		= evans,	$
		dpa		= dpa_in,	$
		pa		= pa,		$
		clockwise	= clockwise)
	1: y=getscan(				$
		image(*,*,img),			$
		index=index(img),		$
		in_xc,				$
		dr		= dr_in,	$
		center		= center,	$
		evans		= evans,	$
		dpa		= dpa_in,	$
		pa		= pa,		$
		clockwise	= clockwise)
	endcase
					; Create output array after 1st image
	if img eq 0 then yy = fltarr(n_elements(y),a(3))
	yy(*,img) = y			; Copy getscan output to output array

    endfor
    return, yy
end
else: message, 'IMAGE must be a 2D image array, or 3D image cube'
endcase

case yohkoh of
0: nparam = n_params()
1: nparam = 2			; accept only image and in_xc arguments
endcase

case n_elements(threepoints) of
0: begin
    case nparam-1 of			; Check arguments
    1: begin				; Only argument in_xc present
	case n_elements(in_xc) of
	1: begin			; in_xc = radius, take center of image
	    case yohkoh of
	    0: begin
		radius = in_xc  &  xc = (nX-1)/2.  &  yc = (nY-1)/2.
	    end
	    1: begin
		radius = in_xc		; Scan radius
;===
;		yoscan, index, yc, rsun, resolution

		resolution = 2^gt_res(index)		; Pixel size in full res pixels
		yc   = sxt_cen(index)/resolution	; Center of Sun (pixels)
		rsun = get_rb0p(index,/r)		; Solar radius (arcsec)
		resolution = pixsiz*gt_pix_size()	; Pixel size (arcsec)
;===
		xc = yc(0) &  yc = yc(1); Center of Sun (pixels)
		rsun_set = 1B  &  reso_set = 1B
	    end
	    endcase
	end
	3: begin			; in_xc = [xc,yc,radius]
	    xc = in_xc(0)  &  yc = in_xc(1)  &  radius = in_xc(2)
	end
	else: message, 'invalid combination of arguments'
	endcase
    end
    2: begin				; Arguments in_xc, in_yc present
	case n_elements(in_xc) of
	2: begin			; in_xc = [xc,yc], in_yc = radius
	    xc = in_xc(0)
	    yc = in_xc(1)
	    radius = in_yc
	end
	else: message, 'invalid combination of arguments
	endcase
    end
    3: begin				; Arguments in_xc, in_yc, in_radius present
	xc = in_xc(0)  &  yc = in_yc(0)  &  radius = in_radius(0)
    end
    else: message, 'invalid combination of arguments'
    endcase
end
6: begin				; 3 points on circle specified
    x = float(threepoints([0,2,4]))	; Calculate center and radius
    y = float(threepoints([1,3,5]))
    r = x^2+y^2
    a = shift(r-shift(r,-1),-1)
    xc = 0.5*total(y*a)/total(y*shift(x-shift(x,-1),-1))
    yc = 0.5*total(x*a)/total(x*shift(y-shift(y,-1),-1))
    radius = sqrt( (xc-x(0))^2+(yc-y(0))^2 )
    print, ' Calculated scan circle from points:'
    print, ' (',x(0),',',y(0),'),(',x(1),',',y(1),'),(',x(2),',',y(2),')'
;   threepoints = [xc,yc,r(0)]		; Return center and radius
end
else: message, 'THREEPOINTS must be a 6 element array [x1,y1,x2,y2,x3,y3]'
endcase

if rsun_set then begin		; rsun converted from arcsec to pixels
    if reso_set then rsun = rsun/resolution
    radius = radius*rsun	; radius converted from rsun to pixels
endif

if evans then begin		; Use aperture width of Sac. Peak Evans coronagraph
    if not reso_set then message, 'If keyword /evans is used, then also specify resolution
    center = 1
    dr  = 1.1*60/resolution	; Radial scan width of 1.1' (converted to pixels)

    pa  = 3*indgen(120)
    dpa = 3
endif else begin

; Defaults: scan from 0 to 360 in steps of 1 degree, using a bin with
; radial size of 0.1*radius and azimuthal size of 1 deg

    if n_elements(dr_in) eq 0 then dr = 0.1 else dr = dr_in
    case rsun_set of
    0: dr = dr*radius		; dr converted from `radius' to pixels
    1: dr = dr*rsun		; dr converted from `rsun' to pixels
    endcase

    if n_elements(pa)     eq 0 then pa = indgen(361)
    if n_elements(dpa_in) eq 0 then dpa = pa(1)-pa(0) else dpa = dpa_in
endelse

rmin = radius-center*0.5*dr  &  rmax = rmin+dr
dpa  = 0.5*abs(dpa)/!radeg

print, ' Image :',strcompress(nX),' x',strcompress(nY)
print, ' Center: (',xc,',',yc,')'
print, ' Input scan radius :',radius,' pixels'
print, ' Inner scan radius :',rmin,' pixels'
print, ' Outer scan radius :',rmax,' pixels'
print, ' Radial scan width :',dr,' pixels'
print, ' Pos. angle width  :',2*dpa*!radeg,' deg'

case n_elements(offset) of
0: offset = [0B,0B]		; Offset is lower left corner of image
2: a = 1			; Dummy statement
else: message, 'OFFSET must be a 2 element array'
endcase

xc = xc-offset(0)  &  yc = yc-offset(1)

if preview then goto, plotscan

x = (indgen(nX)-xc)#replicate(1B,nY)	; Set up x and y arrays in pixel units
y = replicate(1B,nX)#(indgen(nY)-yc)	; Shift origin to center of image

r = sqrt( x*x+y*y )			; Radial distance from center
ok = where( rmin lt r and r lt rmax, nok )
if nok eq 0 then message, 'no pixels inside specified radial range

im = image(ok)				; Extract pixels inside scan
x = x(ok)  &  y = y(ok)  &  r = r(ok)
paxy = -pmclock*atan(x,y)
					; Phase angle in [-pi,pi]
ok = 0					; Not needed anymore

;if bitmap then begin
;    len = fix( (2*!pi)*0.5*(rmin+rmax) )
;    thi = fix( dr )
;    scan = replicate( image(0,0), len, thi )
;    xs = indgen(len)*(360./(len-1))
;    ys = radius+indgen(thi)
;    for i=0,len-1 do begin
;
;    endfor
;endif
r = 0  &  x = 0  &  y = 0		; Not needed anymore

x = pa
x = (x-(x gt 180)*360.)/!radeg		; Stay inside [-pi,pi]

if min(x-dpa lt -!pi) then begin
    a = where( paxy gt !pi-dpa, na )
    if na ne 0 then begin paxy = [paxy,paxy(a)-2*!pi]  &  im = [im,im(a)] & endif
endif

if max(x+dpa gt !pi) then begin
    a = where( paxy lt -!pi+dpa, na )
    if na ne 0 then begin paxy = [paxy,paxy(a)+2*!pi]  &  im = [im,im(a)] & endif
endif

; The do-loop will always work. The second piece of code does the same job
; by manipulating arrays. It creates some scratch arrays which can cause a
; `not enough core' message, or even cause IDL to crash with a `process quota
; exceeded' message. It is also slower than the do-loop (due to the use of the
; temporary function????).

nx = n_elements(x)
y = 0.*x  &  nscan = long(y)  &  minmax = fltarr(2,nx)
for i=0,nx-1 do begin
    a = where( x(i)-dpa lt paxy and paxy lt x(i)+dpa , na)
    if na ne 0 then begin
	a = im(a)
	y(i) = total(a)/na
	nscan(i) = na
	minmax(0,i) = min(a,max=b)  &  minmax(1,i) = b
    endif
endfor
minmax = reform(minmax)

;nx = replicate(1,nx)  &  nok = replicate(1,nok)
;paxy = nx#temporary(paxy)
;x = temporary(x)#nok
;x = abs( temporary(x)-paxy )
;paxy = 0
;x = temporary(x) lt dpa
;im = nx#temporary(im)
;im = temporary(im)*x
;x  = x#nok
;im = im#nok  &  im = float(im)
;y = (x gt 0)*(im/(x > 1))

plotscan:

if preview or plotscan then begin		; Plot results
						; Select Twin windows
    common getscan_save, winopen, wins, winslast
    x = -1  &  twin, winmodes=x
    if n_elements(winopen) eq 0 then begin
	winopen = 0
	winslast = x
    endif
    winopen = winopen and not keyword_set(new) and (where(x ne winslast))(0) eq -1

    if not winopen then begin
	wins = (where(x eq -1))(0:1)
	winopen = 1
    endif
    twin,/show,winnr=wins(0)			; Open first Twin window
    winslast = -1  &  twin, winmodes=winslast

    x = pa/!radeg				; Display image ..
    if flip then begin
	if not oplot then view,in=rotate(image,3); .. unless /oplot is set
	nX = (size(image))(1)
	a = xc  &  xc = yc  &  yc = nX-1-a
    endif else begin
	if not oplot then view,in=image
    endelse

    a = min(x)-dpa+(max(x)-min(x)+2*dpa)*findgen(180+1)/180.
    a = pmclock*a
    plot, xc-rmin*sin(a), yc+rmin*cos(a),  $	; Inner boundary of scan
	xstyle=5, ystyle=5, pos=[0,0,1,1], $
	xrange=[0,!d.x_size], yrange=[0,!d.y_size],$
	/noerase;,linestyle=1
    if rmax ne rmin then begin			; Outer boundary of scan
	oplot, xc-rmax*sin(a), yc+rmax*cos(a),linestyle=1
	for i=0,n_elements(x)-1 do begin	; Mark pa edges of scan elements
	    a = pmclock*x(i)-dpa
	    oplot,xc-[rmin,rmax]*sin(a),yc+[rmin,rmax]*cos(a)
	    a = pmclock*x(i)+dpa
	    oplot,xc-[rmin,rmax]*sin(a),yc+[rmin,rmax]*cos(a)
	endfor
    endif
    oplot, [xc],[yc], psym=1			; Mark center of scan circle
						; Display scan in 2nd window
    if not preview and n_elements(pa) gt 1 then begin

	case hp of
	0: begin
	    twin,/show,winnr=wins(1),xsize=400,ysize=300	; Open 2nd Twin window
	    winslast = -1  &  twin, winmodes=winslast
	end
	1: set_page				; Hardcopy of scan
	endcase

	a = ['counter','']
	plot, pa, y, xstyle=1, /ynozero, $	; Plot intensity vs pa
		xtitle='POS ANGLE ('+a(clockwise)+'clockwise from top)'

	if hp then spitplot			; Send IDL.HP to printer

    endif

endif

if n_elements(y) eq 0 then y = 0	; Make sure to return something
return, y  &  end