in5_sqw_rebin.pro 14.7 KB
Newer Older
1
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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
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
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
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
FUNCTION in5_sqw_rebin, w_in, dQ = dQ, Emin = Emin0, all_angles = all_angles, $
  pos_angles = pos_angles, neg_angles = neg_angles, $
  swap_QE = swap_QE, qb, em, ib, $
  verbose = verbose
  ;--------------------------------------------------------------------------------
  ;********************************************************************************
  ;
  ;
  ; For IN4, IN5, IN6 and D7
  ;
  ;rebins output data from t2e and reb to regular-grid S(Q,w) data using the old
  ;KHA IN6 rebin algorithm. Proper rebionning routine with error analysis (unlike
  ;sqw_interp.pro).
  ;
  ;ARGUMENTS:
  ;	dQ :	Q bin width
  ;	Emin0:	Minimum energy value (meV) - neutron energy gain is negative
  ;
  ;KEYWORDS (- only for D7 data)
  ; /neg_angles	: use only negative angles
  ; /pos_angles	: use only positive angles
  ; /all_angles	: use all angles (default)
  ;		input workspace must be in energy transfer versus scattering angle,
  ;		i.e. only one component or spin phase.
  ; (qb, em and ib are obsolete, kept for backwards compatibility)
  ;
  ;DIMENSIONS:
  ; w_in(nE,nphi) -> w_out(nQs,nEs)
  ;
  ;COMMAND SYNTAX:
  ; w10=spw_rebin(w9,dQ=<dQ>, Emin=<Emin>[,/neg_angles][,/pos_angles][,/all_angles])
  ;
  ; (optional keywords shown in square brackets)
  ;
  ;
  ;
  ; WARNING: Seems to not work with the negative scattering angles of DCSD !!!
  ; First, remove the negative angles and then that's okay.
  ;
  ;
  ;
  ;
  ; Creation:     KHA,JRS 9/02/06
  ; Modification: JO,   Tue Nov 25 18:05:32 CET 2008 : Short-cut of the negative phi (is only for D7)
  ;                                                    IN4, IN5, IN6 only Debye-Scherrer cones Phi > 0
  ;
  ;
  ;--------------------------------------------------------------------------------
  ;********************************************************************************

  COMMON c_lamp_access, inst
  COMMON grid, Qmin, Qmax, Emin, Emax
  FORWARD_FUNCTION string_round, overlap

  iprint = 0	; if iprint>0, show debugging messages


  take_datp, datp

  ibank = 2
  ;  IF N_ELEMENTS(qb) GT 0 THEN dQ    = qb
  ;  IF N_ELEMENTS(em) GT 0 THEN Emin  = em
  ;  IF N_ELEMENTS(ib) GT 0 THEN ibank = ib
  ;
  IF KEYWORD_SET(pos_angles) THEN ibank   = 1
  IF KEYWORD_SET(neg_angles) THEN ibank   = 0
  IF KEYWORD_SET(all_angles) THEN ibank   = 2
  IF KEYWORD_SET(swap_QE)    THEN swap_QE = 1 ELSE swap_QE = 0

  no_small = 0
  IF (datp.y[0] GT 10.) AND (inst EQ 'IN4') THEN BEGIN
    no_small = 1
    PRINT, 'SQW_rebin: IN4 data without small angle bank'
  ENDIF

  ; -------------------------------------------------------------------------------
  ;   Set up starting parameters
  ; -------------------------------------------------------------------------------
  IF N_ELEMENTS(dQ) NE 1 THEN BEGIN
    ii = DIALOG_MESSAGE('SQW_rebin: Error - dQ must be specified', /ERROR)
    RETURN, w_in
  ENDIF

  IF N_ELEMENTS(Emin0) NE 1 THEN Emin0 = -1.E+10

  sw = SIZE(w_in)
  nx = sw[1]
  ny = sw[2]
  IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: nx = ', strtrim(string(nx),2),', ny = ', strtrim(string(ny),2)
  IF keyword_set(verbose) THEN PRINT,'SIZE(w_in) = ', strtrim(string(sw),2)

  IF sw[0] NE 2 THEN BEGIN
    s = 'SQW_REBIN: ERROR: input workspace must be 2-D: E vs. phi'
    ii = DIALOG_MESSAGE(s, /ERROR)
    RETURN, w_in
  ENDIF

  x_in = datp.x	& sx = SIZE(x_in)
  y_in = datp.y	& sy = SIZE(y_in)

;  print, sx
;  print, sy

  IF (nx NE sx[1]) OR (ny NE sy[1]) THEN BEGIN
    s = 'SQW_REBIN: sx = ' + STRTRIM(STRING(sx),2) + $
      ' sy = ' + STRTRIM(STRING(sy),2)
    ii = DIALOG_MESSAGE(s, /ERROR)
    RETURN, w_in
  ENDIF

  e_in = datp.e
  se   = SIZE(e_in)
  IF (se[0] NE sw[0] OR se[1] NE sw[1] OR se[2] NE sw[2]) THEN e_in=w_in*0.

  par=datp.p

  IF keyword_set(verbose) THEN $
    PRINT,'SQW_REBIN: Instrument = ',inst
  IF (inst EQ 'D7') THEN lambda = par[4] ELSE lambda = par[21]

  IF keyword_set(verbose) THEN $
    PRINT, FORMAT = '("SQW_REBIN: lambda     = ",F5.2," A")', lambda

  ; -------------------------------------------------------------------------------------
  ;	Set constants and prepare arrays for rebinning to regular Q-E grid
  ; -------------------------------------------------------------------------------------
  const1 = 5.22697		; E(meV)  = const1*V(m/ms)^2
  const2 = 2.07193571		; E(meV)  = const2*k(A^-1)^2
  const3 = 3.956076		; V(m/ms) = const3/lambda(A)
  const4 = 81.8066		; E(meV)  = const4/lambda(A)^2

  Ei   = const4 / lambda^2
  ki   = SQRT(Ei / const2)
  y_in = y_in*!PI/180.	; convert to radians

  nEps = nx + 1
  Eps  = FLTARR(nEps)
  Eps[0]      =  x_in[0] - (x_in[1] - x_in[0]) / 2.
  Eps[1:nx-1] = (x_in[0:nx - 2] + x_in[1:nx - 1]) / 2.
  Eps[nx]     =  x_in[nx - 1] + (x_in[nx - 1] - x_in[nx - 2]) / 2.

  iEarr = WHERE(x_in GT Emin0)
  Eps   = Eps[iEarr[0]:nx]
  nEps  = nx - iEarr[0] + 1
  Emin  = Eps[0]
  Emax  = Eps[nEps - 1]
  Qmin  = 0
  max_y = MAX(y_in)
  Qmax  = SQRT((2.*Ei - Emin - 2.*SQRT(Ei*(Ei - Emin))*COS(max_y))/const2)
  Qmax  = MAX([Qmax,ki])

  ;;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: iEarr(0) = ',iEarr(0),' nx = ',nx
  IF keyword_set(verbose) THEN PRINT, FORMAT = '("SQW_REBIN: Emin     = ",F7.2," meV")', Emin
  IF keyword_set(verbose) THEN PRINT, FORMAT = '("SQW_REBIN: Emax     = ",F7.2," meV")', Emax

  ;;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: SQW_rebin: E-array prepared'

  ; -----------------------------------------------------------------------------
  ;   angle grid generation
  ; -----------------------------------------------------------------------------
  nQ    = FIX((Qmax - Qmin) / dQ) + 1
  w_out = FLTARR(nQ,nEps)
  e_out = w_out - 1.
  Q     = Qmin + FLOAT(INDGEN(nQ))*dQ

  ;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: Qmin = ', Qmin
  ;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: Qmax = ', Qmax
  IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: nQ   = ', nQ
  IF keyword_set(verbose) THEN PRINT, FORMAT = '("SQW_REBIN: Qmin     = ",F7.2," A-1")', Qmin
  IF keyword_set(verbose) THEN PRINT, FORMAT = '("SQW_REBIN: Qmax     = ",F7.2," A-1")', Qmax
  IF keyword_set(verbose) THEN $
    PRINT, format = '("SQW_REBIN: Phi_in  = [",F7.3,",",F7.3,"]")',min(y_in)*180./!pi,max(y_in)*180./!pi

  IF inst EQ 'D7' THEN BEGIN
    isort = SORT(y_in)
    y_buf1 = y_in[isort]
    w_buf1 = w_in[*,isort]
    e_buf1 = e_in[*,isort]
    i     = WHERE(y_buf1 GT 0.,n)

    IF ibank EQ 2 THEN BEGIN
      twice = 1
      iphi1 = 0		& iphi1next = i[0]
      iphi2 = i[0] - 1	& iphi2next = ny - 1
    ENDIF ELSE BEGIN
      twice = 0
      IF ibank EQ 0 THEN BEGIN
        iphi1 = 0
        iphi2 = i[0] - 1
      ENDIF ELSE IF ibank EQ 1 THEN BEGIN
        iphi1 = i[0]
        iphi2 = ny - 1
      ENDIF
    ENDELSE
  ENDIF ELSE BEGIN
    twice = 0
    iphi1 = 0	& iphi2 = ny - 1
  ENDELSE

  start:
  IF inst EQ 'D7' THEN BEGIN
    nphi          = iphi2 - iphi1 + 2
    phi           = FLTARR(nphi)
    phi[0]        = y_buf1[iphi1] - (y_buf1[iphi1 + 1] - y_buf1[iphi1]) / 2.
    phi[1:nphi-2] = (y_buf1[iphi1:iphi2 - 1] + y_buf1[iphi1 + 1:iphi2]) / 2.
    phi[nphi-1]   = y_buf1[iphi2] + (y_buf1[iphi2] - y_buf1[iphi2 - 1])/2.
    w_buf = w_buf1[iEarr,iphi1:iphi2]
    e_buf = e_buf1[iEarr,iphi1:iphi2]
    y_buf = y_buf1[iphi1:iphi2]
  ENDIF ELSE IF (inst EQ 'IN4') AND (NOT no_small)  THEN BEGIN
    nphi          = ny+2
    phi           = FLTARR(nphi)
    IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: Phi_in=',y_in*180./!pi
    i1=(WHERE(y_in GT 10.*!PI/180.))[0]-1
    phi[0]    = y_in[0] - (y_in[1] - y_in[0]) / 2.
    phi[1:i1] = (y_in[0:i1-1]+y_in[1:i1])/2.
    phi[i1+1] = phi[i1]+(y_in[i1]-y_in[i1-1])
    phi[i1+3:nphi-2] = (y_in[i1+1:ny-2]+y_in[i1+2:ny-1])/2.
    phi[i1+2] = phi[i1+3]-(y_in[i1+2]-y_in[i1+1])
    phi[nphi-1] = phi[nphi-2]+(y_in[ny-1]-y_in[ny-2])
    w_buf = [[w_in[iEarr,0:i1]],[0.*iEarr],[w_in[iEarr,i1+1:ny-1]]]
    e_buf = [[e_in[iEarr,0:i1]],[0.*iEarr-1.],[e_in[iEarr,i1+1:ny-1]]]
    y_buf = (phi[0:nphi-2]+phi[1:nphi-1])/2.
    IF keyword_set(verbose) THEN PRINT,'SQW_REBIN:  phi =',phi*180./!pi
    IF keyword_set(verbose) THEN FOR i=0,50 DO PRINT,'SQW_REBIN: ', i, phi[i]*180./!pi, y_in[i]*180./!pi
  ENDIF ELSE BEGIN
    nphi          = iphi2 - iphi1 + 2
    phi           = FLTARR(nphi)
    phi[0]        = y_in[iphi1] - (y_in[iphi1 + 1] - y_in[iphi1]) / 2.
    phi[1:nphi-2] = (y_in[iphi1:iphi2 - 1] + y_in[iphi1 + 1:iphi2]) / 2.
    phi[nphi-1]   = y_in[iphi2] + (y_in[iphi2] - y_in[iphi2 - 1])/2.
    w_buf         = w_in[iEarr,iphi1:iphi2]
    e_buf         = e_in[iEarr,iphi1:iphi2]
    y_buf         = y_in[iphi1:iphi2]
  ENDELSE
  COSphi = COS(phi)

  ; ----------------------------------------------------------------------------------
  ;    reverse array direction for negative angles
  ; ----------------------------------------------------------------------------------
  IF inst EQ 'D7' THEN BEGIN
    IF phi[0] LT 0. THEN BEGIN
      w_buf =     REVERSE(w_buf,2)
      e_buf =     REVERSE(e_buf,2)
      y_buf = ABS(REVERSE(y_buf))
      phi   = ABS(REVERSE(phi))
      COSphi=     REVERSE(COSphi)
    ENDIF
  ENDIF

  ;;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: phi=',phi*180./!pi
  IF keyword_set(verbose) THEN $
    PRINT, format = '("SQW_REBIN: phi_out  = [",F7.3,",",F7.3,"]")',min(phi)*180./!pi,max(phi)*180./!pi

  ; -------------------------------------------------------------------------------------
  ;   Rebin angles to constant Q grid
  ; -------------------------------------------------------------------------------------
  a       = const2			; E(meV)=a*Q(A**-1)**2
  oldymin = 0.

  FOR iQ = 0,nQ - 1 DO BEGIN
    Qmin = Q[iQ] - (dQ / 2.)
    Qmax = Q[iQ] + (dQ / 2.)
    Q00  = [Qmin, Qmin, Qmax, Qmax]
    FOR iEps = 0, nEps - 2 DO BEGIN
      Emin     = Eps[iEps]
      Emax     = Eps[iEps + 1]
      corrarea = dQ*(Emax - Emin)
      Eps0     = [Emin, Emax, Emax, Emin]
      COSphi0  = (2.*Ei - Eps0 - a*Q00^2) / (2.*SQRT(Ei*(Ei - Eps0)))

      IF MAX(ABS(COSphi0)) GE 1. THEN GOTO, outside
      phi0     = ACOS(COSphi0)
      phimin   = MIN(phi0)
      phimax   = MAX(phi0)
      IF (phimax LT phi[0]) OR (phimin GT phi[nphi-1]) THEN GOTO, outside
      iphi     = WHERE(phi GT phimin AND phi LT phimax, nlines)
      iphi0    = (iphi[0] - 1) > 0
      IF nlines EQ 0 THEN BEGIN
        phimean = (phimin + phimax) / 2.
        ip      = WHERE(phi LT phimean, np)
        iphi0   = ip[np - 1]
      ENDIF
      startrebin:
      Areasum    = 0.
      wsum       = 0.
      e2sum      = 0.
      phiminmeas = 7.
      phimaxmeas = 0.

      FOR iphi = iphi0,(iphi0 + nlines) < (nphi - 2) DO BEGIN
        COSphi1 = COSphi[iphi]
        COSphi2 = COSphi[iphi + 1]
        COSphi0 = [COSphi1, COSphi1, COSphi2, COSphi2]
        Q0      = SQRT((2.*Ei - Eps0 - 2.*SQRT(Ei*(Ei - Eps0))*COSphi0)/a)
        area    = OVERLAP(Q0, Eps0, iprint, oldymin)
        IF area GT 0. THEN BEGIN
          w = w_buf[iEps,iphi]
          e = e_buf[iEps,iphi]
          IF (w NE 0.) OR (e GE 0.) THEN BEGIN
            areasum    = areasum + area
            wsum       = wsum + area*w
            e2sum      = e2sum + (area*e)^2
            phiminmeas = phiminmeas < phi[iphi]
            phimaxmeas = phimaxmeas > phi[iphi + 1]
          ENDIF
307
        ENDIF
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
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
438
439
440
441
442
443
444
445
446
447
      ENDFOR
      IF areasum NE 0. THEN BEGIN
        w_out[iQ,iEps] = wsum / areasum
        e_out[iQ,iEps] = SQRT(e2sum) / areasum
        GOTO, binned
      ENDIF
      outside:	w_out[iQ,iEps] =  0.
      e_out[iQ,iEps] = -1.
      GOTO, nextpoint
      binned:
      p1 = phimin > phiminmeas
      p2 = phimax < phimaxmeas
      IF p2 - p1 LT (phimax - phimin)/2. THEN BEGIN
        w_out[iQ,iEps] =  0.
        e_out[iQ,iEps] = -1.
      ENDIF
      nextpoint:
    ENDFOR
  ENDFOR


  ;;IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: End of rebinning'

  IF twice THEN BEGIN
    IF (iphi1 EQ 0) THEN BEGIN
      w_out1 = w_out
      e_out1 = e_out
      iphi1  = iphi1next
      iphi2  = iphi2next
      GOTO, start
    ENDIF ELSE BEGIN
      ;   take weighted averages of negative and positive banks for D7

      w_out2 = w_out
      e_out2 = e_out
      w_out  = 0.
      e_out  = 0.
      not1   = WHERE(e_out1 LE 0.,n1)
      IF n1 NE 0 THEN e_out1[not1] = 1.
      not2   = WHERE(e_out2 LE 0.,n2)
      IF n2 NE 0 THEN e_out2[not2] = 1.
      w_out = (w_out1/e_out1^2 + w_out2/e_out2^2) / (1./e_out1^2 + 1./e_out2^2)
      e_out = 1./SQRT(1./e_out1^2 + 1./e_out2^2)
      IF n1 NE 0 THEN e_out1[not1] = -1.
      IF n2 NE 0 THEN e_out2[not2] = -1.
      IF n1 NE 0 THEN BEGIN
        w_out[not1] = w_out2[not1]
        e_out[not1] = e_out2[not1]
      ENDIF
      IF n2 NE 0 THEN BEGIN
        w_out[not2] = w_out1[not2]
        e_out[not2] = e_out1[not2]
      ENDIF
    ENDELSE
  ENDIF

  ;;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: End of rebinning section'

  ; -------------------------------------------------------------------------------------
  ;   Chop off superfluous bits
  ; -------------------------------------------------------------------------------------
  i = WHERE(e_out GT -1.,n)
  IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: ', strtrim(n,2), '   non-zeroed points'

  checkQ2:
  iw0 = WHERE(w_out[nQ - 1,*] EQ  0., nw0)
  ie0 = WHERE(e_out[nQ - 1,*] EQ -1., ne0)

  IF (nw0 EQ nEps) AND (ne0 EQ nEps) THEN BEGIN
    nQ    = nQ - 1
    w_out = w_out[0:nQ - 1,*]
    e_out = e_out[0:nQ - 1,*]
    Q     = Q[0:nQ - 1]
    GOTO, checkQ2
  ENDIF

  checkEps1:
  iw0 = WHERE(w_out[*,0] EQ  0., nw0)
  ie0 = WHERE(e_out[*,0] EQ -1., ne0)
  IF (nw0 EQ nQ) AND (ne0 EQ nQ) THEN BEGIN
    nEps  = nEps - 1
    w_out = w_out[*,1:nEps]
    e_out = e_out[*,1:nEps]
    Eps   = Eps[1:nEps]
    GOTO, checkEps1
  ENDIF

  checkEps2:
  iw0 = WHERE(w_out[*,nEps - 1] EQ  0., nw0)
  ie0 = WHERE(e_out[*,nEps - 1] EQ -1., ne0)
  IF (nw0 EQ nQ) AND (ne0 EQ nQ) THEN BEGIN
    nEps  = nEps - 1
    w_out = w_out[*,0:nEps - 1]
    e_out = e_out[*,0:nEps - 1]
    Eps   = Eps[0:nEps - 1]
    GOTO, checkEps2
  ENDIF

  ;;	IF keyword_set(verbose) THEN PRINT,'SQW_REBIN: End of chopping section'

  ;-------------------------------------------------------------------------------------
  ;Return parameters and exit


  IF NOT swap_QE THEN BEGIN
    datp.y_tit = datp.x_tit
    datp.x_tit = 'Wavevector Transfer (A!U-1!N)'
    MOD_DATP, datp, "x", Q
    MOD_DATP, datp, "y", x_in[iEarr]
  ENDIF ELSE BEGIN
    w_out      = TRANSPOSE(w_out)
    e_out      = TRANSPOSE(e_out)
    datp.y_tit = 'Wavevector Transfer (A!U-1!N)'
    MOD_DATP, datp, "x", x_in[iEarr]
    MOD_DATP, datp, "y", Q
  ENDELSE

  MOD_DATP, datp, "e", e_out

  PRINT, FORMAT = '("SQW_REBIN: Rebinned to constant Q-w: dQ = ",F4.2," A-1")', dQ

  s = STRTRIM(STRING(FLOAT(dQ)),2)	& dQ = string_round(s)
  s = ' -sr(dQ=' + dQ
  IF Emin0 GT -9.E+09 THEN BEGIN
    Emin = STRTRIM(STRING(FIX(Emin0)),2)
    s = s + ',emin='+Emin
  ENDIF
  CASE ibank OF
    0: bs = '/neg)'
    1: bs = '/pos)'
    2: bs = '/all)'
  ENDCASE
  IF inst EQ 'D7' THEN s = s + bs ELSE s = s + ')'
  datp.other_tit = datp.other_tit + s

  GIVE_DATP, datp

  finished:
  RETURN, w_out
END
448
449