; file: extcont_vitense_figs.pro = reproduce Erika Vitense "confusograms" ; init: Jul 9 2010 ; last: May 5 2015 Rob Rutten Deil ; note: redoes the selected Vitense figures in RTSA, taken from Novotny's ; book 1973itsa.book.....N which are nicer than the originals by ; Erika Vitense (later Boehm-Vitense) in 1951ZA.....28...81V. ; She had an element mix which differs much from current solar ; abundances, given in Table 3.7 on p134 of Novotny. Much smaller X. ; She included helium and metal edges, I don't. ;; pro,noppes ;RR insert to make IDLWAVE indentation work ; common with atomic data in the Wittmann-Sanchez-Asensio programs common common_atomdata_wsa,at_weight,at_abu,$ at_ionerg1,at_ionerg2,at_ionerg3,at_sym atomdata_wsa,1 ; initialize common with these numbers ; set plot to screen (for testing) or to ps output psplot=1 ; if X then open plot window besides my emacs frame if (psplot eq 0) then begin set_plot,'x' window,xpos=0,ypos=380,xsi=630 endif ; read physics constants (cgs) @constants_cgs.idl ; change atmosphere composition to Vitense values = Novotny Table 3.7 vith=25100. at_abu[2-1]=4570./vith ; more He than modern value at_abu[6-1]=6.31/vith at_abu[7-1]=13.5/vith at_abu[8-1]=24.5/vith at_abu[10-1]=28.8/vith at_abu[11-1]=0.0575/vith at_abu[12-1]=1.55/vith at_abu[13-1]=0.0955/vith at_abu[14-1]=1.41/vith at_abu[16-1]=0.371/vith at_abu[19-1]=0.00794/vith at_abu[20-1]=0.0759/vith at_abu[26-1]=2.34/vith ; Vitense plot details vitensetheta=[0.9,0.05,0.1,0.5,1.3] vitenselogpe=[[-1,0,1,2,3],$ ; each 5 curves at delta Pe = 1 dex [2,3,4,5,6],$ ; these differ from Vitense's [1,2,3,4,5],$ ; a wider range in logpe pe [0,1,2,3,4],$ ; [-1,-1,0,0.5,1]] ; the selection in her graph [-2,-1,0,1,2]] ; my test, see long comment below xranges=[[2.7,5.1],[1.7,4.1],[1.7,4.1],[2.2,4.6],[2.7,5.1]] plotnames=['vitensefig1.ps','vitensefig2.ps','vitensefig3.ps',$ 'vitensefig4.ps','vitensefig5.ps'] ; wavelength grid nwav=400 logwav=indgen(nwav)*0.01+1.5 ; equidistant in log wav=10.^logwav ; Angstrom ; start big loop over the different graphs nplots=size(vitensetheta,/n_elements) for iplot=0,nplots-1 do begin ;; for iplot=4,4 do begin ; for testing the erroneous curve equality theta=vitensetheta[iplot] temp=5040./theta logpeset=vitenselogpe[*,iplot] ; array ncurves=size(logpeset,/n_elements) for icurve=0,ncurves-1 do begin ;; for icurve=0,0 do begin ; for testing ; compute partial hydrogen pressures logpe=float(logpeset[icurve]) pe=10^logpe gaspress_wsa,pe,temp,phneut,phminus,phplus,ph2,ph2plus,pgas,rho eldens=pe/(kcgs*temp) nhneut=phneut/(kcgs*temp) ; compute hydrogen extinctions exthneu=exthneubf_gray(wav,temp)+exthneuff_gray(wav,temp) exthmin=exthminbf_gray(wav,temp,eldens)+exthminff_gray(wav,temp,eldens) extthomson=0.664E-24*eldens/nhneut ; same normalization per neutral H extthomson=wav*0.+extthomson ; convert into array extrayleigh=1.13/2.*0.664E-24*(1216./float(wav))^4 ; only for HI exthtot=exthneu+exthmin+extthomson+extrayleigh ; add some bound-free edges using VALII and VALIII numbers ; Si I edges following VALIII p7ff ; kap=wav ; @@@@@@@ to do..... ; Vitense's extinction coefficients are not per neutral H atom but per gram norm=(nhneut/rho)/pe ; Vitense normalization per gram and per unit Pe ;RR Oct 6 2013 With this normalization I did not get the Hminus curve ; identity in vitensefig5 which she has (Figure RTSA p195 = Figure ; p149 Novotny). I expected such equality since Hminus scales with Pe ; in Saha (eg RTSA (8.2), Gray (8.12)), so where H-minus dominates the ; total extinction should have this scaling. However, at larger P_e ; an increasing fraction of hydrogen sits in H2 molecules. So ; exthminbf/pe should indeed get increasingly reduced, certainly for ; very large Pe, as I get with the above normalization (extending the ; Pe range for more extreme demonstration). Novotny writes (p133) ; that Vitense took H2 into account, but I can reproduce her graph ; when using nhneuterika below which assumes that she neglected the ; hydrogen sitting in H2 in her H-minus evaluation. I then reread ; Vitense's own 1951ZA.....28...81V again and get the impression that ; she indeed did not incude H2 formation (because it would give ; negligible continuous opacity, as she wrote further down). Added ; this conclusion to the caption in RTSA and added a question with ; answer to the RTSA problems. ; let's inspect some numbers ;; print,alog10([pe,phneut,phminus/phneut,phplus/phneut,ph2/phneut,ph2plus/phneut,exthmin[300]/pe,pgas/pe]),format='(20F10.2)' ; H2 formation gets important ; do it Vitense's way (I think) nhneuterika=(phneut+2*ph2)/(kcgs*temp) ;; norm=(nhneuterika/rho)/pe ; insert to get actual Vitense normalization ; make the figure if (icurve eq 0) then begin axrat=0.8 maxy=min(exthtot)*norm*1E3 miny=maxy*1.E-7 psfilename=plotnames[iplot] if (psplot) then openpsplot,psfilename,$ thick=2,fontsize=9,xsize=8.8,ysize=8.8/axrat plot,alog10(wav),exthtot*norm,/ylog,$ xrange=xranges[*,iplot],xstyle=1,$ yrange=[miny,maxy],ystyle=1,$ xtitle=textoidl("log \lambda [")+cgsymbol("angstrom")+']',$ ytitle=textoidl("\alpha_\lambda / P_e ") +$ textoidl(" [cm^2 gram^{-1} (dyne/cm^2)^{-1}]") oplot,alog10(wav),exthneu*norm,linestyle=2 ; without shift as Vitense oplot,alog10(wav),exthmin*norm,linestyle=1 oplot,alog10(wav),extthomson*norm,linestyle=2 oplot,alog10(wav),extrayleigh*norm,linestyle=2 ; add plot labels xyouts,0.75,0.9,/norm,textoidl("\theta = ") $ +ntostr(theta,format='(f5.2)') xyouts,0.75,0.86,/norm,'T = '+ntostr(temp,format='(I6)') endif else begin oplot,alog10(wav),exthtot*norm oplot,alog10(wav),extthomson*norm,linestyle=2 endelse ; add curve identifiers xyouts,alog10(wav[0.4*nwav]),exthtot[0.4*nwav]*norm,$ ntostr(logpe,format='(f5.1)') if (extthomson[0.6*nwav]*norm gt miny) then $ xyouts,alog10(wav[0.6*nwav]),extthomson[0.6*nwav]*norm,$ ntostr(logpe,format='(f5.1)') endfor ; end loop over Pe curves per Vitense plot if (psplot) then closepsplot,psfilename,opengv=0 endfor ; end loop over Vitense graphs end