pro fit_limb, image, x0, y0, r0, r_err, $
		index=index, $
		summed_pixels=summed_pixels, $
		pixel_size=pixel_size, $
                verbose=verbose, $
                oblateness=oblateness,ob_angle=ob_angle, bias=bias,xlimb=xfit,$
                ylimb=yfit, interactive=interactive, quiet=quiet, $
                initial_values=initial,fast=fast, xinitial=hxa_x, $
                yinitial=hxa_y, $
                decompress=decompress, noiterate=noiterate,iterate=iterate, $
                kill=RLOOP,norfit=norfit,date=date, $
                noobcorrect=noobcorrect,ellipse=ellipse,resolution=resolution

;+
;NAME:
;	fit_limb
;PURPOSE:
;	Find the center and radius of an image.  NOTE: For SXT images, the
;	default is to return the result in full resolution pixels.
;CATEGORY:
;CALLING SEQUENCE:
;     fit_limb, image, x, y, r, pixel_size=2.5, date='05-Oct-93'
;     fit_limb, image, x, y, r, pixel_size=2.5, intial=[434,512,403*2.5]
;     fit_limb, data(*,*,nimage), x, y, r, index=index(nimage),/norfit
;INPUTS:
;     image	Full sun image for which the center and radius will 
;		be computed.  It should be a 2-D array (not a 3-D data cube)
;		It should be decompressed, and ideally background subtracted.
;OPTIONAL KEYWORD INPUTS:
;     index	SXT index structure.  If passed, it will use this value to
;		derive the initial guess of the location.
;     summed_pixels - If set, return the pixel coordinate of the sun center in
;		the summed pixels, not the full resolution pixel coordinates
;		which are the default (only relevant if index is set).
;     pixel_size - The pixel size in arcseconds.  If the SXT index is passed
;		it will be derived from that.
;     /verbose = if set, print a bunch of diagnostics
;     /ellipse = if set, fit an ellipse rather than a circle.  May work
;                better for oblate images, but should not be used for circular
;                images since it allows another degree of freedom in the fit.
;     oblateness = return second harmonic of radius
;     ob_angle = return phase of 2nd harmonic
;     /noobcorrect = do not correct final sun center value for oblateness
;     bias = return distortion of limb due higher harmonics
;     xlimb,ylimb = returns the coordinates of the selected limb pixels
;     /interactive = allows interactive selection of initial center
;                    guess and deselection of limb points.
;     /quiet = No images are displayed if set and interactive is not set
;     initial_values = 3-element vector giving initial guess: [x0,y0,r0]
;                      x0 and y0 are in pixels, r0 is in arcseconds !!
;                      If initial values are passed, HXA is NOT used and index 
;                      need not be passed as a parameter.  
;     /fast = If set, do the fast but slightly less accurate computation of
;             the limb position (limb derivative is not smoothed).
;     xinitial = returns  x coordinate of image center for the initial guess
;     yinitial = returns  y coordinate of image center for the initial guess
;     /decompress = If set, use sxt_decomp to decompress the image
;     /iterate = iterate in the circle fitting routine.  In principle this
;                could improve the accuracy. (default)
;     /noiterate = do not iterate.
;     kill = number of iterations eliminating bad limb points.  The 
;            default is 0 for non x-ray images and 2 for x-ray images.
;            X-ray images often need kill of 1 or 2.
;            Smaller numbers (0 or 1) make the algorithm 
;            correspondingly less sensitive to the initial guess.
;            Higher numbers will eliminate more limb points with the potential
;            for a better fit in some cases.
;     /norfit = if set, do not fit the radius of the sun: use the initial guess
;               If you know the pixel size exactly and pass the date or index
;               so that the sun's radius can be computed, this keyword should 
;               be used since it removes one degree of freedom from the circle
;               fit.
;     date = date for computation of solar radius.  Index is used if it is 
;            present.  If date is omitted, a crude value is computed from the
;            image.
;OUTPUTS:
;     x = x position of the center of the image.
;	  For SXT images, the units are FULL res pixels unless /SUMMED_PIXELS
;	  is set.  There is also (0, 1 or 3) pixel offset for FR, HR, and QR.
;     y = y position of the center of the image 
;	  For SXT images, the units are FULL res pixels unless /SUMMED_PIXELS
;	  is set.  There is also (0, 1 or 3) pixel offset for FR, HR, and QR.
;     r = radius of the image in units of (original) pixels
;
;     RSun/RFit = The ratio of the Sun's radius computed for the time of the
;                 image to the fit radius.  This is the size of the SXT pixels
;                 in arcseconds.  However, if the index in NOT suppled, the 
;                 Sun's radius can not be computed and the result is the 
;                 ratio of the initial guess of the radius to the fitted value.
;                 In this case, this is NOT the size of the SXT pixels.
;                 RSun/RFit is printed as the program ends.  It is not returned
;                 in a variable.
;SIDE EFFECTS:
;     Makes plots in the current window.
;RESTRICTIONS:
;     o The algorithm relies on a VERY good initital guess to remove bad limb
;       points.
;     o The correction for oblateness is not done for x-ray images.
;     o If part of the limb is obscurred, you may need to use a large kill 
;       value (around 5).  Such data must be treated on a case by case basis
;       by a smart user!
;PROCEDURE:
;     The initial
;     guess is used to find a rough limb position which is used to compute
;     the points with the maximum derivative in the intensity along the radii. 
;     These maximum derivative points are assumed to delineate the limb and a 
;     circle is fit to these points and the center and the radius are computed.
;     The derivative is smoothed using the deriv_lud procedure.
;MODIFICATION HISTORY:
;     T. Metcalf 4/1992
;     T. Metcalf 10/1992:  Switched from gt_hxa to hxa_suncenter to get the 
;                          initital guess.
;     T. Metcalf 2/1993:   Switched from hxa_suncenter to get_suncenter
;			  *** COPIED FROM SXT_CENTER.PRO TO FIT_LIMB.PRO ***
;	21-Dec-93 (MDM) - Modified to work with non-SXT images
;	 6-Jan-94 (MDM) - Minor changes
;	14-Jan-94 (MDM) - Removed FID_FLAG keyword (was marking if RD_PNT
;			  data had fidutial mark)
;	19-Jan-94 (MDM) - Corrected error not allowing it to work on SXT image
;			- Modified header.
;	 3-Feb-94 (MDM) - Modified to use GT_CORNER with /FROM_SC switch
;			  to determine the initial sun center guess.
;     T. Metcalf 2/1994:   o Default kill for non x-ray images is now 0.
;                          o Added obcorrect keyword.
;                          o Added ellipse keyword
;	15-Feb-94 (MDM) - Removed on_error,2 statement
;			- Corrected error where IF statement was trying to
;			  operate on a one byte array, not a scalar
;			  (the problem was in GET_RB0P)
;	21-Feb-94 (MDM) - Modified the header
;	 9-Dec-94 (MDM) - MOdified the header
;	22-Mar-95 (MDM) - Modified to add a 0.5 pixel offset for
;			  the results when not passing an SXT index.
;-

;on_error,2

CIRCLE_ITER = 50   ; iteration limit for circle fitting

x0 = 0.
y0 = 0.
r0 = 0.
r_err = 0.
ob_angle = 0.
oblateness = 0.
bias = 0.

verbose = keyword_set(verbose)
if keyword_set(noiterate) then noiterate=1 else noiterate=0
;;if keyword_set(iterate) then noiterate=0 else noiterate=1
if keyword_set(nordpnt) then begin
      if n_elements(initial) NE 3 then interactive=1
endif
qindex = keyword_set(index)

save_P_multi = !P.MULTI
!P.MULTI = [0,3,2,0,0]

;---------- Get the image

if (size(image))(0) NE 2 then message,'The image must be a 2-D array'
if qindex then begin
   if n_elements(index) NE 1 then $
      message,'Index must have a single element matching the image'
endif
nx = n_elements(image(*,0))
ny = n_elements(image(0,*))
if keyword_set(decompress) and qindex then image = sxt_decomp(image)

isxray = 0
if (qindex) then begin
    isxray=(gt_filtb(index) GT 1)
    pixel_size = gt_pix_size(index)	;pixel size in arcsec
end else begin
    if (n_elements(pixel_size) eq 0) then begin
	print, 'FIT_LIMB: Need the pixel size'
	return
    end
end

max_value = max(image)		;?? Why did tom set to 255 or 4096

;-------------------- Determine the intial values

if (n_elements(initial) EQ 3) then begin	;user passed them in
    xcen = initial(0)
    ycen = initial(1)
    radius = initial(2)
endif else begin
    if (qindex) then begin
	; ---- Compute the solar radius from the ephemeris
	radius = get_rb0p(index,/radius)
	;;radius = rb0p(0)	;in arcsec

	;;cen = sxt_cen(index)		;MDM removed 3-Feb-94
        ;;xcen = cen(0) / (2^gt_res(index))
        ;;ycen = cen(1) / (2^gt_res(index))

        resolution = 2^gt_res(index)

	corner = gt_corner(index, /from_sc)
	xcen = -corner(0)
	ycen = -corner(1)

        hxa_x = xcen*resolution
        hxa_y = ycen*resolution
    endif else begin
        pp = where(image GT 0.3*max_value)
        img = bytarr(nx,ny)
        img(pp) = 1
        centroidw,img,xcen,ycen
        if keyword_set(date) then radius = get_rb0p(date,/radius) $
        else radius = sqrt(total(img(pp))/!pi)*pixel_size
    end
end

init_radius = radius
radius = radius/pixel_size		;Convert arcsec to pixels

;-------------------- Now get a guess for the center of the circle from the user

	if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then begin
           erase
	   tv,bytscl(image)
	   xyouts,.15,.95,"Initial Guess",align=0.5,/norm,charsize=1.5
	endif

	if keyword_set(interactive) then begin
	    print
	    print,"Guess Disk Center as Accurately as Possible:"
	    print
	    print,"Click left mouse button repeatedly to get the center of the image"
	    print,"Click middle mouse button to redraw"
	    print,"Click right mouse button to continue when you get close"
	endif

	theta = findgen(361) * !dtor
	!ERR = 0
	REPEAT begin
	    if !ERR LT 4 then begin
		if !ERR GE 2 then tv,bytscl(image) else if !ERR EQ 1 then begin
		    xcen = x
		    ycen = y
		endif
		xcir = xcen + radius * cos(theta)
		ycir = ycen + radius * sin(theta)
		if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then plots,xcir,ycir,/device
	    endif
	    if keyword_set(interactive) then cursor,x,y,3,/device
	endrep UNTIL (!ERR GE 4) OR (NOT keyword_set(interactive))

	if keyword_set(interactive) then erase

	if NOT keyword_set(quiet) then begin
	    print
	    print,"Working ..."
	    print
	endif

	;-------------------- Start the work of fitting the limb

	; Number of points on either side of limb to check for the inflection 
	; point in intensity.  Increase if there was a fiducial mark.

   n_points = 10

   ; Compute the positions of the inflection points across the limb

   for t=0,n_elements(theta)-1 do begin  ;Loop through angles, 1 degree per step
      inten = fltarr(n_points*2+1)  ; inten storesthe intensity along a radius
      off_screen = 0
      for r=-n_points,n_points do begin  ; Get the intensity along a radius
         x = fix(xcen + (radius+r)*cos(theta(t)))
         y = fix(ycen + (radius+r)*sin(theta(t)))
         if NOT (x GT (nx-1) OR x LT 0 OR y GT (ny-1) or y LT 0) then begin
            inten(r+n_points) = image(fix(x),fix(y))
         endif $
         else off_screen=1
      endfor 
      if min(inten) EQ max(inten) then off_screen=1 ; Watch for saturated areas
      ; Now we have the intensity along the radius, so compute the point with
      ; the maximum derivative and store its position in the xfit and yfit 
      ; vectors.  
      if NOT off_screen then begin
         rinten = radius-n_points+findgen(n_elements(inten))  ; vector of radii
         if NOT keyword_set(fast) then begin
            ; Compute the *smoothed* derivative
            deriv_lud,rinten,inten,derivative,mid=rmidpoint,reg_param=1.0,/quiet
         endif $
         else begin
            ; Compute the *unsmoothed* derivative
            derivative = shift(inten,1)-inten  
            derivative(n_elements(derivative)-1) = 0.0 & derivative(0)=0.0
            rmidpoint = 0.5*(shift(rinten,1)+rinten)
            rmidpoint(n_elements(rmidpoint)-1)=rinten(n_elements(rinten)-1)
            rmidpoint(0)=rinten(0)
         endelse
         test_function = abs(derivative)
         dmax = max(test_function,inflection_p)  ; Find maximum derivative
         ; Now interpolate to the best value of test_function max
         ; put inflection_p in range
         inflection_p = ((inflection_p > 3) < (n_elements(rmidpoint)-4))
         contrast = abs(inten(inflection_p+3) - inten(inflection_p-3))
         ; line_cent does a parabolic fit to test_function to find the max
         inflection_p = line_cent(abs(test_function),best=inflection_p,n_p=3)
         ; put inflection_p in range
         inflection_p = ((inflection_p > 0) < (n_elements(rmidpoint)-1))
         ; Linear interpolation to get the limb radius
         inflection_p = interpol(rmidpoint,findgen(n_elements(rmidpoint)), $
                                 inflection_p)
         inflection_p = inflection_p(0) ; convert 1-element array to scalar
         if keyword_set(verbose) then begin
            plot,rinten,inten
            oplot,rmidpoint,test_function,lines=2
            oplot,[inflection_p,inflection_p],[-1.e6,1.e6]
         endif
         ; If the maximum derivative is at an endpoint, ignore it,
         ; otherwise add the point to the xfit and yfit vectors.
         ; Also test the significance of the result, ignore if it's not
         ; over 3 sigma or if the contrast over the limb is less than a 
         ; threshold value
         if (inflection_p GT rmidpoint(2)) AND $
            (inflection_p LT rmidpoint(n_elements(rmidpoint)-3)) $
            AND (dmax GE 2.5*stdev(test_function)) AND $
            (contrast GT MAX_VALUE*0.05) AND (contrast LE MAX_VALUE*0.84) $
            AND (max(inten) LT MAX_VALUE*0.86) $
         then begin
            xf = xcen + inflection_p*cos(theta(t))
            yf = ycen + inflection_p*sin(theta(t))
            if n_elements(xfit) LE 0 then xfit = [xf] $
            else xfit=[xfit,xf]
            if n_elements(yfit) LE 0 then yfit = [yf] $
            else yfit=[yfit,yf]
         endif $
         else if keyword_set(verbose) then $
            xyouts,0.75*(!X.CRANGE(1)-!X.CRANGE(0))+!X.CRANGE(0), $
            0.9*(!Y.CRANGE(1)-!Y.CRANGE(0))+!Y.CRANGE(0),align=0.5, $
            'REJECTED',/data
      endif
   endfor

   if n_elements(xfit) LT 3 then begin
      print,"ERROR: fit_limb: Could not get enough limb points" 
      x0 = 0.
      y0 = 0.
      r0 = 0.
      r_err = 0.
      ob_angle = 0.
      oblateness = 0.
      bias = 0.
   endif $
   else begin

      ; Make a cut on the distance from the mean radius.  This assumes that
      ; xcen and ycen are good guesses to the center of the disk!

      ; fit a circle to the limb points
      txfit = xfit  ; temporary variables
      tyfit = yfit
      txcen = xcen
      tycen = ycen
      tradius = radius
      ; Number of iterations killing bad points
      if n_elements(RLOOP) NE 1 then begin
         if (isxray) then RLOOP=2 else RLOOP=0
      endif $
      else RLOOP=fix(abs(RLOOP))
      for loop=0,RLOOP-1 do begin
         if n_elements(xfit) GE 3 then begin
            if keyword_set(ellipse) then begin
               circle_fit,txfit,tyfit, txcen, tycen, tradius, xcrude, ycrude, $
                          rcrude,check,/quiet,/noiterate,norfit=norfit
            endif else begin
               ; tolerance of 0.0 means do not iterate
               fit = fit_circle(txfit,tyfit,[txcen,tycen,tradius], $
                                tolerance=0.,radius_fix=norfit, $
                                limit_iter=CIRCLE_ITER,num_iter=num_iter)
               check = num_iter NE CIRCLE_ITER
               xcrude=fit(0) & ycrude=fit(1) & rcrude=fit(2)
            endelse
            if check then begin
               rfit = sqrt((xfit-xcrude)^2+(yfit-ycrude)^2)
               trfit = sqrt((txfit-xcrude)^2+(tyfit-ycrude)^2)
            endif $
            else begin
               rfit = sqrt((xfit-txcen)^2+(yfit-tycen)^2)
               trfit = sqrt((txfit-txcen)^2+(tyfit-tycen)^2)
               rcrude = total(rfit)/n_elements(rfit)
            endelse
            r_cutoff = (stdev(trfit)>0.25)
            good = where( abs(rfit-rcrude) LT r_cutoff,count)
            if count GT 0 then begin
               txfit = xfit(good)
               tyfit = yfit(good)
            endif
            txcen = xcrude
            tycen = ycrude
            tradius = rcrude
         endif
      endfor
      xfit = txfit
      yfit = tyfit
      xcen = txcen
      ycen = tycen
      radius = tradius
   
      if n_elements(xfit) LT 3 then begin
         print,"ERROR: fit_limb: Could not get enough limb points" 
         x0 = 0.
         y0 = 0.
         r0 = 0.
         r_err = 0.
         ob_angle = 0.
         oblateness = 0.
         bias = 0.
      endif $
      else begin
   
         if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then begin
            erase
            tv,bytscl(image,top=!D.N_COLORS-20)
            xyouts,.15,.95,"Limb Points",align=0.5,/norm,charsize=1.5
            plots,xfit,yfit,psym=1,/device   ; Show the limb points to the user
            wait,2.
         endif
      
         if keyword_set(interactive) then begin
            print
            print
            print,"Click left button near obviously bad limb points to erase them"
            print,"Click right button when finished to continue with the fit"
            REPEAT begin
               cursor,x,y,3,/device    ; Give the user a chance to edit the points
               if !ERR LT 4 then begin
                  distance = (xfit-x)^2+(yfit-y)^2
                  junk = min(distance,closest_p) ; find the closest xfit,yfit point
                  if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then $
                     plots,xfit(closest_p),yfit(closest_p),/device,psym=1,color=0
                  if closest_p NE 0 and closest_p NE n_elements(xfit)-1 then begin
                     xfit = $
                        [xfit(0:closest_p-1),xfit(closest_p+1:n_elements(xfit)-1)]
                     yfit = $
                        [yfit(0:closest_p-1),yfit(closest_p+1:n_elements(yfit)-1)]
                  endif $
                  else if closest_p EQ 0 then begin
                     xfit = xfit(1:n_elements(xfit)-1)
                     yfit = yfit(1:n_elements(yfit)-1)
                  endif $
                  else if closest_p EQ n_elements(xfit)-1 then begin
                     xfit = xfit(0:n_elements(xfit)-2)
                     yfit = yfit(0:n_elements(yfit)-2)
                  endif
               endif
            endrep UNTIL !ERR GE 4
         endif
   
         if keyword_set(interactive) then begin
            erase
            if NOT keyword_set(quiet) then begin
               print
               print,"Working ..."
               print
            endif
         endif
      
         ; fit a circle to the limb points
         !P.MULTI = save_P_multi
         if keyword_set(ellipse) then begin
            circle_fit, xfit, yfit, xcen, ycen, radius, x0, y0, r0, check, $
                        noiterate=noiterate,norfit=norfit
         endif else begin
            if keyword_set(noiterate) then tol = 0.0 else tol = 0.001
            fit = fit_circle(xfit,yfit,[xcen,ycen,radius],radius_fix=norfit, $
                             tolerance=tol,limit_iter=CIRCLE_ITER, $
                             num_iter=num_iter)
            check = (num_iter NE CIRCLE_ITER)
            x0 = fit(0)
            y0 = fit(1)
            r0 = fit(2)
         endelse
   
         if NOT check then begin
            print,"ERROR: fit_limb: No convergence in circle fit" 
            x0 = 0.
            y0 = 0.
            r0 = 0.
            r_err = 0.
            ob_angle = 0.
            oblateness = 0.
            bias = 0.
         endif $
         else begin
      
            if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then begin
               erase
               tv,bytscl(image)
               xyouts,.15,.95,"Best Fit",align=0.5,/norm,charsize=1.5
            endif
            x = x0 + r0 * cos(theta)
            y = y0 + r0 * sin(theta)
            if (NOT keyword_set(quiet)) OR (keyword_set(interactive)) then $
               plots,x,y,/device   ; display the best fit circle on the image
   
            ; Correct for resolution pixel differences  (stolen from find_limb)
            ; Also correct for oblateness of the disk
   
            ;
            ; Analyze limb shape by fitting harmonics (1,2,3,4 * theta).
            ;
            ; Stolen from find_limb
   
            rad = sqrt((xfit-x0)^2 +(yfit-y0)^2) - r0
            theta = atan((yfit-y0),(xfit-x0))
            sin_fun = sin(theta)
            cos_fun = cos(theta)
            sin_1 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
            cos_1 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
            r_err = sigma/sqrt(n_elements(rad))
            ;print, sigma, ' Residuals after fundamental'
            rad = rad - (sin_1(0) + cos_1(0))/2 - sin_1(1)*sin(theta)- cos_1(1)*cos(theta)
            sin_fun = sin(2*theta)
            cos_fun = cos(2*theta)
            sin_2 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
            cos_2 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
            ;print, sigma, ' Residuals after 1st harmonic'
            rad = rad-(sin_2(0)+cos_2(0))/2 - sin_2(1)*sin(2*theta) - cos_2(1)*cos(2*theta)
            sin_fun = sin(3*theta)
            cos_fun = cos(3*theta)
            sin_3 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
            cos_3 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
            ;print, sigma, ' Residuals after 2nd harmonic'
            rad = rad-(sin_3(0)+cos_3(0))/2 - sin_3(1)*sin(3*theta) - cos_3(1)*cos(3*theta)
            sin_fun = sin(4*theta)
            cos_fun = cos(4*theta)
            sin_4 = poly_fit(sin_fun, rad, 1, sin_fit, yband, sigma, a)
            cos_4 = poly_fit(cos_fun, rad, 1, cos_fit, yband, sigma, a)
            ;print, sigma, ' Residuals after 3rd harmonic'
            rad = rad-(sin_4(0)+cos_4(0))/2 - sin_4(1)*sin(4*theta) - cos_4(1)*cos(4*theta)
   
            ob_angle = atan(sin_2(1)/cos_2(1))
            oblateness = sqrt(sin_2(1)^2 + cos_2(1)^2)
            bias = sqrt(sin_3(1)^2 + cos_3(1)^2 + sin_4(1)^2 + cos_4(1)^2)

            ; No oblateness correction for X-rays
            if (NOT keyword_set(noobcorrect)) and (NOT isxray) then begin
               x0 = x0 + sin_1(1) 
               y0 = y0 + cos_1(1)
            endif 

	    if (not keyword_set(summed_pixels) and qindex) then begin
		off_corr = [-1, 0, 1, -1, 3]            ;MDM added 19-Nov-91
		;changing from 1x1 to 1x1,no offset
		;changing from 2x2 to 1x1,offset=1 (1x1 pixels 1&2 are 2x2 pixel 0)
		;changing from 4x4 to 1x1,offset=3 (1x1 pixels 3,4,5&6 are 4x4 pixel 0)
   
		x0 = x0 - 0.5   ; Move the image such that pixel center is defined
		y0 = y0 - 0.5   ; as the pixel

		; -- Output in FULL RESOLUTION pixels
		x0 = x0 * resolution + off_corr(resolution)
		y0 = y0 * resolution
		r0 = r0 * resolution
	    end else begin	; MDM added the ELSE statement 20-Mar-95
		x0 = x0 - 0.5   ; Move the image such that pixel center is defined
		y0 = y0 - 0.5   ; as the pixel
	    end

            print,x0,y0,r0,init_radius/r0, $
               format="('Center: x=',f10.4,' y=',f10.4,' Radius: ',f10.4,' RSun/RFit: ',f10.4)"
         endelse
      endelse
   endelse
   print
   !P.MULTI = save_P_multi

end