FUNCTION in5_gdos, w_in, Temp=temp, Emin=emin, Emax=emax, \$ method=method, cutoff=cutoff, nbbins=nbbins,\$ bgcorr=bgcorr, verbose=verbose, plotting=plotting ;------------------------------------------------------------------------------------ ; ;+ ; ;

; gdos ; ; ;

## PURPOSES:

; ; Calculates the generalized density of states G(w) at a temperature T ; in the case of the incoherent scattering and each time the incoherent ; approximation holds. Does not manage multi-phonons. ; ; ;

## COMMAND SYNTAX:

; w2 = gdos(w1, Emin=<Emin>, Emax=<Emax>, Temp=<T>, method=<method>, cutoff=<cutoff>, nbbins=<nbbins>,/bgcorr,/verbose,/plotting) ; ;

## EXAMPLES:

; w2 = gdos(w1) ; ; With all default parameters ; ; w2 = gdos(w1,emax=50.0) ; ; gdos from E = 0 to 50 meV ; ; 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 ; of the extrapolation fit in 10 groups, with some info on the prompt ; ; ; w2 = gdos(w1, nbbins = 100, emax = 50.0, cutoff = 1.3) ; 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). ; ; ; 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 ; a given cutoff for the extrapolation (default is 1.5 A-1).and with ; the fitlog method instead of the default Egelstaff method. ; ; ; ; The calculus is the following: ; ; p(a,b)=S(2theta,w)*Q^2/w/n(w,T) and ; ; lim p(a,b) = G(w) ; Q^2->0 ; ; ; Workspace data must be as a function of energy transfer. Energies in meV. ; y should contain scattering angles (or averaged scattering angle) ; ;

## ARGUMENTS:

; (all optional) ; ; Emin: Minimum energy transfer to consider ; ; Default: MIN(E) ; ; ; Emax: Maximum energy transfer to consider ; ; Default: MAX(E) ; ; Always >= 0 and Emax>Emin: for IN5 and IN6 the up-scattering only is ; considered. For IN4: downscattering side has to be chosen. ; ; ; Temp: The temperature at which the experiment was performed ; if not in the parameters. ; ; Default: From the parameters or 295K if not in the parameters ; ; ; method: methods holds for the extrapolation method for: ; ; lim p(a,b) = G(w) (a prop. to Q^2) ; Q^2->0 ; ; method = 'egelstaff' is for the (early) Egelstaff approximation ; with the Q_cutoff = threshold (straight line fit of S(a)/a with) ; a cutoff. ; Ref.: Egelstaff PA. & Schofield P. Nucl. Sci. Eng. 12(1962)260 ; ; 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. ; Depending on the estimated multiple scattering let's try one or the other ; method. ; Derived from: JC. Smith et al. J. Chem. Phys. 85(1886)3636 ; ; Default: 'egelstaff' method with a cutoff of 1.5 A-1 ; ; ; cutoff: a threshold for the low-Q cutoff (e.g. 1.0, 1.5 A-1 , ...) ; ; Default: 1.5 for 'egelstaff' and 'fitlog' method ; This value is okay relative to cold neutrons instruments. ; ; nbbins: a slicing allowing to rebin in less than NbMax points ; in energy (for clarity of the presentation only) ; ; Default: no binning ; ; ; /plotting allow one to evaluate by eye the quality of the ; extrapolation on fit samples (for 10 spectra max). ; ; Default: no plotting ; ; ; /verbose write some information about what is doing the routine ; ; Default: no verbose ; ; ;

## VERSION HISTORY:

; ; ToDo: compute the error bar, mean angle by testing the instrument, ; background correction, mutli-phonons treatment. ; ; written by: S. Rols 11/01 rols@gdpc.univ-montp2.fr ; revisions: SR 04.02.02, 07/08/03 ; JOR 29-Oct-2004 (ollivier@ill.fr) New version. ; makes the extrapolation (previous versions above = P(a)/a) ; included lots of new functionalities and keywords: ; methods Egelstaff & "fitlog" ; with Q-cutoff, energy rebinning,... ; JOR 02-Dec-2004: added error bars ; JOR 25-May-2005: added the down-scattering possibility ; in the energy sorting (IN4). ; ; JOR, Mon May 17 13:36:34 CEST 2010: Add prefix in5_ ; ; ; ;- ; ;--------------------------------------------------------------------------------- COMMON c_lamp_access, inst COMMON printing, iprint, outstring ;--------------------------------------------------------------------------------- ; ; Analyses the entries ; ;--------------------------------------------------------------------------------- IF ~keyword_set(method) THEN BEGIN method = 'egelstaff' IF keyword_set(verbose) THEN \$ print,'GDOS: Warning: Method not given, set to ',method ENDIF IF ( ~strcmp(method,'egelstaff',/fold_case) ) && ( ~strcmp(method,'fitlog',/fold_case) ) THEN BEGIN method = 'egelstaff' IF keyword_set(verbose) THEN \$ print,'GDOS: Warning: Undefined method, set to ',method ENDIF IF ~keyword_set(cutoff) THEN BEGIN cutoff = 1.5 IF keyword_set(verbose) THEN \$ print,'GDOS: Warning: Cutoff not given, set to ',strtrim(cutoff,2),' A-1' ENDIF IF keyword_set(bgcorr) THEN BEGIN IF keyword_set(verbose) THEN \$ print,'GDOS: Background correction not yet implemented...' ENDIF take_datp, datp par=datp.p x = datp.x y = datp.y err = datp.e w = w_in lambda=0. CASE strlowcase(inst) OF 'mibemol':BEGIN lambda=FLOAT(par[0]) temp=FLOAT(par[59]) ei=81.805/par[0]^2 + x*0. END 'dcsasc':BEGIN lambda=float(par[8]) ei=81.799/(lambda^2)+x*0. END 'qens_raw':BEGIN ef=3.167 ;average final neutron energy in mev lambda=sqrt(81.799/ef) ei=x+ef END ELSE:BEGIN lambda=par[21] IF (n_elements(temp) EQ 0) THEN temp=par[11] ei=81.799/(lambda^2)+x*0. END ENDCASE IF (temp EQ 0.0) THEN BEGIN temp=293.0 IF keyword_set(verbose) THEN \$ print,'verboseTemperature unknown, T set to: ',strtrim(temp,2),' K' ENDIF IF keyword_set(verbose) THEN BEGIN print, 'GDOS: Instrument is ', strupcase(inst) print, 'GDOS: lambda = ',strtrim(lambda,2),' Angs' print, 'GDOS: Ei = ', strtrim(mean(ei),2),' meV' ENDIF ;------------------------------------------------------ ; Sort energies in ascending order with up-scattering ; positive, i.e, the condition hw = Ef - Ei ; ; Warning: Problem if IN4, i.e. working in down ; scattering mode ;------------------------------------------------------ IF strlowcase(inst) EQ 'in4' THEN BEGIN ENDIF ELSE BEGIN a1 = x[0] & a2 = x[n_elements(x)-1] IF abs(a1) GT abs(a2) THEN BEGIN IF a1 LT a2 THEN BEGIN ind = sort(-1.0*x) x = -1.0*x[ind] w = w[ind,*] err = err[ind,*] ENDIF ELSE BEGIN x = reverse(x) err = reverse(err) w = reverse(w) ENDELSE ENDIF ELSE BEGIN IF a1 GT a2 THEN BEGIN ind = sort(x) x = x[ind] w = w[ind,*] err = err[ind,*] ENDIF ENDELSE ENDELSE ;---------------------------------------------------- ; reduce the number of points if any ;---------------------------------------------------- IF (n_elements(emax) EQ 0) THEN emax=MAX(x) IF (n_elements(emin) EQ 0) THEN emin=0.0 IF emin LT 0. OR emax LT 0. THEN BEGIN print,'Error , Emin and Emax must be >= 0.0' return,0 ENDIF IF emin GT emax THEN BEGIN print,'Error , Emax must be > Emin' return,0 ENDIF ind = where(x GE emin AND x LE emax) x = x[ind] err = err[ind,*] w = w[ind,*] ;------------------------------------------------------- IF keyword_set(verbose) THEN BEGIN print,'GDOS: Temperature=',strtrim(temp,2),' K' print,'GDOS: emin=',strtrim(emin,2),' meV' print,'GDOS: emax=',strtrim(emax,2),' meV' ENDIF ;------------------------------------------------------- ;Calcul de Q^2 ;------------------------------------------------------- q2=w*0.0 x=float(x) 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 q2[*,i]=2*(2*!PI/lambda)*(1+0.5*x[*]/ei[*] \$ - sqrt(1+x[*]/ei[*])*cos(y[i]*!PI/180)) ENDFOR alpha = 2.0721*q2/temp ;------------------------------------------------------------- ; Calcul de w/n(w) ou n(w) est le facteur de temperature dans ; le cas Stokes et anti Stokes ;------------------------------------------------------------- beta = 11.605*x/temp bosex= beta*(exp(beta)-1) ;----------------------------------------------------- ; Corrections ;----------------------------------------------------- indnul=where(abs(bosex) LE 1.e-12) IF n_elements(indnul) GT 1 THEN BEGIN w[indnul,0:n_elements(y)-1]=0.000001 bosex[indnul]=1. END ;------------------------------------------------------ ; The main transform: p(a,b) ;------------------------------------------------------ FOR i=0,n_elements(y)-1 DO BEGIN err[*,i]=err[*,i]/w[*,i] w[*,i]=w[*,i]*bosex[*]/alpha[*,i] err[*,i]=err[*,i]*w[*,i] ENDFOR ;------------------------------------------------------- ; slicing in energy if requested ;------------------------------------------------------- IF keyword_set(nbbins) THEN BEGIN 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' xq = fltarr(nbbins) dxq = fltarr(nbbins) ww = fltarr(nbbins,n_elements(y)) qq2 = fltarr(nbbins,n_elements(y)) dbins = float(nchannels)/float(nbbins) dxq = x[1:nchannels-1] - x[0:nchannels-2] dxq[0] = dxq[1] dxq = dbins*dxq[uint(dbins*indgen(nbbins))] en = emin FOR i=0,nbbins-1 DO BEGIN FOR j=0,n_elements(y)-1 DO BEGIN ind = where((x GT en) AND (x LT en+dxq[i])) IF (n_elements(ind) GT 0) && (ind[0] GE 0) THEN BEGIN xq[i] = mean(x[ind],/double) ww[i,j] = mean(w[ind,j],/double) qq2[i,j]= mean(q2[ind,j],/double) ENDIF ENDFOR en = en + dxq[i] ENDFOR w = ww x = xq q2 = qq2 ENDIF ;--------------------------------------------------- ; extrapolate the p(a,b) to the GDOS ;--------------------------------------------------- nchannels = n_elements(x) nspectra = n_elements(y) w2 = fltarr(nchannels) e2 = fltarr(nchannels) ;--------------------------------------------------- ; 'Egelstaff' method (with threshold) ;--------------------------------------------------- IF strcmp(method,'egelstaff',/fold_case) THEN BEGIN IF keyword_set(verbose) THEN \$ print, 'GDOS: Extrapolate with the ',method,' method and with Q_cutoff = ',strtrim(cutoff,2),' A-1' nmod = round(nchannels/10.0) cutoff2 = cutoff*cutoff FOR ichan=0,nchannels-1 DO BEGIN ind = where(q2[ichan,*] GE cutoff2) IF (ind[0] LT 0) THEN BEGIN cutoff2 = q2[ichan,n_elements(y)-3] ind = where(q2[ichan,*] GE cutoff2) IF keyword_set(verbose) THEN \$ print, 'GDOS: Warning: ',strtrim(ichan,2),' no Q values for extrapolation Q ->',strtrim(sqrt(cutoff2),2),' A-1' ENDIF qq2 = fltarr(n_elements(ind)) & qq2[*] = q2[ichan,ind] ww = fltarr(n_elements(ind)) & ww[*] = w[ichan,ind] ; linear regression (linear fit ) ; result = regress(qq2,ww,const=dos,sigma=err,/double) result = linfit(qq2,ww,sigma=err,/double) w2[ichan] = result[0] ; e2(ichan) = err(0)*w2(ichan) ; Semble etre bon cette fois (2007-01-23) e2[ichan] = err[0] ; Semble etre bon cette fois (2007-01-23) IF keyword_set(plotting) THEN BEGIN IF (ichan MOD nmod) EQ 0 THEN BEGIN qfit = indgen(101)*max(qq2)/100 fit=qfit fit=result[0]+result[1]*qfit plot, q2[ichan,*], w[ichan,*], \$ color=1324565, psym=4, \$ xrange = [0.0,max(q2[ichan,*])], \$ title='Method '+method+', spectrum '+string(ichan)+'/'+string(nchannels) oplot,qfit,fit wait,1.0 ENDIF ENDIF ENDFOR ENDIF ;--------------------------------------------------- ; 'fitlog' method (with/without threshold) ;--------------------------------------------------- IF strcmp(method,'fitlog',/fold_case) THEN BEGIN IF keyword_set(verbose) THEN \$ print, 'GDOS: Extrapolate with the ',method,' method and with Q_cutoff = ',strtrim(cutoff,2),' A-1' logw=alog(w) nmod = round(nchannels/10.0) cutoff2 = cutoff*cutoff FOR ichan=0,nchannels-1 DO BEGIN ind = where(q2[ichan,*] GE cutoff2) IF ( ind[0] LT 0 ) THEN BEGIN cutoff2 = q2[ichan,n_elements(y)-3] ind = where(q2[ichan,*] GE cutoff2) IF keyword_set(verbose) THEN \$ print, 'GDOS: Warning: no Q values for extrapolation Q ->',strtrim(sqrt(cutoff2),2),' A-1' ENDIF qq2 = fltarr(n_elements(ind)) & qq2[*] = q2[ichan,ind] ww = fltarr(n_elements(ind)) & ww[*] = logw[ichan,ind] ; linear regression (linear fit ) result =linfit(qq2,ww,sigma=err,/double) w2[ichan] = exp(result[0]) e2[ichan] = err[0]*w2[ichan] ; Semble etre bon cette fois (2007-01-23) ; e2(ichan) = err(0) ; Semble etre bon cette fois (2007-01-23) IF keyword_set(plotting) THEN BEGIN IF (ichan MOD nmod) EQ 0 THEN BEGIN qfit = indgen(101)*max(qq2)/100 fit=qfit fit=result[0]+result[1]*qfit plot, q2[ichan,*], w[ichan,*], \$ color=1324565, psym=4, \$ xrange = [0.0,max(q2[ichan,*])], \$ title='Method '+method+'spectrum '+string(ichan)+'/'+string(nchannels) oplot,qfit,exp(fit) wait,1.0 ENDIF ENDIF ENDFOR ENDIF ;---------------------------- ; output ;---------------------------- x_out = x e_out = e2 y_out = 0 w_out = w2 output:mod_datp, datp, "e", e_out mod_datp, datp, "x", x_out mod_datp, datp, "y", y_out give_datp, datp return, w2 END