;+
; Project     :	STEREO - SECCHI
;
; Name        :	COR1_QUICKPOL
;
; Purpose     :	Return COR1 total and polarized brightness, mu angle
;
; Category    :	STEREO, SECCHI, COR1, Analysis
;
; Explanation :	This routine performs the same function as COR_POLARIZ, but can
;               be up to 20 times faster.  It returns the total brightness,
;               polarized brightness, and angle of polarization (MU).  It can
;               also handle multiple sequences.
;
; Syntax      :	COR1_QUICKPOL, IMAGE, TOTB  [, PB  [, MU ]]
;
; Examples    :	CAT = COR1_PBSERIES('1-Apr-2007', 'Ahead')
;               SECCHI_PREP, CAT.FILENAME, HEADER, IMAGE
;               COR1_QUICKPOL, IMAGE, TOTB, PB, MU
;
; Inputs      :	IMAGE   = Array containing the pB images to process, with
;                         polarizer angles of 0,120,240.  This must have either
;                         the dimensions (Nx,Ny,3*Ns) or (Nx,Ny,3,Ns), where Ns
;                         is the number of sequences.
;
; Opt. Inputs :	None.
;
; Outputs     :	TOTB    = The total brightness image(s)
;
; Opt. Outputs:	PB      = The polarized brightness image(s)
;               MU      = The polarization angle image(s)
;
; Keywords    :	DOUBLE  = If set, then the output images are returned as double
;                         precision.  The calculation is performed in double
;                         precision, but returned as floating point unless
;                         IMAGE is double precision or /DOUBLE is passed.
;
;               HEADER  = Array containing the image header structures.
;                         Passing in the headers allows the MU values to be
;                         corrected for the use of other polarization
;                         sequences, such as [10,130,250].
;
;               TANGENTIAL = If set, then return only the tangential
;                            component of pB.  Ignored if HEADER not passed.
;                            The pB values are multiplied by the cosine of the
;                            angle to the tangential direction.  This tends to
;                            sppress noisy regions, while not significantly
;                            effecting other regions.
;
;               RADIAL  = If set, then return only the radial component of pB.
;                         The solar corona should not be polarized radially,
;                         but comet tails can be polarized radially or
;                         tangentially, depending on the phase angle.
;
;               ERRMSG  = If defined and passed, then any error messages will
;                         be returned to the user in this parameter rather than
;                         depending on the MESSAGE routine in IDL.  If no
;                         errors are encountered, then the input string is
;                         returned unchanged.  In order to use this feature,
;                         ERRMSG must be defined first, e.g.
;
;                               ERRMSG = ''
;                               COR1_QUICKPOL, ERRMSG=ERRMSG, ...
;                               IF ERRMSG NE '' THEN ...
;
; Calls       :	None.
;
; Common      :	None.
;
; Restrictions:	The calculation assumes that the images are presented with
;               polarizer angles of 0,120,240.  The calculations of TOTB and
;               PB will be correct so long as the three polarizer angles are
;               separated by either 60 or 120 degrees.  For cases other than
;               [0,120,240], the correct MU will be returned when the optional
;               HEADER keyword is passed.
;
; Side effects:	None.
;
; Prev. Hist. :	None.
;
; History     :	Version 1, 6-Apr-2007, William Thompson, GSFC
;               Version 2, 2-Jul-2007, WTT, make more memory friendly
;               Version 3, 25-Sep-2007, WTT, fix double precision bug
;               Version 4, 15-Feb-2008, WTT, added /TANGENTIAL keyword
;
; Contact     :	WTHOMPSON
;-
;
pro cor1_quickpol, image, totb, pb, mu, double=k_double, header=header, $
                   tangential=k_tangential, radial=k_radial, errmsg=errmsg
on_error, 2
;
;  Check the number of input parameters.
;
if n_params() lt 2 then begin
    message = 'Syntax: COR1_QUICKPOL, IMAGE, TOTB  [, PB  [, MU ]]'
    goto, handle_error
endif
;
tangential = keyword_set(k_tangential) and (n_elements(header) gt 0)
radial = keyword_set(k_radial) and (n_elements(header) gt 0)
;
;  Make sure that IMAGE has the right dimensions.
;
sz = size(image)
case sz[0] of
    3: begin
        if (sz[3] mod 3) ne 0 then begin
            message = 'Number of images must be divisible by 3'
            goto, handle_error
        endif
        save_dim = sz[1:3]
        nseq = sz[3]/3
        image = reform(image, sz[1], sz[2], 3, nseq, /overwrite)
    end
    4: begin
        if sz[3] ne 3 then begin
            message = 'Third IMAGE dimension must be 3'
            goto, handle_error
        endif
        nseq = sz[4]
        save_dim = sz[1:4]
    end
    else: begin
        message = 'IMAGE must have dimensions (Nx,Ny,3*Ns) or (Nx,Ny,3,Ns)'
        goto, handle_error
    end
endcase
;
;  If IMAGE is not double precision, then save the original image and make a
;  double precision version for processing.
;
test_dbl = datatype(image) eq 'DOU'
if test_dbl then begin
    delvarx, save_image
    a = temporary(image)
end else begin
    save_image = temporary(image)
    a = double(save_image)
endelse
;
;  Calculate the total brightness.  The 2/3 factor will be applied later.
;
totb = total(a,3)
;
;  If requested, calculate the polarized brightness.
;
if n_params() ge 3 then begin
    pb = totb
    for iseq = 0,nseq-1 do pb[*,*,iseq] = totb[*,*,iseq]^2 - $
      3.d0*reform(a[*,*,0,iseq]*a[*,*,1,iseq] + $
                  a[*,*,0,iseq]*a[*,*,2,iseq] + $
                  a[*,*,1,iseq]*a[*,*,2,iseq])
    pb = (4.d0/3.d0) * sqrt(0 > temporary(pb))
end else pb = 0                 ;Place-holder
;
;  Apply the 2/3 factor to total brightness.
;
totb = (2.d0/3.d0) * temporary(totb)
;
;  If requested, calculate the angle of polarization.
;
if (n_params() eq 4) or tangential or radial then begin
    w = where(pb gt 0, n0)
    if n0 eq 0 then mu = pb else begin
        pmin = min(pb[w])
        mu = pb
        for iseq = 0,nseq-1 do begin
            mu0 = acos(sqrt(0 > ((reform(a[*,*,0,iseq]) - $
                                  0.5d0*(totb[*,*,iseq]-pb[*,*,iseq])) / $
                                 (pb[*,*,iseq]>pmin))) < 1)
            w = where(a[*,*,2,iseq] lt a[*,*,1,iseq], count)
            if count gt 0 then mu0[w] = -mu0[w]
            mu[*,*,iseq] = mu0
        endfor
    endelse
;
;  If the header was passed, then modify MU based on the polarizer angles in
;  the header.
;
    if n_elements(header) gt 0 then begin
        mu0 = header[0].polar * !dpi / 180.d0
        halfpi = !dpi / 2.d0
        mu = temporary(mu) + mu0
        w = where(mu gt halfpi, count)
        if count gt 0 then mu[w] = mu[w] - !dpi
        w = where(mu lt -halfpi, count)
        if count gt 0 then mu[w] = mu[w] + !dpi
        dpolar = header[1].polar - header[0].polar
        if (dpolar eq 60) or (dpolar eq 240) or (dpolar eq -120) then $
          mu = -temporary(mu)
    endif
end else mu = 0                 ;Place-holder
;
;  Apply the /tangential or /radial keyword.
;
if tangential or radial then begin
    wcs = fitshead2wcs(header[0])
    cen = wcs_get_pixel(wcs, [0,0])
    x = dindgen(sz[1]) - cen[0]
    y = dindgen(sz[2]) - cen[1]
    x = rebin(reform(x,sz[1],1), sz[1],sz[2])
    y = rebin(reform(y,1,sz[2]), sz[1],sz[2])
    theta = atan(y,x)
    for iseq = 0,nseq-1 do begin
        dmu = mu[*,*,iseq] - theta
        corr = abs(cos(dmu))
        if radial then corr = abs(sin(dmu))
        pb[*,*,iseq] = pb[*,*,iseq] * corr
    endfor
endif
;
;  Restore the original IMAGE array.
;
if n_elements(save_image) gt 0 then image = temporary(save_image) else $
  image = temporary(a)
image = reform(image, save_dim, /overwrite)
;
;  If the original IMAGE was not double precision, then convert the output to
;  floating point.
;
if (not test_dbl) and (not keyword_set(k_double)) then begin
    totb = float(temporary(totb))
    pb   = float(temporary(pb))
    mu   = float(temporary(mu))
endif
return
;
;  Error handling point.
;
handle_error:
if n_elements(errmsg) eq 0 then message, message else errmsg=message
return
end