FUNCTION  EMISS_CALC, IZ, ION, TEMP=TEMP, DENS=DENS, RADT=RADT, DIL=DIL, $
                           PATH=PATH, NO_DE=NO_DE, PROTON=PROTON, $
                           QUIET=QUIET, PRESSURE=PRESSURE, DIEL=DIEL

;+
;
; PROJECT:  CHIANTI
;
;       CHIANTI is an atomic database package for the calculation of
;       astrophysical emission line spectra.  It is a collaborative project
;       involving Ken Dere (Naval Research Laboratory, Washington DC), 
;       Brunella Monsignori-Fossi and Enrico Landi (Arcetri Observatory, 
;       Florence), and Helen Mason and Peter Young (DAMTP, Cambridge Univ.).
;
;
; NAME: EMISS_CALC()
;       
; PURPOSE:
;
;       To compute the emissivities of all lines of a specified ion over
;	given ranges of temperature and density.
;
; CATEGORY:
;
;       Scientific analysis
;
; EXPLANATION:
;
; 	This routine calculates:
;
;		hc
;       	--  * N_j  * A_ji
;      	       lamb
;
; 	where hc = 1.986 * 10^-8 erg AA, lamb is in angstroms, N_j is the 
;	fraction of ions in the upper emitting level j, and A_ji is the 
;	radiative decay rate for the transition.
;
;	The emissivities are stored in a structure called EMISS that also
;	holds the wavelength of the transition, the level numbers i and j 
;	and also a 'flag', which is set to -1 if the wavelength is negative.
;
;	The temperature and density ranges can be specified directly using
;	the TEMP and DENS keywords. Setting TMAX to the log T_max of the 
;	ion, gives emissivities for 3 temperatures: log T_max +- 0.15. 
;	If DENS is not set, then it is set to 8 to 12 in 0.5 dex intervals. 
;	STDENS allows the start density (of 8) to be changed to some other 
;	value; ND allows the number of densities to be varied (default 9); 
;	DINT allows the density interval to be varied (default 0.5).
;
; CALLING SEQUENCE:
;
;       EMISS=EMISS_CALC (IZ, ION, [ TEMP=TEMP, DENS=DENS, RADT=RADT, $
;				DIL=DIL, PATH=PATH, /NO_DE, /PROTON, $
;				QUIET, PRESSURE=PRESSURE)
;
; EXAMPLES:
;
;	EMISS=EMISS_CALC(26,13)
;	EMISS=EMISS_CALC(26,13,temp=[6.2],dens=findgen(5)+8)
;	EMISS=EMISS_CALC(26,13,temp=findgen(11)/100.+5.5,press=10.^15)
;
; INPUTS:
;
;	IZ:	The atomic number of the ion
;	ION:	The spectroscopic number of the ion (e.g., 12 = XII)
;
; OPTIONAL INPUTS:
;
;	TEMP:	Direct specification of the temperature range (log T)
;	DENS:	Direct specification of the density range (log Ne)
;	RADT:	Specify background radiation temperature (default: 6000 K)
;	DIL:	Specify radiation dilution factor (default: 0.0)
;	PATH:	If specified, the routine will look for the atomic data in 
;		the PATH directory, rather than in the CHIANTI database
;	PRESSURE:	If a temperature array is given, and PRESSURE set, 
;			then the emissivities will be evaluated at the 
;			specified temperatures, but for densities = 
;			pressure/temperature. If DENS is set, then it will 
;			be ignored. The pressure is assumed to be in units 
;			K * cm^-3.
;
;
; KEYWORDS:
;
;	NO_DE:	Drops the hc/lambda factor in the computation of the 
;		emissivities. Useful for emission measure analyses involving 
;		photon fluxes
;	PROTON: Allows proton rates to be added to level balance if the 
;		.prot file exists.
;	QUIET:	If set, don't list the temperatures and densities at which 
;		the emissivities are caculated.
;       DIEL:   If the dielectronic recombination files exist for the ion, 
;               then these are used to derive the emissivities.
;
; OUTPUT:
;
;	The structure that is output has the following tags:
;
;	.ion_name	string; contains ion name, e.g., 'Fe XIII'
;	.lambda		float; contains wavelength
;	.level1		integer; contains lower level of transition
;	.level2		integer; contains upper level of transition
;	.flag		integer; a flag to mark particular transitions
;	.em		fltarr(nt,nd); contains emissivities at particular 
;			  temperatures and densities.
;
; PROGRAMMING NOTES:
;
;	Transitions where only theoretical energies are available for at 
;	least one of the two levels are assigned negative wavelengths in 
;	the .wgfa file. With emiss_calc, the wavelength is set to be 
;	positive, but emiss.flag is set to -1 for that transition, so 
;	that other routines can keep track of the theoretical wavelengths.
;
; COMMON BLOCKS:
;
;	ELVLC, WGFA, UPSILON, RADIATIVE, PROTON
;
; CALLS:
;
;	SETUP_ION, POP_SOLVER
;	ZION2NAME, ZION2FILENAME, ZION2SPECTROSCOPIC, CHOOSE_IONEQ, 
;	READ_PROT
;
; HISTORY:
;
;	Ver 1, PRY 28-Jun-97
;	Ver 2, PRY 26-Jul-97  - corrected problem with size of emiss
;	Ver 3, PRY 22-Sep-97  - allowed photo-excitation to be included
;	Ver 4, PRY 6-Jul-98   - added PATH
;	Ver 5, PRY 5-Sep-98   - added call to choose_ioneq
;       Ver 6, PRY 3-Dec-98   - dosen't crash if no params given
;	Ver 7, PRY 9-Jan-99   - allowed proton rates to be added through 
;				/PROTON keyword.
;	Ver 8, PRY 10-Feb-99  - added /QUIET keyword
;	Ver 9, PRY 8-Oct-99   - for H-like ions, there's a 2-photon 
;				transition with a non-zero A-value that is 
;				assigned a zero wavelength (as it does not 
;				produce an emission line). This caused 
;				problems for dens_plotter, so emiss_calc 
;				now removes this transition from the emiss 
;				structure.
;	Ver 10, PRY 15-Dec-99 -	added deu to the upsilon common block in 
;				order to be consistent with the main Chianti 
;				routines.
;	Ver 11, PRY 8-May-00  - added PRESSURE
;       Ver 12, PRY 17-Aug-00 - changed elvlc common block to match new 
;                               version of pop_solver
;       Ver 13, PRY 10-Oct-00 - now calls setup_ion to read ion data
;
;
; CONTACT:
;
;	Peter Young, CfA, pyoung@cfa.harvard.edu
;-


COMMON elvlc,l1,term,conf,ss,ll,jj,ecm,eryd,ecmth,erydth,eref
COMMON wgfa, wvl,gf,a_value
COMMON upsilon,t_type,deu,c_ups,splups
COMMON radiative, radtemp,dilute
COMMON proton, prot_struc


;-------------------------------[]
; There was a problem when using this routine in the CDS environment as 
; !xuvtop was not set, and an error would be generated. The following 
; prevents this problem
;
chianti_var='!xuvtop'
defsysv,chianti_var,exist=exist
;
IF exist EQ 0 THEN xuvtop_help,chianti_var
;-------------------------------[]


IF N_PARAMS() LT 2 THEN BEGIN
  PRINT,'Use:  IDL> emiss=emiss_calc(iz,ion [, temp=temp, dens=dens, '
  PRINT,'                              radt=radt, dil=dil, /no_de,   '
  PRINT,'                              path=path, /proton, /quiet, $ '
  PRINT,'                              pressure=pressure ] )         '
  RETURN,0
ENDIF

;
; The following extracts the names of the files to be read. I want to allow 
; a different path to be chosen, and I extract only the information I need 
; from filename.
;
IF keyword_set(diel) THEN diel = 1 ELSE diel = 0
zion2name,iz,ion,name,diel=diel
zion2filename,iz,ion,filename,diel=diel
IF N_ELEMENTS(path) NE 0 THEN BEGIN
  filename=path+'/'+name
ENDIF


setup_ion,name,-1,-1,wvltst,lvl1,lvl2,wvl1,gf1,a_value1,path=path

pname=filename+'.prot'


;-----------------------<>
; Read proton data, if required
;
IF KEYWORD_SET(proton) THEN read_prot,pname,prot_struc
;-----------------------<>


IF N_ELEMENTS(radt) NE 0 THEN radtemp=radt ELSE radtemp=6000.
IF N_ELEMENTS(dil) NE 0 THEN dilute=dil ELSE dilute=0.


;------------------------------------+
; If TEMP not specified, then set TEMP=[log(T_max),log(T_max)+-0.15]
; using the ion balance data in CHIANTI. 
;
IF N_ELEMENTS(temp) EQ 0 THEN BEGIN
  choose_ioneq,iz,ion,temp_all,ioneq
 ;
  f_all=ioneq(*,iz-1,ion-1)
  ind=where(f_all EQ max(f_all))
  tmax=temp_all(ind) & tmax=tmax(0)
  temp=[tmax-0.15,tmax,tmax+0.15]
ENDIF
nt=N_ELEMENTS(temp)
;------------------------------------+


;;-----------------<>
; Default density range is 8 to 12 in .5 intervals:
;
IF NOT KEYWORD_SET(dens) THEN dens=findgen(9)/2. + 8.
nd=n_elements(dens)
;;-----------------<>


IF N_ELEMENTS(pressure) NE 0 THEN BEGIN
  nd=1
  dens=ALOG10(pressure/10.^temp)
ENDIF

;----------X
; Find the size of the population array
;
pop_solver,10.^6,10.^9,pop
pop=pop(0,0,*)
nlvls2=n_elements(pop)
;----------X


;--------------------O
; Adjust a_value and wvl, and work out ntrans2
;
a_value2 = a_value(0:nlvls2-1 , 0:nlvls2-1)
wvl2     = wvl(0:nlvls2-1 , 0:nlvls2-1)
;
; 8-Oct-99, I've added the wvl2 bit below to prevent problems with H-like ions
;
ind_A=WHERE((a_value2 NE 0.) AND (wvl2 NE 0.))
ntrans2=n_elements(ind_A)
;--------------------O

;----------------------------------o
; Create the EMISS structure, using ntrans2 to set size
;
str={   ion_name: '', $
	lambda: 0., $
	level1: 0, $
	level2: 0, $
	flag:   0, $
	em: dblarr(nt,nd) $
        }
;
emiss=replicate(str,ntrans2)
;
zion2spectroscopic,iz,ion,species
emiss[0].ion_name=species
;----------------------------------o

upp=fltarr(nlvls2,nlvls2) & low=upp
FOR j=0,nlvls2-1 DO upp(*,j)=j
FOR i=0,nlvls2-1 DO low(i,*)=i



ind2=SORT(ABS(wvl2[ind_A]))


av=a_value2[ind_A[ind2]]
wa=wvl2[ind_A[ind2]]
lo=low[ind_A[ind2]]
up=upp[ind_A[ind2]]


emiss.lambda=ABS(wa)
emiss.level1=lo+1
emiss.level2=up+1
ind=WHERE(wa LT 0.)
IF ind(0) NE -1 THEN emiss(ind).flag=-1

; Note that emiss.flag is set if the wavelength is negative, while emiss.lambda
; is made positive
;
n=N_ELEMENTS(av)

IF NOT KEYWORD_SET(quiet) THEN BEGIN
  PRINT,''
  PRINT,'Log_10 temperatures...'
  PRINT,FORMAT='("   ",9f6.2)',temp
  PRINT,''
  PRINT,'Log_10 densities...'
  PRINT,FORMAT='("   ",9f6.2)',dens
ENDIF

IF N_ELEMENTS(pressure) NE 0 THEN BEGIN
  FOR j=0,nt-1 DO BEGIN
    pop_solver,10.^temp[j],10.^dens[j],pop,proton=proton
    FOR i=0,n-1 DO BEGIN
      po=pop(*,*,up(i))
      IF KEYWORD_SET(no_de) THEN $
            emiss(i).em[j]=po*av(i) ELSE $
            emiss(i).em[j]=1.986*10.^(-8)/ABS(wa(i)) * po*av(i)
    ENDFOR
  ENDFOR
ENDIF ELSE BEGIN
  pop_solver,10.^temp,10.^dens,pop,proton=proton
  FOR i=0,n-1 DO BEGIN
    po=pop(*,*,up(i))
    IF KEYWORD_SET(no_de) THEN $
          emiss(i).em=po*av(i) ELSE $
        emiss(i).em=1.986*10.^(-8)/ABS(wa(i)) * po*av(i)
  ENDFOR
ENDELSE

RETURN,emiss

END