pro xttx_plot_s_usage ;+ ; ; procedure: xttx_plot_s ; ; purpose: compute & plot spectra fitted by xttx.f ; ; routines: xttx_plot_s_usage xttx_s_click xttx_plot_duo xttx_plot_s ; ; author: paul@ncar 11/95 ; ;============================================================================== ; ; Check number of parameters. ; if 1 then begin print print, "usage: xttx_plot_s, ifile, tfile, xfile" print print, " Compute & plot spectra fitted by xttx.f print print, " Arguments:" print, " ifile - path to spectra input to xttx" print, " tfile - path to t matrix file output by xttx" print, " xfile - path to x matrix file output by xttx" print, " ex:" print print, " ifile = '~paul/stokes/src/xttx/test/tinv_cal'" print, " run = '~paul/stokes/src/xttx/test/test1'" print, " xttx_plot_s, ifile, run+'_o_t', run+'_o_x'" print return endif ;- end ;============================================================================== ; ; procedure: xttx_s_click ; ; purpose: wait for user to click window. ; ;============================================================================== pro xttx_s_click, ww, xx, yy ;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 exit. if down eq 4 then ww.exit=1 ;Wait for button up. repeat cursor, xx,yy,/device,/nowait until !err eq 0 end ;============================================================================== ; ; procedure: xttx_plot_duo ; ; purpose: plot p_vectors & s_vectors in duo windows. ; ; ;============================================================================== pro xttx_plot_duo, ww, file, azimuth, p_dat, s_dat, p_mod ;Number of stokes vectors. xdim = n_elements(p_dat)/4 xxxx = indgen(xdim) ;Stokes parameter subscripts. I = 0 Q = 1 U = 2 V = 3 ;Model vectors normalized to ;polarimeter vectors. factor = p_dat(I,*)/p_mod(I,*) normP = replicate( 1., 4, xdim ) normP(Q,*) = p_mod(Q,*)*factor normP(U,*) = p_mod(U,*)*factor normP(V,*) = p_mod(V,*)*factor ;Plot I profile. wset,ww.w(0) !p.region = [ 0.0, 0.5, 0.5, 1.0 ] mx = max( p_dat(I,*) ) plot, [0,xdim-1], [0,mx], /nodata, background=!d.n_colors-1, color=0 plots, xxxx, p_dat(I,*), thick=2., color=0 xyouts, /data, align=1., .8*xdim, .1*mx, 'p_0' $ , charsize=2., charthick=2., color=0 ;Plot Q profiles & residuals. !p.region = [ 0.5, 0.5, 1.0, 1.0 ] mx = max(abs( normP(Q:U,*) )) > max(abs( p_dat(Q:U,*) )) plot, [0,xdim-1], [-mx,mx], /nodata, /noerase $ , background=!d.n_colors-1, color=0 plots, xxxx, p_dat(Q,*), thick=2., color=0 plots, xxxx, normP(Q,*), color=0 plots, xxxx, p_dat(Q,*)-normP(Q,*), color=0 xyouts, /data, align=1., .8*xdim, -mx, 'p_1' $ , charsize=2., charthick=2., color=0 ;Plot U profiles & residuals. !p.region = [ 0.0, 0.0, 0.5, 0.5 ] plot, [0,xdim-1], [-mx,mx], /nodata, /noerase $ , xtitle = 'wavelength pixel number' $ , background=!d.n_colors-1, color=0 plots, xxxx, p_dat(U,*), thick=2., color=0 plots, xxxx, normP(U,*), color=0 plots, xxxx, p_dat(U,*)-normP(U,*), color=0 xyouts, /data, align=1., .8*xdim, -mx, 'p_2' $ , charsize=2., charthick=2., color=0 ;Plot V profiles & residuals. !p.region = [ 0.5, 0.0, 1.0, 0.5 ] mx = max(abs( normP(V,*) )) > max(abs( p_dat(V,*) )) plot, [0,xdim-1], [-mx,mx], /nodata, /noerase $ , background=!d.n_colors-1, color=0 plots, xxxx, p_dat(V,*), thick=2., color=0 plots, xxxx, normP(V,*), color=0 plots, xxxx, 20.*(p_dat(V,*)-normP(V,*)), color=0 xyouts, /data, align=1., .8*xdim, -mx, 'p_3' $ , charsize=2., charthick=2., color=0 xyouts, /data, align=1., .97*xdim, .8*mx $ , '20x residual', color=0 ;Print some info. xyouts, /device, align=.5, !d.x_size/2, !d.y_size-12, file+' p-vec', color=0 xyouts, /device, align=.5, !d.x_size/2, !d.y_size/2, color=0 $ , 'azimuth = '+string(azimuth,format='(f6.1)') ;Plot calibrated I profile. wset,ww.w(1) !p.region = [ 0.0, 0.5, 0.5, 1.0 ] mx = max( s_dat(I,*) ) plot, [0,xdim-1], [0,mx], /nodata, background=!d.n_colors-1, color=0 plots, xxxx, s_dat(I,*), thick=2., color=0 xyouts, /data, align=1., .8*xdim, .1*mx, 's_0' $ , charsize=2., charthick=2., color=0 ;Plot calibrated Q profile. !p.region = [ 0.5, 0.5, 1.0, 1.0 ] mx = max( abs( s_dat(Q:U,*) ) ) plot, [0,xdim-1], [-mx,mx], /nodata, /noerase $ , background=!d.n_colors-1, color=0 plots, xxxx, s_dat(Q,*), thick=2., color=0 plots, xxxx, s_dat(Q,*)*0., color=0 xyouts, /data, align=1., .8*xdim, -mx, 's_1' $ , charsize=2., charthick=2., color=0 ;Plot calibrated U profile. !p.region = [ 0.0, 0.0, 0.5, 0.5 ] plot, [0,xdim-1], [-mx,mx], /nodata, /noerase $ , xtitle = 'wavelength pixel number' $ , background=!d.n_colors-1, color=0 plots, xxxx, s_dat(U,*), thick=2., color=0 plots, xxxx, s_dat(U,*)*0., color=0 xyouts, /data, align=1., .8*xdim, -mx, 's_2' $ , charsize=2., charthick=2., color=0 ;Plot calibrated V profile. !p.region = [ 0.5, 0.0, 1.0, 0.5 ] mx = max( abs( s_dat(V,*) ) ) plot, [0,xdim-1], [-mx,mx], /nodata, /noerase $ , background=!d.n_colors-1, color=0 plots, xxxx, s_dat(V,*), thick=2., color=0 xyouts, /data, align=1., .8*xdim, -mx, 's_3' $ , charsize=2., charthick=2., color=0 ;Print some info. xyouts, /device, align=.5, !d.x_size/2, !d.y_size-12, file+' s-vec', color=0 xyouts, /device, align=.5, !d.x_size/2, !d.y_size/2, color=0 $ , 'azimuth = '+string(azimuth,format='(f6.1)') end ;============================================================================== ; ; procedure: xttx_plot_s ; ; purpose: compute & plot spectra fitted by xttx.f ; ; ;============================================================================== pro xttx_plot_s, ifile, tfile, xfile ;Check number of parameters. if n_params() eq 0 then begin xttx_plot_s_usage return end ;Read spectra data file input to xttx.f. xttx_read, ifile, dd if n_elements(dd) eq 0 then return print, 'p-vec_{data} file: ', ifile ;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 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] $ , click: 'click:(continue)(continue)(exit)' $ , exit: 0 $ } ;Set iquv array indices. I = 0 Q = 1 U = 2 V = 3 ;Initialize arrays for storage. ; pvecs: polarimeter vectors. ; tvecs: telescope vectors. ; svecs: calibrated solar vectors. ; tmodel: model telescope vectors. ; pmodel: model polarimeter vecotrs. pvecs = fltarr( 4, dd.dim ) tvecs = fltarr( 4, dd.dim ) svecs = fltarr( 4, dd.dim ) tmodel = fltarr( 4, dd.dim ) pmodel = fltarr( 4, dd.dim ) ;Set polarimeter vectors. pvecs(I,*) = dd.iiii pvecs(Q,*) = dd.qqqq pvecs(U,*) = dd.uuuu pvecs(V,*) = dd.vvvv ;Read X matrix file. openr, /get_lun, unit, xfile str = '' for ix=1,4 do readf, unit, str xmtx = fltarr(4,4) readf, unit, xmtx free_lun, unit xmtx = transpose(xmtx) ;Calculate telescope vectors. tvecs = invert( xmtx ) # pvecs ;First pixel index of next pixel set. next = 0L ;ASP op count. opdim = 0L ;Loop through profile sets. while next lt dd.dim do begin opdim = opdim+1 first = next ;Scan elevation to find end of profile. while next lt dd.dim $1 and dd.elevation(next < (dd.dim-1)) eq dd.elevation(first) $ do next=next+1 ;Index of last vector in band pass. last = next-1 ;Number of vectors in band pass. ;Index ramp. xdim = last-first+1 xxxx = lindgen(xdim) ;Flag zero if band pass changes. if n_elements(save_xdim) eq 0 then save_xdim=xdim if xdim ne save_xdim then save_xdim = 0 ;Get t matrix. tmtx = get_t( 0 $ ,dd.azimuth(first),dd.elevation(first),dd.table_angle(first) $ ,tfile ) ;Isolate polarimeter & telescope vector ;profiles. vecsP = pvecs(*,first:last) vecsT = tvecs(*,first:last) ;Compute calibrated solar vectors. vecsS = invert( tmtx ) # vecsT svecs(*,first:last) = vecsS ;Assume zero linear polarization model. modelS = vecsS modelS(Q:U,*) = 0. ;Compute model telescope vectors. modelT = tmtx # modelS tmodel(*,first:last) = modelT ;Compute model polarimeter vectors. modelP = xmtx # modelT pmodel(*,first:last) = modelP end ;Return if band pass changed. if save_xdim eq 0 then begin stop, 'Band pass changed; enter .c to continue' goto, clean_up end ;Where arrays for 0th wavelength pixel. whr0 = lindgen(opdim)*xdim ;Sort where array on azimuth. ;;;;;;;;if opdim gt 1 then whr0 = whr0( sort( dd.azimuth(whr0) ) ) ;Initialize structure to hold sub arrays. tmp = fltarr(opdim,3) zz = $ { file: dd.file $ , azimuth: dd.azimuth(whr0) $ , elevation: dd.elevation(whr0) $ , table_angle: dd.table_angle(whr0) $ , pvQ: tmp $ , pvU: tmp $ , pvV: tmp $ , tvQ: tmp $ , tvU: tmp $ , tvV: tmp $ , svQ: tmp $ , svU: tmp $ , svV: tmp $ , tmQ: tmp $ , tmU: tmp $ , tmV: tmp $ , pmQ: tmp $ , pmU: tmp $ , pmV: tmp $ , tsQ: tmp $ , tsU: tmp $ , tsV: tmp $ , tgQ: tmp $ , tgU: tmp $ , tgV: tmp $ } ;Arrays for continuum plot. cctI = fltarr(opdim) cctQ = fltarr(opdim) cctU = fltarr(opdim) cctV = fltarr(opdim) ;Arrays for Q U solar residual plot. rsdQ = fltarr(dd.dim) rsdU = fltarr(dd.dim) ;Loop through ops to get sub arrays. for iop=0,opdim-1 do begin ;Wavelength pixels for op. w0 = whr0(iop) w1 = w0+xdim-1 ;Number of points to average for continuum. npix = 9 hlf = npix/2 ;Where for extrema. tmp = min( pvecs(V,w0:w1), pmn ) & pmn = w0+pmn tmp = max( pvecs(V,w0:w1), pmx ) & pmx = w0+pmx tmp = max( pvecs(I,w0:w1), cmx ) & cmx = w0+cmx whr = [pmn,pmx,cmx] ;Average solar signal in continuum. c0 = w0 > (cmx-hlf) < (w1-npix+1) c1 = c0+npix-1 cctI(iop) = total(svecs(I,c0:c1))/npix cctQ(iop) = total(svecs(Q,c0:c1))/svecs(I,cmx)/npix cctU(iop) = total(svecs(U,c0:c1))/svecs(I,cmx)/npix cctV(iop) = total(svecs(V,c0:c1))/svecs(I,cmx)/npix ;Solar Q U residuals. r0 = iop*xdim r1 = r0+xdim-1 rsdQ(r0:r1) = svecs(Q,w0:w1)/svecs(I,cmx) rsdU(r0:r1) = svecs(U,w0:w1)/svecs(I,cmx) zz.pvQ(iop,*) = pvecs(Q,whr)/pvecs(I,cmx) zz.pvU(iop,*) = pvecs(U,whr)/pvecs(I,cmx) zz.pvV(iop,*) = pvecs(V,whr)/pvecs(I,cmx) zz.tvQ(iop,*) = tvecs(Q,whr)/tvecs(I,cmx) zz.tvU(iop,*) = tvecs(U,whr)/tvecs(I,cmx) zz.tvV(iop,*) = tvecs(V,whr)/tvecs(I,cmx) zz.svQ(iop,*) = svecs(Q,whr)/svecs(I,cmx) zz.svU(iop,*) = svecs(U,whr)/svecs(I,cmx) zz.svV(iop,*) = svecs(V,whr)/svecs(I,cmx) zz.tmQ(iop,*) = tmodel(Q,whr)/tmodel(I,cmx) zz.tmU(iop,*) = tmodel(U,whr)/tmodel(I,cmx) zz.tmV(iop,*) = tmodel(V,whr)/tmodel(I,cmx) zz.pmQ(iop,*) = pmodel(Q,whr)/pmodel(I,cmx) zz.pmU(iop,*) = pmodel(U,whr)/pmodel(I,cmx) zz.pmV(iop,*) = pmodel(V,whr)/pmodel(I,cmx) end ;Open workstation window. xsz=660 & ysz=round(xsz*11./8.5) window, /free, xsize=xsz, ysize=ysz $ , xpos=0, ypos=900-ysz $ , title=ww.click ww.w(2) = !d.window ;Set margin sizes as if they were inches. lft = 7./8. rgt = 5./16. btm = 5./8. top = 3./4. xgp = 1./8. ygp = 1./4. xfr = (8.5-lft-xgp-rgt)/2. yfr = (11.-btm-ygp-ygp-top)/3. ;Loop over wave length pixels. for ipx=0L,2L do begin pvQ = zz.pvQ(*,ipx) & pmQ = zz.pmQ(*,ipx) pvU = zz.pvU(*,ipx) & pmU = zz.pmU(*,ipx) pvV = zz.pvV(*,ipx) & pmV = zz.pmV(*,ipx) tvQ = zz.tvQ(*,ipx) & tmQ = zz.tmQ(*,ipx) tvU = zz.tvU(*,ipx) & tmU = zz.tmU(*,ipx) tvV = zz.tvV(*,ipx) & tmV = zz.tmV(*,ipx) svQ = zz.svQ(*,ipx) svU = zz.svU(*,ipx) svV = zz.svV(*,ipx) ;Average solar vector for smooth telescope ;plot. std_q = stdev( svQ, avg_q ) std_u = stdev( svU, avg_u ) std_v = stdev( svV, avg_v ) print print, 'average solar vector: 1.0', avg_q, avg_u, avg_v $ , format = '(a,3f14.9)' print, ' standard deviation: 0.0', std_q, std_u, std_v $ , format = '(a,3f14.9)' chr_v='"data"{symbols} vs T*(1,0,0,'+stringit(avg_v)+'){line}' ;Initialize to average solar vectors. iquv = fltarr(4,opdim) iquv(I,*) = 1. iquv(V,*) = avg_v ;Loop to compute instrument vectors. for idx=0,opdim-1 do begin ;Get t matrix. tmtx = get_t( 0 $ , zz.azimuth(idx), zz.elevation(idx), zz.table_angle(idx) $ , tfile ) ;Compute telescope vector. iquv(*,idx) = tmtx # iquv(*,idx) end ;Normalize to intensity. zz.tsQ(*,ipx) = iquv(Q,*)/iquv(I,*) & tsQ = zz.tsQ(*,ipx) zz.tsU(*,ipx) = iquv(U,*)/iquv(I,*) & tsU = zz.tsU(*,ipx) zz.tsV(*,ipx) = iquv(V,*)/iquv(I,*) & tsV = zz.tsV(*,ipx) ;Set curves with goodness of fit. zz.tgQ(*,ipx) = tsQ+tvQ-tmQ & tgQ = zz.tgQ(*,ipx) zz.tgU(*,ipx) = tsU+tvU-tmU & tgU = zz.tgU(*,ipx) zz.tgV(*,ipx) = tsV+tvV-tmV & tgV = zz.tgV(*,ipx) ;Set y scale. if ipx eq 2 then kk=2 else kk=[0,1] mx = max(abs( $ [ zz.pvQ(*,kk), zz.pvU(*,kk), zz.pvV(*,kk) $ , zz.pmQ(*,kk), zz.pmU(*,kk), zz.pmV(*,kk) $ , zz.tsQ(*,kk), zz.tsU(*,kk), zz.tsV(*,kk) $ , zz.tgQ(*,kk), zz.tgU(*,kk), zz.tgV(*,kk) $ ] )) yrange = 1.05*[-mx,mx] ;Set xscale. xrange = [65.,295.] ;Set x anointation. if ipx eq 0 then begin xtitle='azimuth (degrees)' xcharsize = 1. end else begin xtitle='' xcharsize = .00001 end if ipx eq 2 $ then title='p-vec: data{symbols}, model{line}' $ else title='' ;Draw plot frame for polarimeter vectors. !p.region = 0. plot, /nodata, noerase=ipx<1, background=!d.n_colors-1, color=0 $ , xrange, xstyle=1, xtitle=xtitle, xcharsize=xcharsize $ , yrange, ystyle=1, ytitle='polarization' $ , position = $ [ lft, btm+ipx*(yfr+ygp), lft+xfr, btm+ipx*(yfr+ygp)+yfr ] $ / [8.5,11.,8.5,11.] $ , title = title ;Plot observed data. oplot, zz.azimuth, pvQ, psym=1, color=0 oplot, zz.azimuth, pvU, psym=7, color=0 oplot, zz.azimuth, pvV, psym=8, color=0 ;Plot model fit. plots, zz.azimuth, pmQ, color=0 plots, zz.azimuth, pmU, color=0 plots, zz.azimuth, pmV, color=0 ;Plot magnitude of residual. ;plots, zz.azimuth $ ;, sqrt((pvQ-pmQ)^2+(pvU-pmU)^2+(pvV-pmV)^2) $ ;, color=0 ;Draw plot frame for telescope vectors. plot, /nodata, /noerase, background=!d.n_colors-1, color=0 $ , xrange, xstyle=1, xtitle=xtitle, xcharsize=xcharsize $ , yrange, ystyle=1, ycharsize=.00001 $ , position = $ [ lft+xfr+xgp, btm+ipx*(yfr+ygp) $ , lft+xfr+xgp+xfr, btm+ipx*(yfr+ygp)+yfr ] $ / [8.5,11.,8.5,11.] $ , title=chr_v ;Plot observed data. oplot, zz.azimuth, tgQ, psym=1, color=0 oplot, zz.azimuth, tgU, psym=7, color=0 oplot, zz.azimuth, tgV, psym=8, color=0 ;Plot model fit. plots, zz.azimuth, tsQ, color=0 plots, zz.azimuth, tsU, color=0 plots, zz.azimuth, tsV, color=0 plots, zz.azimuth, sqrt(tsQ^2+tsU^2), color=0 ;Plot magnitude of residual. ;plots, zz.azimuth, sqrt((tgQ-tsQ)^2+(tgU-tsU)^2+(tgV-tsV)^2), color=0 end ;Print spectra file path. xyouts, (8.5-rgt)/8.5, (11.-top+1./4.)/11., zz.file $ , color=0, /normal, align=1. ;Open window for I_c Q_c U_c V_c plot. xsize=320 & ysize=round(xsize*11./8.5) window, /free, xsize=2*xsize, ysize=2*ysize, xpos=1130-2*xsize, ypos=40 $ , title = 'Continuum' ww.w(3) = !d.window ;Draw plot frame. !p.region = [ 0.0, 0.5, 0.5, 1.0 ] plot, [65.,295.], [0.,max(cctI)], xstyle=1, /nodata $ , background=!d.n_colors-1, color=0 $ , xtitle='azimuth (degrees)' $ , ytitle='I _ c (continuum)' $ , title='Solar Continuum Signal' oplot, zz.azimuth, cctI, color=0, thick=2. oplot, zz.azimuth, cctI, psym=8, color=0, thick=2. mx = 1.3*max( [ abs(cctQ), abs(cctU), abs(cctV) ] ) ;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 _ c / I _ c' $ , title=zz.file oplot, zz.azimuth, cctQ, psym=1, color=0 oplot, [65.,295.], [0.,0.], color=0 stdQ = stdev( cctQ, 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 _ c / I _ c' oplot, zz.azimuth, cctU, psym=7, color=0 oplot, [65.,295.], [0.,0.], color=0 stdU = stdev( cctU, 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 _ c / I _ c' oplot, zz.azimuth, cctV, psym=8, color=0 oplot, [65.,295.], [0.,0.], color=0 stdV = stdev( cctV, 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. ;X size of residual plot. rsize = 720 ;Position in window for risduals plot. rposition= [.08,.12,1.-.02,1.-.06] ;Open window for residual plot. if 0 then begin window, /free, xsize=rsize, ysize=rsize/2, xpos=100, ypos=0 $ , title = 'click (profile set, profile set, exit)' ww.w(4) = !d.window ;Solar residuals. rsd = sqrt((svecs(Q,*)/svecs(I,*))^2+(svecs(U,*)/svecs(I,*))^2) mx = 1.3*max(rsd) ;Draw plot frame. plot, [0,dd.dim-1], [0.,mx], xstyle=1, ystyle=1., /nodata $ , position = rposition $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' (dashed lines separate operations)' $ , xtitle='vector number', ytitle='residual = | | s >_{d} - | s >_{m} |' ;Plot residuals. plots, indgen(dd.dim), rsd, color=0 ;Loop over operations. for xop=0,dd.dim-xdim,xdim do begin ;Dashed lines separate operations. if xop ne 0 then $ oplot, [xop,xop], [mx/2.,mx], linestyle=2, color=0 ;Anotate operation with azimuth. azstr = string(dd.azimuth(xop),format='(f5.1)') if (xop mod 4) eq 1 then azstr='azimuth '+azstr xyouts, xop+.5*xdim, .97*mx, azstr, color=0, align=1., orient=90. end end ;Open window for residual plot. if 1 then begin window, /free, xsize=rsize, ysize=rsize/2, xpos=100, ypos=0 $ , title = 'click (profile set, profile set, exit)' ww.w(4) = !d.window mx = ( max(abs(rsdQ)) > max(abs(rsdU)) ) ;Positions in windows for risdual plot. qposition= [ .08, (1.-.06)/2, 1.-.02, 1.-.06 ] uposition= [ .08, 0., 1.-.02, (1.-.06)/2 ] ;Draw plot frame. plot, [0,dd.dim-1], [-mx,1.3*mx], xstyle=1, ystyle=1., /nodata $ , position = qposition $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' (dashed lines separate operations)' $ , xcharsize = .000001, xticklen=.000001 $ , ytitle='residual s_1' ;Plot residuals. plots, indgen(dd.dim), rsdQ, color=0 oplot, [0,dd.dim-1], [0,0], color=0 ;Loop over operations. for xop=0,dd.dim-xdim,xdim do begin ;Dashed lines separate operations. if xop ne 0 then $ oplot, [xop,xop], [mx/3.,1.3*mx], linestyle=2, color=0 ;Anotate operation with azimuth. azstr = string(dd.azimuth(xop),format='(f5.1)') if (xop mod 4) eq 1 then azstr='AZ '+azstr xyouts, xop+.5*xdim, .97*1.3*mx, azstr, color=0, align=1., orient=90. end ;Draw plot frame. plot, [0,dd.dim-1], [-1.3*mx,mx], xstyle=1, ystyle=1., /nodata, /noerase $ , position = uposition $ , background=!d.n_colors-1, color=0 $ , xcharsize = .000001, xticklen=.000001 $ , ytitle='residual s_2' ;Plot residuals. plots, indgen(dd.dim), rsdU, color=0 oplot, [0,dd.dim-1], [0,0], color=0 ;Loop over operations. for xop=0,dd.dim-xdim,xdim do begin ;Dashed lines separate operations. if xop ne 0 then $ oplot, [xop,xop], [-.97*1.3*mx,-mx/3], linestyle=2, color=0 ;Anotate operation with azimuth. azstr = string(dd.azimuth(xop),format='(f5.1)') if (xop mod 4) eq 1 then azstr=azstr + ' AZ' xyouts, xop+.5*xdim, - .97*1.3*mx, azstr, color=0, align=0., orient=90. end end ;Open window for residual plot. if 0 then begin window, /free, xsize=720, ysize=360, xpos=200, ypos=200 $ , title = 'residuals' ww.w(4) = !d.window ;Polarimeter residuals. rsd = sqrt( $ (pvecs(Q,*)/pvecs(I,*)-pmodel(Q,*)/pmodel(I,*))^2 $ + (pvecs(U,*)/pvecs(I,*)-pmodel(U,*)/pmodel(I,*))^2 $ + (pvecs(V,*)/pvecs(I,*)-pmodel(V,*)/pmodel(I,*))^2 ) mx = 1.3*max(rsd) ;Draw plot frame. plot, [0,dd.dim-1], [0.,mx], xstyle=1, ystyle=1., /nodata $ , position = rposition $ , background=!d.n_colors-1, color=0 $ , title=dd.file+' (dashed lines separate operations)' $ , xtitle='vector number', ytitle='residual = | | p >_{d} - | p >_{m} |' ;Plot residuals. plots, indgen(n_elements(rsd)), rsd, color=0 ;Loop over operations. for xop=0,dd.dim-xdim,xdim do begin ;Dashed lines separate operations. if xop ne 0 then $ oplot, [xop,xop], [mx/2.,mx], linestyle=2, color=0 ;Anotate operation with azimuth. azstr = string(dd.azimuth(xop),format='(f5.1)') if (xop mod 4) eq 1 then azstr='azimuth '+azstr xyouts, xop+.5*xdim, .97*mx, azstr, color=0, align=1., orient=90. end end ;Loop to plot profile sets. xpfl = -1 while 1 do begin ;Wait for window click. xttx_s_click, ww, xx, yy w_clicked = !d.window ;Open windows for profile set plots. if xpfl eq -1 then begin xsize=256 & ysize=xsize window, /free, xsize=2*xsize, ysize=2*ysize $ , xpos=0, ypos=900-2*ysize $ , title=ww.click ww.w(0) = !d.window window, /free, xsize=2*xsize, ysize=2*ysize $ , xpos=2*xsize+20, ypos=900-2*ysize $ , title=ww.click ww.w(1) = !d.window end ;Step profile number. xpfl = xpfl+1 ;Set profile set to plot corresponding ;to residual window click. if w_clicked eq ww.w(4) then begin xpfl = fix( opdim*(xx-rposition(0)*rsize) $ / ((rposition(2)-rposition(0))*rsize) ) xpfl = 0 > xpfl < (opdim-1) end ;Exit if right button clicked or stepped ;off end of profile set. if ww.exit or xpfl ge opdim then goto, clean_up ;Plot profile set xop = xpfl*xdim xttx_plot_duo, ww, dd.file, dd.azimuth(xop) $ , pvecs(*,xop:xop+xdim-1), svecs(*,xop:xop+xdim-1) $ , pmodel(*,xop:xop+xdim-1) end ;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