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