Otsuji's homepage

Beckers' Cloud model

知ってて得するIDLプログラム

思いついた/気づいた順番に書いていきます。
上の方が新しい。

epsでいろんなフォント

ps環境内でのフォント変更にいつも悩む人へ
基本的にps環境内では!p.font=0に設定。
IDL> set_plot,'ps'
IDL> !p.font=0
どんなフォントが使えるのかをリストアップ
IDL> ps_show_fonts
例えば、Helveticaフォントの斜体を使いたければ、
IDL> device,/HELVETICA,/OBLIQU,font_index=20
これで、'!20'と'!X'の間にはさまれた文字列は指定されたフォントで出力される。
IDL> plot,indgen(10),title='test !20test!X test'
IDL> device,/close
ODL> set_plot,'x'
IDL> $display idl.ps

一度に複数の変数を同じ値に初期化する方法

IDL> a=(b=(c=indgen(10)))
IDL> help,a,b,c
A               INT       = Array[10]
B               INT       = Array[10]
C               INT       = Array[10]

これは何気にすごく便利。例えば、ある配列をその最小値と最大値で0～1に規格化したい場合、

IDL> array=findgen(10)+5    ;最小値5、最大値14の配列を用意
IDL> norm=(array-(minv=min(array,max=maxv)))/(maxv-minv)

のたった一行で規格化できる(しかもSSW非依存)。

structureの要素の指定方法

IDL> print,structure.tagname
または、目的の要素がstructure内の何番目か分かっているなら、
IDL> print,structure.(番号)
で指定可能。

今更ながらarrayの:の用法

IDL> x=indgen(10)
x[2:8]はxの二番目から八番目をひとつずつ。(2,3,4,5,6,7,8)
x[2:8:1]は上と同じ。(2,3,4,5,6,7,8)
x[2:8:2]はxの二番目から八番目をひとつ飛ばしに。(2,4,6,8)
x[8:2:-1]はxの八番目から二番目をひとつずつ。(8,7,6,5,4,3,2)

reverseとか使わずに配列をひっくりかえせるのは便利かも。

for,while,repeatなどloopに関して知っておくと便利なこと

if <条件文> then continue ; それ以下のloop行を飛ばして次のloopへ。
if <条件文> then break ; それ以下のloop行を飛ばしてloopの外へ。

gotoとかでやるよりスマート。

早い配列比較
AとBの配列があって、完全に一致するか否かを調べたい場合、
```
IDL> print,min(A eq B)
```
だと、最後までA eq Bを判定してからminを取るので、配列が大きい場合には顕著に遅くなる。そこで、
```
IDL> print,array_equal(A,B)
```
だと、配列の最初から比較していって不一致な要素を検出した段階で0 (false)を返すのでかなり早い。
文字列配列をrebinしたい場合
```
IDL> str1=['abc','def']
IDL> tmp=byte(str1)
IDL> str2=string(rebin(tmp,3,2,2))
```
と、一旦バイトに変換してからrebinして、stringで文字に戻してやるといい。
ps環境でplot_mapするといちいち
```
% LOADCT: Loading table B-W LINEAR
```
とか表示されて鬱陶しい時は、plot_mapのキーワード指定で
```
IDL> plot_map,map,mask_value=255
```
としてやるとおとなしくなる。
ps環境ではmargin部分をデフォルトでは白く塗るが、この時呼び出されるplot_map_whiteの中でloadctしているのが原因。mask_valueで値を与えてやればloadctしなくなる。
MapのB angleとか
index2mapで出来たmapのB0やRSUNの値は、もとのindex中に以下のどれかのキーワードが入ってないとB0=0、RSUN=960で固定される。
```
'spacec'
'obsrvtry'
'observa'
'telesc'
'instrume'
'detect'
'wavelnth'
'wave_len'
```
各キーワードは前方一致でok(TELESCOPとかでも認識される)。
詳細はget_map_angles.pro及びget_fits_par.proを参照。
EPSで記号とか
illustratorとかでちょこちょこ記号を入力するのが面倒！
でもidlのps環境でどうやったらいいか…。
そんな時は
```
showfont,'symbol','symbol',/TT
```
とすると記号一覧が現れる。
欲しい記号のoctal値をチェック(一番左の列と一番上の行の表示が24xと07なら24*10+7=247)。
PS環境内で
```
str='!M'+string("247b)
xyouts,0.5,0.5,str,/norm
```
とすればOK。
Soho mappingの際の注意点
SohoのMDIなどのデータをindex2mapでmapにする際、元のデータが180度回転している場合はhave_soho_roll.proというプログラムが判断して勝手に回転補正する。
すなわち、
```
mreadfits,mdifits,index,data
index2map,index,data,map
```
とした際にdataとmap.dataが一致しない場合があるので注意すること。
index中のCRVAL系とXCEN系
CRVAL系：CRVAL/SC_ATTX/SC_ATTY/->CRPIXでの位置
XCEN系：XCEN->画像中心での位置
CRPIXが画像中心であれば上二つは一致
Array matching (配列の要素の比較)
要素の個数が異なる2つの配列AとBの内、共通している要素のindexを全て求める。
```
match2,A,B,Asub,Bsub
```
とすれば、
```
A(Asub)=B(Bsub)
```
となる。
ただし、負の数は無視される。
階乗
```
a=factorial(10)
```
で10!。
plotなどでのオングストローム記号(Å)の出力方法
以下の4つの方法がある。
```
angst=string("305B)
angst=string(197B)
angst='!6!sA!r!u!9 %!6!n'
angst='!Z(00c5)!3'
```
自分としては2番めがお薦めかな。

分岐

if flg eq 1 then a=4 else a=6

は実は以下の記述でもできる。

a=flg eq 1 ? 4 : 6

utplotでバグったときの対処法
光度曲線などを時間プロットできるutplotはとても便利だが、 !xなどの環境変数をいじるため、下手にプログラム内でエラーで止まると、プロットの軸がおかしくなったままになることが良くある。
これを防止するためには、あらかじめ.idl_startup内に、
```
plot,indgen(10) ;環境変数の初期化
store_plotvar
```
の一行を加えておくと良い。
もしutplotでバグったときは、
```
restore_plotvar
```
で復活できる。
plotのx軸の数字を消す
```
plot,indgen(10),xtickformat='(a1)'
```
FITS index中のCROTA系とP_ANGLE系の関係
CROT/CROTA/CROTA1/CROTA2の値はindex2mapでMAPのROLL_ANGLEの値にそのまま渡される。
SC_ROLL/P_ANGLE/SOLAR_P0/ANGLEの値はMAP変換の際に符号が逆転するので注意
太陽のp0(=P_ANGLE=-CROTA)は、地上から見て太陽の最上端から、太陽の自転軸(北の方向)が反時計回り(北→東が正)に何度傾いているかを示す値。
つまり、MAPのROLL_ANGLEはデータのy軸から太陽の自転軸が何度傾いているかを"時計回り"に測った値。
なので、FITS indexを作成する際には、
```
time='2009-02-26 12:41:45'
pa=get_rb0p(time,/pa,/deg)
index.CROTA1=-pa
index.CROTA2=-pa
index.P_ANGLE=pa
```
としておくと、どのタグを使ってもmap変換の際にROLL_ANGLE=-paとなり、
```
rmap=rot_map(map,-map.ROLL_ANGLE)
```
としたときに回転角度がpa("時計回り"が正)となるため、正しく補正される。
ちなみに、太陽の回転関連のタグが複数ある場合、CROT->CROTA->CROTA1->CROTA2->SC_ROLL->P_ANGLE->SOLAR_P0->ANGLEの順番で値が採用される(詳しくはget_fits_roll.pro参照)。
配列の差分
```
def=deriv(a)
```
":"は使わない方が得？
既存の配列の特定の範囲に別の小配列を入れたい時に、
```
a=fltarr(200,200)
b=dist(50)
a(100:149,100:149)=b
tvscl,a
```
とかしなくても、
```
a=fltarr(200,200)
b=dist(50)
a(100,100)=b
tvscl,a
```
とした方が早いし簡単。

readせずにファイルのサイズなどを知る

Result = FSTAT(Unit)
Result = QUERY_TIFF ( Filename [, Info] [, IMAGE_INDEX=index] )
spawn,'identify '+file,tmp & ss=strsplit((strsplit(tmp,/ex))(2),'x',/ex)

配列の結合
```
result=append_arr(input,append)
```
文字列の置き換え
```
new=str_replace(strings,'old','new')
```
"not"
```
~
```
ファイルの終わり
```
EOF(unit)
```
直線とシンボルでplot plot,indgen(10),psym=-2
単位行列 identity(m)でmxmの単位行列
対角行列 diag_matrix([a11,a22,...,ann])

mxnのarrayの各点に原点からの距離を入れたい場合

x=indgen(m,n) mod m
y=indgen(m,n)  /  m
r=sqrt(x^2.+y^2.)

でも、普通に

x=rebin(indgen(m),m,n)
y=rebin(indgen(1,n),m,n)
r=sqrt(x^2.+y^2.)

とした方がいいかも。
さらにSSWなら、

dist_circle,r,[m,n],0,0

で一発。

b(i)=∑_k=0ⁱa(k) という配列bを作りたい場合

n=n_elements(a)
c=(reform([indgen(n*(2*n-1)) mod (2*n)]-n lt 0,2*n-1,n))(0:n-1,*)
b=c#a

最近知ったが、

b=total(a,/cumlative)

で出来た。。。

roll_xy,xin,yin,angle,xout,yout
(xin,yin)をangle度回転させて(xout,yout)に。
u=limb_dark(lambda)
lambda[nm]の波長に対する周辺減光定数wを返す。
因みに周辺減光の式は、Ioをdisk centerのIntensity、qを知りたい領域のRadiusと置くと、
I(q)=Io(1-u(1-cos(q)))
plot,y,xtickformat='(a1)'
x軸の数値を表示しない。
BLAS_AXPY
2D配列や、3D配列の全体や一部に1Dや2Dの配列を足したいときに、いちいちfor文で回さなくてよいプログラム。

配列を用意。

Y=intarr(5,5)
print,Y
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0

0行目だけにX=[0,2,4,6,8]を加えたい場合、

Y=intarr(5,5)
X=[0,2,4,6,8]
blas_axpy,Y,1,X,1,[0,0]
print,Y
       0       2       4       6       8
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0

0行目だけにXを3倍してから加えたい場合、

Y=intarr(5,5)
blas_axpy,Y,3,X,1,[0,0]
print,Y
       0       6      12      18      24
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0

1行目だけにXを加えたい場合、

Y=intarr(5,5)
blas_axpy,Y,1,X,1,[0,1],2,0
print,Y
       0       0       0       0       0
       0       2       4       6       8
       0       0       0       0       0
       0       0       0       0       0
       0       0       0       0       0

1行目と3行目だけにXを加えたい場合、

Y=intarr(5,5)
blas_axpy,Y,1,X,1,[0,1],2,[0,2]
;blas_axpy,Y,1,X,1,[0,0],2,[1,3]だとうまくいかない。
print,Y
       0       0       0       0       0
       0       2       4       6       8
       0       0       0       0       0
       0       2       4       6       8
       0       0       0       0       0

各行にXを加えたい場合、

Y=intarr(5,5)
blas_axpy,Y,1,X,1,[0,0],2,[0,1,2,3,4]
print,Y
       0       2       4       6       8
       0       2       4       6       8
       0       2       4       6       8
       0       2       4       6       8
       0       2       4       6       8

0列目にXを加えたい場合、

Y=intarr(5,5)
blas_axpy,Y,1,X,2,[0,0]
print,Y
       0       0       0       0       0
       2       0       0       0       0
       4       0       0       0       0
       6       0       0       0       0
       8       0       0       0       0

3D配列への使いかたはidl help参照。

REPLICATE_INPLACE
既存の配列を変更する。速い、らしい。
centroid
そのまんま。質量中心を返す。
簡単！色付き画像ファイル作成方法
色んなカラーテーブルを使って絵を描く。
write_image,'test.gif','gif',tvrd(/true)
これでok。ただしzバッファだと無理。

zバッファでも使える色々な色を使った画像ファイル作成方法

set_plot,'z'
device,set_pixel_depth=24,decomposed=0

;(描画開始)

wdef,0,400,400
loadct,2
tvscl,dist(200),  0,  0
loadct,3
tvscl,dist(200),200,  0
loadct,4
tvscl,dist(200),  0,200
loadct,5
tvscl,dist(200),200,200

;(描画終了)

write_bmp,'test.bmp',reverse(tvrd(true=1),1)
write_jpeg,'test.jpg',tvrd(true=1),true=1,quality=100

set_plot,'x'

test.bmp

なんと鮮やか!!

いろんな色でplot

device,decomposed=1
wdef,0,200,100
plots,[0,400],[ 0,  0],/dev,color='0000ff'xl
plots,[0,400],[10, 10],/dev,color='00ccff'xl
plots,[0,400],[20, 20],/dev,color='ffcc22'xl
plots,[0,400],[30, 30],/dev,color='3355aa'xl
plots,[0,400],[40, 40],/dev,color='012345'xl
plots,[0,400],[50, 50],/dev,color='6789ab'xl
plots,[0,400],[60, 60],/dev,color='cdefff'xl
plots,[0,400],[70, 70],/dev,color='123f4b'xl
plots,[0,400],[80, 80],/dev,color='888888'xl
plots,[0,400],[90, 90],/dev,color='ffffff'xl

write_bmp,'lines.bmp',reverse(tvrd(true=1),1)

lines.bmp

MDI index keyword 解説ページ

HOME