這份程式的概念很簡單，但也沒想到會寫成這麼長的程式碼(273行)…
學長說還可以再精簡，而且還有一些小問題。
不過我目前不想再改動它了(論文code不能再拖了…)
初始條件1，用免費軟體AstroArt讀出來的星表，把三色影像都讀出。
初始條件2，用眼睛對好星表幾號和幾號的星是同一顆，然後把號碼寫成一個txt檔。(這步驟應該可以自動化，但我目前還沒有概念)
有以上兩個條件，我先把參考星迴歸，得到三幅影像平移旋轉的參數，
然後把三份星表改成同座標系統，位置相同(小於某畫素)就當作同一顆星，
製作新的星表，附上三色星等，繪製星圖、赫羅圖、雙色圖等。
不過我至今仍無法解決把星圖的星星依RGB比例上色的問題… 我以為很簡單，但我不會用IDL寫。
以下就是程式碼了，智慧財產權聲明同本BLOG。
;This is a program try to read an initial data from AstroArt4.0 .sta file.
;And in this program, we can get a Color-Magnitude Diagram and Two-Color Diagram.
;This NGC6802 images are taken by LOT, 2008/05/24, with B=400s, V=140s, R=70s.
;And this program is written by Hsieh Hsiang-Yu, 2011/10/12 update.
pro HRD6802
;==================enviromental setting=====================
close,/all
Q=' '
cir = findgen(17) * (!pi*2/16.)
print,'=====Start!!!!!========'
set_plot, 'ps'
device, filename='D:\NGC6802\6802CMD.ps'
xx=fltarr(1)
yy=fltarr(1)
;===========================reference star reading=============================
openr, 5, 'D:\NGC6802\RefStar.prn'
readf, 5, Q
ref=0
while not (EOF(5)) do begin
readf, 5, format='(x)'
ref=ref+1
endwhile
close,5
print,' total number of reference stars = ',ref
Vnum=fltarr(ref)
Bnum=fltarr(ref)
Rnum=fltarr(ref)
PXV =fltarr(ref)
PXB =fltarr(ref)
PXR =fltarr(ref)
PYV =fltarr(ref)
PYB =fltarr(ref)
PYR =fltarr(ref)
CMV =fltarr(ref)
CMB =fltarr(ref)
CMR =fltarr(ref)
;=========================V img. ref. & catalog writting============================
plot,xx,yy,psym=3,xrange=[0,1340],yrange=[1300,0],/xstyle,/ystyle,/isotropic,title='V band image'
openr, 5, 'D:\NGC6802\RefStar.prn'
readf, 5, Q
readf, 5, format='(3(4x,i4))',NB,NV,NR
rf=0
openr,1, 'D:\NGC6802\6802V.sta'
openw,2, 'D:\NGC6802\6802VM.prn'
v=0
for i=0,1 do readf,1, Q
while not (EOF(1)) do begin
readf,1, format='(f7.2,x,f7.2,x,i8,94x)',xc,yc,aduv
xx(0)=xc
yy(0)=yc
usersym, cos(cir), sin(cir), /fill
oplot,xx,yy,psym=8,symsize=(-2.5*alog10(ADUV)+31)*(-0.063)+1.5 ;fake mag. diameter!!
printf,2,format='(i4,3(x,f7.2))',v+1,xc,yc,-2.5*alog10(ADUV)+31
v=v+1
if v eq NV then begin
Vnum(rf)=NV
PXV(rf)=xc
PYV(rf)=yc
CMV(rf)=-2.5*alog10(ADUV)+31
xx(0)=xc
yy(0)=yc
usersym, cos(cir), sin(cir)
oplot,xx,yy,psym=8,symsize=1,color=0
rf=rf+1
if rf ge ref then goto, jumpV
readf, 5, format='(3(4x,i4))',NB,NV,NR
jumpV:
endif
endwhile
close,/all
print,' total number of stars in V frame = ',v
;==========================B img. ref. & catalog writting=============================
usersym, cos(cir), sin(cir), /fill
plot,xx,yy,psym=3,xrange=[0,1340],yrange=[1300,0],/xstyle,/ystyle,/isotropic,title='B band image'
openr, 5, 'D:\NGC6802\RefStar.prn'
readf, 5, Q
readf, 5, format='(3(4x,i4))',NB,NV,NR
rf=0
openr,1, 'D:\NGC6802\6802B.sta'
openw,2, 'D:\NGC6802\6802BM.prn'
b=0
for i=0,1 do readf,1, Q
while not (EOF(1)) do begin
readf,1, format='(f7.2,x,f7.2,x,i8,94x)',xc,yc,adub
xx(0)=xc
yy(0)=yc
usersym, cos(cir), sin(cir), /fill
oplot,xx,yy,psym=8,symsize=(-2.5*alog10(ADUB)+31)*(-0.063)+1.5,color=0 ;fake mag. diameter!!
printf,2,format='(i4,3(x,f7.2))',b+1,xc,yc,-2.5*alog10(ADUB)+31
b=b+1
if b eq NB then begin
Bnum(rf)=NB
PXB(rf)=xc
PYB(rf)=yc
CMB(rf)=-2.5*alog10(ADUB)+31
xx(0)=xc
yy(0)=yc
usersym, cos(cir), sin(cir)
oplot,xx,yy,psym=8,symsize=1,color=0
rf=rf+1
if rf ge ref then goto, jumpB
readf, 5, format='(3(4x,i4))',NB,NV,NR
jumpB:
endif
endwhile
close,/all
print,' total number of stars in B frame = ',b
;========================R img. ref. & catalog writting==============================
plot,xx,yy,psym=3,xrange=[0,1340],yrange=[1300,0],/xstyle,/ystyle,/isotropic,title='R band image'
openr, 5, 'D:\NGC6802\RefStar.prn'
readf, 5, Q
readf, 5, format='(3(4x,i4))',NB,NV,NR
rf=0
openr,1, 'D:\NGC6802\6802R.sta'
openw,2, 'D:\NGC6802\6802RM.prn'
r=0
for i=0,1 do readf,1, Q
while not (EOF(1)) do begin
readf,1, format='(f7.2,x,f7.2,x,i8,94x)',xc,yc,adur
xx(0)=xc
yy(0)=yc
usersym, cos(cir), sin(cir), /fill
oplot,xx,yy,psym=8,symsize=(-2.5*alog10(ADUR)+31)*(-0.063)+1.5,color=0 ;fake mag. diameter!!
printf,2,format='(i4,3(x,f7.2))',r+1,xc,yc,-2.5*alog10(ADUr)+31
r=r+1
if r eq NR then begin
Rnum(rf)=NR
PXR(rf)=xc
PYR(rf)=yc
CMR(rf)=-2.5*alog10(ADUr)+31
xx(0)=xc
yy(0)=yc
usersym, cos(cir), sin(cir)
oplot,xx,yy,psym=8,symsize=1,color=0
rf=rf+1
if rf ge ref then goto, jumpR
readf, 5, format='(3(4x,i4))',NB,NV,NR
jumpR:
endif
endwhile
close,/all
print,' total number of stars in R frame = ',r
;======================regressing B&R img. to V img.============================
BXY=[transpose(PXB),transpose(PYB)]
BXV=regress(BXY,PXV,SIGMA=BXVS,CONST=BXVC,MEASURE_ERRORS=BXVM)
BYV=regress(BXY,PYV,SIGMA=BYVS,CONST=BYVC,MEASURE_ERRORS=BYVM)
print,' Regress Sigma in BtoV X-axis = ',BXVS
print,' Regress Sigma in BtoV Y-axis = ',BYVS
RXY=[transpose(PXR),transpose(PYR)]
RXV=regress(RXY,PXV,SIGMA=RXVS,CONST=RXVC,MEASURE_ERRORS=RXVM)
RYV=regress(RXY,PYV,SIGMA=RYVS,CONST=RYVC,MEASURE_ERRORS=RYVM)
print,' Regress Sigma in RtoV X-axis = ',RXVS
print,' Regress Sigma in RtoV Y-axis = ',RYVS
openr, 1, 'D:\NGC6802\6802VM.prn'
openr, 2, 'D:\NGC6802\6802BM.prn'
openr, 3, 'D:\NGC6802\6802RM.prn'
openw, 4, 'D:\NGC6802\6802C.prn'
iv=fltarr(v)
xv=fltarr(v)
yv=fltarr(v)
mv=fltarr(v)
for i=0,v-1 do begin
readf, 1, format='(i4,3(x,f7.2))',n,x,y,m
iv(i)=n
xv(i)=x
yv(i)=y
mv(i)=m
endfor
close, 1
ib=fltarr(b)
xb=fltarr(b)
yb=fltarr(b)
mb=fltarr(b)
for i=0,b-1 do begin
readf, 2, format='(i4,3(x,f7.2))',n,x,y,m
ib(i)=n
xb(i)=x
yb(i)=y
mb(i)=m
endfor
close, 2
xbv=xb*BXV[0]+yb*BXV[1]+BXVC
ybv=xb*BYV[0]+yb*BYV[1]+BYVC
ir=fltarr(r)
xr=fltarr(r)
yr=fltarr(r)
mr=fltarr(r)
for i=0,r-1 do begin
readf, 3, format='(i4,3(x,f7.2))',n,x,y,m
ir(i)=n
xr(i)=x
yr(i)=y
mr(i)=m
endfor
close, 3
xrv=xr*RXV[0]+yr*RXV[1]+RXVC
yrv=xr*RYV[0]+yr*RYV[1]+RYVC
;=========================Match stars and make catalog===============================
printf,4, ' ck NB X(V) Y(V) Bmag Vmag Rmag B-V V-R'
plot,xx,yy,psym=3,xrange=[0,1340],yrange=[1300,0],/xstyle,/ystyle,/isotropic,title='Match stars image (size in V mag.)'
;=============match circle setting=================
xcent=650 ;match center in x
ycent=650 ;match center in y
matchD=300 ;match circle diameter
StarD=1. ;Stars pos. less then 1 pixel = same star.
;================end of setting=====================
ck=0
for i=0,b-1 do begin
for j=0,v-1 do begin
if ((xv(j)-xbv(i))^2+(yv(j)-ybv(i))^2)^0.5 lt StarD then begin
for k=0,r-1 do begin
if ((xv(j)-xrv(k))^2+(yv(j)-yrv(k))^2)^0.5 lt StarD then begin
if ((xv(j)-xcent)^2+(yv(j)-ycent)^2)^0.5 lt matchD then begin
printf,4, format='(2(i4,x),7(x,f7.2))',ck,ib(i),xv(j),yv(j),mb(i),mv(j),mr(k),mb(i)-mv(j),mv(j)-mr(k)
xx(0)=xv(j)
yy(0)=yv(j)
usersym, cos(cir), sin(cir),/fill
oplot,xx,yy,psym=8,symsize=(mv(j))*(-0.063)+1.5
ck=ck+1
endif
endif
endfor
endif
endfor
endfor
close, /all
print,' total number of picked stars = ',ck
;=========================Plot the CMD and TCD===============================
Rmag=fltarr(ck)
Vmag=fltarr(ck)
Bmag=fltarr(ck)
openr, 4, 'D:\NGC6802\6802C.prn'
readf, 4, Q
for i=0,ck-1 do begin
readf, 4, format='(2(i4,x),7(x,f7.2))',nni,nnb,nnx,nny,nmb,nmv,nmr,nmbv,nmvr
Rmag(i)=nmr
Vmag(i)=nmv
Bmag(i)=nmb
endfor
usersym, cos(cir), sin(cir), /fill
plot, Bmag-Vmag, Vmag, psym=8, symsize=0.2, xtitle='B-V', ytitle='Vmag', $
title='Color-Magnitude Diagram', xrange=[-0.8,1.2], yrange=[23,16], /xstyle, /ystyle
plot, xx, yy, psym=3, xtitle='V-R', ytitle='B-V',$
title='Two-Color Diagram (size in V mag.)', xrange=[-0.5,1.5], yrange=[1.0,-0.8], /xstyle, /ystyle, /isotropic
for i=0,ck-1 do begin
xx(0)=Vmag(i)-Rmag(i)
yy(0)=Bmag(i)-Vmag(i)
oplot, xx, yy, psym=8, symsize=Vmag(i)*(-0.063)+1.5
endfor
close, /all
device, /close
print,'=====The End!!!!!====='
end
以上，希望日後還會再改善。
沒有留言:
張貼留言