pro EIS_JITTER, index, off, a_only= a_only, b_only= b_only, model= model, $ xrt_pix= xrt_pix, read_hk2= read_hk2, dir= dir, verbose=verbose, $ ref= ref, tdir= tdir ;+ ; NAME: ; EIS_JITTER ; INPUTS: ; index - eis_windata structure (or any structure with ; time_ccsds and exposure_time elements) ; ; OPTIONAL INPUTS: ; model - the file name of the XRT orbital variation model ; such as 'alignment_model_sine_070325.txt' ; If not specified, automatically search the file closest to the time of ; the dataset. ; dir - the directory where HK2 or alignment FITS files are located ; ref - the image frame which should be referenced for UFSS bias ; correction. Default is the 1st image (ref=0). If the 1st ; image was acquired within 5 minutes after the satellite ; repointing, it may be suitable to specify the frame which ; was exposured 15 minutes or later after the repointing. ; tdir - (for test purpose) directrory where the model/table ; files are stored. ; OUTPUTS: ; off - fltarr[2, n] n= number of obs points ; this array contains shift values in arcsec ; off[0,*] - E-W direction ; off[1,*] - N-S direction ; For removing jitter and orbital variation, perform ; data_ca= shift_img(data_x, off) ; ; OPTIONAL KEYWORDS: ; a_only, b_only - Specify a_only if UFSS_A only used. b_only if UFSS_B ; only. Default is the average of UFSS-A result and UFSS-B result. ; xrt_pix - Unit of the output is in XRT pixel scale instead of arcsec. ; read_hk2 - read original Hinode HK2 FITS files instead of alignment ; FITS files ; verbose - if set, plot AOL/ATT/UFSS time behaviors. ; ; NOTES: ; 0) See Shimizu et al. PASJ 2007 Hinode special issue for alignment details. ; 1) If the observation was made with satellite attitude tracking mode, ; you may see a very small drift in the long-term XRT cube after ; correcting with this routine. This is because the tracking curve ; used in satellite attitude is identical to the propoer motion of ; the sunspots seen in the field of view. ; 2) This routine requires that the user has the ; $SSWDB/hinode/gen/alignment directory. This is not part of the ; standard EIS Solarsoft distribution. ; ; HISTORY: ; 2007/05/02 T.Shimizu Initial release (under evaluation) ; 2007/09/25 T.Shimizu add read_hk2 keyword ; use alignment FITS as the default ; use UFSS angle with polynominal formula ; correct UFSS bias with AOL=0 value -> incorrect ; add verbose keyword ; 2008/04/04 T.Shimizu correct XRT vs UFSS scale conversion ; correct UFSS bias correction ; add tdir for test purpose ; 2008/07/18 T.Shimizu correct the header information. No ; program change ; 2008-Dec-17 J. Mariska: an EIS version that uses times of EIS obs. Also ; cleaned our a lot of commented out code. ; 2012-Sep-18 P. Young added check to see if user has $SSWDB directory. ;- print,'% EIS_JITTER: Satellite jitter and orbital variation for Hinode EIS data' ;print,'% EIS_JITTER: Questions and bug reports to shimizu@solar.isas.jaxa.jp' ; number of EIS observation points num = n_elements(index.time_ccsds) aol_x= fltarr(num) uf_bx_x= fltarr(num) uf_by_x= fltarr(num) uf_ax_x= fltarr(num) uf_ay_x= fltarr(num) att_x_x= fltarr(num) att_y_x= fltarr(num) ; ; Reading HK data from Hinode alignment FITS database ; per day, can be reduced if read more limited time ; ; reading HK2 FITS files if keyword_set(read_hk2) then begin if keyword_set(dir) then dir0= dir else dir0= '/soda/hinode/cmn/HK2/' ; this is the environment at ISAS/JAXA ds= anytim2utc(index.time_ccsds[0]) de= anytim2utc(index.time_ccsds[num-1]) days= de.mjd-ds.mjd ; for next days in= index.time_ccsds[0] yy= strmid(in, 0, 4) mm= strmid(in, 5, 2) dd= strmid(in, 8, 2) mon_dir= yy + '_' + mm +'/' date_f= 'HK2-'+ yy + mm + dd date0= strmid(yy,2,2) + mm + dd files= file_search(dir0 + mon_dir + date_f + '*.fits.gz', count= count) if count eq 0 then begin print, 'No HK2 FITS found: ' + date_f + '. Returned.' return endif else hk2= mrdfits(files[0], 1, head) if count gt 1 then $ for i=1, count-1 do hk2= [hk2, mrdfits(files[i], 1, head)] ; ### Set first yy/mm/dd as **0 ### yy0 = yy mm0 = mm dd0 = dd ; ########## if days ge 1 then begin for j=1, days do begin ;ds= ymd2jd(yy,mm,dd) ds= ymd2jd(yy0,mm0,dd0) jd2ymd, ds+j, y, m, d yy= strcompress(y, /remove_all) mm= strcompress(m, /remove_all) if m lt 10 then mm= '0'+ mm dd= strcompress(d, /remove_all) if d lt 10 then dd= '0'+ dd mon_dir= yy + '_' + mm +'/' date_f= 'HK2-'+ yy + mm + dd files= file_search(dir0 + mon_dir + date_f + '*.fits.gz', count=count) if count eq 0 then begin print, 'No HK2 FITS found: ' + date_f + '. Just warning.' endif else begin for i=0, n_elements(files)-1 do hk2= [hk2, mrdfits(files[i], 1, head)] endelse endfor endif aol= hk2.hk2_aol att_x= hk2.hk2_att_x * 60 * 60 uf_ax= hk2.HK2_UFSS_A_XANG * 60 * 60 uf_bx= hk2.HK2_UFSS_B_XANG * 60 * 60 att_y= hk2.hk2_att_y * 60 * 60 uf_ay= hk2.HK2_UFSS_A_YANG * 60 * 60 uf_by= hk2.HK2_UFSS_B_YANG * 60 * 60 buttime= '1-jan-2000 0:00:00' ti= hk2.PACKET_EDITION_TIME ; < 2 msec deviation may exit ; better to use hk2.hk2_ti, but need to have a logic for ti reset ; end of HK2 reading endif else begin ; start of the default alignment FITS files in SSWDB if keyword_set(dir) then dir0= dir else begin ;dir0= '/soda/hinode/cmn/alignment/' dir0= getenv('SSWDB') + '/hinode/gen/alignment/' tst=file_test(dir0,/dir) IF tst NE 1 THEN BEGIN print,'' print,'%EIS_JITTER: the directory $SSWDB/hinode/gen/alignment was not found' print,' on your computer. Please check your Solarsoft installation.' print,' Returning...' return ENDIF endelse ds= anytim2utc(index.time_ccsds[0]) de= anytim2utc(index.time_ccsds[num-1]) days= de.mjd-ds.mjd ; for next days in= index.time_ccsds[0] yy= strmid(in, 0, 4) mm= strmid(in, 5, 2) dd= strmid(in, 8, 2) mon_dir= yy + '_' + mm +'/' date_f= 'alignment-'+ yy + mm + dd date0= strmid(yy,2,2) + mm + dd ; date0 is used for getting the propoer orbital variation file files= file_search(dir0 + mon_dir + date_f + '.fits.gz', count=count) if count eq 0 then begin print, 'No alignment FITS found: ' + date_f +'. Returned.' return endif else begin print, 'Reading ', files[0] hk2= mrdfits(files[0], 1, head) endelse ; ### Set first yy/mm/dd as **0 ### yy0 = yy mm0 = mm dd0 = dd ; ########## if days ge 1 then begin for j=1, days do begin ;ds= ymd2jd(yy,mm,dd) ds= ymd2jd(yy0,mm0,dd0) jd2ymd, ds+j, y, m, d yy= strcompress(y, /remove_all) mm= strcompress(m, /remove_all) if m lt 10 then mm= '0'+ mm dd= strcompress(d, /remove_all) if d lt 10 then dd= '0'+ dd mon_dir= yy + '_' + mm +'/' date_f= 'alignment-'+ yy + mm + dd files= file_search(dir0 + mon_dir + date_f + '*.fits.gz', count=count) if count eq 0 then begin print, 'No alignment FITS found: ' + date_f +'. Returned.' return endif else begin for i=0, n_elements(files)-1 do begin print, 'Reading ', files[i] hk2= [hk2, mrdfits(files[i], 1, head)] endfor endelse endfor endif aol= hk2.hk2_aol att_x= hk2.hk2_att_x * 60 * 60 uf_ax= hk2.MRG_UFSS_A_XANG * 60 * 60 uf_bx= hk2.MRG_UFSS_B_XANG * 60 * 60 att_y= hk2.hk2_att_y * 60 * 60 uf_ay= hk2.MRG_UFSS_A_YANG * 60 * 60 uf_by= hk2.MRG_UFSS_B_YANG * 60 * 60 buttime= '1-jan-2000 0:00:00' ti= hk2.hk2_time2000 ; HK2_TIME2000 is seconds since 1-Jan-2000 0:00:00UT. ; This is generated from HK2_TI ; endelse ; end of alignment FITS reading ; * * * EIS changes here * * * ; Provide data as seconds since first index entry ; and everything should work. ; ; (I think the second statemnt is incorrect JTM) ;mdpclk_x= ulong(index.obt_time)/512. - ulong(index[0].obt_time)/512. ;obtend_x= ulong(index.obt_end)/512. - ulong(index[0].obt_end)/512. mdpclk_x = anytim2tai(index.time_ccsds) - anytim2tai(index.time_ccsds[0]) obtend_x = anytim2tai(index.time_ccsds) + index.exposure_time $ - anytim2tai(index.time_ccsds[0]) basetime= index.time_ccsds[0] a= anytim2jd(basetime) b= anytim2jd(buttime) d_t= (a.int + a.frac - b.int - b.frac)*24*60*60 corr_ti= mdpclk_x + d_t mdpclk_x= corr_ti obtend_x= obtend_x + d_t ; ; UFSS and AOL at each of XRT observation time ; for i=0, num-1 do begin x= where(ti ge mdpclk_x(i) and ti le obtend_x(i), count) if count le 0 then begin x= where(ti ge mdpclk_x(i)) aol_x(i)= aol(x[0]) uf_bx_x[i]= average(uf_bx(x[0])) uf_by_x[i]= average(uf_by(x[0])) uf_ax_x[i]= average(uf_ax(x[0])) uf_ay_x[i]= average(uf_ay(x[0])) att_x_x[i]= average(att_x(x[0])) att_y_x[i]= average(att_y(x[0])) dt= ti(x[0]) - mdpclk_x[i] if dt gt 2 then $ print, 'Warning: AOL and UFSS data at the time deviated 10 sec or more ',$ 'No= ', i, ' dT(sec)= ', dt endif else begin aol_x(i)= average(aol[x]) uf_bx_x[i]= average(uf_bx[x]) uf_by_x[i]= average(uf_by[x]) uf_ax_x[i]= average(uf_ax[x]) uf_ay_x[i]= average(uf_ay[x]) att_x_x[i]= average(att_x[x]) att_y_x[i]= average(att_y[x]) endelse endfor ; ---------------------------------------------------- ; UFSS angle calibration ; The default UFSS angle is the result of the convertion ; with the simple formula. ; For alignment, it is more suitable to use the conversion ; with the polynominal formula, which is used on board ; the attitude control system. ; ---------------------------------------------------- ang= ufss_cal_ang(uf_ax_x/60./60., uf_ay_x/60./60., 'A') uf_ax_x= reform(ang[0,*]) uf_ay_x= reform(ang[1,*]) ang= ufss_cal_ang(uf_bx_x/60./60., uf_by_x/60./60., 'B') uf_bx_x= reform(ang[0,*]) uf_by_x= reform(ang[1,*]) ; ; Coordinate information ; ATT_X ; N-S direction large positive value when looking at Solar South limb ; increase value when the sun move north-ward ; UFSS_*_XANG ; N-S direction large positive value when looking at Solar South limb ; increase value when the sun move north-ward ; XRT ; Y center address ; increase value when the sun move north-ward ; SOT ; North limb moves upward (increase value) when the sun move north-ward ; ; ATT_Y ; E-W direction large positive value when looking at Solar East limb ; increase value when the sun move west-ward ; UFSS_*_YANG ; E-W direction large positive value when looking at Solar East limb ; increase value when the sun move west-ward ; XRT ; X center address ; increase value when the sun move west-ward ; SOT ; East limb moves rightward (increase value) when the sun move west-ward ; ; Give a caution to users, who want to co-align the data acquired ; in early and last period of eclipse seasons. ; mmdd= long(strmid(date0, 2, 4)) if mmdd ge 501 and mmdd lt 516 then begin print, '%%%% WARNING %%%%' print, ' Coalignment property changed raipdly in 5/1 - 5/16.' endif if mmdd ge 725 and mmdd lt 805 then begin print, '%%%% WARNING %%%%' print, ' Coalignment property changed raipdly in 7/25 - 8/5.' endif if keyword_set(model) then begin ; -------------------------------------------- ; Reading orbital variation model files ; -------------------------------------------- xrt_coalign_model, date0, c_xrt_n_a, c_xrt_n_b, c_xrt_e_a, c_xrt_e_b, $ file= model, dir0= tdir ; -------------------------------------------- ; Offset from Time variation model (in XRT pixel) ; -------------------------------------------- ; N-S direction c= c_xrt_n_a xx= aol_x fourier_s3,xx*!dtor,c,xrt_fit_ax fourier_s3,[0]*!dtor,c,off xrt_fit_ax= xrt_fit_ax -off[0] c= c_xrt_n_b xx= aol_x fourier_s3,xx*!dtor,c,xrt_fit_bx fourier_s3,[0]*!dtor,c,off xrt_fit_bx= xrt_fit_bx -off[0] ; E-W direction c= c_xrt_e_a xx= aol_x fourier_s3,xx*!dtor,c,xrt_fit_ay fourier_s3,[0]*!dtor,c,off xrt_fit_ay= xrt_fit_ay -off[0] c= c_xrt_e_b xx= aol_x fourier_s3,xx*!dtor,c,xrt_fit_by fourier_s3,[0]*!dtor,c,off xrt_fit_by= xrt_fit_by -off[0] endif else begin ; -------------------------------------------- ; Reading orbital variation table files ; -------------------------------------------- xrt_coalign_table, date0, fit_aol, $ xrt_fit_ax0, xrt_fit_bx0, xrt_fit_ay0, xrt_fit_by0, $ dir0= tdir ; -------------------------------------------- ; Offset from Time variation table (in XRT pixel) ; -------------------------------------------- ; N-S direction xx= round(aol_x) xrt_fit_ax= xrt_fit_ax0[xx] xrt_fit_bx= xrt_fit_bx0[xx] ; E-W direction xrt_fit_ay= xrt_fit_ay0[xx] xrt_fit_by= xrt_fit_by0[xx] end ; ----------------------------------------------------- ; Satellite jitter ; Removing the solar rotation term from UFSS signal ; UFSS signal is used to remove satellite jitter, but it also ; contain the drift due to the solar rotation ; ----------------------------------------------------- ; UFSS offset value (values when watching the disk center address) ;uf_ax_of= -92.0 uf_ax_of= -82.0 uf_bx_of= 78.0 ;uf_ay_of= 70.0 uf_ay_of= 67.0 uf_by_of= -23.0 uf_ax_xc= uf_ax_x - att_x_x - uf_ax_of uf_bx_xc= uf_bx_x - att_x_x - uf_bx_of uf_ay_xc= uf_ay_x - att_y_x - uf_ay_of uf_by_xc= uf_by_x - att_y_x - uf_by_of ; ------------------------------------------------------ ; Correct the UFSS bias ; Zero value at AOL=0 deg ; ------------------------------------------------------ ; ; case 1) aligned to the reference frame (at n) ; This case should be used if the sunspot feature position ; in the reference frame is used to coaligned to SOT data ; ; ; UFSS values at n if n_elements(ref) eq 0 then n=0 else n= ref uf_ax_bi= uf_ax_xc[n] uf_ay_bi= uf_ay_xc[n] uf_bx_bi= uf_bx_xc[n] uf_by_bi= uf_by_xc[n] ; shift due to orbital drift at reference image aol_bi= round(aol_x[n]) if keyword_set(model) then begin fourier_s3,[0,aol_bi]*!dtor,c_xrt_n_a,xrt_fit_ax0 fourier_s3,[0,aol_bi]*!dtor,c_xrt_n_b,xrt_fit_ay0 fourier_s3,[0,aol_bi]*!dtor,c_xrt_e_a,xrt_fit_bx0 fourier_s3,[0,aol_bi]*!dtor,c_xrt_e_b,xrt_fit_by0 fit_ax_bi= - (xrt_fit_ax0[1]-xrt_fit_ax0[0])*1.031 fit_ay_bi= - (xrt_fit_ay0[1]-xrt_fit_ay0[0])*1.031 fit_bx_bi= - (xrt_fit_bx0[1]-xrt_fit_bx0[0])*1.031 fit_by_bi= - (xrt_fit_by0[1]-xrt_fit_by0[0])*1.031 endif else begin ;2008/04/04 xrt_fit_ax0 in XRT pixel ; fit_ax_bi should be in arcsec fit_ax_bi= - (xrt_fit_ax0[aol_bi])*1.031 fit_ay_bi= - (xrt_fit_ay0[aol_bi])*1.031 fit_bx_bi= - (xrt_fit_bx0[aol_bi])*1.031 fit_by_bi= - (xrt_fit_by0[aol_bi])*1.031 end uf_ax_bi= uf_ax_bi - fit_ax_bi uf_ay_bi= uf_ay_bi - fit_ay_bi uf_bx_bi= uf_bx_bi - fit_bx_bi uf_by_bi= uf_by_bi - fit_by_bi if keyword_set(verbose) then begin print, '% UFSS Bias used is below (arcsec):' a= anytim2ints(buttime, offset=mdpclk_x[n]) a= anytim(a, /vms) if not keyword_set(read_hk2) then print, ' '+ a print, ' UFSS_A_X= ', uf_ax_bi print, ' UFSS_B_X= ', uf_bx_bi print, ' UFSS_A_Y= ', uf_ay_bi print, ' UFSS_B_Y= ', uf_by_bi endif ; Correct UFSS bias ; uf_ax_xc= uf_ax_xc - uf_ax_bi uf_bx_xc= uf_bx_xc - uf_bx_bi uf_ay_xc= uf_ay_xc - uf_ay_bi uf_by_xc= uf_by_xc - uf_by_bi ; -------------------------------------------------------- ; plotting the offset history ; -------------------------------------------------------- if keyword_set(verbose) then begin ; UTplot on the screen ts= anytim2ints(buttime, offset=mdpclk_x[0]) te= anytim2ints(buttime, offset=mdpclk_x[n_elements(mdpclk_x)-1]) x= where(ti ge mdpclk_x(0)) xs= x[0] ; start x= where(ti ge mdpclk_x[n_elements(mdpclk_x)-1]) xe= x[0] time = anytim2ints(buttime, offset=ti) window, xs=600, ys=800, /free !p.multi=[0, 2, 4] !x.charsize=1.4 !y.charsize=1.4 !p.charsize=1.4 utplot, time[xs:xe], att_x[xs:xe], /xst, /yst, $ title='Raw ATT_X in the period', $ ytitle='arcsec', timerange=[ts,te] utplot, time[xs:xe], att_y[xs:xe], /xst, /yst, $ title='Raw ATT_Y in the period', $ ytitle='arcsec', timerange=[ts,te] utplot, time[xs:xe], uf_ax[xs:xe], /xst, /yst, $ title='Raw UFSS_A_XANG in the period', $ ytitle='arcsec', timerange=[ts,te] utplot, time[xs:xe], uf_ay[xs:xe], /xst, /yst, $ title='Raw UFSS_A_YANG in the period', $ ytitle='arcsec', timerange=[ts,te] utplot, time[xs:xe], uf_bx[xs:xe], /xst, /yst, $ title='Raw UFSS_B_XANG in the period', $ ytitle='arcsec', timerange=[ts,te] utplot, time[xs:xe], uf_by[xs:xe], /xst, /yst, $ title='Raw UFSS_B_YANG in the period', $ ytitle='arcsec', timerange=[ts,te] utplot, time[xs:xe], aol[xs:xe], /xst, /yst, $ title='AOL in the period', $ ytitle='deg', timerange=[ts,te] ;utplot, time, xx, timerange=['13-may-2007 16:20','13-may-2007 16:50'] !p.multi=0 endif ; -------------------------------------------- ; Correct (xrt_cen_x, xrt_cen_y) address with orbital variation model ; -------------------------------------------- ; E-W in arcsec dev_ay= xrt_fit_ay *1.031 + uf_ay_xc dev_by= xrt_fit_by *1.031 + uf_by_xc ;N-S in arcsec dev_ax= xrt_fit_ax *1.031 + uf_ax_xc dev_bx= xrt_fit_bx *1.031 + uf_bx_xc ; ; Jitter correction ; off[0,*] - E-W direction ; off[1,*] - N-S direction ; For jitter correction, perform ; data_ca= image_translate(data_x, off) ; /interp) ; off= fltarr(2, num) off[0,*]= (dev_ay*(-1) + dev_by*(-1))/2. off[1,*]= (dev_ax*(-1) + dev_bx*(-1))/2. if keyword_set(a_only) then begin off[0,*]= dev_ay*(-1) off[1,*]= dev_ax*(-1) endif if keyword_set(b_only) then begin off[0,*]= dev_by*(-1) off[1,*]= dev_bx*(-1) endif ; if keyword_set(xrt_pix) then begin off= off/1.031 endif if keyword_set(verbose) then begin ; UTplot the offset data on the screen ts= anytim2ints(buttime, offset=mdpclk_x[0]) te= anytim2ints(buttime, offset=mdpclk_x[n_elements(mdpclk_x)-1]) window, xs=600, ys=600, /free !p.multi=[0, 2, 3] !x.charsize=1.4 !y.charsize=1.4 !p.charsize=1.4 utplot, mdpclk_x, dev_ax*(-1), buttime, /xst, /yst, psym=10, $ title='N-S Offset derived with UFSS-A', $ ytitle='arcsec', timerange=[ts,te] utplot, mdpclk_x, dev_ay*(-1), buttime, /xst, /yst, psym=10, $ title='E-W Offset derived with UFSS-A', $ ytitle='arcsec', timerange=[ts,te] utplot, mdpclk_x, dev_bx*(-1), buttime, /xst, /yst, psym=10, $ title='N-S Offset derived with UFSS-B', $ ytitle='arcsec', timerange=[ts,te] utplot, mdpclk_x, dev_by*(-1), buttime, /xst, /yst, psym=10, $ title='E-W Offset derived with UFSS-B', $ ytitle='arcsec', timerange=[ts,te] utplot, mdpclk_x, off[1,*], buttime, /xst, /yst, psym=10, $ title='N-S Values in returned offset', $ ytitle=' ', timerange=[ts,te] utplot, mdpclk_x, off[0,*], buttime, /xst, /yst, psym=10, $ title='E-W Values in returned offset', $ ytitle=' ', timerange=[ts,te] !p.multi=0 !x.charsize=1. !y.charsize=1. !p.charsize=1. endif return end