; file: makefigSBJ.pro = main to produce SBJ figure ; init: Feb 7 2012 ; note: = fig-S2level without the 2-level S ; note: parameter-free except xrange, yrange axis tuning etc ; compute J_bar = R_lu/B_lu (holds also for PRD) @constants_si.idl Rlu=*thisatom.transition[itrans].Rij_ptr gl=thisatom.g[ilow] gu=thisatom.g[jupp] Aul=thisatom.transition[itrans].strength Blu=(gu/gl)*Aul*(wavlc*nm_to_m)^3/(2.*hplanck*clight) Jbar=Rlu/Blu ; compute various source functions readj,iwavplot readopacity,iwavplot Sline=eta_as/chi_as Stotal=(eta_as + eta_c + J*scatt) / (chi_c + chi_as) if (PRD) then Jplot=J else Jplot=Jbar ; monochromatic or freq average ; make SBJ figure ; --------------- if (opengv) then spawn,'pkill -s 0 gv' psfilename='fig-SBJ-'+modelname+'-'+filelinename+'.ps' openpsplot,psfilename,thick=1.5 xtitle='height [km]' ytitle='log S, J, B [J m!U-2!N s!U-1!N Hz!U-1!N sr!U-1!N]' bnu=planck(atmos.T,wavlc,/hz) logbnu=alog10(bnu) if (n_elements(maxplotheight) eq 0) then maxplotheight=max(heightkm) if (maxplotheight eq 0) then maxplotheight=max(heightkm) xrange=[min(heightkm),maxplotheight] ;;xrange=[-100,2100] ymin=min(alog10(Sline)) ymax=max(logbnu) yrange=[-8.5,-6.3] yrange=[ymin-0.1*(ymax-ymin),ymax+0.1*(ymax-ymin)] plot,heightkm,logbnu,$ position=[0.2,0.2,0.95,0.95],/normal,$ ; set margins around plot xticklen=0.03,yticklen=0.03/1.6,$ ; same-length ticks xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1, yrange=yrange,ystyle=1 ; ,ytickinterval=1 oplot,heightkm,alog10(Stotal),linestyle=0 ; solid, total source function oplot,heightkm,alog10(Sline),linestyle=1 ; dotted, S^l ;; oplot,heightkm,alog10(Jplot),linestyle=2 ; dashed, J_bar or J ;; if (PRD) then oplot,heightkm,alog10(Jbar),linestyle=3 ; dotdashed, J_bar ; add tau ticks tau=fltarr(nh) tau=gettau(height,chi_as+chi_c) ; getTau in readopacity.pro tabinv,tau,[0.3,1.0,3.0],tau_eff xtau=linear(heightkm,tau_eff) ytau=linear(alog10(Sline),tau_eff) ydash=(yrange[1]-yrange[0])/20. for itau=0,2 do begin xtick=[xtau[itau],xtau[itau]] ytick=[ytau[itau]-0.10,ytau[itau]+0.10] ytick=[ytau[itau]-ydash,ytau[itau]+ydash] oplot,xtick,ytick endfor ; add model name xyouts,0.25,0.85,/norm,modelname ; add line name xyouts,0.25,0.78,/norm,pslinelabel ; add run dir if (addrundir eq 1) then begin cd,c=thedir lastdir=strsplit(thedir,'/',/extract) sizelastdir=size(lastdir) xyouts,0.01,0.01,/norm,lastdir[sizelastdir[1]-1] endif ; close plot closepsplot,psfilename,opengv=opengv end