; file: saha_one.pro = obtain Saha fraction for one stage ; init: Mar 28 2015 Rob Rutten Deil ; last: Apr 2 2015 Rob Rutten Deil function saha_one,elemnr,temp,eldens,ionstage,$ range=range,verbose=verbose ;+ ; NAME: ; saha_one ; PURPOSE: ; compute LTE ion population fraction for one stage ; DESCRIPTION: ; uses Chianti or WSA data ; CALL: ; saha=saha_one(elemnr,temp,eldens,ionstage,range=range,verbose=verbose) ; INPUTS: ; elemnr = element number (hydrogen=1) ; temp = temperature [K] (scalr or array) ; eldens = electrion density [cm^-3] (scalar or same-size array) ; ionstage = ion stage number (0 for neutral, one less than "spectrum") ; OPTIONAL KEYWORD INPUTS: ; range: do Saha for ions [ionstage-range:ionstage+range] (default 2) ; verbose 1/0: print intermediate results (default 0) ; OUTPUTS: ; result = Saha population fraction, array if temp and/or eldens array ; RESTRICTIONS: ; needs SSW Chianti and other Rob Rutten ltelib programs ; maybe some Chianti data availability issues ; HISTORY: ; Apr 1 2015 RR: start ;- ; answer no-parameter query if (n_params(0) lt 4) then begin print,'xx=saha_one(elemnr,temp,eldens,ionstage,$' print,' range=range,verbose=verbose)' return,-1 ;RR return,-1 for a function endif ; defaults for keywords if (not keyword_set(verbose)) then verbose=0 if (not keyword_set(range)) then range=2 ; set actual range to not exceed the stage limits ionmin=ionstage-range if (ionmin lt 0) then ionmin=0 ionmax=ionstage+range if (ionmax gt elemnr) then ionmax=elemnr ; various checks if (elemnr gt 30) then begin print,' #### ABORT saha_one: element not in Chianti (only 1-30)' return,-2 endif ntemp=n_elements(temp) neldens=n_elements(eldens) if (ntemp gt 1 and neldens gt 1 and ntemp ne neldens) then begin print,' #### ABORT saha_one: temp and eldens inequal array sizes' return,-3 endif ; get Saha fractions (1D or 2D array) sahafrac=saha_all(elemnr,temp,eldens,$ ionmin=ionmin,ionmax=ionmax,verbose=verbose) ;; result=sahafrac[ionstage,*] result=sahafrac[ionstage-ionmin,*] if (n_elements(result) eq 1) then result=result[0] return,result end ; =============== main for testing per IDLWAVE H-c ====================== elemnr=56 ; Ba (error mssage test, not in Chianti) elemnr=20 ; Ca elemnr=26 ; Fe elemnr=14 ; Si ionstage=3 verbose=0 eldens=[1.E14,1.E13,1.E12] ; check array input eldens=1.E14 temp=[5000.,4500,4000.] ; check array input temp=25000 sahafrac=saha_one(elemnr,temp,eldens,ionstage,verbose=verbose) help,sahafrac print,' ==== saha_one: ',ntostr(sahafrac,format='(G15.3)') ; check against saha_rr.pro (first three stages only) atomdata_wsa,elemnr,atomwgt,abund,ion1,ion2,ion3,elemsymbol u0=partfunc_rr(temp,elemnr,0) u1=partfunc_rr(temp,elemnr,1) u2=partfunc_rr(temp,elemnr,2) saha_rr,temp,eldens,ion1,ion2,u0,u1,u2,n0_ntot,n1_ntot,n2_ntot print,' ==== saha_rr: ',ntostr([n0_ntot,n1_ntot,n2_ntot],format='(G15.3)') end