;+ IDL PROCEDURE
;
; NAME		ftscan_rpb_ne.pro
;
; PURPOSE	Perform a Mk4 rpb azimuthal scan
;		using a FITS image [cartesian coord].
;               Compute the electron density from the pB intensity.
;		Generates azimuthal scan plots for pB [Bsun] and Ne/cm2.
;
; SYNTAX	ftscan_rpb_ne, 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.
;		nedpb.pro	computes Ne/pb[Bsun]*cm2
;
; OUTPUTS	GIF file showing the position of the azimuthal scan.
;		POSTSCRIPT file containing the azimuthal scan plots.
;		TEXT file containing the scan values [ASCII].
;
; NOTE 		This procedure works best in private colormap mode (256 levels).
;
; AUTHOR	Andrew L. Stanger   HAO/NCAR   18 Jun 2004
;-

PRO ftscan_rpb_ne, fits_file, radius, thmin, thmax, thinc, $
		   ymin=ymin, ymax=ymax

text = 1
SET_PLOT, 'X'

;--- Load color table.

stars = '***'
date_obs = ''
time_obs = ''
bscale  =  1.0e-10
bzero   = -1.0e-7

;--- dispmin & dispmax are set explicitly, since "readfits" fails to
;    read these values properly in the FITS header of Mk4 images,
;    apparently because they do NOT contain a decimal point (-5e-08, 5e-07).

dispmin = -5.0e-08
dispmax =  5.0e-07

PRINT, 'dispmin/dispmax: ', dispmin, dispmax
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.

ftspos = STRPOS (fits_file, '.fts')
basename = STRMID (fits_file, 0, ftspos)
;PRINT, 'basename: ', basename

;--- Read FITS image into memory.

;img = readfits (fits_file, hdu, /noscale)
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
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
type_obs = fxpar (hdu, 'TYPE-OBS')	; Observation type (e.g., pb)
roll     = fxpar (hdu, 'CROTA1')	; rotation angle (degrees).
crradius = fxpar (hdu, 'CRRADIUS')	; pixels/Rsun
pixrs    = crradius
bunit    = fxpar (hdu, 'BUNIT')		; Brightness unit  (e.g., Bsun)

bscale   = fxpar (hdu, 'BSCALE')	; physical = data * bscale + bzero
bzero    = fxpar (hdu, 'BZERO')
datamin  = fxpar (hdu, 'DATAMIN')	; Data    minimum intensity.
datamax  = fxpar (hdu, 'DATAMAX')	; Data    maximum intensity.
;dispmin  = fxpar (hdu, 'DISPMIN')	; Display minimum intensity.
;dispmax  = fxpar (hdu, 'DISPMAX')	; Display maximum intensity.

PRINT, 'date_obs: ', date_obs
PRINT, 'time_obs: ', time_obs
PRINT, 'xdim/ydim: ', xdim, ydim
PRINT, 'xcen/ycen: ', xcen, ycen
PRINT, 'crradius: ', crradius
;PRINT, 'bscale/bzero: ', bscale, bzero
;PRINT, 'datamin/datamax: ', datamin, datamax
;PRINT, 'dispmin/dispmax: ', dispmin, dispmax

;dispmin = (dispmin - bzero) / bscale
;dispmax = (dispmax - bzero) / bscale
;dispmin = dispmin * bscale + bzero
;dispmax = dispmax * bscale + bzero
;PRINT, 'dispmin/dispmax: ', dispmin, dispmax

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

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

;--- Verify that image is NOT vignetted.

type_obs = STRTRIM (type_obs, 2)
vigpos = -1
vigpos = STRPOS (type_obs, 'VIG')
IF (vigpos NE -1) THEN		$
BEGIN
  PRINT, 'type_obs: ', type_obs
  PRINT, 'This program will NOT produce valid results for VIGNETTED data."
  RETURN
END

ylab = type_obs + ' [' + bunit + ']'

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

nescan = fltarr (ns)

;--- Compute the electron density/pB[Bsun] at the scan height.

necm2 = nedpb (radius)
PRINT, radius, ' Rsun  ', necm2, ' Ne/cm2'

;--- Convert pB [Bsun] to Ne/cm2.

FOR i = 0, ns - 1 DO	$
BEGIN
  nescan [i] = scan [i] * necm2
END

IF (KEYWORD_SET (text)) THEN	$
BEGIN
  pfile = basename + '.tscan.ne.dat'
  PRINT, 'text file: ', pfile
  CLOSE,  11
  OPENW,  11, pfile
  PRINTF, 11, 'ftscan_rpb_ne.pro'
  PRINTF, 11, fits_file, '   Azimuthal scan @ ', srad, ' Rsun'

  FOR i = 0, ns - 1 DO		$
  BEGIN
    PRINTF, 11, scandx [i], ' degrees   ', scan [i], ' pB [B/Bsun]   ', nescan [i], ' Ne/cm2'
  END

  CLOSE, 11
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)
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, 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 + '.ne.gif'
PRINT, 'gif  file: ', 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 + '.ne.ps'
PRINT, 'plot file: ', 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]

PLOT, scandx, nescan,				$
      title = 'Theta Scan ' + fits_file + ' ' + srad + ' Rsun',		$
      xtitle = 'Position Angle',		$
      ytitle = 'Electron Density [Ne/cm2]',	$
      yrange = [ymin, ymax]

END						$
ELSE						$
BEGIN
  PLOT, scandx, scan,				$
      title = fits_file + ' Theta Scan @' + srad + ' Rsun',	$
      xtitle = 'Position Angle',		$
      ytitle = ylab,				$
      yrange = [0.0, max (scan)]

  PLOT, scandx, nescan,				$
      title = fits_file + ' Theta Scan @' + srad + ' Rsun',	$
      xtitle = 'Position Angle',		$
      ytitle = 'Electron Density [Ne/cm2]' ,	$
      yrange = [0.0, max (nescan)]
END

DEVICE, /close
SET_PLOT, 'X'

END