function cor_polariz,images,hdrs,hdr, pB=pB, mu=mu, PERCENT=percent, date_on=date_on,calibrate_off=calibrate_off, PIPELINE=pipeline, update_hdr_off=update_hdr_off,silent=silent,_extra=ex
;+
;$Id: cor_polariz.pro,v 1.24 2008/01/17 21:43:02 thompson Exp $
;
; Project   : STEREO SECCHI
;                   
; Name      : cor_polariz
;               
; Purpose   : This program takes three polarized images and calculates
;             the total brightness image.
;               
; Explanation: Uses the Mueller matrix elements returned by
;              cor_mueller.pro to calculate the total brightness image, the
;              polarized brightness and polarization angle.
;             
;               
; Use       : IDL> totB_image = cor_polariz(images,hdrs)
;    
; Inputs    : images - cube of 3 images from the same polarization sequence
;             hdrs - array of 3 SECCHI header structures
;
; Outputs   :  totB - total brightness image
;
; Keywords  :	pB - set to change output to polarized brightness
;		mu - set to change output to polarization angle
;		/PIPELINE - define pphdr
;               UPDATE_HDR_OFF- if set, don't call SCC_UPDATE_HDR (saves time)
;
; Common    : SECCHI_DAILY defined in daily/?.pro
;
; Calls     : SCC_FITSHDR2STRUCT
;               
; Category  : Polarization
;
; $Log: cor_polariz.pro,v $
; Revision 1.24  2008/01/17 21:43:02  thompson
; Use scc_update_history
;
; Revision 1.23  2007/10/15 17:48:12  nathan
; defined ex for not-secchi_prep case (Bug 240)
;
; Revision 1.22  2007/06/28 20:59:37  colaninn
; added on_error,2
;
; Revision 1.21  2007/06/15 21:47:25  colaninn
; corrected last update
;
; Revision 1.20  2007/06/15 21:36:22  colaninn
; added filename update
;
; Revision 1.19  2007/05/24 17:07:47  colaninn
; another typo error
;
; Revision 1.18  2007/05/23 14:21:29  thernis
; Add info in message concerning polarizers to avoid stop
;
; Revision 1.17  2007/05/23 14:05:57  colaninn
; corrected typo bug
;
; Revision 1.16  2007/05/21 20:16:06  colaninn
; added silent keyword
;
; Revision 1.15  2007/04/10 17:52:30  thernis
; Avoid program stop when obs_id is not the same for the image triplet
;
; Revision 1.14  2007/04/09 15:45:58  colaninn
; removed check of cadence keyword
;
; Revision 1.13  2007/03/30 22:11:35  thompson
; *** empty log message ***
;
; Revision 1.12  2007/03/30 22:02:50  thompson
; Add keyword update_hdr_off
;
; Revision 1.11  2007/02/16 18:41:54  nathan
; Replace PPIMAGE,PPHDR keywords with /PIPELINE and common block
;
; Revision 1.10  2007/02/16 16:32:34  avourlid
;
; Fixed problem with %p image and corrected typos - AV
;
; Revision 1.9  2007/02/02 00:17:36  nathan
; Shorten filename to standard length to be compat with scc_read_summary.pro;
; EXPTIME=-1; update CROTA, PCi_j
;
; Revision 1.8  2007/01/31 22:04:03  nathan
; header updates; set BUNIT; use correct Level digit in filename
;
; Revision 1.7  2007/01/29 22:43:31  nathan
; avg exptime in hdr; in filename use 2 if calibrated else 0
;
; Revision 1.6  2007/01/27 00:16:41  nathan
; major changes for pipeline
;
; Revision 1.5  2007/01/17 19:14:08  colaninn
; added keyword pipeline
;
; Revision 1.4  2007/01/17 15:31:30  colaninn
; add date stamp for polarized images
;
; Revision 1.3  2006/12/07 21:16:04  colaninn
; fixed typo
;
; Revision 1.2  2006/12/06 20:12:46  colaninn
; added output messages
;
; Revision 1.1  2006/10/03 15:29:44  nathan
; moved from dev/analysis/cor
;
; Revision 1.3  2006/09/07 14:17:16  colaninn
; check that image angles are different
;
; Revision 1.2  2006/09/06 21:33:52  colaninn
; check cadence images are within the cadence
;
; Revision 1.1  2006/09/06 19:48:29  colaninn
; first entry
;
;-    
COMMON secchi_daily, ppimage, pphdr
ON_ERROR,2

info="$Id: cor_polariz.pro,v 1.24 2008/01/17 21:43:02 thompson Exp $"
len=strlen(info)
histinfo=strmid(info,1,len-2)

taiobs=anytim2tai(hdrs.DATE_OBS)
srtd=sort(taiobs)
taiobs=taiobs(srtd)

IF ~keyword_set(silent) THEN $
 message, 'checking angles, cadence, obs_id...',/inform

;--Check that the three angles are different
xxx = where(hdrs[0].POLAR EQ hdrs.POLAR, n_angle)
IF n_angle GT 1 THEN $
	message,'IMAGES MUST HAVE 3 DIFFERENT POLARIZATION ANGLES',/info

;--Check all three image were taken before the following sequence
;  do not check if CADENCE=0
;dt=taiobs[2]-taiobs[0]
;xxx = where(hdrs.CADENCE LE 0,nzcadence)
;print,'Cadences:',hdrs.cadence
;IF (nzcadence EQ 0) AND (dt GE min(hdrs.CADENCE)) THEN $
;	message, 'IMAGES SEQUENCE DURATION MUST BE LT CADENCE'
;********CADENCE IS NOT RELIABLE***************

;--If OBS_ID is not equal then check other pararmeters
xxx = where(hdrs[0].OBS_ID NE hdrs.OBS_ID,n_obsid)
IF n_obsid NE 0 THEN $
	message, 'IMAGES ARE NOT FROM THE SAME OBS SEQUENCE: OBS_ID='+strjoin(hdrs.obs_id),/info

;--Get Muller Matrix
mueller = cor_mueller(hdrs[0], hdrs[1], hdrs[2])

im1 = images[*,*,0]
im2 = images[*,*,1]
im3 = images[*,*,2]

;--Calculate Stokes Parameters for Incident Intensity
I = im1*mueller[0,0] + im2*mueller[1,0] + im3*mueller[2,0];total brigtness
Q = im1*mueller[0,1] + im2*mueller[1,1] + im3*mueller[2,1] 
U = im1*mueller[0,2] + im2*mueller[1,2] + im3*mueller[2,2]


pbim = sqrt(Q^2+U^2)

CASE 1 OF
  
  ;--Calculate incident polarized intensity
  keyword_set(pB): BEGIN
    IF ~keyword_set(silent) THEN $
     message,'RETURNING POLARIZED BRIGHTNESS',/inform
    im=pbim
    polval = 1002
    fnstr='P'
    unit='DN/s'
  END
  
  ;--Calculate polarization angle
  keyword_set(mu): BEGIN
    IF ~keyword_set(silent) THEN $
     message,'RETURNING POLARIZATION ANGLE',/inform
    im = 0.5 *atan(U,Q)
    polval = 1004
    fnstr='A'
    unit='radians'
  END
  
  ELSE: BEGIN
    IF ~keyword_set(silent) THEN $
     message,'RETURNING TOTAL BRIGHTNESS',/inform
     im = I
     polval = 1001
     fnstr='B'
     unit='DN/s'
  END    
ENDCASE

;--Percent Polarized image (in COMMON block)
;nzero = where(I ne 0)
ppimage = 100.*RATIO(pbim,I)

IF keyword_set(PERCENT) THEN BEGIN
  IF ~keyword_set(silent) THEN $
   message,'RETURNING PERCENT POLARIZED',/inform
  im = ppimage
  polval = 1003
  fnstr='p'
  unit='percent'
ENDIF



;--Define base header
hdr=hdrs[srtd[0]]
;-Update filename
filename = hdr.filename
strput, filename,'2',16 
hdr.filename = filename

hdr.expcmd	=avg(hdrs.expcmd)
hdr.exptime	=-1.			; N/A for polarized product images
hdr.biasmean	=avg(hdrs.biasmean)
hdr.biassdev	=max(hdrs.biassdev)
hdr.ceb_t	=avg(hdrs.ceb_t)
hdr.temp_ccd	=avg(hdrs.temp_ccd)
hdr.readtime	=avg(hdrs.readtime)
hdr.cleartim	=avg(hdrs.cleartim)
hdr.ip_time	=total(hdrs.ip_time)
hdr.compfact	=avg(hdrs.compfact)
hdr.nmissing	=total(hdrs.nmissing)
hdr.tempaft1	=avg(hdrs.tempaft1)
hdr.tempaft2	=avg(hdrs.tempaft2)
hdr.tempmid1	=avg(hdrs.tempmid1)
hdr.tempmid2	=avg(hdrs.tempmid2)
hdr.tempfwd1	=avg(hdrs.tempfwd1)
hdr.tempfwd2	=avg(hdrs.tempfwd2)
hdr.tempthrm	=avg(hdrs.tempthrm)
hdr.temp_ceb	=avg(hdrs.temp_ceb)
taiend=taiobs[2]+hdrs[srtd[2]].exptime
hdr.date_end	=utc2str(tai2utc(taiend),/noz)
hdr.date_avg	=utc2str(tai2utc(taiobs[0]+(taiend-taiobs[0])/2),/noz)
hdr.N_IMAGES=3
hdr.ENCODERP=-1
hdr.CROTA	=avg(hdrs.crota)
hdr.PC1_1	=cos(hdr.CROTA*!DPI/180.d0)
hdr.PC1_2	=-sin(hdr.CROTA*!DPI/180.d0) ;*(scch.CDELT1/scch.CDELT2) ;may be added later
hdr.PC2_1	=sin(hdr.CROTA*!DPI/180.d0) ;*(scch.CDELT1/scch.CDELT2)
hdr.PC2_2	=cos(hdr.CROTA*!DPI/180.d0)

;--Create header for polarized image
IF datatype(ex) EQ 'UND' THEN ex = create_struct('blank',-1)
 	;   for non-secchi_prep case
if ~keyword_set(update_hdr_off) THEN hdr = SCC_UPDATE_HDR(im,hdr,/level_2)
;--add rev info and 3 filenames to HISTORY
hdr = scc_update_history(hdr,histinfo)

;--add comment for combined values
c_dex= 20-total(strmatch(hdr.COMMENT,''),1)
hdr.COMMENT[c_dex]='EXPTIMEs: '	+trim(string(hdrs[srtd[0]].exptime))+', ' $
				+trim(string(hdrs[srtd[1]].exptime))+', ' $
				+trim(string(hdrs[srtd[2]].exptime))
c_dex=c_dex+1
hdr.COMMENT[c_dex]='CROTAs: ' +trim(string(hdrs[srtd[0]].crota))+', ' $
				+trim(string(hdrs[srtd[1]].crota))+', ' $
				+trim(string(hdrs[srtd[2]].crota))
c_dex=c_dex+1
hdr.POLAR=polval
hdr.bunit=unit

;--add component filenames to comments
for i=1,3 do hdr.comment[c_dex+i]='FILEORIG'+trim(string(i))+': '+hdrs[srtd[i-1]].filename
c_dex=c_dex+3
filename = hdr.filename
strput, filename,fnstr,16
;stop
IF ~keyword_set(CALIBRATE_OFF) THEN BEGIN
	; if calibrated, but the 1 back in that was inserted in scc_update_hdr.pro
	prt1=strmid(filename,0,16)
	prt2=strmid(filename,16,9)
	filename=prt1+'1'+prt2
ENDIF	
strput, filename,fnstr,17

; if not calibrated, put 0 over 2 that was inserted in scc_update_hdr.pro
status = tag_exist(ex,'exptime_off')+tag_exist(ex,'normal_off') $
	+tag_exist(ex,'bias_off')   +tag_exist(ex,'calimg_off') $
	+tag_exist(ex,'calimg_off') +tag_exist(ex,'calibrate_off')

IF status GE 1 THEN strput, filename,'0',16
IF ~keyword_set(silent) THEN message,'NEW filename :'+filename,/inform
hdr.filename = filename

IF keyword_set(PIPELINE) THEN BEGIN
	message,'RETURNING PERCENT POLARIZED in COMMON secchi_daily',/inform
	pphdr=scc_update_hdr(ppimage,hdr,/level_2)
	pphdr.POLAR=1003
	strput, filename, 'p',17
	help,filename
	wait,2
	pphdr.filename=filename
	pphdr.bunit='percent'
ENDIF


;--Date and Logo Stamp (OFF)
IF keyword_set(date_on) THEN im = scc_add_datetime(im,hdr)

RETURN,im
END