; -----------------------------------
;   
;   @mupho_example
;   
;  w1 contains a reduced data set. 
;
;
; -----------------------------------
RDSET,inst="IN4"
rdset, path = '../../rawdata/'

; data numbers for file naming
c = '189138-189167'

; -------------------------------------------------------------
; gdos (Bredov & Ostowskii) through the MUPHOCOR routine
; -------------------------------------------------------------
; -------------------------------------------------------------
; Create the MUPHOCOR input with various options.
; WARNING: Sample parameters to be adjusted for the UO2 case:
; amass: atomic mass = 238.03 + 2*16.0 = 270.03 g            (per formula unit) (90.01 g/at)
; sig  : total scattering cross section: 8.908+ 2* 4.232 = 13.14 barns (per formula unit) (4.38 barns/at)
; conc : concentration if one component: = 1
; unt  : flat background
; 
; -------------------------------------------------------------
s = 'input_mupho.txt'

print, '--------------------------------------------------------------'
print, ' UNT = ', strtrim(string(g),2)
print, '--------------------------------------------------------------'

; -------------------------------------------------------------
;  the gdos sensible to   emax and unt
; -------------------------------------------------------------
w59 = write_mupho(w1, file=s, emax=70.0, abs=0.0, phimin=10.0, phimax=115.5, iemp=0, amass=270.03, sig=13.14, unt=g)

; For check
;  w21 = total(w1,2)
;  w22 = read_mupho(file=s) & see, w22, /below


; -------------------------------------------------------------
; Compute the GDOS with the MUPHOCOR routine
; -------------------------------------------------------------

; without multi-phonon corrections (/first_guess)
; -------------------------------------------------------------
w10 = muphocor(s, /first_guess) & see, w10, /below

; with multi-phonon corrections
; -------------------------------------------------------------
w11 = muphocor(s) & see, w11, /below
y_tit(10) = 'GDOS [a.u.]'
y_tit(11) = 'GDOS [a.u.]'
      
; -------------------------------------------------------------
; plotting  (better thsan the basics above)     
; -------------------------------------------------------------
p = plot(x10, w10, '-db', xrange=[0,140], yrange=[0, 0.12], name='wo. multi-ph corr.')
q = plot(x11, w11, '-tur', xrange=[0,140], yrange=[0, 0.12], name='with multi-ph corr.', /overplot)

p.xtitle='Energy [meV]'
p.ytitle='G(w) [a.u.]'
p.XTickLen=1.0
p.YTickLen=1.0
p.XGridStyle=1
p.YGridStyle=1 & p.xminor=0 & p.yminor=0
p.title= 'IN4 GDOS - UO2, 1.7 $\AA$ ,  T ~ 446-346 K'
l= LEGEND(TARGET= [p,q], POSITION=[140, 0.115], /DATA, /AUTO_TEXT_COLOR)
p.save, "Figures/gdos_xxx.png", BORDER=10, RESOLUTION=300 ;;, /TRANSPARENT


      
; -------------------------------------------------------------
;  write in column format
; -------------------------------------------------------------
WRITE_LAMP,"gdos_muphocor_"+c, w=10, format="column"     
WRITE_LAMP,"gdos_muphocor_multi_"+c, w=11, format="column"     

; -------------------------------------------------------------
;  write in hdf
; -------------------------------------------------------------
s = "gdos_muphocor_"+c
WRITE_LAMP, s, w=10, format="hdf"     
spawn,"mv -f "+s+"_LAMP.hdf "+s+".hdf"
s = "gdos_muphocor_multi_"+c
WRITE_LAMP, s, w=11, format="hdf"     
spawn,"mv -f "+s+"_LAMP.hdf "+s+".hdf"