; file: makefigbetaslinefrac.pro = main; adds curves population fraction ; init: Apr 10 2015 Rob Rutten Deil from makefigbetasline,.pro ; last: Jul 15 2022 Rob Rutten Deil ; note: first runrh on the desired XXX.atmos and XXX.atom RH setup ; adapt RH run dir doinitrh.pro ; then in new IDL session first @doinitrh, @doinitline or progran ; parameter-free except xrange, yrange axis tuning etc ; get populations, compute departure coefficients nlow=(*thisatom.n_ptr)[*,ilow] nstarlow=(*thisatom.nstar_ptr)[*,ilow] nupp=(*thisatom.n_ptr)[*,jupp] nstarupp=(*thisatom.nstar_ptr)[*,jupp] blow=nlow/nstarlow bupp=nupp/nstarupp ; make betas figure ; ----------------- if (opengv) then spawn,'pkill -s 0 gv' psfilename='fig-betas-'+modelname+'-'+filelinename+'.ps' openpsplot,psfilename,thick=2,xsize=8.,ysize=8./plotaspect xtitle='height [km]' ytitle='log b' logbl=alog10(blow) logbu=alog10(bupp) xrange=heightrange yrange=betarange if (yrange[1]-yrange[0] lt 5) then ytickinterval=1 $ else if (yrange[1]-yrange[0] lt 12) then ytickinterval=2 $ else ytickinterval=5 plot,heightkm,logbl,$ position=[0.2,0.2,0.85,0.95],/normal,$ ; set margins around plot xticklen=0.03,yticklen=0.03/plotaspect,$ ; same-length ticks xtitle=xtitle,ytitle=ytitle,ytickinterval=ytickinterval,$ xrange=xrange,xstyle=1,yrange=yrange,ystyle=9,yminor=2 oplot,heightkm,logbu ; also solid ; add level names upper,lower labelheight=heightrange[1]-300. ;; labelheight=1500. ; Jul 15 2022 out for smaller height FeI6302-FeI6173 hindex=where(min(abs(heightkm-labelheight)) eq abs(heightkm-labelheight)) if (logbl[hindex] gt yrange[0]+0.2 and logbl[hindex] lt yrange[1]-0.5) then $ xyouts,heightkm[hindex],logbl[hindex]+0.1,'lower' if (logbu[hindex] gt yrange[0]+0.2 and logbu[hindex] lt yrange[1]-0.5) then $ xyouts,heightkm[hindex],logbu[hindex]-0.2,'upper' ; add line tau ticks on b_low curve tau=fltarr(nh) tau=gettau(height,chi_as+chi_c) ; getTau in readopacity.pro ;; tabinv,tau,[0.3,1.0,3.0],tau_eff tabinv,tau,[1E10,1E8,1E6,1E4,1E3,1E2,1E1,1,1E-1,1E-2,1E-3,1E-4],tau_eff tauname=['10','8','6','4','3','2','1','0','-1','-2','-3','-4'] xtau=linear(heightkm,tau_eff) ytau=linear(logbl,tau_eff) ydash=(yrange[1]-yrange[0])/20. sizetau=size(tau_eff) for itau=0,sizetau[1]-1 do begin xtick=[xtau[itau],xtau[itau]] ytick=[ytau[itau]-ydash,ytau[itau]+ydash] if (xtick[0] gt xrange[0]+50 and xtick[0] lt xrange[1]-50 and $ ytick[0] gt yrange[0]+0.01*(yrange[1]-yrange[0]) and $ ytick[1] lt yrange[1]-0.1*(yrange[1]-yrange[0])) then begin oplot,xtick,ytick xyouts,xtau[itau],ytau[itau]+ydash*1.5,tauname[itau],$ alignment=0.5,charsize=0.8 endif endfor ; Apr 10 2015 add dotted curve population fraction ; ------------------------------------------------ ; sum all pops to find total population (= density/abundance or such) ntotlevels=thisatom.nlevel sumpop=(*thisatom.n_ptr)[*,0] for ilevel=1,ntotlevels-1 do sumpop=sumpop+(*thisatom.n_ptr)[*,ilevel] ; get fraction fracpop=nlow/sumpop ; find index for h=0 indh0=where(abs(heightkm) eq min(abs(heightkm))) normzero=fracpop[indh0[0]] ; bloody IDL array > scalar norm=normzero[0] offset=1.5 normfrac=alog10(fracpop/norm)+offset oplot,heightkm,normfrac,linestyle=2 ; add LTE curve fracpopstar=nstarlow/sumpop normfracstar=alog10(fracpopstar/norm)+offset oplot,heightkm,normfracstar,linestyle=1 ; add label? ;; xyouts,heightkm[hindex],normfracstar[hindex]+0.01,'LTE' ; Apr 29 2015 add fraction axis on right (Jorrit suggestion) fracrange=[(alog10(norm)-offset)+yrange[0],(alog10(norm)-offset)+yrange[1]] axis,yaxis=1,yrange=fracrange,ystyle=1,ytickinterval=ytickinterval,yminor=2,$ ytitle=textoidl("log (fraction)") ; add model name xyouts,0.23,0.23,/norm,modelname ; add line name xyouts,0.62,0.23,/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