; file: extline_wsa.pro = compute Fraunhofer line extinction ; init: Sep 25 2013 Rob Rutten ; last: Nov 5 2015 Rob Rutten Deil ; note: Jun 16 2014 Natasha Shchukina correction to gam6 applied (@check) ; Apr 1 2015 was extline.pro before ;=========================================================================== function extline_wsa,$ temp,eldens,nhtot,vmicro,$ delwav,linewav,$ elemnr,ionstage,$ gfval,excev,gup,$ blow=blow,bup=bup,logstarkc4=logstarkc4,totalext=totalext ;+ ;; ext=extline_wsa($ ;; temp,eldens,nhtot,vmicro,$ ; arrays(h), cgs, total H, km/s ;; delwav,linewav,$ ; delwav array dl AA, airwav AA ;; elemnr,ionstage,$ ; neutral=0 ;; gfval,excev,gup,$ ; excitation energy in eV ;; blow=blow,bup=bup,$ ; Zwaan beta arrays(h) ;; logstarkc4=logstarkc4,$ ; quadratic Stark log(c4) value ;; totalext=totalext ; return total extinction, not monochrom ;; returns float array ext(delwav,h) or ext(delwav) ;- ; common in Wittmann programs common common_atomdata_wsa,at_weight,at_abu,$ at_ionerg1,at_ionerg2,at_ionerg3,at_sym if (n_params() lt 11) then begin print,' ext=extline_wsa($' print,' temp,eldens,nhtot,vmicro,$ ; arrays(h), cgs, total H, km/s' print,' delwav,linewav,$ ; delwav array dl AA, airwav AA' print,' elemnr,ionstage,$ ; neutral=0' print,' gfval,excev,gup,$ ; excev in eV' print,' blow=blow,bup=bup,$ ; Zwaan beta arrays(h)' print,' logstarkc4=logstarkc4,$ ; quadratic Stark log(c4)' print,' totalext=totalext ; return total extinction' print,' returns float array(delwav,h) or (delwav)' return,-1 endif ; keyword defaults if (n_elements(blow) eq 0) then blow=1. if (n_elements(bup) eq 0) then bup=1. if (n_elements(totalext) eq 0) then totalext=0. if (n_elements(logstarkc4) eq 0) then logstarkc4=0 ; physics constants hcgs=6.62607D-27 ; Planck constant (erg s) ccgs=2.99792458D10 ; velocity of light (cm/s) kcgs=1.38062D-16 ; B;; oltzmann constant in erg/deg kev=8.61734D-5 ; Boltzmann constant (eV/deg) hckcgs=hcgs*ccgs/kcgs ; hc/k in cgs matom=1.6605D-24 ; atomic mass constant melectron=9.1094E-28 ; electron mass in gram eelectron=4.803204E-10 ; electron charge in cgs ; get and check array sizes (height or z in 1-D atmosphere) ntemp=n_elements(temp) neldens=n_elements(eldens) nnhtot=n_elements(nhtot) if (neldens ne ntemp) then print,' #### temp, eldens not same dim' if (nnhtot ne ntemp) then print,' #### temp, nhtot not same dim' natmos=max([ntemp,neldens,nnhtot]) ; get _wsa atom parameters atomdata_wsa,elemnr,atomwgt,abund,ion1,ion2,ion3,elemsymbol ionerg=[ion1,ion2] ; compute gas pressure including all elements from common_atomdata_wsa totabu=total(at_abu) pgas=(totabu*nhtot+eldens)*kcgs*temp ; partition function (array) parfun0=partfunc_rr(temp,elemnr,0) parfun1=partfunc_rr(temp,elemnr,1) parfun2=partfunc_rr(temp,elemnr,2) if (ionstage eq 0) then parfun=parfun0 if (ionstage eq 1) then parfun=parfun1 ; Saha (three stages) saha_rr,temp,eldens,ionerg[0],ionerg[1],parfun0,parfun1,parfun2,$ n0_ntot,n1_ntot,n2_ntot if (ionstage eq 0) then ionfrac=n0_ntot if (ionstage eq 1) then ionfrac=n1_ntot ; ========= line broadening helpwav=linewav airtovac,helpwav ;RR it changes! old function was better vacwavaa=helpwav vacwavcm=vacwavaa*1.E-8 ; radiation damping (number) gamrad=6.67E13*(gfval/gup)/((vacwavaa/10.)^2) ; van der Waals = Gray 3rd p245 ff, has N_H factored out already ; correction: RH has Mihalas II table 9-1 p.286 = 8.08 ; while Gray has prefactor 17 in his eq.11.28 on p.245 ; Mihalas value fits RH results (see cdrhidl testgammas.idl) c6=0.3*1E-30*(1/(ionerg[ionstage]-excev-1.2398E4/vacwavaa)^2 $ -1/(ionerg[ionstage]-excev)^2) loggam6=20.+0.4*alog10(c6)+alog10(pgas)-0.7*alog10(temp) ;; gam6=10^loggam6*8.08/17. ; array ;RR Jun 15 2014 Natasha said this is wrong ;; @RR check, also again against RH gam6=10^loggam6 ; quadratic Stark with specified logc4 ;RR @ ?? todo: replace C4 for =0 by Traving estimate (RH) ; or by Bates-Damgaard in 1978ApJ...224.1079F for neutrals, Griem for ions if (logstarkc4 eq 0) then gam4=0 else begin ; Gray III 11.27 p244 loggam4=19.+0.666667*logstarkc4 $ +alog10(eldens*kcgs*temp)-0.833333*alog10(temp) gam4=10^loggam4 endelse ; Doppler width (array) doppwidthcm=(vacwavcm/ccgs)*$ sqrt(2*kcgs*temp/(atomwgt*matom)+(vmicro*1.E5)^2) doppwidthaa=doppwidthcm*1.E8 ; Voigt a (array) voigta=(vacwavcm)^2/(4*!pi*ccgs)*(gamrad+gam6+gam4)/doppwidthcm ; get total extinction extlinetot= $ 1.E8*!pi*eelectron^2/(melectron*ccgs) $ ; units total = cm^-1 AA *vacwavcm^2/ccgs $ ; wavelength units *ionfrac*abund*nhtot*blow $ ; lower level density *(gfval/parfun)*exp(-excev/(kev*temp)) $ ; gf * Boltzmann *(1.-(bup/blow)*exp(-hckcgs/(vacwavcm*temp))) ; stimulated emission ; get monochromatic extinction in loop over wav (to keep atmosphere arrays) if (totalext) then extline=extlinetot else begin sizedelwav=size(delwav) if (sizedelwav[0] eq 0) then nwav=1 else nwav=sizedelwav[1] extline=fltarr(nwav,natmos) for iwav=0,nwav-1 do extline[iwav,*]=extlinetot $ *voigt(voigta,delwav[iwav]/doppwidthaa)/(sqrt(!pi)*doppwidthaa) endelse return,extline end ; ============= main to test per IDLWAVE Hyper-c ======================== ; put plot window besides emacs frame window,xpos=0,ypos=450,xsi=630,title='WSA' ; read FALC falcfile='/home/rutten/rr/web/rjr-edu/exercises/ssb/falc.dat' readcol,falcfile,height,tau5,colmass,temp,vmicro,nhtot,nprot,nel,ptot,$ pgas_ptot,rho,skipline=4 ; set wavelength sample grid ;; delwav=indgen(1000)*0.01-5. ; fine spacing ;; delwav=[-5.,0.,+5.] ; test spacing delwav=0 ; only line center ; get MgIb2 line parameters read_line_params,5173,$ linewav,elemnr,ionstage,excev,$ glow,gup,llow,lup,gfval,logc4 ; compute extinction extmgb2=extline_wsa(temp,nel,nhtot,vmicro,$ delwav,linewav,elemnr,ionstage,$ gfval,excev,gup,logstarkc4=logc4) plot,height,alog10(extmgb2) ;; ; make plot to compare with Fig 13 of 2011A&A...531A..17R ;; ; that is NLTE with increasing overpopulation; compare deep LTE only! ;; plot,height,extmgb2[0,*]/rho*1E-3 ;; ; get NaID1 line parameters ;; read_line_params,5896,$ ;; linewav,elemnr,ionstage,excev,$ ;; glow,gup,llow,lup,gfval,logc4 ;; ; compute NaD1 extinction ;; extnad1=extline_wsa(temp,nel,nhtot,vmicro,$ ;; delwav,linewav,elemnr,ionstage,$ ;; gfval,excev,gup,logstarkc4=logc4) ;; oplot, height,extnad1[0,*]/rho*1E-3 ;; NaI looks similar, MgI same shape but higher here (max 2.8, Jorrit 2.0) ;; differences with Jorrit's plot may be gfval differences ;; Jorrit's plot has initial turnup, I don't, why?? Some damping? end