[IDL] 黑體輻射曲線程式

又是IDL文… 這次只是個基本的作業，比較特別的就是希臘字母在IDL的寫法。
基本架構在一年半以前就做好了，這次只是refine罷了。



pro PlanckFunction

; This program is trying to draw the Planck Function
; Made by Hsiang yu,Hsieh, 2011/11/10.

h=6.626068765d-27      ;erg*s
c=2.99792458d10        ;cm/s
k=1.380650324d-16      ;erg/K

TVLCT,[0,255,0,0,255],[0,0,255,0,255],[0,0,0,255,0] ;0=black,1=red,2=green,3=blue,4=yellow
set_plot,'ps'
device, filename="D:\PlanckFunction\PlanckFunction.ps",/color

plot,[0],[0],xrange=[0d0,2d-4], yrange=[0d0,1.5d15], xtitle='WaveLength (cm)', ytitle='Energy (erg/cm^2/sec)', title='Planck Function',psym=3

for i=0,2 do begin
  T=6d3+i*1d3 ; Temperature, K
  L_step=1d6
  L=findgen(L_step)/(L_step)
  B=(2d0*h*c^(2d0)/(L^(5d0)))/(exp((h*c)/(L*k*T))-1d0)  ;unit in erg/cm^2/sec
oplot,L,B, color=i+1
endfor

iv=10000  ;K max
WIENL=findgen(iv) ;K
BB=0.2897768551  ;cm*K
LM=dblarr(iv)
BLM=dblarr(iv)
  for i=0,iv-1 do begin
    LM(i)=BB/WIENL(i)
    BLM(i)=(2d0*h*c^(2d0)/(LM(i)^(5d0)))/(exp((h*c)/(LM(i)*k*i))-1d0)
  if i eq 6000 then xyouts,[9d-5,11d-5],[4d14,4d14],['6000K','!4k!X!IMax!X!N='],color=1
    if i eq 6000 then xyouts,12d-5,4d14,LM(i),color=1
  if i eq 7000 then xyouts,[9d-5,11d-5],[8d14,8d14],['7000K','!4k!X!IMax!X!N='],color=2
    if i eq 7000 then xyouts,12d-5,8d14,LM(i),color=2
  if i eq 8000 then xyouts,[9d-5,11d-5],[12d14,12d14],['8000K','!4k!X!IMax!X!N='],color=3
    if i eq 8000 then xyouts,12d-5,12d14,LM(i),color=3
  endfor
 
oplot,LM,BLM,color=0

device,/close
close,/all
end