;+ IDL PROCEDURE
;
; NAME		frscan_rpb.pro
;
; PURPOSE	Perform a Mk4 rpb radial scan using a FITS image [cartesian].
;
; SYNTAX	frscan_rpb, fits_file, angle, radmin, radmax, radinc
;
; PARAMETERS	fits_file	name of FITS file.
;		angle		position angle [degrees CCW from solar North]
;		radmin		beginning radius [Rsun]
;		radmax		ending    radius [Rsun]
;		radinc		angle increment  [Rsun]
; OPTIONAL
;		ymin		Y-axis minimum value
;		ymax		Y-axis maximum value
;
; EXT.ROUTINES	rscan.pro	Performs radial 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   03 Aug 2001
;-

PRO frscan_rpb, fits_file, angle, radmin, radmax, radinc, ymin=ymin, ymax=ymax

SET_PLOT, 'X'

;--- Load color table.

stars = '***'
date_obs = ''
time_obs = ''
dispmin = -5.0e-08
dispmax =  5.0e-07
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
dataform = fxpar (hdu, 'DATAFORM')	; data format
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

;dispmin  = fxpar (hdu, 'DISPMIN')
;dispmax  = fxpar (hdu, 'DISPMAX')

;--- Verify that dataform is CARTESIAN.

dataform = STRTRIM (dataform, 2)
rectpos = -1
cartpos = -1
rectpos = STRPOS (dataform, 'Rect')
cartpos = STRPOS (dataform, 'Cart')
PRINT, 'rectpos: ', rectpos
PRINT, 'cartpos: ', cartpos

IF (rectpos EQ -1 AND cartpos EQ -1)  THEN	$
BEGIN
   PRINT, 'DATAFORM: ', dataform, ' is incorrect.'
   PRINT, 'It should be RECT or CARTESIAN.'
   PRINT, 'ONLY Cartesian coordinate data is valid with this program.'
   RETURN
END

;IF (dataform NE '')  THEN	$
IF (cartpos NE 0) THEN		$
BEGIN
    PRINT, 'Must be cartesian coordinate image.'
    RETURN
END

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.


sangle  = STRING  (angle, FORMAT='(F7.2)')
sangle  = STRTRIM (sangle, 2)
sradmin = STRING  (radmin, FORMAT='(F7.2)')
sradmin = STRTRIM (radmin, 2)
sradmax = STRING  (radmax, FORMAT='(F7.2)')
sradmax = STRTRIM (sradmax, 2)

print, 'xdim/ydim: ', xdim, ydim

;--- Do radial scan.

rscan, namimg, img, crradius, roll, xcen, ycen,		$
       radmin, radmax, radinc, angle,			$
       scan, scandx, ns

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

;--- 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)
imgb = BYTSCL (img, min=dispmin, max=dispmax, 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, sangle + ' deg', /device, color=250, charsize=0.6
XYOUTS,        2,   5, sradmin + ' - ' + sradmax + ' Rsun', 	$
	/device, color=250, charsize=0.6

;--- Plot scan on image.

radius = radmin - radinc
FOR i = 0, ns - 1 DO	$
   BEGIN ;{
   radius = radius + radinc
   ierr = rcoord (radius, angle, 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) + '.pa' + sangle + '.gif'
;PRINT, 'gif_name: ', gif_name

;--- Write GIF image to disk.

;WRITE_GIF, gif_name, imgnew, red, green, blue

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

ps_name  = STRMID (fits_file, 0, extpos) + '.pa' + sangle + '.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 = fits_file + ' Radial Scan @' + sangle + ' degrees',	$
      xtitle = 'Radius [Rsun]',		$
      ytitle = 'Pixel Magnitude',		$
      yrange = [ymin, ymax]

END						$
ELSE						$
PLOT, scandx, scan,				$
      title = fits_file + ' Radial Scan @' + sangle + ' degrees',	$
      xtitle = 'Radius [Rsun]',		$
      ytitle = 'Pixel Magnitude'

DEVICE, /close
SET_PLOT, 'X'

END