function popsol,temp,pres,chianti,n_e=n_e,delec=delec,z=z,jon=jon,$
radtemp=radtemp,dilute=dilute,pe_ratio=pe_ratio,$
n_levels=n_levels,data_str=data_str
;+
;function popsol
; returns the level populations at given temperature and
; pressure for a set of ions of a given element as an array
; of size [N(TEMP),N(PRES or DENSITY),N(N_LEVELS)]
;
;usage
; pop=popsol(temp,pres,chianti,n_e=n_e,/delec,$
; radtemp=radtemp,dilute=dilute,pe_ratio=pe_ratio,$
; n_levels=n_levels,data_str=data_str)
;
;parameters
; temp [INPUT; required] Temperature [K]
; pres [INPUT; required] pressure [cm^-3 K]
; chianti [INPUT; required] structure that contains all the relevant
; information gleaned from the CHIANTI database (type
; DBHELP=RD_CHIANTI() & HELP,DBHELP,/STRUCTURE
; for a description of what this structure contains)
;
;keywords
; n_e [INPUT] electron density in cm^-3
; * if set, overrides PRES/TEMP
; delec [INPUT] if set, assumes these are DR lines
; z [INPUT] atomic number of element
; jon [INPUT] ionic state of species
; * Z and JON are used only if DELEC is set
; radtemp [INPUT] The blackbody radiation field temperature
; * default is 6000 K
; dilute [INPUT] a dilution factor for photoionization, based
; on the distance above a star's surface -- see CHIANTI
; routine R2W()
; * if not set, assumed to be 0
; pe_ratio [INPUT] proton-to-electron ratio
; * if not given, assumed to be 1.0
; n_levels [INPUT] number of levels in the model
; * (from CHIANTI/idl/low_level/pop_solver.pro)
; This allows the number of levels in the model to be reduced.
; E.g., if the full model contains 100 levels, one could set
; N_LEVELS=50. This can be useful if one is interested in
; looking at the effects of cascading from higher levels.
; data_str [OUTPUT] some potentially useful data
; * (from CHIANTI/idl/low_level/pop_solver.pro)
; If POP_SOLVER is called for just 1 temperature and density,
; then the individual data arrays for each of the physical
; processes can be output through DATA_STR. This allows the
; user to check for the dominant processes affecting the
; population of a given level. DATA_STR is a structure with
; the following tags:
; .aa A-values
; .aax Photoexcitation/stimulated emission
; .cc Electron rate coefficients
; .ccp Proton rate coefficients
; Each tag is a 2D array, and is such that, e.g.,
; aa[0,20] corresponds to an excitation, while aa[20,0]
; is a de-excitation.
;
;A note about the default value of PE_RATIO:
; The original version of this program,
; CHIANTI/dbase/idl/low_level/pop_solver.pro,
; set the default to be 0.83, which is the limit for a high-T
; plasma with cosmic abundances. Later versions include a
; calculated value passed in via a common block. Here we
; exclude this factor because we want to store the emissivities
; without reference to their ion balance and abundance, on which
; this factor depends. The intensity of a line transition from
; upper level j to lower level i is in general
; dI(lam_ij) = (1/4pi)(hc/lam_ij) N_j A_ji dV [ergs/s/sr]
; where lam_ij is the wavelength of the transition, A_ji is the
; Einstein coefficient (the spontaneous transition probability),
; dV is the emitting volume, and N_j is the number density of the
; ions in level j. We can write N_j as the product
; N_j = (N_j/N_XIon) (N_XIon/N_X) (N_X/N_H) (N_H/N_e) N_e
; where N_XIon is the number density of the ions,
; N_j/N_XIon is the population fraction of the ions at level j,
; N_XIon/N_X is the ion fraction,
; N_X/N_H is the relative abundance of element X wrt H, and
; N_H/N_e is the H density relative to the free-electron density.
;
; In PINTofALE, we calculate and store the emissivities as
; e_ji = (hc/lam_ij) A_ji (1/N_e) (N_j/N_XIon)
; and write the line intensity as
; f(lam_ij) = e_ji (N_XIon/N_X) (N_X/N_H) EM (N_H/N_e)
; where EM = N_e^2 V is the emission measure, and the factor
; (1/4pi) is always ignored and the factor (hc/lam_ij) is usually
; taken out again. This calculation is carried out in LINEFLX(),
; which provides for the inclusion of the ion balance, abundance,
; and the proton-to-electron ratio. Hence the use of a default of
; 1.0 for PE_RATIO here.
;
;subroutines
; CHIANTI routine DESCALE_ALL
; NENH
;
;history
; POP_SOLVER.PRO written by Peter Young
;{..
;.. Ver 1, PRY 29-Mar-99
;.. Ver 2, PRY 30-Apr-99, added call to get_prot_rates
;.. Ver 3, PRY 15-Dec-99, added deu to upsilon common block in order
;.. to be consistent with the main Chianti routines.
;.. Ver 4, PRY 9-May-00, corrected problem with threshold when dealing
;.. with sparse matrices. Basically values less than 1.e-30 in
;.. the c-matrix were being set to zero and giving rise to
;.. NaN's in certain circumstances.
;.. Ver.5, PRY 14-Jul-00, changed elvl common block to the elvlc common
;.. block which is now the Chianti standard. Also, when
;.. descaling upsilons, the routine now uses the Delta-E from
;.. the .splups file.
;.. Ver.6, PRY 9-Aug-00, changed routine to deal better with the
;.. dielectronic recombination files
;.. Ver.7, PRY 17-Aug-00, routine does not call LINBCG now if radtemp
;.. is non-zero.
;.. Ver.8, PRY 29-Aug-00, the sparse matrix section has been disabled.
;.. Ver.9, PRY 12-Nov-01, calls routine proton_dens() to calculate the
;.. proton to electron ratio.
;.. Ver.10, PRY, 6-Dec-01, corrected bug when there are more levels
;.. in .splups file than in .elvlc file (ZnXXV).
;.. Ver.11, PRY, 11-Jul-02, removed ION keyword
;.. Ver.12, PRY, 9-Aug-02, within the equation solving section, I've set
;.. the population of the ground level (rather the n_level level)
;.. to 1, and this seems to stop negative populations appearing
;.. in extreme conditions.
;.. Ver.12, PRY, 21-Aug-02, changed exp(-1/1/a) to exp(-a) in electron
;.. excitation section which caused a hang-up in some
;.. circumstances. Also, the routine now uses vector ECMC
;.. (combined experimental and theoretical energies) in
;.. determining if a level lies above or below another level.
;.. Previously only used the observed energy vector. Also, the
;.. exponential in the electron excitation section now uses the
;.. (accurate) .elvlc energy separation rather than the .splups
;.. energy separation, which can cause significant (~20-30%)
;.. differences in level populations of high-lying levels at
;.. low temperatures.
;.. Ver.13, PRY, 10-Sep-02, corrected bug for proton rates. The excitation
;.. and de-excitation rates were being swapped.
;..
;.. V. 14 4-Oct-2003 Giulio Del Zanna (GDZ).
;.. -removed all COMMON blocks (note that only proton_dens.pro has
;.. one: COMMON elements,abund,abund_ref,ioneq,ioneq_logt,ioneq_ref)
;.. -only the essential information input is passed to the routine
;.. via a new input structure.
;.. -fixed a bug, that affected all the satellite lines, and was
;.. introduced in v.12, included in CHIANTI v.4.0.
;.. basically the ionization potential was not subtracted when
;.. calculating the Delta E in the exponential.
;..
;.. VERSION : 14, 4-Oct-2003
;..}
; modified to eliminate common blocks and the call to SETUP_ION,
; and altered the I/O interface to mesh with PINTofALE, changed
; call to PROTON_DENS() to call NENH() and changed default PE_RATIO
; to 1.0 (Vinay Kashyap, Apr03)
; modified to account for v14, for CHIANTI v4.2 (VK; Apr04)
;-
; usage
ok='ok' & np=n_params(0)
nt=n_elements(temp) & npr=n_elements(pres) & nch=n_tags(chianti)
if np lt 3 then ok='Insufficient parameters' else $
if nt eq 0 then ok='Temperature: not defined' else $
if npr eq 0 and not keyword_set(n_e) then ok='Pressure: not defined' else $
if nch eq 0 then ok='CHIANTI structure is undefined'
if ok ne 'ok' then begin
print,'Usage: pop=popsol(temperature,pressure,chianti_structure,$'
print,' n_e=edensity,radtemp=radiation_temperature,$'
print,' dilute=dilution_factor,pe_ratio=proton_to_electron_ratio,$'
print,' n_levels=n_levels,data_str=data_str)'
print,' returns level populations at given temperature and pressure'
if np ne 0 then message,ok,/info
return,-1L
endif
;------------------------------------------------------
;{VK the following bit is how POP_SOLVER is set up
;------------------------------------------------------
;;these are required:
;;-------------------
;
; gname = input.gname
; jj= input.jj
; ecm= input.ecm
; ecmth= input.ecmth
; wvl= input.wvl
; a_value= input.a_value
; splstr= input.splstr
;
;;optional ones:
;;-------------------
;
;IF tag_exist(input, 'radtemp') THEN radtemp= input.radtemp
;IF tag_exist(input, 'dilute') THEN dilute= input.dilute
;IF tag_exist(input, 'prot_struc') THEN prot_struc= input.prot_struc
;IF tag_exist(input, 'pe_ratio') THEN pe_ratio= input.pe_ratio
;
;
;;COMMON elvlc,l1a,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref
;;COMMON wgfa, wvl,gf,a_value
;;COMMON upsilon, splstr
;;COMMON radiative, radtemp,dilute
;;COMMON proton, prot_struc, pe_ratio
;
;IF n_params(0) LT 4 THEN BEGIN
; print,' use> pop_solver,input, temperature,density,populations, $'
; print,' [n_levels=n_levels, data_str=data_str] '
; return
;ENDIF
;
;convertname,gname,iz,ion
;ion2spectroscopic,gname,snote, dielectronic=dielectronic
;
;if dielectronic then begin
;read_ip,concat_dir(concat_dir(!xuvtop, 'ip'), 'chianti.ip'),ionpot,ipref
;ip=ionpot(iz-1,ion-1)
;endif else ip=0
;
;
;xne = DOUBLE(xne)
;t = DOUBLE(t)
;;
;; need the following to turn t into an array if it only has 1 element
;;
;IF n_elements(t) EQ 1 THEN BEGIN
; t0=t
; t=dblarr(1)
; t[0]=t0
;ENDIF
;------------------------------------------------------
;VK: thus ends POP_SOLVER setup}{below is PoA translation
;------------------------------------------------------
; inputs
t=[double(temp)]
if keyword_set(n_e) then xne=double(n_e) else xne=double(pres/t)
; decode CHIANTI variables
l1=chianti.lev1 & l2=chianti.lev2
w1=chianti.wvl & g1=chianti.gf & a1=chianti.a
nlev=max([l1,l2]) & wvl=fltarr(nlev,nlev) & gf=wvl & a_value=wvl
;
o1=where(w1 eq 0,mo1) & o2=where(w1 ne 0,mo2)
wvl(l1-1,l2-1)=abs(w1) & gf(l1-1,l2-1)=g1
if mo1 gt 0 then begin
a_value(l1(o1)-1,l2(o1)-1) = a1(o1)
if mo2 gt 0 then a_value(l1(o2)-1,l2(o2)-1)=a_value(l1(o2)-1,l2(o2)-1)+a1(o2)
endif else begin
if mo2 gt 0 then a_value(l1(o2)-1,l2(o2)-1)=a1(o2)
endelse
;
ecm=chianti.ecmobs & ecmth=chianti.ecmthr & jj=chianti.j & mult=2*jj+1
splstr=chianti.splstr
prot_struc=chianti.psplstr
;
deu=chianti.deryd
t_type=chianti.trans & c_ups=chianti.c_ups & splups=chianti.spups
ionpot=chianti.ip
if not keyword_set(delec) then ip=0 else begin
iz=1 & if keyword_set(Z) then iz=fix(Z[0])
ion=1 & if keyword_set(jon) then ion=fix(jon[0])
ip=ionpot(iz-1,ion-1)
endelse
; special keywords
; RADTEMP
if not keyword_set(radtemp) then radtemp=6000.D
; DILUTE
if not keyword_set(dilute) then dilute=0.
; PE_RATIO
if not keyword_set(pe_ratio) then pe_ratio=dblarr(nt)+1.0
;------------------------------------------------------
;VK: thus ends PoA setup}
;------------------------------------------------------
;------------------------------------------------------
;{VK: from here on down is uninterrupted POP_SOLVER.PRO
; {except for the bit about calling PROTON_DENS()
;------------------------------------------------------
ecmc=ecm
ind=where(ecm EQ 0.)
IF ind[0] NE -1 THEN ecmc[ind]=ecmth[ind]
mult=2.*jj+1.
;
hck=1.98648d-16/1.38062d-16
ryd2cm=109737.31534d
;
n_elvl=n_elements(ecm)
wsize=size(wvl)
n_wgfa1=wsize(1)
n_wgfa2=wsize(2)
usize=max([splstr.lvl1,splstr.lvl2])
;
IF N_ELEMENTS(n_levels) EQ 0 THEN n_levels=min([n_elvl,n_wgfa1,n_wgfa2,usize])
;
;
c=dblarr(n_levels,n_levels)
d=dblarr(n_levels,n_levels)
b=dblarr(n_levels)
diag=dblarr(n_levels)
nt=N_ELEMENTS(t) ; no. of temperatures
nxne=N_ELEMENTS(xne) ; no. of densities
pop=dblarr(nt,nxne,n_levels)
;;------------------------------[]
; The arrays
;
; e.g., aa(0,19) will be zero (no 0 -> 19 A value)
; aa(19,0) will be non-zero
;
; qq(0,19) electron excitation
; qq(19,0) electron de-excitation
;;------------------------------[]
ident=make_array(n_levels,val=1.) ; use for making de arrays
ecmn=ecmc(0:n_levels-1) ; n_levels
den=ecmn#ident & dem=ident#ecmn
aat=a_value(0:n_levels-1,0:n_levels-1) ; transpose of aa
aa=DOUBLE(TRANSPOSE(aat))
aax=c
IF N_ELEMENTS(dilute) EQ 0 THEN dilute=0.
;;------------------------------------------------------------------[-]
; The following loads up the photoexcitation (pexc) and stimulated
; emission (stem) arrays)
;
IF dilute NE 0. THEN BEGIN
stem=c & pexc=c & ede=c
;
multn=mult[0:n_levels-1] ; in case mult and ecm are bigger than
;
mm=TRANSPOSE(multn#(1/multn)) ; needed for photoexcitation
;
dd=ABS(den-dem)*hck/radtemp
ede=exp(dd) - 1.
;
ind=WHERE( (aat NE 0.) AND (ede NE 0.) )
pexc[ind]=aat[ind]*dilute*mm[ind]/ede[ind]
;
ind=where( (aa NE 0.) AND (ede NE 0.) )
stem[ind]=aa[ind]*dilute/ede[ind]
;
aax=pexc+stem
ENDIF
;;------------------------------------------------------------------[-]
;________________
; create a ppr array for the proton rates
;
ppr=MAKE_ARRAY(n_levels,n_levels,nt,/double)
IF n_tags(prot_struc) NE 0 THEN BEGIN
;
;------------------------------------------------------
;{VK: replace call to PROTON_DENS() by NENH()
;------------------------------------------------------
;
IF (n_elements(pe_ratio) NE nt) THEN BEGIN
print,'%POP_SOLVER: WARNING, pe_ratio size does not match temp'
print,n_elements(pe_ratio),nt
; pe_ratio=proton_dens(alog10(t))
pe_ratio=1./(nenh(alog10(t),/nproton)>1.)
ENDIF
;
;------------------------------------------------------
;VK: replaced call to PROTON_DENS() by NENH()}
; and from here on out, back to uninterrupted POP_SOLVER}
;------------------------------------------------------
;
;
FOR i=0,n_elements(prot_struc)-1 DO BEGIN
l1=prot_struc[i].lvl1-1
l2=prot_struc[i].lvl2-1
de=ABS(prot_struc[i].de)
descale_all,t,prot_struc,i,prate
IF ecmc(l1) LT ecmc(l2) THEN BEGIN
ppr[l1,l2,*]=prate*pe_ratio
ppr[l2,l1,*]=prate*pe_ratio*mult[l1]/mult[l2]* $
exp(de*13.61/8.617/10.^(-5)/t)
ENDIF ELSE BEGIN
ppr[l2,l1,*]=prate*pe_ratio
ppr[l1,l2,*]=prate*pe_ratio*mult[l2]/mult[l1]* $
exp(de*13.61/8.617/10.^(-5)/t)
ENDELSE
ENDFOR
ENDIF
;________________
;______________
; Create a qq array for electron rates
;
qq=MAKE_ARRAY(n_levels,n_levels,nt,/double)
;
l1=splstr.lvl1-1
l2=splstr.lvl2-1
;GDZ- added ip
kte=(hck*abs(ecmc[l1]-(ecmc[l2]-ip))) # (1d0/t)
;******************************************
; kte=(hck*abs(ecmc[l1]-ecmc[l2])) # (1d0/t)
ind_pos=where(ecmc[l2] GT ecmc[l1])
ind_neg=where(ecmc[l2] LT ecmc[l1])
xx=dblarr(n_elements(splstr),nt)
yy=xx
;
; xx and yy contain all factors in the expression for the rate coefficient,
; except for the upsilon. They can be generated using array operations - the
; upsilons need a for loop.
;
IF ind_neg[0] NE -1 THEN BEGIN
xx[ind_neg,*]=(8.63d-6/(mult[l1[ind_neg]])#(1./sqrt(t)))
yy[ind_neg,*]=8.63d-6* exp(-kte[ind_neg,*]) * $
1./( mult[l2[ind_neg]] # sqrt(t) )
ENDIF
IF ind_pos[0] NE -1 THEN BEGIN
yy[ind_pos,*]=(8.63e-6/mult[l2[ind_pos]]) # (1./sqrt(t))
xx[ind_pos,*]=8.63e-6* exp(-kte[ind_pos,*]) * $
1./(mult[l1[ind_pos]] # sqrt(t))
ENDIF
;
; this is the for loop for the upsilons
;
FOR i=0,n_elements(splstr)-1 DO BEGIN
IF (l1[i] LE n_levels-1) AND (l2[i] LE n_levels-1) THEN BEGIN
descale_all,t,splstr,i,ups
qq[l1[i],l2[i],*]=xx[i,*]*ups
qq[l2[i],l1[i],*]=yy[i,*]*ups
ENDIF
ENDFOR
;______________
FOR j=0,nt-1 DO BEGIN ;; loop over temperatures
;
; We are solving the equation Ax=0 where A is the matrix c, and x is a
; vector containing the level populations. There are an infinity of
; solutions to Ax=0 of the form kx',
; where k is real. This is the case here, and so we throw away one of
; the linear equations, and replace it with an equation setting the
; ground level population to 1. Once the linear equations have been
; solved, the populations are renormalised such that the sum of the level
; populations is 1.
; The equation we throw away is the first row in A.
; Note with the normalisation, we are solving Ax=b, where all elements
; of b are zero, except the first which is one.
;
; Most of the c arrays we are dealing with in CHIANTI are sparse, with
; many zero elements. IDL has special routines for dealing with sparse
; matrices, and we use this below. We define frac_not_zero to denote the
; number of non-zero elements in c. If this is less than a specified number,
; then we switch to using the sparse matrix routine SPRSIN.
;
FOR i=0,nxne-1 DO BEGIN ;; loop over densities
;
c=aa + xne(i)*qq[*,*,j] + xne(i)*ppr[*,*,j] + aax ; creates the c matrix
;
ind=where(c NE 0.)
frac_not_zero=float(n_elements(ind))/float(n_elements(c))
;
diag = -total(c,2)
c[findgen(n_levels),findgen(n_levels)] = diag
;
c(*,0)=0d0 ; set this row to zeros...
c(0,0)=1d0 ; ...except for ground level
;
b(0)=1d0 ; b is zero except for last element
;
;
; Currently (v4 of CHIANTI) the sparse matrix solver (SPRSIN) is not used
; in CHIANTI so the IF statement below is set to always use the INVERT
; routine.
;
IF frac_NOT_zero GT 0. THEN BEGIN
pp1=invert(transpose(temporary(c)),status)#b
ENDIF ELSE BEGIN
pp1=LINBCG(SPRSIN(temporary(c),thresh=1d-40),b,b)
ENDELSE
;
pp=pp1/total(pp1) ; force total population to be 1
pop[j,i,*]=pp ; temp first, then dens
ENDFOR
ENDFOR
IF nt EQ 1 AND nxne EQ 1 THEN $
data_str={aa: aa, aax: aax, cc: xne[0]*qq[*,*,0], ccp: xne[0]*ppr[*,*,0]}
;-------------------------------------------------------------
;VK: above was uninterrupted POP_SOLVER.PRO}
;-------------------------------------------------------------
return,pop
end