とりあえず、分光データのReductionメモ。

必要なファイル

手順

簡単に以下の順番で解析(一意ではない)。

  1. overscan 処理
  2. (bias 処理、dark 処理)
  3. flat 補正
  4. 1次元化
  5. 波長較正
  6. 連続光割り
  7. 太陽中心に補正
以下簡単に説明。
  1. overscan 処理
    overscan 補正用のスクリプト:overscan.coverscan.mk

    補正前(左)と後(右)では以下のようになる。
     

  2. (bias 処理、dark 処理)
    bias 処理はoverscan 処理をしていれば必要ないのだが、充分なS/Nがあれば、各ピクセルのbiasということで行うのはあり。
    dark 処理は充分低くなる様に設定されていれば行わなくてもよい。

    各bias、darkは1枚にまとめる。
    IRAFはファイルリストを作って一括処理ができるので、初めにbiasファイルリスト、targetファイルリスト、等を作っておくとよい。

    IRAFで、input とoutput のファイル名を同じにしておくと上書き処理される。また、IRAF はdefault では現存するファイルを上書き処理してくれない。
    ※因みに、IRAFでは拡張子の省略は可能(output ファイルをデフォルトで"fits"としておけばoutput だけ"XXX.imh"になることもない。)。
    尚、個々に処理するときは以下のよう。

  3. flat 補正
    flatはここで行うか"apflatten"コマンドで行うか2通りあるので、自分が良いと思う方で行う。apflattenコマンドを使う場合は、1次元化をしてから行う。

    flat 画像を作る。

    flat処理をする

  4. 1次元化
    >
    以降1次元画像に関するtask は、以下のpackage 下にあるので、そこに移動する。
    ※一般的な分光観測データの場合(noao-onespec package)
    尚、そのpackege から出たいときには、
    とする。
    ※echelle分光観測データの場合(noao-imred-echelle package)
    とする。
    また、それぞれ、動かすパラメータ(epar task で編集できる。interactive=yes とすれば処理をしながらパラメータの変更もできる。)も初めにある程度自分の好きなように変えておくと便利。

    1次元化は、明るい画像(flat 画像がお勧め)を用いて基準となる道を作って、その道をなぞる様に行う。

    まず、flat 画像の一つ(補足1)を用いて、道を作る

    道に沿って各画像を1次元ファイルにする

    readnoise とgain については、OAO/HIDES の場合は次の値を用いる(仕様のページから)。但し値は左右の平均を用いることにした。

    CCDReadnoise[ADU]Gain[e-/ADU]
    1(Red)1.882.618
    2(Green)1.502.985
    0(Blue)1.742.697

    その他のそれぞれのtaskのパラメータは例えば次のように設定。

    epar aptrace
    	input   =                       List of input images to trace
    	(apertur=                     ) Apertures
    	(referen=                     ) List of reference images
    	(interac=                  yes) Run task interactively?
    	(find   =                   no) Find apertures?
    	(recente=                   no) Recenter apertures?
    	(resize =                   no) Resize apertures?
    	(edit   =                  yes) Edit apertures?
    	(trace  =                  yes) Trace apertures?
    	(fittrac=                  yes) Fit the traced points interactively?
    	(line   =                INDEF) Starting dispersion line
    	(nsum   =                   10) Number of dispersion lines to sum
    	(step   =                   10) Tracing step
    	(nlost  =                  100) Number of consecutive times profile is lost befo
    	(functio=             legendre) Trace fitting function
    	(order  =                    4) Trace fitting function order
    	(sample =                    *) Trace sample regions
    	(naverag=                    1) Trace average or median
    	(niterat=                    3) Trace rejection iterations
    	(low_rej=                   3.) Trace lower rejection sigma
    	(high_re=                   3.) Trace upper rejection sigma
    	(grow   =                   0.) Trace rejection growing radius
    	(mode   =                   ql)
    
    epar apsum
    	input   =                       List of input images
    	(output =                     ) List of output spectra
    	(apertur=                     ) Apertures
    	(format =              echelle) Extracted spectra format
    	(referen=                     ) List of aperture reference images
    	(profile=                     ) List of aperture profile images
    	(interac=                  yes) Run task interactively?
    	(find   =                   no) Find apertures?
    	(recente=                   no) Recenter apertures?
    	(resize =                   no) Resize apertures?
    	(edit   =                  yes) Edit apertures?
    	(trace  =                   no) Trace apertures?
    	(fittrac=                   no) Fit the traced points interactively?
    	(extract=                  yes) Extract apertures?
    	(extras =                  yes) Extract sky, sigma, etc.?
    	(review =                  yes) Review extractions?
    	(line   =                INDEF) Dispersion line
    	(nsum   =                   10) Number of dispersion lines to sum or median
    	(backgro=                 none) Background to subtract (none|average|fit)
    	(weights=             variance) Extraction weights (none|variance)
    	(pfit   =                fit1d) Profile fitting type (fit1d|fit2d)
    	(clean  =                  yes) Detect and replace bad pixels?
    	(skybox =                    1) Box car smoothing length for sky
    	(saturat=                50000) Saturation level
    	(readnoi=                  5.0) Read out noise sigma (photons)
    	(gain   =                  4.2) Photon gain (photons/data number)
    	(lsigma =                   3.) Lower rejection threshold
    	(usigma =                   3.) Upper rejection threshold
    	(nsubaps=                    1) Number of subapertures per aperture
    	(mode   =                   ql)
    

  5. 波長較正
    1次元化されたデータは座標軸がpixel になっている。これを波長に変換する。
    まず、comparison 画像で、pixel-波長を対応させ、これをobject 画像に参照させる。

    comarison画像を用いて、波長を割付ける。
    Th-Arランプのアトラスがあるので、それを見ながら、1つのaperture 辺り、4、5本(万遍なく)割付けるのが望ましい。
    じっくりするのは1枚だけでよく、他は、これを基準に割付けることができる。

    (上のコマンドが割付け、下のコマンドが他のcomparison に割付ける為のもの。)

    この2つのパラメータは例えば次のように設定。

    epar ecidentify
    	images  =                       Images containing features to be identified
    	(databas=             database) Database in which to record feature data
    	(coordli=   linelists$thar.dat) User coordinate list
    	(units  =                     ) Coordinate units
    	(match  =                   1.) Coordinate list matching limit in user units
    	(maxfeat=                  100) Maximum number of features for automatic identif
    	(zwidth =                  10.) Zoom graph width in user units
    	(ftype  =             emission) Feature type
    	(fwidth =                   4.) Feature width in pixels
    	(cradius=                   5.) Centering radius in pixels
    	(thresho=                  10.) Feature threshold for centering
    	(minsep =                   2.) Minimum pixel separation
    	(functio=             legendre) Coordinate function
    	(xorder =                    4) Order of coordinate function along dispersion
    	(yorder =                    4) Order of coordinate function across dispersion
    	(niterat=                    3) Rejection iterations
    	(lowreje=                   3.) Lower rejection sigma
    	(highrej=                   3.) Upper rejection sigma
    	(autowri=                   no) Automatically write to database?
    	(graphic=             stdgraph) Graphics output device
    	(cursor =                     ) Graphics cursor input
    	(mode   =                   ql)
    
    epar ecreidentify
    	images  =                       Spectra to be reidentified
    	referenc=                       Reference spectrum
    	(shift  =                   0.) Shift to add to reference features
    	(cradius=                   5.) Centering radius
    	(thresho=                   10) Feature threshold for centering
    	(refit  =                  yes) Refit coordinate function?
    	(databas=             database) Database
    	(logfile=       STDOUT,logfile) List of log files
    	(mode   =                   ql)
    

    pixel-波長の対応付けが終わったら、これをobject画像に参照させる。objectの前後にconparisonをとって、挟み込むのが良い。なければ1毎の画像で"nearest"にするといい。

    上2つのパラメータの設定は以下の通り

    epar refsp
    	input   =                       List of input spectra
    	(referen=                *.imh) List of reference spectra
    	(apertur=                     ) Input aperture selection list
    	(refaps =                     ) Reference aperture selection list
    	(ignorea=                  yes) Ignore input and reference apertures?
    	(select =              average) Selection method for reference spectra ***This need 2 references. If only 1 reference, select = "nearest"
    	(sort   =                     ) Sort key
    	(group  =                     ) Group key
    	(time   =                   no) Is sort key a time?
    	(timewra=                  17.) Time wrap point for time sorting
    	(overrid=                   no) Override previous assignments?
    	(confirm=                  yes) Confirm reference spectrum assignments?
    	(assign =                  yes) Assign the reference spectra to the input spectr
    	(logfile=       STDOUT,logfile) List of logfiles
    	(verbose=                   no) Verbose log output?
    	answer  =                       Accept assignment?
    	(mode   =                   ql)
    
    epar dispco
    	input   =                       List of input spectra
    	output  =                       List of output spectra
    	(lineari=                   no) Linearize (interpolate) spectra?
    	(databas=             database) Dispersion solution database
    	(table  =                     ) Wavelength table for apertures
    	(w1     =                INDEF) Starting wavelength
    	(w2     =                INDEF) Ending wavelength
    	(dw     =                INDEF) Wavelength interval per pixel
    	(nw     =                INDEF) Number of output pixels
    	(log    =                   no) Logarithmic wavelength scale?
    	(flux   =                  yes) Conserve flux?
    	(samedis=                   no) Same dispersion in all apertures?
    	(global =                   no) Apply global defaults?
    	(ignorea=                   no) Ignore apertures?
    	(confirm=                   no) Confirm dispersion coordinates?
    	(listonl=                   no) List the dispersion coordinates only?
    	(verbose=                  yes) Print linear dispersion assignments?
    	(logfile=                     ) Log file
    	(mode   =                   ql)
    

  6. 連続光割り
    連続光を割ります。輝線。吸収線を、間違えてfit しないように注意。また、fit させる関数の次数は3〜5(多くても6とか、7位まで)で調整するのがいいかも。

    パラメータの設定は以下の通り。

    epar cont
    	input   =                       Input images
    	output  =                       Output images
    	(lines  =                    *) Image lines to be fit
    	(bands  =                    1) Image bands to be fit
    	(type   =                ratio) Type of output
    	(replace=                   no) Replace rejected points by fit?
    	(wavesca=                  yes) Scale the X axis with wavelength?
    	(logscal=                   no) Take the log (base 10) of both axes?
    	(overrid=                   no) Override previously fit lines?
    	(listonl=                   no) List fit but don't modify any images?
    	(logfile=              logfile) List of log files
    	(interac=                  yes) Set fitting parameters interactively?
    	(sample =                  *:*) Sample points to use in fit
    	(naverag=                    1) Number of points in sample averaging
    	(functio=             legendre) Fitting function
    	(order  =                    5) Order of fitting function
    	(low_rej=                   3.) Low rejection in sigma of fit
    	(high_re=                   3.) High rejection in sigma of fit
    	(niterat=                   10) Number of rejection iterations
    	(grow   =                   1.) Rejection growing radius in pixels
    	(markrej=                  yes) Mark rejected points?
    	(graphic=             stdgraph) Graphics output device
    	(cursor =                     ) Graphics cursor input
    	ask     =                  yes                                       
    	(mode   =                   ql)
    

  7. 太陽中心に補正
    星と、地球との位置関係が日々変化するので、太陽に座標を統一させる。
    まず、fits header から、日付、UT、天体座標(R.A.-DEC.)情報を以下のフォーマットで読み取る(このファイル名をrvcor.obs とする)。
    20xx xx xx 00:00:00.0 00:00:00.0 +00:00:00
    (Date,UT,R.A.,DEC.の順)
    Fits header から上の情報を読み取るスクリプト:rvcor.c

    このファイルを使って、補正する速度を求める(別の方法)。
    このコマンドはnoao-astutil packeage にあるので、そこに移動すること。

    設定パラメータは以下の通り。observatoryも編集する必要がある。

    epar rvcor
    	(files  =                     ) List of files containing observation data
    	(images =                     ) List of images containing observation data
    	(header =                  yes) Print header?
    	(input  =                   no) Print input data?
    	(imupdat=                   no) Update image header with corrections?
    	(epoch  =                INDEF) Epoch of observation coordinates (years)
    	(observa=       )_.observatory) Observatory
    	(vsun   =                  20.) Solar velocity (km/s)
    	(ra_vsun=                  18.) Right ascension of solar velocity (hours)
    	(dec_vsu=                  30.) Declination of solar velocity (degrees)
    	(epoch_v=                1900.) Epoch of solar coordinates (years)
    	(year   =                     ) Year of observation
    	(month  =                     ) Month of observation (1-12)
    	(day    =                     ) Day of observation
    	(ut     =                     ) UT of observation (hours)
    
    epar obser
    	command =                 list  Command (set|list|images)
    	obsid   =                    ?  Observatory to set, list, or image default
    	images  =                       List of images
    	(verbose=                   no) Verbose output?
    	(observa=               obspar) Observatory identification
    	(name   =              Okayama) Observatory name
    	(longitu=      226.40361111111) Observatory longitude (degrees)
    	(latitud=      34.573888888889) Observatory latitude (degrees)
    	(altitud=                 372.) Observatory altitude (meters)
    	(timezon=                  -9.) Observatory time zone
    	(mode   =                   ql)
    
    	override=                       Observatory identification
    
    ※ぐんま天文台の場合、4ヵ所変更:
    	(name   =                Gunma) Observatory name
    	(longitu=     -138.97277777778) Observatory longitude (degrees)
    	(latitud=      36.596388888889) Observatory latitude (degrees)
    	(altitud=                 885.) Observatory altitude (meters)
    
    次のような出力がされる。
    # RVCORRECT: Observatory parameters for Okayama
    #       latitude = 34.573888888889
    #       longitude = 226.40361111111
    #       altitude = 372.
    ##   HJD          VOBS   VHELIO     VLSR   VDIURNAL   VLUNAR  VANNUAL   VSOLAR
    2454844.33368     0.00    27.87    24.27     -0.034   -0.009   27.917   -3.606
    2454844.34702     0.00    25.84    37.82      0.251   -0.012   25.605   11.976
    2454844.21767     0.00     2.60   -15.36     -0.178    0.002    2.780  -17.968
    
    必要なのはVHELIO。LSR(局所静止系)に対する太陽の運動に関してのパラメータ(vsun,ra_vsun.dec_vsun,epoch_v)には問題があるがVHELIOの値は変わらないようだ(補足3)。

    VHELIOを用いて、太陽中心に補正する。但し、正負を逆にして値を入れること。

    設定パラメータは以下の通り。
    epar dopcor
    	input   =                       List of input spectra
    	output  =                       List of output spectra
    	redshift=                       Redshift or velocity (Km/s)
    	(isveloc=                  yes) Is the redshift parameter a velocity?
    	(add    =                   no) Add to previous dispersion correction?
    	(dispers=                  yes) Apply dispersion correction?
    	(flux   =                   no) Apply flux correction?
    	(factor =                   3.) Flux correction factor (power of 1+z)
    	(apertur=                     ) List of apertures to correct
    	(verbose=                   no) Print corrections performed?
    	(mode   =                   ql)
    

補足的なこと

aptrace の時に使う画像

flatと星の光の入り方は違う(右図参照)。
ので、traceするときにflat画像で真中を選ぶのは難しい。そこで、なるべくmaskを小さくして(光が入る幅を小さくして)flat 画像を1枚とってtraceする。
とはいえ、露出によってスリットのその部分に天体が来るかはまちまちなので、apsum の時に中心を決め直す必要はありますが…。

rvcor に読み込ませるファイルの作り方2

IRAFのhselect コマンドを用いればrvcor に読ませるべきファイルを作ってくれる。

hselect @INfile.lst obs-DATE,UT,RA,DEC yes

太陽向点の値について

太陽向点は色々不定性があるらしいが、Gupta 2005 "Observers Handbook" によれば、これらはJ2000で、
だそうだ。太陽中心補正をする場合にはVHELIOしか使わないので古い分点のものを使っても影響はないが、他の値は変化するので注意が必要。