function cor_mueller, hdr1, hdr2, hdr3
;+
;$Id: cor_mueller.pro,v 1.2 2007/05/21 20:16:06 colaninn Exp $
;
; Project   : STEREO SECCHI
;                   
; Name      : cor_mueller
;               
; Purpose   : This function returns the matrix of elements needed to
;             calculate the Stokes parameters of the input irradiance.
;               
; Explanation: Self-Study Manual on Optical Radiation Measurements:
;              Part1-Concepts, Chapter 6 Distribution of Optical Radiation
;              with Respect to Polarization, John B. Shumaker 
;              To calculate the Stokes parameters of the input irradiance
;              from three polarized irradiance measurements, we only need
;              the first row of the Mueller matrix of each polarizer inverted.
;               
; Use       : IDL> matrix = cor_mueller(hdr1,hdr2,hdr3)
;    
; Inputs    : hdr1 -image header, either fits or SECCHI structure
;             hdr2 -image header, either fits or SECCHI structure
;             hdr3 -image header, either fits or SECCHI structure
;             
;
; Outputs   :  matrix - 3x3 matrix
;
; Keywords  :
;
; Common    : 
;
; Calls     : SCC_FITSHDR2STRUCT
;               
; Category  : Polarization
;
; $Log: cor_mueller.pro,v $
; Revision 1.2  2007/05/21 20:16:06  colaninn
; added silent keyword
;
; Revision 1.1  2006/10/03 15:29:44  nathan
; moved from dev/analysis/cor
;
; Revision 1.2  2006/09/01 19:42:16  colaninn
; add check of level
;
; Revision 1.1  2006/09/01 18:57:06  colaninn
; first entry
;
;-    

caller = scope_traceback( )
IF total(strmatch(caller, '*COR_POLARIZ*',/fold_case)) LT 1 THEN BEGIN
  message,/inform, $
    'DOES NOT CHECK IF HEADERS ARE FROM THE SAME POLARIZATION SEQUENCE'
  
  IF(DATATYPE(hdr1) NE 'STC') THEN hdr=SCC_FITSHDR2STRUCT(hdr)
  IF(DATATYPE(hdr2) NE 'STC') THEN hdr=SCC_FITSHDR2STRUCT(hdr)
  IF(DATATYPE(hdr3) NE 'STC') THEN hdr=SCC_FITSHDR2STRUCT(hdr)
  
  IF (TAG_NAMES(hdr1, /STRUCTURE_NAME) NE 'SECCHI_HDR_STRUCT') OR $
  (TAG_NAMES(hdr2, /STRUCTURE_NAME) NE 'SECCHI_HDR_STRUCT') OR $
  (TAG_NAMES(hdr3, /STRUCTURE_NAME) NE 'SECCHI_HDR_STRUCT') THEN $
  MESSAGE,'ONLY SECCHI HEADER STRUCTURE ALLOWED'
ENDIF

angle1 = hdr1.polar
angle2 = hdr2.polar
angle3 = hdr3.polar

;--Calculate Mueller matrix output
;----for ideal polarizers k1 = 1 and k2 = 0
;----we will measure s and d after launch

s = .5 ; s = .5*(K1 + k2) 
d = .5 ; d = .5*(k1 - k2)

;First row of Mueller matrix for each polarization angle
;  [im1,im2,im3] = x*[I,Q,U]
x = [[s, d*cos(2*angle1*!DtoR ), d*sin(2*angle1*!DtoR)],$
     [s, d*cos(2*angle2*!DtoR ), d*sin(2*angle2*!DtoR)],$
     [s, d*cos(2*angle3*!DtoR ), d*sin(2*angle3*!DtoR)]]

;We then invert the matrix to find [I,Q,U]
x = invert(x)

RETURN, x
END