dome = 1

; The difference between the two sets of centers is that
; the first two were determined from the moon for medium exposures
; The last two were determined by looking at minimizing Ut which
; is centered on the Sun, not the moon
; These are for 1/3 third sized images
xcen=356.5
ycen=323.
xcen=356.5 - 5.
ycen=323. - 4.
rsun = 110.

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

close,1 & close,13

loadct,1

; These two philes must be in phase

lfilenm = 'list'
cfilenm = 'poise_140_centers'
filenm='         '
close,13 & openr,13,lfilenm
close,14 & openr,14,cfilenm

a=0.
b=0.

; Read in FTS image
    readf,14,time,imgx,imgy
    print,imgx,imgy
    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)
; Rather than using a dark, take bias from a bottom row
;    sbias=total(img(*,ysz-1)/n_elements(img(*,ysz-1)))
;    img = img - sbias
; 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

;;;;;;;;;;;temporary to see if offset correction works better first;;;
    sbias=total(img(*,ysz-1)/n_elements(img(*,ysz-1)))
    img = img - sbias

    simg1 = rebin(img,xsz/3,ysz/3)
    tv,simg1
    img=img+rimg/100.
    img=rot(img,0.,1.,imgx,imgy,/pivot,/interp)
; Here is where opal is divided into data
;
;
;
; Scaling from here on will be changed
    simg1 = rebin(img,xsz/3,ysz/3 )
    tv,simg1

; Read in FTS image
    readf,14,time,imgx,imgy
    print,imgx,imgy
    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(*,ysz-1)/n_elements(img(*,ysz-1)))
	img = img - sbias
    simg2 = rebin(img,xsz/3,ysz/3 )
    tv,simg2
    img=img+rimg/100.
    img=rot(img,0.,1.,imgx,imgy,/pivot,/interp)
; Here is where opal is divided into data
; Scaling from here on may be changed
    simg2 = rebin(img,xsz/3,ysz/3 )
	tv,simg2

; Read in FTS image
    readf,14,time,imgx,imgy
    print,imgx,imgy
    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(*,ysz-1)/n_elements(img(*,ysz-1)))
	img = img - sbias
    simg3 = rebin(img,xsz/3,ysz/3 )
    tv,simg3
    img=img+rimg/100.
    img=rot(img,0.,1.,imgx,imgy,/pivot,/interp)
; Here is where opal is divided into data
; Scaling from here on may be changed
    simg3 = rebin(img,xsz/3,ysz/3 )
	tv,simg3

; Read in FTS image
    readf,14,time,imgx,imgy
    print,imgx,imgy
    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(*,ysz-1)/n_elements(img(*,ysz-1)))
	img = img - sbias
    simg4 = rebin(img,xsz/3,ysz/3 )
    tv,simg4
    img=img+rimg/100.
    img=rot(img,0.,1.,imgx,imgy,/pivot,/interp)
; Here is where opal is divided into data
; Scaling from here on may be changed
    simg4 = rebin(img,xsz/3,ysz/3 )
	tv,simg4

close,13

;goto,noplots

	intensity=simg1+simg2+simg3+simg4
	tv,intensity/1024.

	q = simg1-simg3
	u = simg2-simg4
	pol = sqrt(q*q+u*u)
	tv,pol/64.
	stop

;mask = makemask(xsz/3,ysz/3,xcen,ycen,rsun)

stop

	sob=sobel(pol)
	tv,sob/32.
	stop

	mad=sob/(pol+250.)
	tv,mad*128.
	stop

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

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

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

noplots:

goto, noaz
;Performing azimuthal analysis
print,'azimuthal analysis'

az1=fltarr(360)
az2=fltarr(360)
az3=fltarr(360)
az4=fltarr(360)
toti = fltarr(360)
angle=fltarr(360)

for ii=0,359 do begin
	rang = ii*.01745
	angle(ii)=ii
	px = fix(xcen + 120. * cos(rang))
	py = fix(ycen + 120. * sin(rang))
	az1(ii)=simg1(px,py)
	az2(ii)=simg2(px,py)
	az3(ii)=simg3(px,py)
	az4(ii)=simg4(px,py)
	toti(ii)=(intensity(px,py))/4.
end
!y.range=[0.5,1.5]
plot,angle,az1/toti,line=0
oplot,angle,az2/toti,line=1
oplot,angle,az3/toti,line=2
oplot,angle,az4/toti,line=3

stop

noaz:

;
; Compute matrix form of Q and U, and invert to see how well we did
;
; Definition of X matrix from evaluation of shifts above
x = fltarr(3,3)
x = [[1.,0.04,0.045],[0.,0.627,0.315],[0.,0.32,-0.632]]
xinv = invert(x)

qmeas = simg1-simg2-simg3+simg4
umeas = simg1+simg2-simg3-simg4
q = xinv(1,1)*qmeas+xinv(2,1)*umeas
u = xinv(1,2)*qmeas+xinv(2,2)*umeas

; use these lines to get qt and ut without crosstalk correction from
; X matrix
;q = 1.414 * qmeas
;u = -1.414 * umeas

tv,q
tv,u
		
qt = q
ut = u
ca=q
sa=u
for ii=0,xsz/3-1 do begin
	px = float(ii - xcen)
	for jj=0,ysz/3-1 do begin
		py = float(jj - ycen)
		ang = atan(py/px)
		ca(ii,jj) = cos(2.*ang)
		sa(ii,jj) = sin(2.*ang)
	end
print,ii
end

qt = q*ca - u*sa
ut = q*sa + u*ca

tv,qt
tv,ut
tv,ut/16. + 128.

end