;+ IDL PROCEDURE
;
; NAME		ftscan.pro
;
; PURPOSE	Perform an azimuthal scan using a FITS image [cartesian coord].
;
; SYNTAX	ftscan, fits_file, radius, thmin, thmax, thinc
;
; PARAMETERS	REQUIRED:
;		fits_file	name of FITS file.
;		radius		radius value [Rsun units]
;		thmin		beginning angle [degrees]
;		thmax		ending    angle [degrees]
;		thinc		angle increment [degrees]
;
;		OPTIONAL:
;		ymin		Y-axis minimum value
;		ymax		Y-axis maximum value
;
; EXT.ROUTINES	tscan.pro	Performs azimuthal (theta) scan.
; 		lct.pro		loads a color table from an ASCII file.
;		rcoord.pro	converts between [r,th] and [x,y] coordinates.
;
; NOTE 		This procedure works best in private colormap mode (256 levels).
;
; AUTHOR	Andrew L. Stanger   HAO/NCAR   08 Feb 2001
;
; HISTORY	7 Jun 2001 Adapt for SOHO/LASCO C2 (French version).
;-

PRO ftscan, fits_file, radius, thmin, thmax, thinc, ymin=ymin, ymax=ymax

SET_PLOT, 'X'

;--- Load color table.

stars = '***'
date_obs = ''
time_obs = ''
lct, '/home/stanger/color/bwcp.lut'
lct, '/home/stanger/color/bwy.lut'

red   = bytarr (256)
green = bytarr (256)
blue  = bytarr (256)

tvlct, red, green, blue, /get		; Read color table into arrays.

;--- Read FITS image into memory.

img = readfits (fits_file, hdu)

imin = min (img, max=imax)
print, 'imin/imax: ', imin, imax

;--- Get parameters from FITS header.

telescop = fxpar (hdu, 'TELESCOP')	; Telescope name
instrume = fxpar (hdu, 'INSTRUME')	; Instrument name
detector = fxpar (hdu, 'DETECTOR')	; Detector name
date_obs = fxpar (hdu, 'DATE-OBS')	; observation date
time_obs = fxpar (hdu, 'TIME-OBS')	; observation time
xdim     = fxpar (hdu, 'NAXIS1')	; X dimension
ydim     = fxpar (hdu, 'NAXIS2')	; Y dimension
xcen     = fxpar (hdu, 'CRPIX1')	; X center
ycen     = fxpar (hdu, 'CRPIX2')	; Y center
type_obs = fxpar (hdu, 'TYPE-OBS')	; Observation type (e.g., pb)
bunit    = fxpar (hdu, 'BUNIT')		; Brightness unit  (e.g., Bsun)
roll     = fxpar (hdu, 'CROTA1')	; rotation angle (degrees).
crradius = fxpar (hdu, 'CRRADIUS')	; pixels/Rsun
pixrs    = crradius

IF (xcen LT 0.0) THEN xcen = (xdim - 1) / 2.0
IF (ycen LT 0.0) THEN ycen = (ydim - 1) / 2.0

telescop = STRLOWCASE (STRCOMPRESS (telescop, /remove_all))
instrume = STRLOWCASE (STRCOMPRESS (instrume, /remove_all))
detector = STRLOWCASE (STRCOMPRESS (detector, /remove_all))

;--- LASCO C2 uses different keywords.

IF (detector EQ 'c2') THEN	$	
BEGIN
   PRINT, 'detector: ', detector
   date_obs = fxpar (hdu, 'DATE_OBS')
   time_obs = fxpar (hdu, 'TIME_OBS')
   xcen  = fxpar (hdu, 'XSUN')
   ycen  = fxpar (hdu, 'YSUN')
   roll  = 0.0
   crradius = fxpar (hdu, 'RSUN_PIX')
   pixrs = crradius
   img = img * 1.0e-10			; French B/Bsun scaling factor.
   type_obs = ''
   bunit    = 'Pixel Intensity'
END

ylab = type_obs + ' [' + bunit + ']'
PRINT, 'date_obs: ', date_obs
PRINT, 'time_obs: ', time_obs
PRINT, 'xdim/ydim: ', xdim, ydim
PRINT, 'xcen/ycen: ', xcen, ycen
PRINT, 'crradius: ', crradius

;--- Convert numerical values to strings.

srad   = STRING (radius, FORMAT='(F5.2)')
srad   = STRTRIM (srad, 2)
sthmin = STRING (thmin, FORMAT='(F7.2)')
sthmin = STRTRIM (sthmin, 2)
sthmax = STRING (thmax, FORMAT='(F7.2)')
sthmax = STRTRIM (sthmax, 2)

;--- Do theta scan.

tscan, namimg, img, crradius, roll, xcen, ycen,	$
       thmin, thmax, thinc, radius,			$
       scan, scandx, ns

;FOR i = 0, ns - 1 DO		$
;BEGIN
;    print, 'scandx: ', scandx [i], '  scan: ', scan [i]
;END

;--- Reduce image size for display (if needed).

sizimg = SIZE (img)

WHILE (sizimg [1] GT 1000 OR sizimg [2] GT 1000) DO	$
    BEGIN
    img = REBIN (img, xdim/2, ydim/2)
    xdim = xdim / 2
    ydim = ydim / 2
    xcen = xcen / 2.0
    ycen = ycen / 2.0
    pixrs = pixrs / 2.0
    sizimg = SIZE (img)
    END

;--- Display image.

WINDOW, xsize=xdim, ysize=ydim
imin = MIN (img, max=imax)
imgb = BYTSCL (img, min=imin, max=imax, top=249)
TV, imgb

;--- Label image.

XYOUTS,        2, ydim-10, fits_file, /device, color=250, charsize=0.6
XYOUTS,        2, ydim-25, telescop,  /device, color=250, charsize=0.6
XYOUTS, xdim-110, ydim-10, date_obs,  /device, color=250, charsize=0.6
XYOUTS, xdim-110, ydim-25, time_obs + ' UT',  /device, color=250, charsize=0.6
XYOUTS,        2,  20, srad + ' Rsun', /device, color=250, charsize=0.6
XYOUTS,        2,   5, sthmin + ' - ' + sthmax + ' deg.', 	$
	/device, color=250, charsize=0.6

;--- Plot scan on image.

th = thmin - thinc
FOR i = 0, ns - 1 DO	$
   BEGIN ;{
   th = th + thinc
   ierr = rcoord (radius, th, x, y, 1, roll, xcen, ycen, pixrs)
   ixg = FIX (x + 0.5)
   iyg = FIX (y + 0.5)
   PLOTS, [ixg, ixg], [iyg, iyg], /device, color=254
   END ;}

;--- Draw a dotted circle (10 degree increments) at 1.0 Rsun.

th = -10.0
r  =   1.0
FOR i = 0, 360, 10 DO	$
   BEGIN
   th = th + 10.0
   ierr = rcoord (r, th, x, y, 1, roll, xcen, ycen, pixrs)
   ixg = FIX (x + 0.5)
   iyg = FIX (y + 0.5)
   PLOTS, [ixg, ixg], [iyg, iyg], /device, color=251
   END

;--- Draw a dotted line (0.2 Rsun increments) at   0.0 degrees.

r  =  -0.2
th =   0.0
FOR i = 0, 5 DO		$
   BEGIN
   r = r + 0.2
   ierr = rcoord (r, th, x, y, 1, roll, xcen, ycen, pixrs)
   ixg = FIX (x + 0.5)
   iyg = FIX (y + 0.5)
   PLOTS, [ixg, ixg], [iyg, iyg], /device, color=251
   END

;--- Draw a dotted line (0.2 Rsun increments) at  90.0 degrees.

r  =  -0.2
th =  90.0
FOR i = 0, 5 DO		$
   BEGIN
   r = r + 0.2
   ierr = rcoord (r, th, x, y, 1, roll, xcen, ycen, pixrs)
   ixg = FIX (x + 0.5)
   iyg = FIX (y + 0.5)
   PLOTS, [ixg, ixg], [iyg, iyg], /device, color=251
   END

;--- Draw a dotted line (0.2 Rsun increments) at 180.0 degrees.

r  =  -0.2
th = 180.0
FOR i = 0, 5 DO		$
   BEGIN
   r = r + 0.2
   ierr = rcoord (r, th, x, y, 1, roll, xcen, ycen, pixrs)
   ixg = FIX (x + 0.5)
   iyg = FIX (y + 0.5)
   PLOTS, [ixg, ixg], [iyg, iyg], /device, color=251
   END

;--- Draw a dotted line (0.2 Rsun increments) at 270.0 degrees.

r  =  -0.2
th = 270.0
FOR i = 0, 5 DO		$
   BEGIN
   r = r + 0.2
   ierr = rcoord (r, th, x, y, 1, roll, xcen, ycen, pixrs)
   ixg = FIX (x + 0.5)
   iyg = FIX (y + 0.5)
   PLOTS, [ixg, ixg], [iyg, iyg], /device, color=251
   END

;--- Read displayed image into 2D array.

imgnew = TVRD ()
extpos = STRPOS (fits_file, ".fts")
gif_name = STRMID (fits_file, 0, extpos) + '.r' + srad + '.gif'
PRINT, 'gif_name: ', gif_name

;--- Write GIF image to disk.

WRITE_GIF, gif_name, imgnew, red, green, blue

;--- Plot theta scan to a postscript file.

ps_name  = STRMID (fits_file, 0, extpos) + '.r' + srad + '.ps'
PRINT, 'plot name: ', ps_name
SET_PLOT, 'PS'
DEVICE, filename=ps_name


yminset = 0
ymaxset = 0
IF (KEYWORD_SET (ymin)) THEN yminset = 1  ELSE ymin = 0.0
IF (KEYWORD_SET (ymax)) THEN ymaxset = 1  ELSE ymax = 0.0

print, 'ymin/ymax: ', ymin, ymax

IF (yminset OR ymaxset) THEN			$
BEGIN
print, 'yrange: ', ymin, ymax
PLOT, scandx, scan,				$
      title = 'Theta Scan ' + fits_file + ' ' + srad + ' Rsun',		$
      xtitle = 'Position Angle',		$
      ytitle = ylab,		$
      yrange = [ymin, ymax]

END						$
ELSE						$
PLOT, scandx, scan,				$
      title = fits_file + ' Theta Scan @' + srad + ' Rsun',	$
      xtitle = 'Position Angle',		$
      ytitle = ylab 

DEVICE, /close
SET_PLOT, 'X'

END