Commit 85085fdd authored by Miguel Angel Gonzalez's avatar Miguel Angel Gonzalez
Browse files

Update some IN5 files from /home/cs/lambda/macros/IN5/JOR_LIBRARY

parent a9afb575
FUNCTION in5_gdos, w_in, Temp=temp, Emin=emin, Emax=emax, $ FUNCTION in5_gdos, w_in, Temp=temp, Emin=emin, Emax=emax, $
method=method, cutoff=cutoff, nbbins=nbbins,$ method=method, cutoff=cutoff, nbbins=nbbins,$
bgcorr=bgcorr, verbose=verbose, plotting=plotting bgcorr=bgcorr, verbose=verbose, plotting=plotting
;------------------------------------------------------------------------------------ ;------------------------------------------------------------------------------------
; ;
;+ ;+
; ;
;<h2>NAME:</h2> ;<h2>NAME:</h2>
; gdos ; gdos
; ;
; ;
; <h2>PURPOSES:</h2> ; <h2>PURPOSES:</h2>
; ;
; Calculates the generalized density of states G(w) at a temperature T ; Calculates the generalized density of states G(w) at a temperature T
; in the case of the incoherent scattering and each time the incoherent ; in the case of the incoherent scattering and each time the incoherent
; approximation holds. Does not manage multi-phonons. ; approximation holds. Does not manage multi-phonons.
; ;
; ;
;<h2>COMMAND SYNTAX:</h2> ;<h2>COMMAND SYNTAX:</h2>
; w2 = gdos(w1, Emin=&lt;Emin&gt;, Emax=&lt;Emax&gt;, Temp=&lt;T&gt;, method=&lt;method&gt;, cutoff=&lt;cutoff&gt;, nbbins=&lt;nbbins&gt;,/bgcorr,/verbose,/plotting) ; w2 = gdos(w1, Emin=&lt;Emin&gt;, Emax=&lt;Emax&gt;, Temp=&lt;T&gt;, method=&lt;method&gt;, cutoff=&lt;cutoff&gt;, nbbins=&lt;nbbins&gt;,/bgcorr,/verbose,/plotting)
; ;
;<h2>EXAMPLES:</h2> ;<h2>EXAMPLES:</h2>
; w2 = gdos(w1) ; w2 = gdos(w1)
; ;
; With all default parameters ; With all default parameters
; ;
; w2 = gdos(w1,emax=50.0) ; w2 = gdos(w1,emax=50.0)
; ;
; gdos from E = 0 to 50 meV ; gdos from E = 0 to 50 meV
; ;
; w2 = gdos(w1, nbbins = 100, emax = 50.0, /plotting, /verbose) ; w2 = gdos(w1, nbbins = 100, emax = 50.0, /plotting, /verbose)
; With a given number of points in the gdos, an enery max, with a plotting ; With a given number of points in the gdos, an enery max, with a plotting
; of the extrapolation fit in 10 groups, with some info on the prompt ; of the extrapolation fit in 10 groups, with some info on the prompt
; ;
; ;
; w2 = gdos(w1, nbbins = 100, emax = 50.0, cutoff = 1.3) ; w2 = gdos(w1, nbbins = 100, emax = 50.0, cutoff = 1.3)
; With a given number of points in the gdos, a enery max, with ; With a given number of points in the gdos, a enery max, with
; a given cutoff for the extrapolation (default is 1.5 A-1). ; a given cutoff for the extrapolation (default is 1.5 A-1).
; ;
; ;
; w2 = gdos(w1, nbbins = 100, emax = 50.0, cutoff = 1.1, method = 'fitlog', /plotting) ; w2 = gdos(w1, nbbins = 100, emax = 50.0, cutoff = 1.1, method = 'fitlog', /plotting)
; With a given number of points in the gdos, a enery max, with ; With a given number of points in the gdos, a enery max, with
; a given cutoff for the extrapolation (default is 1.5 A-1).and with ; a given cutoff for the extrapolation (default is 1.5 A-1).and with
; the fitlog method instead of the default Egelstaff method. ; the fitlog method instead of the default Egelstaff method.
; ;
; ;
; ;
; The calculus is the following: ; The calculus is the following:
; ;
; p(a,b)=S(2theta,w)*Q^2/w/n(w,T) and ; p(a,b)=S(2theta,w)*Q^2/w/n(w,T) and
; ;
; lim p(a,b) = G(w) ; lim p(a,b) = G(w)
; Q^2->0 ; Q^2->0
; ;
; ;
; Workspace data must be as a function of energy transfer. Energies in meV. ; Workspace data must be as a function of energy transfer. Energies in meV.
; y should contain scattering angles (or averaged scattering angle) ; y should contain scattering angles (or averaged scattering angle)
; ;
; <h2>ARGUMENTS:</h2> ; <h2>ARGUMENTS:</h2>
; (all optional) ; (all optional)
; ;
; <b>Emin</b>: Minimum energy transfer to consider ; <b>Emin</b>: Minimum energy transfer to consider
; ;
; Default: MIN(E) ; Default: MIN(E)
; ;
; ;
; <b>Emax</b>: Maximum energy transfer to consider ; <b>Emax</b>: Maximum energy transfer to consider
; ;
; Default: MAX(E) ; Default: MAX(E)
; ;
; Always >= 0 and Emax>Emin: for IN5 and IN6 the up-scattering only is ; Always >= 0 and Emax>Emin: for IN5 and IN6 the up-scattering only is
; considered. For IN4: downscattering side has to be chosen. ; considered. For IN4: downscattering side has to be chosen.
; ;
; ;
; <b>Temp</b>: The temperature at which the experiment was performed ; <b>Temp</b>: The temperature at which the experiment was performed
; if not in the parameters. ; if not in the parameters.
; ;
; Default: From the parameters or 295K if not in the parameters ; Default: From the parameters or 295K if not in the parameters
; ;
; ;
; <b>method</b>: methods holds for the extrapolation method for: ; <b>method</b>: methods holds for the extrapolation method for:
; ;
; lim p(a,b) = G(w) (a prop. to Q^2) ; lim p(a,b) = G(w) (a prop. to Q^2)
; Q^2->0 ; Q^2->0
; ;
; method = 'egelstaff' is for the (early) Egelstaff approximation ; method = 'egelstaff' is for the (early) Egelstaff approximation
; with the Q_cutoff = threshold (straight line fit of S(a)/a with) ; with the Q_cutoff = threshold (straight line fit of S(a)/a with)
; a cutoff. ; a cutoff.
; Ref.: Egelstaff PA. & Schofield P. Nucl. Sci. Eng. 12(1962)260 ; Ref.: Egelstaff PA. & Schofield P. Nucl. Sci. Eng. 12(1962)260
; ;
; method = 'fitlog' is an "improvement" of the above method by fitting ; method = 'fitlog' is an "improvement" of the above method by fitting
; a straight line to the natural log of S(a)/a with or without a cutoff. ; a straight line to the natural log of S(a)/a with or without a cutoff.
; Depending on the estimated multiple scattering let's try one or the other ; Depending on the estimated multiple scattering let's try one or the other
; method. ; method.
; Derived from: JC. Smith et al. J. Chem. Phys. 85(1886)3636 ; Derived from: JC. Smith et al. J. Chem. Phys. 85(1886)3636
; ;
; Default: 'egelstaff' method with a cutoff of 1.5 A-1 ; Default: 'egelstaff' method with a cutoff of 1.5 A-1
; ;
; ;
; <b>cutoff</b>: a threshold for the low-Q cutoff (e.g. 1.0, 1.5 A-1 , ...) ; <b>cutoff</b>: a threshold for the low-Q cutoff (e.g. 1.0, 1.5 A-1 , ...)
; ;
; Default: 1.5 for 'egelstaff' and 'fitlog' method ; Default: 1.5 for 'egelstaff' and 'fitlog' method
; This value is okay relative to cold neutrons instruments. ; This value is okay relative to cold neutrons instruments.
; ;
; <b>nbbins</b>: a slicing allowing to rebin in less than NbMax points ; <b>nbbins</b>: a slicing allowing to rebin in less than NbMax points
; in energy (for clarity of the presentation only) ; in energy (for clarity of the presentation only)
; ;
; Default: no binning ; Default: no binning
; ;
; ;
; <b>/plotting</b> allow one to evaluate by eye the quality of the ; <b>/plotting</b> allow one to evaluate by eye the quality of the
; extrapolation on fit samples (for 10 spectra max). ; extrapolation on fit samples (for 10 spectra max).
; ;
; Default: no plotting ; Default: no plotting
; ;
; ;
; <b>/verbose</b> write some information about what is doing the routine ; <b>/verbose</b> write some information about what is doing the routine
; ;
; Default: no verbose ; Default: no verbose
; ;
; ;
;<h2>VERSION HISTORY:</h2> ;<h2>VERSION HISTORY:</h2>
; ;
; ToDo: compute the error bar, mean angle by testing the instrument, ; ToDo: compute the error bar, mean angle by testing the instrument,
; background correction, mutli-phonons treatment. ; background correction, mutli-phonons treatment.
; ;
; written by: S. Rols 11/01 rols@gdpc.univ-montp2.fr ; written by: S. Rols 11/01 rols@gdpc.univ-montp2.fr
; revisions: SR 04.02.02, 07/08/03 ; revisions: SR 04.02.02, 07/08/03
; JOR 29-Oct-2004 (ollivier@ill.fr) New version. ; JOR 29-Oct-2004 (ollivier@ill.fr) New version.
; makes the extrapolation (previous versions above = P(a)/a) ; makes the extrapolation (previous versions above = P(a)/a)
; included lots of new functionalities and keywords: ; included lots of new functionalities and keywords:
; methods Egelstaff & "fitlog" ; methods Egelstaff & "fitlog"
; with Q-cutoff, energy rebinning,... ; with Q-cutoff, energy rebinning,...
; JOR 02-Dec-2004: added error bars ; JOR 02-Dec-2004: added error bars
; JOR 25-May-2005: added the down-scattering possibility ; JOR 25-May-2005: added the down-scattering possibility
; in the energy sorting (IN4). ; in the energy sorting (IN4).
; ;
; JOR, Mon May 17 13:36:34 CEST 2010: Add prefix in5_ ; JOR, Mon May 17 13:36:34 CEST 2010: Add prefix in5_
; ;
; ;
; ;
;- ;-
; ;
;--------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------
COMMON c_lamp_access, inst COMMON c_lamp_access, inst
COMMON printing, iprint, outstring COMMON printing, iprint, outstring
;--------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------
; ;
; Analyses the entries ; Analyses the entries
; ;
;--------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------
if ~keyword_set(method) then begin IF ~keyword_set(method) THEN BEGIN
method = 'egelstaff' method = 'egelstaff'
if keyword_set(verbose) then $ IF keyword_set(verbose) THEN $
print,'GDOS: Warning: Method not given, set to ',method print,'GDOS: Warning: Method not given, set to ',method
endif ENDIF
if ( ~strcmp(method,'egelstaff',/fold_case) ) && ( ~strcmp(method,'fitlog',/fold_case) ) then begin IF ( ~strcmp(method,'egelstaff',/fold_case) ) && ( ~strcmp(method,'fitlog',/fold_case) ) THEN BEGIN
method = 'egelstaff' method = 'egelstaff'
if keyword_set(verbose) then $ IF keyword_set(verbose) THEN $
print,'GDOS: Warning: Undefined method, set to ',method print,'GDOS: Warning: Undefined method, set to ',method
endif ENDIF
if ~keyword_set(cutoff) then begin IF ~keyword_set(cutoff) THEN BEGIN
cutoff = 1.5 cutoff = 1.5
if keyword_set(verbose) then $ IF keyword_set(verbose) THEN $
print,'GDOS: Warning: Cutoff not given, set to ',strtrim(cutoff,2),' A-1' print,'GDOS: Warning: Cutoff not given, set to ',strtrim(cutoff,2),' A-1'
endif ENDIF
if keyword_set(bgcorr) then begin IF keyword_set(bgcorr) THEN BEGIN
if keyword_set(verbose) then $ IF keyword_set(verbose) THEN $
print,'GDOS: Background correction not yet implemented...' print,'GDOS: Background correction not yet implemented...'
endif ENDIF
take_datp, datp take_datp, datp
par=datp.p par=datp.p
x = datp.x x = datp.x
y = datp.y y = datp.y
err = datp.e err = datp.e
w = w_in w = w_in
lambda=0. lambda=0.
CASE strlowcase(inst) OF CASE strlowcase(inst) OF
'mibemol':begin 'mibemol':BEGIN
lambda=FLOAT(par[0]) lambda=FLOAT(par[0])
temp=FLOAT(par[59]) temp=FLOAT(par[59])
ei=81.805/par[0]^2 + x*0. ei=81.805/par[0]^2 + x*0.
end END
'dcsasc':begin 'dcsasc':BEGIN
lambda=float(par[8]) lambda=float(par[8])
ei=81.799/(lambda^2)+x*0. ei=81.799/(lambda^2)+x*0.
end END
'qens_raw':begin 'qens_raw':BEGIN
ef=3.167 ;average final neutron energy in mev ef=3.167 ;average final neutron energy in mev
lambda=sqrt(81.799/ef) lambda=sqrt(81.799/ef)
ei=x+ef ei=x+ef
end END
else:begin ELSE:BEGIN
lambda=par[21] lambda=par[21]
if (n_elements(temp) eq 0) then temp=par[11] IF (n_elements(temp) EQ 0) THEN temp=par[11]
ei=81.799/(lambda^2)+x*0. ei=81.799/(lambda^2)+x*0.
end END
endcase ENDCASE
if (temp eq 0.0) then begin IF (temp EQ 0.0) THEN BEGIN
temp=293.0 temp=293.0
if keyword_set(verbose) then $ IF keyword_set(verbose) THEN $
print,'verboseTemperature unknown, T set to: ',strtrim(temp,2),' K' print,'verboseTemperature unknown, T set to: ',strtrim(temp,2),' K'
endif ENDIF
if keyword_set(verbose) then begin IF keyword_set(verbose) THEN BEGIN
print, 'GDOS: Instrument is ', strupcase(inst) print, 'GDOS: Instrument is ', strupcase(inst)
print, 'GDOS: lambda = ',strtrim(lambda,2),' Angs' print, 'GDOS: lambda = ',strtrim(lambda,2),' Angs'
print, 'GDOS: Ei = ', strtrim(mean(ei),2),' meV' print, 'GDOS: Ei = ', strtrim(mean(ei),2),' meV'
endif ENDIF
;------------------------------------------------------ ;------------------------------------------------------
; Sort energies in ascending order with up-scattering ; Sort energies in ascending order with up-scattering
; positive, i.e, the condition hw = Ef - Ei ; positive, i.e, the condition hw = Ef - Ei
; ;
; Warning: Problem if IN4, i.e. working in down ; Warning: Problem if IN4, i.e. working in down
; scattering mode ; scattering mode
;------------------------------------------------------ ;------------------------------------------------------
if strlowcase(inst) eq 'in4' then begin IF strlowcase(inst) EQ 'in4' THEN BEGIN
endif else begin ENDIF ELSE BEGIN
a1 = x[0] & a2 = x[n_elements(x)-1] a1 = x[0] & a2 = x[n_elements(x)-1]
if abs(a1) gt abs(a2) then begin IF abs(a1) GT abs(a2) THEN BEGIN
if a1 lt a2 then begin IF a1 LT a2 THEN BEGIN
ind = sort(-1.0*x) ind = sort(-1.0*x)
x = -1.0*x[ind] x = -1.0*x[ind]
w = w[ind,*] w = w[ind,*]
err = err[ind,*] err = err[ind,*]
endif else begin ENDIF ELSE BEGIN
x = reverse(x) x = reverse(x)
err = reverse(err) err = reverse(err)
w = reverse(w) w = reverse(w)
endelse ENDELSE
endif else begin ENDIF ELSE BEGIN
if a1 gt a2 then begin IF a1 GT a2 THEN BEGIN
ind = sort(x) ind = sort(x)
x = x[ind] x = x[ind]
w = w[ind,*] w = w[ind,*]
err = err[ind,*] err = err[ind,*]
endif ENDIF
endelse ENDELSE
endelse ENDELSE
;---------------------------------------------------- ;----------------------------------------------------
; reduce the number of points if any ; reduce the number of points if any
;---------------------------------------------------- ;----------------------------------------------------
if (n_elements(emax) EQ 0) then emax=MAX(x) IF (n_elements(emax) EQ 0) THEN emax=MAX(x)
if (n_elements(emin) EQ 0) then emin=0.0 IF (n_elements(emin) EQ 0) THEN emin=0.0
if emin lt 0. or emax lt 0. then begin IF emin LT 0. OR emax LT 0. THEN BEGIN
print,'Error , Emin and Emax must be >= 0.0' print,'Error , Emin and Emax must be >= 0.0'
return,0 return,0
endif ENDIF
if emin gt emax then begin IF emin GT emax THEN BEGIN
print,'Error , Emax must be > Emin' print,'Error , Emax must be > Emin'
return,0 return,0
endif ENDIF
ind = where(x ge emin and x le emax) ind = where(x GE emin AND x LE emax)
x = x[ind] x = x[ind]
err = err[ind,*] err = err[ind,*]
w = w[ind,*] w = w[ind,*]
;------------------------------------------------------- ;-------------------------------------------------------
if keyword_set(verbose) then begin IF keyword_set(verbose) THEN BEGIN
print,'GDOS: Temperature=',strtrim(temp,2),' K' print,'GDOS: Temperature=',strtrim(temp,2),' K'
print,'GDOS: emin=',strtrim(emin,2),' meV' print,'GDOS: emin=',strtrim(emin,2),' meV'
print,'GDOS: emax=',strtrim(emax,2),' meV' print,'GDOS: emax=',strtrim(emax,2),' meV'
endif ENDIF
;------------------------------------------------------- ;-------------------------------------------------------
;Calcul de Q^2 ;Calcul de Q^2
;------------------------------------------------------- ;-------------------------------------------------------
q2=w*0.0 q2=w*0.0
x=float(x) x=float(x)
FOR i=0,n_elements(y)-1 DO begin FOR i=0,n_elements(y)-1 DO BEGIN
if y[i] eq 0.0 then y[i]= 61.96 ; IN6 only, 'instr' must be tested IF y[i] EQ 0.0 THEN y[i]= 61.96 ; IN6 only, 'instr' must be tested
q2[*,i]=2*(2*!PI/lambda)*(1+0.5*x[*]/ei[*] $ q2[*,i]=2*(2*!PI/lambda)*(1+0.5*x[*]/ei[*] $
- sqrt(1+x[*]/ei[*])*cos(y[i]*!PI/180)) - sqrt(1+x[*]/ei[*])*cos(y[i]*!PI/180))
endfor ENDFOR
alpha = 2.0721*q2/temp alpha = 2.0721*q2/temp
;------------------------------------------------------------- ;-------------------------------------------------------------
; Calcul de w/n(w) ou n(w) est le facteur de temperature dans ; Calcul de w/n(w) ou n(w) est le facteur de temperature dans
; le cas Stokes et anti Stokes ; le cas Stokes et anti Stokes
;------------------------------------------------------------- ;-------------------------------------------------------------
beta = 11.605*x/temp beta = 11.605*x/temp
bosex= beta*(exp(beta)-1) bosex= beta*(exp(beta)-1)
;----------------------------------------------------- ;-----------------------------------------------------
; Corrections ; Corrections
;----------------------------------------------------- ;-----------------------------------------------------
indnul=where(abs(bosex) LE 1.e-12) indnul=where(abs(bosex) LE 1.e-12)
if n_elements(indnul) gt 1 then begin IF n_elements(indnul) GT 1 THEN BEGIN
w[indnul,0:n_elements(y)-1]=0.000001 w[indnul,0:n_elements(y)-1]=0.000001
bosex[indnul]=1. bosex[indnul]=1.
end END
;------------------------------------------------------ ;------------------------------------------------------
; The main transform: p(a,b) ; The main transform: p(a,b)
;------------------------------------------------------ ;------------------------------------------------------
for i=0,n_elements(y)-1 do begin FOR i=0,n_elements(y)-1 DO BEGIN
err[*,i]=err[*,i]/w[*,i] err[*,i]=err[*,i]/w[*,i]
w[*,i]=w[*,i]*bosex[*]/alpha[*,i] w[*,i]=w[*,i]*bosex[*]/alpha[*,i]
err[*,i]=err[*,i]*w[*,i] err[*,i]=err[*,i]*w[*,i]
endfor ENDFOR
;-------------------------------------------------------
; slicing in energy if requested
;------------------------------------------------------- ;-------------------------------------------------------
if keyword_set(nbbins) then begin ; slicing in energy if requested
nchannels = n_elements(x) ;-------------------------------------------------------
if nbbins gt nchannels then nbbins = nchannels IF keyword_set(nbbins) THEN BEGIN
if keyword_set(verbose) then $ nchannels = n_elements(x)
IF nbbins GT nchannels THEN nbbins = nchannels
IF keyword_set(verbose) THEN $
print,'GDOS: Slicing energies in ',strtrim(nbbins,2),' values' print,'GDOS: Slicing energies in ',strtrim(nbbins,2),' values'
xq = fltarr(nbbins) xq = fltarr(nbbins)
dxq = fltarr(nbbins) dxq = fltarr(nbbins)
ww = fltarr(nbbins,n_elements(y)) ww = fltarr(nbbins,n_elements(y))
qq2 = fltarr(nbbins,n_elements(y)) qq2 = fltarr(nbbins,n_elements(y))
dbins = float(nchannels)/float(nbbins) dbins = float(nchannels)/float(nbbins)
dxq = x[1:nchannels-1] - x[0:nchannels-2] dxq = x[1:nchannels-1] - x[0:nchannels-2]
dxq[0] = dxq[1] dxq[0] = dxq[1]
dxq = dbins*dxq[uint(dbins*indgen(nbbins))] dxq = dbins*dxq[uint(dbins*indgen(nbbins))]
en = emin en = emin
for i=0,nbbins-1 do begin FOR i=0,nbbins-1 DO BEGIN
for j=0,n_elements(y)-1 do begin FOR j=0,n_elements(y)-1 DO BEGIN
ind = where((x gt en) and (x lt en+dxq[i])) ind = where((x GT en) AND (x LT en+dxq[i]))
if (n_elements(ind) gt 0) && (ind[0] ge 0) then begin IF (n_elements(ind) GT 0) && (ind[0] GE 0) THEN BEGIN
xq[i] = mean(x[ind],/double) xq[i] = mean(x[ind],/double)
ww[i,j] = mean(w[ind,j],/double) ww[i,j] = mean(w[ind,j],/double)
qq2[i,j]= mean(q2[ind,j],/double) qq2[i,j]= mean(q2[ind,j],/double)
endif ENDIF
endfor ENDFOR
en = en + dxq[i] en = en + dxq[i]
endfor ENDFOR
w = ww w = ww
x = xq