;+ IDL PROCEDURE
;
; NAME		ftscans.pro
;
; PURPOSE	Perform an azimuthal scan using a FITS image.
;
; SYNTAX	ftscan, fits_file, radius, thmin, thmax, thinc
;
; PARAMETERS	fits_file	name of FITS file.
;		radius		radius value [Rsun units]
;		thmin		beginning angle [degrees]
;		thmax		ending    angle [degrees]
;		thinc		angle increment [degrees]
;		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 ftscans, 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)

;--- 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
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.
   END

PRINT, 'date_obs: ', date_obs
PRINT, 'time_obs: ', time_obs

;--- 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

;--- 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 = 'Pixel Magnitude',		$
      yrange = [ymin, ymax]

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

DEVICE, /close
SET_PLOT, 'X'

END