; file: makefigsepstaulambda.pro = main ; init: Jan 24 2016 Rob Rutten train to Bern ; last: Jan 26 2016 Rob Rutten Bern ; note: in IDL session ; @doinitrh ; @doinitline after selecting Lyalpha ; .r makefigS2leveltemp.pro ; compute eps(h) in two-level approx (OK where eta S^D negligible) Aul=thisatom.transition[itrans].strength Cul=(*thisatom.Cij_ptr)[*,jupp,ilow] @constants_si.idl epsacc=(1.-exp(-hplanck*clight/(wavlc*nm_to_m*kboltzmann*atmos.t)))*Cul/Aul eps2lev=epsacc/(1.+epsacc) ; compute photon mean free path in cm tau=fltarr(nh) tau=gettau(height,chi_as+chi_c) ; getTau in readopacity.pro path=1./(-deriv(heightkm,tau)*1.E-5) ; compute Lamba = 1/eps for Gaussian profile in cm Lambda=(1./eps2lev)/(-deriv(heightkm,tau)) yfigs=[[alog10(eps2lev)],[alog10(path)],[alog10(Lambda)]] ynames=['epsilon','path','Lambda'] ytitles=['log (epsilon)','log (mean free path) [cm]','log (Lambda) [km]'] ; loop over figures for ifig=0,2 do begin y=yfigs[*,ifig] yname=ynames[ifig] ytitle=ytitles[ifig] ; make figure ; ---------------------------- psfilename='fig-'+yname+'-'+modelname+'-'+filelinename+'.ps' openpsplot,psfilename,thick=2 xtitle='height [km]' if (n_elements(maxplotheight) eq 0) then maxplotheight=max(heightkm) if (maxplotheight eq 0) then maxplotheight=max(heightkm) xrange=[min(heightkm),maxplotheight] ymin=min(y) ymax=max(y) yrange=[ymin-0.1*(ymax-ymin),ymax+0.1*(ymax-ymin)] plot,heightkm,y,$ position=[0.2,0.2,0.95,0.95],/normal,$ ; ample margins xticklen=0.03,yticklen=0.03/plotaspect,$ ; same-length ticks xtitle=xtitle,ytitle=ytitle,$ xrange=xrange,xstyle=1, yrange=yrange,ystyle=1,$ thick=2 ; 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 endfor ; end loop over figures end