;
; Polarization calibration routine for a single exposure
; time of poise images
; Opal normalization is not yet included
; X matrix is included
;
dome = 1

if ( dome ) then begin ;{
    window,xsize=678,ysize=678
endif ;}

close,1 & close,13

loadct,1

lfilenm = 'list'
filenm='         '
close,13 & openr,13,lfilenm

a=0.
b=0.

; Read in FTS image
    readf,13,filenm
    img = readfits(filenm,hdr)
    xsz = fix(sxpar (hdr , 'NAXIS1'))
    ysz = fix(sxpar (hdr , 'NAXIS2'))
; Define these arrays now that dimensions are known (right)
    rimg=fltarr(xsz,ysz)
    row=fltarr(ysz)
; Fix the byte order
    byteorder,img
; Correct for integer wraparound
    img = long(img)
    img = (img lt 0 ) * ( img + 32768 + 32768)  + (img ge 0)*(img)
; Flip top to bottom for IDL
    img = rotate(img,7)
; Record where intensity is too high for future use
    a1 = where(img gt 40000)
; Compute an average value for each row and turn that into an image
; Add a fraction of the average row image to eliminate horizontal
; offset in the CCD which is proportional to total signal of row
; Factor of 1/100. is just a guess
    for jj=0,ysz-2 do begin
	row(jj)=total(img(*,jj+1))
	row(jj)=row(jj)/(1.*xsz)
	rimg(*,jj)=row(jj)
    end
; Rather than using a dark, take bias from a bottom row
    sbias=total(img(4,*))/(1.*ysz)
    img = img - sbias
    simg1 = rebin(img,xsz/3,ysz/3 )
    tv,simg1
    img1=img+rimg/100.
; Here is where opal is divided into data
;
;
;
; Scaling from here on will be changed
    simg1 = rebin(img1,xsz/3,ysz/3 )
    tv,simg1

; Read in FTS image
    readf,13,filenm
    img = readfits(filenm,hdr)
; Fix the byte order
    byteorder,img
; Correct for integer wraparound
    img = long(img)
    img = (img lt 0 ) * ( img + 32768 + 32768)  + (img ge 0)*(img)
; Flip top to bottom for IDL
    img = rotate(img,7)
; Record where intensity is too high for future use
	a2 = where(img gt 40000)
; Compute an average value for each row and turn that into an image
; Add a fraction of the average row image to eliminate horizontal
; offset in the CCD which is proportional to total signal of row
; Factor of 1/100. is just a guess
	for jj=0,ysz-2 do begin
		row(jj)=total(img(*,jj+1))
		row(jj)=row(jj)/(1.*xsz)
		rimg(*,jj)=row(jj)
	end
; Rather than using a dark, take bias from a bottom row
    sbias=total(img(4,*))/(1.*ysz)
	img = img - sbias
    simg2 = rebin(img,xsz/3,ysz/3 )
    tv,simg2
    img2=img+rimg/100.
; Here is where opal is divided into data
; Scaling from here on may be changed
    simg2 = rebin(img2,xsz/3,ysz/3 )
    tv,simg2

; Read in FTS image
    readf,13,filenm
    img = readfits(filenm,hdr)
; Fix the byte order
    byteorder,img
; Correct for integer wraparound
    img = long(img)
    img = (img lt 0 ) * ( img + 32768 + 32768)  + (img ge 0)*(img)
; Flip top to bottom for IDL
    img = rotate(img,7)
; Record where intensity is too high for future use
	a3 = where(img gt 40000)
; Compute an average value for each row and turn that into an image
; Add a fraction of the average row image to eliminate horizontal
; offset in the CCD which is proportional to total signal of row
; Factor of 1/100. is just a guess
	for jj=0,ysz-2 do begin
		row(jj)=total(img(*,jj+1))
		row(jj)=row(jj)/(1.*xsz)
		rimg(*,jj)=row(jj)
	end
; Rather than using a dark, take bias from a bottom row
    sbias=total(img(4,*))/(1.*ysz)
	img = img - sbias
    simg3 = rebin(img,xsz/3,ysz/3 )
    tv,simg3
    img3=img+rimg/100.
; Here is where opal is divided into data
; Scaling from here on may be changed
    simg3 = rebin(img3,xsz/3,ysz/3 )
    tv,simg3

; Read in FTS image
    readf,13,filenm
    img = readfits(filenm,hdr)
; Fix the byte order
    byteorder,img
; Correct for integer wraparound
    img = long(img)
    img = (img lt 0 ) * ( img + 32768 + 32768)  + (img ge 0)*(img)
; Flip top to bottom for IDL
    img = rotate(img,7)
; Record where intensity is too high for future use
	a4 = where(img gt 40000)
; Compute an average value for each row and turn that into an image
; Add a fraction of the average row image to eliminate horizontal
; offset in the CCD which is proportional to total signal of row
; Factor of 1/100. is just a guess
	for jj=0,ysz-2 do begin
		row(jj)=total(img(*,jj+1))
		row(jj)=row(jj)/(1.*xsz)
		rimg(*,jj)=row(jj)
	end
; Rather than using a dark, take bias from a bottom row
    sbias=total(img(4,*))/(1.*ysz)
	img = img - sbias
    simg4 = rebin(img,xsz/3,ysz/3 )
    tv,simg4
    img4=img+rimg/100.
; Here is where opal is divided into data
; Scaling from here on may be changed
    simg4 = rebin(img4,xsz/3,ysz/3 )
    tv,simg4

close,13

    intensity=img1+img2+img3+img4
    sint = rebin(intensity,xsz/3,ysz/3 )

;
; Compute matrix form of Q and U, and invert to see how well we did

	x = fltarr(3,3)
	x = [[1.,0.04,0.045],[0.,0.627,0.315],[0.,0.32,-0.632]]
	xinv = invert(x)

	qmeas = img1-img2-img3+img4
	umeas = img1+img2-img3-img4
	q = xinv(1,1)*qmeas+xinv(2,1)*umeas
	u = xinv(1,2)*qmeas+xinv(2,2)*umeas

	pol = sqrt(q*q+u*u)
    spol = rebin(pol,xsz/3,ysz/3 )
	tv,spol/64.
	stop

print, 'masked pb'
; Read in 1/3 scale mask.  See pb.pro for centers and call to make it
	mask=spol
	openr,1,'/thisbe/d/elmore/Eclipse/mask'
	readu,1,mask
	close,1
	tv,spol*mask/256.
	stop

	sob=sobel(pol)
    ssob = rebin(sob,xsz/3,ysz/3 )
	tv,ssob/32.
	stop

	mad=sob/(pol+250.)
    smad = rebin(mad,xsz/3,ysz/3 )
	tv,smad*128.
	stop

print,"log intensity"
	tv,bytscl(alog10(sint),2,6)
	stop

print,"log pb"
	tv,bytscl(alog10(spol),1,5)
	stop

print,"log pb/intensity"
; next good for long exposure
;	tv,bytscl(alog10(spol/sint),-1.8,-.1)  
; next good for medium
	tv,bytscl(alog10(spol/sint),-3.4,.4)
	stop

end