in5_gdos.pro 14.6 KB
Newer Older
1
FUNCTION in5_gdos, w_in, Temp=temp, Emin=emin, Emax=emax, $
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
  method=method, cutoff=cutoff, nbbins=nbbins,$
  bgcorr=bgcorr, verbose=verbose, plotting=plotting
  ;------------------------------------------------------------------------------------
  ;
  ;+
  ;
  ;<h2>NAME:</h2>
  ; 	gdos
  ;
  ;
  ; <h2>PURPOSES:</h2>
  ;
  ; 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.
  ;
  ;
  ;<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)
  ;
  ;<h2>EXAMPLES:</h2>
  ;    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)
  ;
  ; <h2>ARGUMENTS:</h2>
  ;     (all optional)
  ;
  ;  <b>Emin</b>:   Minimum energy transfer to consider
  ;
  ;             Default: MIN(E)
  ;
  ;
  ;  <b>Emax</b>:   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.
  ;
  ;
  ;  <b>Temp</b>:   The temperature at which the experiment was performed
  ;             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:
  ;
  ;                     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
  ;
  ;
  ; <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
  ;                     This value is okay relative to cold neutrons instruments.
  ;
  ; <b>nbbins</b>: a slicing allowing to rebin in less than NbMax points
  ;                in energy (for clarity of the presentation only)
  ;
  ;                     Default: no binning
  ;
  ;
  ;  <b>/plotting</b> allow one to evaluate by eye the quality of the
  ;                       extrapolation on fit samples (for 10 spectra max).
  ;
  ;                       Default: no plotting
  ;
  ;
  ;  <b>/verbose</b> write some information about what is doing the routine
  ;
  ;                       Default: no verbose
  ;
  ;
  ;<h2>VERSION HISTORY:</h2>
  ;
  ;        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_
  ;
  ;
  ;
  ;-
  ;
  ;---------------------------------------------------------------------------------
145

146
147
  COMMON c_lamp_access, inst
  COMMON printing, iprint, outstring
148

149
150
151
152
153
154
155
156
  ;---------------------------------------------------------------------------------
  ;
  ;  Analyses the entries
  ;
  ;---------------------------------------------------------------------------------
  IF ~keyword_set(method) THEN BEGIN
    method = 'egelstaff'
    IF keyword_set(verbose) THEN $
157
      print,'GDOS: Warning: Method not given,  set to ',method
158
159
160
161
  ENDIF
  IF ( ~strcmp(method,'egelstaff',/fold_case) ) && ( ~strcmp(method,'fitlog',/fold_case) ) THEN BEGIN
    method = 'egelstaff'
    IF keyword_set(verbose) THEN $
162
      print,'GDOS: Warning: Undefined method,  set to ',method
163
164
165
166
  ENDIF
  IF ~keyword_set(cutoff) THEN BEGIN
    cutoff = 1.5
    IF keyword_set(verbose) THEN $
167
      print,'GDOS: Warning: Cutoff not given,  set to ',strtrim(cutoff,2),' A-1'
168
169
170
  ENDIF
  IF keyword_set(bgcorr) THEN BEGIN
    IF keyword_set(verbose) THEN $
171
      print,'GDOS: Background correction not yet implemented...'
172
  ENDIF
173
174


175
176
177
178
179
180
  take_datp, datp
  par=datp.p
  x = datp.x
  y = datp.y
  err = datp.e
  w = w_in
181

182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  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
204

205
  IF (temp EQ 0.0) THEN BEGIN
206
    temp=293.0
207
208
209
    IF keyword_set(verbose) THEN $
      print,'verboseTemperature unknown, T set to: ',strtrim(temp,2),' K'
  ENDIF
210

211
212
213
214
215
  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
216

217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
  ;------------------------------------------------------
  ;  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
247

248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
  ;----------------------------------------------------
  ;  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
282

283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
  ;-------------------------------------------------------------
  ; 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 $
314
      print,'GDOS: Slicing energies in ',strtrim(nbbins,2),' values'
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
    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
333
      en = en + dxq[i]
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
    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 $
351
      print, 'GDOS: Extrapolate with the ',method,' method and with Q_cutoff = ',strtrim(cutoff,2),' A-1'
352
353
354
355
356
357
358
359
360
361
    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
362
363
364
      qq2 = fltarr(n_elements(ind)) & qq2[*] = q2[ichan,ind]
      ww  = fltarr(n_elements(ind)) & ww[*]  = w[ichan,ind]
      ; linear regression (linear fit )
365
      ;     result = regress(qq2,ww,const=dos,sigma=err,/double)
366
367
      result = linfit(qq2,ww,sigma=err,/double)
      w2[ichan] = result[0]
368
      ;      e2(ichan) = err(0)*w2(ichan) ; Semble etre bon cette fois (2007-01-23)
369
      e2[ichan] = err[0] ; Semble etre bon cette fois (2007-01-23)
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
      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
385

386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
  ;---------------------------------------------------
  ;   '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