pro cal_pol
;+
;NAME:
;       cal_pol
;PURPOSE:
;       analyze polarization test data: for each paw position do a sine-fit
;	to resulting intensity (cal and obs) as input polarization is rotated
;	- for lcp/rcp use rotating input quarterwave (behind fixed polarizer)
;	- for s/p use rotating input halfwave (behind fixed polarizer)
;	store resulting phase, amplitude, offset and sigma as fits images
;	to be read later and displayed by 'cal_fig_pol'
;INPUT:
;OPTIONAL INPUT:
;OPTIONAL KEYWORD INPUT:
;OUTPUT:
;	four fits files per set, 32 total in directory 'outdir'
;HISTORY:
;       Written  9-Sep-94 by I. Zayer (Ver. 1.0)
;	Ver 1.1 12-Sep-94 IZ: fixed bug when switching to lin. pol.
;	Ver 2.0 12-Sep-94 IZ: add fit for fixed input & rotating PAW, and also
;		add fit fo rotating input for all PAW positions.
;-
;

progver='CAL_POL  Ver. 2.0'
serno = 619
outdir=	'$MDI_CAL_DATA'			; = /mdisw/dbase/cal/results I guess..

d=indgen(44)*9+8			;index of all darks
indx=intarr(44,2,4)			;array of indices (inputs,modes,pawpos)
indx(*,0,3)=d-1					;inputs: stim waveplate pos.
for i=2,0,-1 do indx(*,0,i)=indx(*,0,i+1)-1	;modes: 0=obs, 1=cal
indx(*,1,*)=indx(*,0,*)-4			;pawpos: 0=17,1=53,2=89,3=125

inpol=intarr(19,2)		;indices (within indx) of input polarization:
inpol(*,0)=indgen(19)+22	;0=quarterwave, 1=halfwave
inpol(*,1)=indgen(19)+1	
inpol(7:*,0)=inpol(7:*,0)+1	;qwave @140deg is repeated in data, fix it...

fixin=intarr(9,4)		;indices (within datafile) of MDI polarization:
imset=[0,2,1,3,4,6,5,7,8]	;for input pol: 0=p-wave,1=s-wave,2=lcp,3=rcp
fixin(*,0)=imset+180		;MDI pol. being: 17,53,89,125 --> need to make 
fixin(*,1)=imset+189		;sequence 17,89,53,125 for sine fit
fixin(*,2)=imset+378
fixin(*,3)=imset+387	

paw=['_p017','_p053','_p089','_p125','_pall']	;use p to subsrcibe
intit=['_inqwave','_inhwave','_inpwave','_inswave','_inlcp','_inrcp']
						;use inp to subsrcibe
comtit=['quarterwave','halfwave','p-wave','s-wave','lcp','rcp']
						;use inp to subsrcibe
mode=['_obs','_cal']				;use m to subsrcibe
restit=['_ph','_amp','_off','_sig']		;use r to subsrcibe
rescom=['Phase','Amplitude','Offset','Sigma']	;use r to subsrcibe

angle=indgen(19)*40	;good for quarterwave, will need to double for halfwave
res=fltarr(256,256,4)

read_mdi,serno,d,dim,hdr & hdr=hdr(*,0)
SXADDPAR, hdr, 'BITPIX', -32
SXADDPAR, hdr, 'PROGVER', progver
SXADDPAR, hdr, 'PROG_RUN', !stime
SXADDPAR, hdr, 'DATASET', serno

; fit sine wave to fixed PAW position and rotating input
ang=angle
for inp=0,1 do begin
  for p=0,3 do begin
    if (inp ge 1) or (p ge 2) then ang=angle*2 else ang=angle
    for m=0,1 do begin
      read_mdi,serno,indx(inpol(*,inp),m,p),im & im=im-dim(*,*,inpol(*,inp))
      sinfit,ang,im,ph,amp,off,sig,rad,dummy,fwhm=-1
      res(*,*,0)=ph & res(*,*,1)=amp & res(*,*,2)=off & res(*,*,3)=sig
      for r=0,3 do begin
        head=hdr
        fnam=outdir+'/s'+strtrim(serno,2)+paw(p)+mode(m)+intit(inp)+restit(r)+'.fits'
        SXADDPAR, head, 'COMMENT',$
	  rescom(r)+' of sine fit to rotating input '+comtit(inp)
        wrt_fits, fnam, head, reform(res(*,*,r))
      endfor
    endfor
  endfor
endfor

; now fit sine wave to rotating PAW and fixed input
ang=indgen(4)*90
for inp=2,5 do begin
  read_mdi,serno,fixin(*,inp-2),im & for i=0,7 do im(*,*,i)=im(*,*,i)-im(*,*,8)
  for m=0,1 do begin
    sinfit,ang,im(*,*,4-m*4:7-m*4),ph,amp,off,sig,rad,dummy,fwhm=-1
    res(*,*,0)=ph & res(*,*,1)=amp & res(*,*,2)=off & res(*,*,3)=sig
    for r=0,3 do begin
      head=hdr
      fnam=outdir+'/s'+strtrim(serno,2)+paw(4)+mode(m)+intit(inp)+restit(r)+'.fits'
      SXADDPAR, head, 'COMMENT',$
	rescom(r)+' of sine fit to fixed input '+comtit(inp)
      wrt_fits, fnam, head, reform(res(*,*,r))
    endfor
  endfor
endfor


end