pro xttx_plot_cal_usage, ifile, ofile ;+ ; ; procedure: xttx_plot_cal ; ; purpose: plot xttx.f fit to calibration optics or clear port. ; ; routines: xttx_plot_cal_usage xttx_cal_click xttx_plot_discrim ; xttx_plot_sku xttx_text_in xttx_plot_cal ; ; author: paul@ncar 11/95 ; ;============================================================================== ; if 1 then begin print print, "usage: xttx_plot_cal, ifile, ofile" print print, " Plot xttx.f fit to calibration optics or clear port." print print, " Arguments:" print, " ifile - path to data file input to xttx" print, " ofile - path to fitted file output by xttx" print print, " Keyword:" print, " xfile - path to xfile as output by xttx.f;" print, " used to inverse X & D from the data." print, " tfile - path to tfile as output by xttx.f;" print, " used to inverse T from the data." print, " date - if date is set to vector" print, " [year,month,day], final plot" print, " will have disk center t_vectors" print, " plotted as smooth curves" print, " (no effect if linear polarizer" print, " present)." print, " ex:" print print, " ifile = '~paul/stokes/src/xttx/test/92.03.24.clr'" print, " ofile = '~paul/stokes/src/xttx/test/test0_o_1'" print, " xttx_plot_cal, ifile, ofile" print print, " xttx_plot_cal, ifile, ofile $ print, " , xfile=xfile, tfile=tfile, date=[1998,5,28]" print return endif ;- end ;============================================================================== ; ; procedure: xttx_cal_click ; ; purpose: wait for user to click window. ; ;============================================================================== pro xttx_cal_click, ww, message ;Don't wait if right button clicked. if ww.finish eq 0 then begin ;Wait for button down. num = -1 repeat begin num = (num+1) mod n_elements(ww.w) down = 0 if ww.w(num) ge 0 then begin wset, ww.w(num) cursor, xx,yy,/device,/nowait down = !err if down eq 0 then wait, .002 end end until down ne 0 ;Check if right button clicked to finish. if down eq 2 then ww.finish=1 if down eq 4 then ww.exit=1 ;Wait for button up. repeat cursor, xx,yy,/device,/nowait until !err eq 0 end end ;============================================================================== ; ; procedure: xttx_plot_discrim ; ; purpose: do xttx_plot_cal discrimanator plots. ; ;============================================================================== pro xttx_plot_discrim, dvecs, mvecs, dd, mm, ww ;Stokes component indices. I = 0 Q = 1 U = 2 V = 3 ;Normalized vectors. dnorm=dvecs & mnorm=mvecs dnorm(I,*)=dvecs(I,*)/dvecs(I,*) & mnorm(I,*)=mvecs(I,*)/mvecs(I,*) dnorm(Q,*)=dvecs(Q,*)/dvecs(I,*) & mnorm(Q,*)=mvecs(Q,*)/mvecs(I,*) dnorm(U,*)=dvecs(U,*)/dvecs(I,*) & mnorm(U,*)=mvecs(U,*)/mvecs(I,*) dnorm(V,*)=dvecs(V,*)/dvecs(I,*) & mnorm(V,*)=mvecs(V,*)/mvecs(I,*) ;Compute residual polarization fraction. rrrr = sqrt( $ (dnorm(Q,*)-mnorm(Q,*))^2 $ + (dnorm(U,*)-mnorm(U,*))^2 $ + (dnorm(V,*)-mnorm(V,*))^2 ) ;Set discriminator for table angle. tbl = (dd.table_angle-dd.azimuth-30.+(360.+360.+360.)) mod 360. dtable = ( round(tbl/45.) mod 8 )*45L ;Set discriminator array. case 1 of ;Clear port case. dd.ix_linear(0) eq 0 and dd.ix_retarder(0) eq 0: $ begin ;Use table discriminator on clear port. discrim = lonarr(dd.dim) ;Check for old spectra file case. if max(dd.ix_dark) eq 1 then dtable(*)=0 end ;Retarder only case. dd.ix_linear(0) eq 0 and dd.ix_retarder(0) ne 0: $ begin discrim = (dd.ix_retarder-1)*45L end ;Linear polarizer only case. dd.ix_linear(0) ne 0 and dd.ix_retarder(0) eq 0: $ begin discrim = (dd.ix_linear-1)*45L end ;Linear polarizer with retarder case. else: begin discrim = (dd.ix_linear-1)*45L whr = where( ( (dd.ix_retarder-dd.ix_linear+8) mod 8 ) $ eq 1, nwhr ) if nwhr ne 0 then discrim(whr)=discrim(whr)-360 end end ;Find number of negative discriminators. whr = where( discrim lt 0, n_negative ) ;Size of windows to be opened. xsize=512 & ysize=round(xsize*11./8.5) ;Open windows. window, /free, xsize=xsize, ysize=ysize $ , xpos=0, ypos=900-ysize $ , title=ww.click ww.w(0) = !d.window window, /free, xsize=xsize, ysize=ysize $ , xpos=xsize+30, ypos=900-ysize $ , title=ww.click ww.w(1) = !d.window if n_negative ne 0 then begin window, /free, xsize=xsize, ysize=ysize $ , xpos=30, ypos=900-ysize-30 $ , title=ww.click ww.w(2) = !d.window window, /free, xsize=xsize, ysize=ysize $ , xpos=30+xsize+30, ypos=900-ysize-30 $ , title=ww.click ww.w(3) = !d.window end window, /free, xsize=xsize, ysize=ysize $ , xpos=60, ypos=900-ysize-60 $ , title=ww.click ww.w(4) = !d.window window, /free, xsize=xsize, ysize=ysize $ , xpos=60+xsize+30, ypos=900-ysize-60 $ , title=ww.click ww.w(5) = !d.window ;Maximum color index. mx_col = !d.n_colors-1 ;Maximums for plot scale. mx_pol = max(abs(dnorm(Q:V,*))) > max(abs(mnorm(Q:V,*))) ;Draw plot frmes in accumulation windows. wset, ww.w(0) !p.region = 0. plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , title=dd.file+' '+dd.vtype $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype wset, ww.w(1) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , title=mm.file+' '+dd.vtype $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype if n_negative ne 0 then begin wset, ww.w(2) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , title=dd.file+' '+dd.vtype $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype wset, ww.w(3) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , title=mm.file+' '+dd.vtype $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype end ;Vector plot separated with discriminator. for tab= 0,360,45 do begin for dis=-360,315,45 do begin whr = where( dtable eq tab and discrim eq dis, nwhr ) if nwhr gt 0 then begin if nwhr eq 1 then begin ;Set where array sorted on azimuth. awhr = whr aaaa = dd.azimuth(awhr) ;Do plots. if dis ge 0. then wset,ww.w(0) else wset,ww.w(2) xyouts, aaaa, dnorm(Q,awhr), 'q', /data, color=0 xyouts, aaaa, dnorm(U,awhr), 'u', /data, color=0 xyouts, aaaa, dnorm(V,awhr), 'v', /data, color=0 xyouts, aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts, aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts, aaaa, mnorm(V,awhr), 'V', /data, color=0 if dis ge 0. then wset,ww.w(1) else wset,ww.w(3) xyouts, aaaa, dnorm(Q,awhr), 'q', /data, color=0 xyouts, aaaa, dnorm(U,awhr), 'u', /data, color=0 xyouts, aaaa, dnorm(V,awhr), 'v', /data, color=0 xyouts, aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts, aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts, aaaa, mnorm(V,awhr), 'V', /data, color=0 wset, ww.w(4) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' '+dd.vtype $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype xyouts, aaaa, dnorm(Q,awhr), 'q', /data, color=0 xyouts, aaaa, dnorm(U,awhr), 'u', /data, color=0 xyouts, aaaa, dnorm(V,awhr), 'v', /data, color=0 xyouts, aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts, aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts, aaaa, mnorm(V,awhr), 'V', /data, color=0 wset, ww.w(5) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' '+dd.vtype $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype xyouts, aaaa, dnorm(Q,awhr), 'q', /data, color=0 xyouts, aaaa, dnorm(U,awhr), 'u', /data, color=0 xyouts, aaaa, dnorm(V,awhr), 'v', /data, color=0 xyouts, aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts, aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts, aaaa, mnorm(V,awhr), 'V', /data, color=0 end else begin ;Set where array sorted on azimuth. awhr = whr( sort( dd.azimuth(whr) ) ) aaaa = dd.azimuth(awhr) ;Do plots. if dis ge 0. then wset,ww.w(0) else wset,ww.w(2) plots, aaaa, dnorm(Q,awhr), color=0 plots, aaaa, dnorm(U,awhr), color=0 plots, aaaa, dnorm(V,awhr), color=0 xyouts,aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts,aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts,aaaa, mnorm(V,awhr), 'V', /data, color=0 if dis ge 0. then wset,ww.w(1) else wset,ww.w(3) xyouts,aaaa, dnorm(Q,awhr), 'q', /data, color=0 xyouts,aaaa, dnorm(U,awhr), 'u', /data, color=0 xyouts,aaaa, dnorm(V,awhr), 'v', /data, color=0 plots, aaaa, mnorm(Q,awhr), color=0 plots, aaaa, mnorm(U,awhr), color=0 plots, aaaa, mnorm(V,awhr), color=0 wset, ww.w(4) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' '+dd.vtype $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype plots, aaaa, dnorm(Q,awhr), color=0 plots, aaaa, dnorm(U,awhr), color=0 plots, aaaa, dnorm(V,awhr), color=0 xyouts,aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts,aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts,aaaa, mnorm(V,awhr), 'V', /data, color=0 wset, ww.w(5) plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' '+dd.vtype $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of '+dd.vtype xyouts,aaaa, dnorm(Q,awhr), 'q', /data, color=0 xyouts,aaaa, dnorm(U,awhr), 'u', /data, color=0 xyouts,aaaa, dnorm(V,awhr), 'v', /data, color=0 plots, aaaa, mnorm(Q,awhr), color=0 plots, aaaa, mnorm(U,awhr), color=0 plots, aaaa, mnorm(V,awhr), color=0 end ;Print vector numbers. print print, ' vec# optics ' $ +'q(data) ' $ +'u(data) ' $ +'v(data) ' $ +' PHI_ma L_ang D_ang' print, ' ' $ +'Q(model) ' $ +'U(model) ' $ +'V(model) ' $ +' residual' for ipt=0,n_elements(awhr)-1 do begin pt = awhr(ipt) tbl = round( dd.table_angle(pt)-dd.azimuth(pt)-30. ) tbl = (tbl+360+360+360) mod 360 lang = ' ' if dd.ix_linear(pt) ne 0 then $ lang = string((dd.ix_linear(pt)-1)*45,format='(i7)') dang = ' ' if dd.ix_retarder(pt) ne 0 then $ dang = string(((dd.ix_retarder(pt)+1)*45) mod 360 $ ,format='(i7)') print, format='(i6,2x,4i1,3f11.6,i19,a,a)' $ , dd.whr(pt) $ , dd.ix_dark(pt), dd.ix_bulb(pt) $ , dd.ix_linear(pt), dd.ix_retarder(pt) $ , dnorm(Q,pt), dnorm(U,pt), dnorm(V,pt) $ , tbl, lang, dang print, format='(12x,3f11.6,f12.6)' $ , mnorm(Q,pt), mnorm(U,pt), mnorm(V,pt), rrrr(pt) end ;Get click to continue. xttx_cal_click, ww if ww.exit then return end end end ;Delete windows with partial data sets. wdelete, ww.w(4), ww.w(5) ww.w(4:5) = -1 end ;============================================================================== ; ; procedure: xttx_plot_sku ; ; purpose: do sku residual plots in solar frame ; ;============================================================================== pro xttx_plot_sku, tvecs, tmodel, dvecs, mvecs, dd, mm, ww, date, tfile ;Stokes component indices. I = 0 Q = 1 U = 2 V = 3 ;Normalized telescope vectors. dnorm=tvecs & mnorm=tmodel dnorm(I,*)=tvecs(I,*)/tvecs(I,*) & mnorm(I,*)=tmodel(I,*)/tmodel(I,*) dnorm(Q,*)=tvecs(Q,*)/tvecs(I,*) & mnorm(Q,*)=tmodel(Q,*)/tmodel(I,*) dnorm(U,*)=tvecs(U,*)/tvecs(I,*) & mnorm(U,*)=tmodel(U,*)/tmodel(I,*) dnorm(V,*)=tvecs(V,*)/tvecs(I,*) & mnorm(V,*)=tmodel(V,*)/tmodel(I,*) ;Compute residual polarization fraction. rrrr = sqrt( $ (dnorm(Q,*)-mnorm(Q,*))^2 $ + (dnorm(U,*)-mnorm(U,*))^2 $ + (dnorm(V,*)-mnorm(V,*))^2 ) ;Set discriminator for table angle. tbl = (dd.table_angle-dd.azimuth-30.+(360.+360.+360.)) mod 360. dtable = ( round(tbl/45.) mod 8 )*45L ;Set discriminator array. case 1 of ;Clear port case. dd.ix_linear(0) eq 0 and dd.ix_retarder(0) eq 0: $ begin ;Use table discriminator on clear port. discrim = lonarr(dd.dim) ;Check for old spectra file case. if max(dd.ix_dark) eq 1 then dtable(*)=0 end ;Retarder only case. dd.ix_linear(0) eq 0 and dd.ix_retarder(0) ne 0: $ begin discrim = (dd.ix_retarder-1)*45L end ;Linear polarizer only case. dd.ix_linear(0) ne 0 and dd.ix_retarder(0) eq 0: $ begin discrim = (dd.ix_linear-1)*45L end ;Linear polarizer with retarder case. else: begin discrim = (dd.ix_linear-1)*45L whr = where( ( (dd.ix_retarder-dd.ix_linear+8) mod 8 ) $ eq 1, nwhr ) if nwhr ne 0 then discrim(whr)=discrim(whr)-360 end end ;Find number of negative discriminators. whr = where( discrim lt 0, n_negative ) ;Size of windows to be opened. xsize=320 & ysize=round(xsize*11./8.5) ;Open windows. window, /free, xsize=2*xsize, ysize=2*ysize, xpos=1130-2*xsize, ypos=40 $ , title = ww.click ww.w(7) = !d.window ;Maximum color index. mx_col = !d.n_colors-1 ;Maximums for plot scale. mx_pol = max(abs(dnorm(Q:V,*))) > max(abs(mnorm(Q:V,*))) ;Get disk center t_vectors if date is ;present. if n_elements(date) eq 3 then begin t_plot, tfile, date(0), date(1), date(2), /noplot $ , azimuth=az, elevation=el, tvecs=tvs sqrtqquu = sqrt( tvs(Q,*)^2 + tvs(U,*)^2 ) mx_pol = mx_pol > max(sqrtqquu) end ;Draw plot frame. !p.region = [ 0.0, 0.5, 0.5, 1.0 ] plot, [65.,295.], [-mx_pol,mx_pol], xstyle=1, /nodata $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth (degrees)' $ , ytitle='polarization components of | t >' ;Vector plot separated with discriminator. for tab= 0,360,45 do begin for dis=-360,315,45 do begin whr = where( dtable eq tab and discrim eq dis, nwhr ) if nwhr gt 0 then begin if nwhr eq 1 then begin ;Set where array sorted on azimuth. awhr = whr aaaa = dd.azimuth(awhr) ;Do plots. oplot, aaaa, dnorm(Q,awhr), psym=1, color=0 oplot, aaaa, dnorm(U,awhr), psym=7, color=0 oplot, aaaa, dnorm(V,awhr), psym=8, color=0 xyouts, aaaa, mnorm(Q,awhr), 'Q', /data, color=0 xyouts, aaaa, mnorm(U,awhr), 'U', /data, color=0 xyouts, aaaa, mnorm(V,awhr), 'V', /data, color=0 end else begin ;Set where array sorted on azimuth. awhr = whr( sort( dd.azimuth(whr) ) ) aaaa = dd.azimuth(awhr) ;Do plots. oplot,aaaa, dnorm(Q,awhr), psym=1, color=0 oplot,aaaa, dnorm(U,awhr), psym=7, color=0 oplot,aaaa, dnorm(V,awhr), psym=8, color=0 if ( tab eq 0 and dis eq 0 ) $ or ( tab eq 270 and dis eq 0 ) then begin plots, aaaa, mnorm(Q,awhr), color=0 plots, aaaa, mnorm(U,awhr), color=0 if tab eq 0 then $ plots, aaaa, mnorm(V,awhr), color=0 end end ;Plot disk center t_vectors if date given. if n_elements(date) eq 3 then begin ewhr = where( el gt 0. ) plots, az(ewhr), tvs(Q,ewhr) $ , linestyle=2, color=0 plots, az(ewhr), tvs(U,ewhr) $ , linestyle=2, color=0 plots, az(ewhr), tvs(V,ewhr) $ , linestyle=2, color=0 ;plots, az(ewhr), sqrtqquu(ewhr) $ ;, linestyle=2, color=2 end end end end ;Normalized solar vectors. dnorm(I,*)=dvecs(I,*)/dvecs(I,*) & mnorm(I,*)=mvecs(I,*)/mvecs(I,*) dnorm(Q,*)=dvecs(Q,*)/dvecs(I,*) & mnorm(Q,*)=mvecs(Q,*)/mvecs(I,*) dnorm(U,*)=dvecs(U,*)/dvecs(I,*) & mnorm(U,*)=mvecs(U,*)/mvecs(I,*) dnorm(V,*)=dvecs(V,*)/dvecs(I,*) & mnorm(V,*)=mvecs(V,*)/mvecs(I,*) ;Set where array sorted on azimuth. awhr = sort( dd.azimuth ) aaaa = dd.azimuth(awhr) aamn = min(aaaa,max=aamx) ;Plot bounds. mx = 1.3*( max(abs(dnorm(Q:V,*))) > max(abs(mnorm(Q:V,*))) ) ;Draw Q residual plot. !p.region = [ 0.5, 0.5, 1.0, 1.0 ] plot, [65.,295.], [-mx,mx], xstyle=1, ystyle=1, /nodata, /noerase $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth', ytitle='Q solar residual' plots, aaaa, 0.*mnorm(Q,awhr), color=0 oplot, aaaa, dnorm(Q,awhr), psym=1, color=0 oplot, [ 65.,aamn], [0.,0.], color=0 oplot, [295.,aamx], [0.,0.], color=0 stdQ = stdev( dnorm(Q,awhr), meanQ ) oplot, [ 65.,295.], [meanQ,meanQ], linestyle=2, color=0 xyouts, 285., .9*mx, 'mean'+string(meanQ,format='(f10.6)'), color=0, align=1. xyouts, 285., .83*mx, 'stdev'+string(stdQ,format='(f10.6)'), color=0, align=1. ;Draw U residual plot. !p.region = [ 0.0, 0.0, 0.5, 0.5 ] plot, [65.,295.], [-mx,mx], xstyle=1, ystyle=1, /nodata $ , background=!d.n_colors-1, color=0, /noerase $ , xtitle='azimuth', ytitle='U solar residual' plots, aaaa, 0.*mnorm(U,awhr), color=0 oplot, aaaa, dnorm(U,awhr), psym=7, color=0 oplot, [ 65.,aamn], [0.,0.], color=0 oplot, [295.,aamx], [0.,0.], color=0 stdU = stdev( dnorm(U,awhr), meanU ) oplot, [ 65.,295.], [meanU,meanU], linestyle=2, color=0 xyouts, 285., .9*mx, 'mean'+string(meanU,format='(f10.6)'), color=0, align=1. xyouts, 285., .83*mx, 'stdev'+string(stdU,format='(f10.6)'), color=0, align=1. ;Draw V residual plot. !p.region = [ 0.5, 0.0, 1.0, 0.5 ] plot, [65.,295.], [-mx,mx], xstyle=1, ystyle=1, /nodata, /noerase $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth', ytitle='V solar residual' plots, aaaa, 0.*mnorm(V,awhr), color=0 oplot, aaaa, dnorm(V,awhr), psym=8, color=0 oplot, [ 65.,aamn], [0.,0.], color=0 oplot, [295.,aamx], [0.,0.], color=0 stdV = stdev( dnorm(V,awhr), meanV ) oplot, [ 65.,295.], [meanV,meanV], linestyle=2, color=0 xyouts, 285., .9*mx, 'mean'+string(meanV,format='(f10.6)'), color=0, align=1. xyouts, 285., .83*mx, 'stdev'+string(stdV,format='(f10.6)'), color=0, align=1. ;Print a title. xyouts, xsize, 2*ysize-14, mm.file, align=0.5, color=0, /device end ;============================================================================== ; ; procedure: xttx_plot_cal ; ; purpose: plot xttx.f fit to calibration optics or clear port. ; ;============================================================================== pro xttx_plot_cal, ifile, ofile, xfile=xfile, tfile=tfile, date=date ;Check number of parameters. if n_params() ne 2 then begin xttx_plot_cal_usage return end ;Read xttx format data files. xttx_read, ifile, dd & if n_elements(dd) eq 0 then return xttx_read, ofile, mm & if n_elements(mm) eq 0 then return ;Check that i/o data sets are similar. if dd.dim ne mm.dim then begin print print,'xttx_plot_cal: ifile & ofile files differ in length.' return end ;Save device type & system variables. sav_n=!d.name & sav_w=!d.window sav_o=!order & sav_p=!p & sav_x=!x & sav_y=!y & sav_z=!z tvlct, sav_r, sav_g, sav_b, /get ;Set output device type to X windows. set_plot, 'X' !p.font = -1 ;Set special psym=8, filled dot. theta = findgen(16)*(!pi*2/16.) usersym, .8*cos(theta),.8*sin(theta),/fill ;Load black & white color table. loadct, 0 ;Structure to handle windows. ww = $ { w: [-1,-1,-1,-1,-1,-1,-1,-1] $ , click: 'click:(continue)(finish)(exit)' $ , finish: 0 $ , exit: 0 $ } ;Open window for residual plot. xsize=512 & ysize=round(xsize*8.5/11.) window, /free, xsize=xsize, ysize=ysize, xpos=0, ypos=40 $ , title='type in reject fraction (or RETURN)' ww.w(6) = !d.window ;Stokes component array indices. I = 0 Q = 1 U = 2 V = 3 ;Polarization residuals. rsid = sqrt( $ (dd.qqqq/dd.iiii-mm.qqqq/mm.iiii)^2 $ + (dd.uuuu/dd.iiii-mm.uuuu/mm.iiii)^2 $ + (dd.vvvv/dd.iiii-mm.vvvv/mm.iiii)^2 ) ;No initially reject fraction. rjct = 1. dmiss = 0 ;Loop to set reject fraction. repeat begin ;Where residual is satisfactory. gdwhr = where( rsid lt rjct ) print print, 'p-vec_{data} file: '+dd.file print, 'p-vec_{model} file: '+mm.file if gdwhr(0) lt 0 then begin print, 'All vectors rejected' print, 'Reject residual polarization reset to 1.' rjct = 1. gdwhr = where( rsid lt rjct ) endif ;Plot polarization residuals. !p.region = 0. plot, (lindgen(dd.dim))(gdwhr), rsid(gdwhr) $ , background=!d.n_colors-1, color=0 $ , xtitle = 'vector number from 0' $ , ytitle = 'residual = | | p >_{d} - | p >_{m} |' $ , title = mm.file ;Print vector numbers. if n_elements(gdwhr) eq dd.dim then begin print, 'No reject vectors' end else begin print print, 'Reject vectors:' print, ' vec# optics ' $ +'q(data) ' $ +'u(data) ' $ +'v(data) ' $ +' PHI_ma L_ang D_ang' print, ' ' $ +'Q(model) ' $ +'U(model) ' $ +'V(model) ' $ +' residual' awhr = where( rsid ge rjct ) for ipt=0,n_elements(awhr)-1 do begin pt = awhr(ipt) tbl = round( dd.table_angle(pt)-dd.azimuth(pt)-30. ) tbl = (tbl+360+360+360) mod 360 lang = ' ' if dd.ix_linear(pt) ne 0 then $ lang = string((dd.ix_linear(pt)-1)*45,format='(i7)') dang = ' ' if dd.ix_retarder(pt) ne 0 then $ dang = string(((dd.ix_retarder(pt)+1)*45) mod 360 $ ,format='(i7)') print, format='(i6,2x,4i1,3f11.6,i19,a,a)' $ , pt $ , dd.ix_dark(pt), dd.ix_bulb(pt) $ , dd.ix_linear(pt), dd.ix_retarder(pt) $ , [dd.qqqq(pt),dd.uuuu(pt),dd.vvvv(pt)]/dd.iiii(pt) $ , tbl, lang, dang print, format='(12x,3f11.6,f12.6)' $ , [mm.qqqq(pt),mm.uuuu(pt),mm.vvvv(pt)]/mm.iiii(pt) $ , rsid(pt) end end ;Print instructions. print, 'Reject residual polarization is presently ', rjct ;Read reject fraction. tmp = '' read, 'Enter max rejection residual polarization (or RETURN): ', tmp ;Interpret input string. tmp = strcompress(tmp) if tmp eq '' then dmiss=1 else reads,tmp,rjct end until dmiss ;Reset data structures. dd = $ { dim: n_elements(gdwhr) $ , ix_dark: dd.ix_dark(gdwhr) $ , ix_bulb: dd.ix_bulb(gdwhr) $ , ix_linear: dd.ix_linear(gdwhr) $ , ix_retarder: dd.ix_retarder(gdwhr) $ , elevation: dd.elevation(gdwhr) $ , azimuth: dd.azimuth(gdwhr) $ , table_angle: dd.table_angle(gdwhr) $ , iiii: dd.iiii(gdwhr) $ , qqqq: dd.qqqq(gdwhr) $ , uuuu: dd.uuuu(gdwhr) $ , vvvv: dd.vvvv(gdwhr) $ , file: dd.file $ , whr: gdwhr $ , vtype: '| s >' $ } mm = $ { dim: n_elements(gdwhr) $ , ix_dark: mm.ix_dark(gdwhr) $ , ix_bulb: mm.ix_bulb(gdwhr) $ , ix_linear: mm.ix_linear(gdwhr) $ , ix_retarder: mm.ix_retarder(gdwhr) $ , elevation: mm.elevation(gdwhr) $ , azimuth: mm.azimuth(gdwhr) $ , table_angle: mm.table_angle(gdwhr) $ , iiii: mm.iiii(gdwhr) $ , qqqq: mm.qqqq(gdwhr) $ , uuuu: mm.uuuu(gdwhr) $ , vvvv: mm.vvvv(gdwhr) $ , file: mm.file $ } ;Initialize stokes vector storage. ; pvecs: polarimeter vectors. ; pmodel: model polarimeter vectors. ; dvecs: data vectors transformed. ; mvecs: model vectors transformed. pvecs = fltarr( 4, dd.dim ) pmodel = fltarr( 4, dd.dim ) dvecs = fltarr( 4, dd.dim ) mvecs = fltarr( 4, dd.dim ) ;Set polarimeter vectors. pvecs(I,*) = dd.iiii pvecs(Q,*) = dd.qqqq pvecs(U,*) = dd.uuuu pvecs(V,*) = dd.vvvv ;Set model vectors normalized to ;polarimeter intensity. pmodel(I,*) = mm.iiii*dd.iiii/mm.iiii pmodel(Q,*) = mm.qqqq*dd.iiii/mm.iiii pmodel(U,*) = mm.uuuu*dd.iiii/mm.iiii pmodel(V,*) = mm.vvvv*dd.iiii/mm.iiii ;Initial untransform data & model vectors. dd.vtype = '| p >' dvecs = pvecs mvecs = pmodel ;On keyword transform vectors with ;inverse X matrix. if n_elements(xfile) ne 0 then begin ;Read X matrix file. openr, /get_lun, unit, xfile str='' cd=0. & yd=0. & dlc=0. & erar=0. xmtx=fltarr(4,4) readf, unit, str readf, unit, str readf, unit, txx, tyy, format='(f11.0, 4x,f11.0)' readf, unit, dlc, erar, format='(f11.0,11x,f15.0)' readf, unit, xmtx free_lun, unit xmtx = transpose(xmtx) ;Transform data with inverse X matrix. dvecs = invert( xmtx ) # dvecs mvecs = invert( xmtx ) # mvecs print print, 'xttx_plot_cal: p_vectors transformed with inverse X.' print, xfile ;Check if linear polarizer is present. if max(dd.ix_linear) ne 0 then dd.vtype='| c >' else dd.vtype='| t >' ;Apply inverse D matrix if applicable. if max(dd.ix_retarder) ne 0 and dd.vtype eq '| t >' then begin ;Retardence matrix. ang = dlc*!pi/180. sn = sin(ang) cs = cos(ang) rmtx = setmtx( $ 1., 0., 0., 0. $ , 0., 1., 0., 0. $ , 0., 0., cs, sn $ , 0., 0.,-sn, cs ) ;Residual linear polarization matrix. lmtx = setmtx( $ .5*(txx+tyy), .5*(txx-tyy), 0., 0. $ , .5*(txx-tyy), .5*(txx+tyy), 0., 0. $ , 0., 0., sqrt(txx*tyy), 0. $ , 0., 0., 0., sqrt(txx*tyy) ) ;Compute D matrix. dmtx = lmtx # rmtx ;Loop through vectors. for kk=0L,dd.dim-1 do begin if dd.ix_retarder(kk) gt 0 then begin ;Angle CW to fast axis. ang = ((dd.ix_retarder(kk)-1)*45.+90.+erar) * !pi/180. sn = sin(2.*ang) cs = cos(2.*ang) ;Rotation matrix to fast axis. rot_fast = setmtx( $ 1., 0., 0., 0. $ , 0., cs, sn, 0. $ , 0.,-sn, cs, 0. $ , 0., 0., 0., 1. ) ;Inverse rotation matrix from fast axis. rot_inv = setmtx( $ 1., 0., 0., 0. $ , 0., cs,-sn, 0. $ , 0., sn, cs, 0. $ , 0., 0., 0., 1. ) ;Inverse calibration D matrix. dinv = invert( rot_inv # dmtx # rot_fast ) ;Apply inverse calibration D matrix. dvecs(*,kk) = dinv # dvecs(*,kk) mvecs(*,kk) = dinv # mvecs(*,kk) end end print print, 'xttx_plot_cal: cal_vectors transformed with inverse D.' print, xfile end ;Save t vectors for sku residual plot. tvecs = dvecs tmodel = mvecs ;On keyword transform vectors with ;inverse T matrix if applicable. if n_elements(tfile) ne 0 and dd.vtype eq '| t >' then begin ;Loop through vectors. for kk=0L,dd.dim-1 do begin ;Get T matrix. tmtx = get_t( 0 $ , dd.azimuth(kk) $ , dd.elevation(kk) $ , dd.table_angle(kk) $ , tfile ) ;Apply inverse T matrix. dvecs(*,kk) = invert( tmtx ) # dvecs(*,kk) mvecs(*,kk) = invert( tmtx ) # mvecs(*,kk) end dd.vtype = '| s >' print print, 'xttx_plot_cal: t_vectors transformed with inverse T.' print, tfile end end ;Plot data in windows. xttx_plot_discrim, dvecs, mvecs, dd, mm, ww if ww.exit then goto, clean_up ;Plot sku residual plot in solar frame. if dd.vtype eq '| s >' then $ xttx_plot_sku, tvecs, tmodel, dvecs, mvecs, dd, mm, ww, date, tfile ;Wait for click. ww.finish=0 xttx_cal_click, ww ;Clean up clean_up: for wn=0,n_elements(ww.w)-1 do if ww.w(wn) ge 0 then wdelete,ww.w(wn) set_plot, sav_n if sav_n eq 'X' and sav_w ge 0 then wset,sav_w tvlct, sav_r, sav_g, sav_b !order=sav_o & !p=sav_p & !x=sav_x & !y=sav_y & !z=sav_z end