pro mipump

readcol, 'XSED.dat', wav, flx, f = 'd,d'
readcol, 'CI.dat', x1, x2, x3, x4, f = 'd,d,d,d'

gl = x2   ; An array; occupation number of level lambda.
g1 = 1.
g2 = 3.
g3 = 5.

c = 2.99792458D10  ; cm/s
h = 6.6260755D-27  ; erg.s
cte = (c^2.)/(2. * h)

Al1 = x4(where(x3 eq 1.0))
Al2 = x4(where(x3 eq 2.0))
Al3 = x4(where(x3 eq 3.0))

x1t = x1(UNIQ(x1, sort(x1)))+0.0D

n = n_elements(x1)
m = n_elements(x1t)

Alm = fltarr(m, 4) ; Alm = A^lambda_m
Blm = fltarr(m, 4)
Bml = fltarr(m, 4)
Ecm = fltarr(m, 4) ; energy in cm^-1.
Ean = fltarr(m, 4) ; energy in Angstrom.
Klm = fltarr(m, 4)
Kml = fltarr(m, 4)

for k = 0 , m - 1 do begin
for i = 0 , n - 1 do begin

Alm[k, 0] = x1t[k]; Energy [cm^-1]
Blm[k, 0] = x1t[k]; Energy [cm^-1]
Bml[k, 0] = x1t[k]; Energy [cm^-1]
Ecm[k, 0] = x1t[k]; Energy [cm^-1]
Ean[k, 0] = x1t[k]; Energy [cm^-1]
Klm[k, 0] = cmtoangs(x1t[k])
Kml[k, 0] = cmtoangs(x1t[k])

if (x1t[k] eq x1[i]) and (x3[i] eq 1) then begin
Alm[k,  1] = x4[i]
Blm[k,  1] = (cte / (cmtonu(x1t[k]))^3.0) * x4[i]
Bml[k,  1] = (gl[i]/g1) * Blm[k,  1]
Ecm[k,  1] = (x1t[k] - 0.0)          ; energy in cm^-1.
Ean[k,  1] = cmtoangs(x1t[k] - 0.0)  ; energy in Angstrom.
ntmp       = (closest((Ean[k,  1]), wav))(0)
Klm[k,  1] = flx[ntmp] * Blm[k,  1]
Kml[k,  1] = flx[ntmp] * Bml[k,  1]
endif

if (x1t[k] eq x1[i]) and (x3[i] eq 2) then begin
Alm[k,  2] = x4[i]
Blm[k,  2] = (cte / (cmtonu(x1t[k]))^3.0) * x4[i]
Bml[k,  2] = (gl[i]/g2) * Blm[k,  2]
Ecm[k,  2] = (x1t[k] - 16.40)
Ean[k,  2] = cmtoangs(x1t[k] - 16.40)
ntmp       = (closest((Ean[k,  2]), wav))(0)
Klm[k,  2] = flx[ntmp] * Blm[k,  2]
Kml[k,  2] = flx[ntmp] * Bml[k,  2]
endif

if (x1t[k] eq x1[i]) and (x3[i] eq 3) then begin
Alm[k,  3] = x4[i]
Blm[k,  3] = (cte / (cmtonu(x1t[k]))^3.0) * x4[i]
Bml[k,  3] = (gl[i]/g3) * Blm[k,  3]
Ecm[k,  3] = (x1t[k] - 43.40)          
Ean[k,  3] = cmtoangs(x1t[k] - 43.40)  
ntmp       = (closest((Ean[k,  3]), wav))(0)
Klm[k,  3] = flx[ntmp] * Blm[k,  3]
Kml[k,  3] = flx[ntmp] * Bml[k,  3]
endif

endfor
endfor  ; This is for k.

;%%%%%%%%%%%% calculating Gamma (Excitation rates by fluorescence) %%%%%%%%%%%%%%
gam12 = fltarr(m)
gam13 = fltarr(m)
gam23 = fltarr(m)
for l = 0 , m - 1 do begin
gamS = Kml[l, 1] * (Alm[l, 2] + Klm[l, 2])
gamM = (Alm[l, 1] + Klm[l, 1]) + (Alm[l, 2] + Klm[l, 2]) + (Alm[l, 3] + Klm[l, 3])
gam12[l] = gamS / gamM
endfor
for l = 0 , m - 1 do begin
gamS = Kml[l, 1] * (Alm[l, 3] + Klm[l, 3])
gamM = (Alm[l, 1] + Klm[l, 1]) + (Alm[l, 2] + Klm[l, 2]) + (Alm[l, 3] + Klm[l, 3])
gam13[l] = gamS / gamM
endfor
for l = 0 , m - 1 do begin
gamS = Kml[l, 2] * (Alm[l, 3] + Klm[l, 3])
gamM = (Alm[l, 1] + Klm[l, 1]) + (Alm[l, 2] + Klm[l, 2]) + (Alm[l, 3] + Klm[l, 3])
gam23[l] = gamS / gamM
endfor
print, 'total(gam12) = ', total(gam12)
print, 'total(gam13) = ', total(gam13)
print, 'total(gam23) = ', total(gam23)

end
