in5_pab.pro 3.56 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
FUNCTION in5_pab,w_in,emin=emin,emax=emax, dw=dw, verbose = verbose
  ;** *********************************************************
  ;
  ; Examples:
  ;
  ;  w2 = in5_pab(w1, /verbose)
  ;
  ; w2 =  in5_pab(w1, emin = -150, emax = 0.0, /verbose)
  ;
  ;
  ; Calculates (P(a,b))
  ;
  ;    P(a,b)=S(2theta,w)*Q^2/w/n(w,T)
  ;
  ;  w1 should be a horizontal S(Q,w) (after sqw_rebin and transpose)
  ;
  ;
  ;----------------------------------------------------
  ; empty cell 300K-500 5.1A
  ;----------------------------------------------------
  ;  a = '180632>180662'
  ;  w1 = rdopr(a) & w1 = remove_spectra(w1,M) & w1 = normalise(w1) & w1 = sumbank(w1) & w1 = in5_vnorm(w1, w60, /verbose)
  ;  w2 = in5_t2e(w1, w60, /verbose)
  ;  w3 = in5_pab(w2, /verbose) ; axes in E, Phi for sqw_rebin 
  ;  w4 = in5_sqw_rebin(w3, emin=-100., dq=0.02, /verbose)
  ;  w5 = transpose(w4)
  ;
  ;
  ;
  ;
  ; New: JO, Fri Jul  7 15:47:16 2017 Pab for nuclear data
  ;
  ;
  ;** *********************************************************

  COMMON c_lamp_access, inst
  COMMON printing, iprint, outstring

  take_datp, datp
  par  = datp.p
  x_in = datp.x
  y_in = datp.y
  e_in = datp.e

  lambda = par(21)
  temp   = par(11)
  ei = 81.799/(lambda^2)+x_in*0.
  sz = size(w_in)

  IF (temp EQ 0.0) THEN BEGIN
    print, 'WARNING: Temperature not found, set to 300 K'
    temp=300.0
  ENDIF
  IF NOT keyword_set(emax) THEN emax = MAX(x_in)
  IF NOT keyword_set(emin) THEN emin = MIN(x_in)
  IF NOT keyword_set(dw)   THEN dw = 0.0

  IF keyword_set(verbose) THEN BEGIN
    print, 'IN5_PAB: T      =',strtrim(string(temp),2),' K'
    print, 'IN5_PAB: lambda =',strtrim(string(lambda),2),' A'
    print, 'IN5_PAB: emin=',strtrim(string(emin),2),' meV, emax=',strtrim(string(emax),2),' meV'
  ENDIF

  points=where(x_in GE emin AND x_in LE emax)

  w_buf=w_in & e_buf=e_in & y_out=y_in

  ; -------------------------
  ; Compute Alpha and Beta
  ; -------------------------
  ;
  ; a = hbar^2 Q^2/2mkT and DW factor
  q2=w_in*0.0 & deb_wal=q2
  x_in=FLOAT(x_in)
  FOR i=0,n_elements(y_in)-1 DO q2[*,i]=2*ei[*]-x_in[*]-2*sqrt(ei[*]*ABS(ei[*]-x_in[*]))*cos(y_in[i]*!PI/180)
  q2=q2/2.072
  deb_wal=exp(-dw*q2)
  a = q2*deb_wal

  ; b = hbar omega/kT
  ;  b = x_in/(1-exp(-1.*x_in*11.6045/temp))
  b = x_in*11.6045/temp




  indnul=WHERE(ABS(b) LE 1.e-12)
  IF n_elements(indnul) GT 1 THEN BEGIN
    w_buf(indnul,0:n_elements(y_in)-1)=0.
    b(indnul)=1.
  END

  ; -------------------------------------
  ; Compute S(a,b)
  ; -------------------------------------
  w_buf = w_buf/a
  e_buf = e_buf/a
  ;  FOR i=0,n_elements(x_in)-1 DO BEGIN
  ;    w_buf[i,*] = w_buf[i,*]/a
  ;    e_buf[i,*] = e_buf[i,*]/a
  ;ENDFOR

  ; simplification: 2*b*sinh(b/2)*exp(-b/2) = x (1 - Cosh[x] + Sinh[x])

  FOR i=0,n_elements(y_in)-1 DO BEGIN
;    w_buf[*,i] = 2*b*sinh(b/2)*w_buf[*,i] *exp(-b/2)
;    e_buf[*,i] = 2*b*sinh(b/2)*e_buf[*,i] *exp(-b/2)
    w_buf[*,i] = 2*b*sinh(b/2)*w_buf[*,i] *exp(-b/2)
    e_buf[*,i] = 2*b*sinh(b/2)*e_buf[*,i] *exp(-b/2)
;    w_buf[*,i] = b*(exp(b)-1)*w_buf[*,i]
;    e_buf[*,i] = b*(exp(b)-1)*e_buf[*,i]
  ENDFOR



  ;************************
  ;output
  ;************************
  ;  y_out = a[points]   ; alpha is x Q^2
  ;  x_out = b[points]   ; beta is x hw

  ; in Phi and energy as for Sqw_rebin:
  y_out = y_in           ; in 2theta
  x_out = x_in[points]   ;
  e_out = e_buf[points,*]
  w_out = w_buf[points,*]


  output:
  mod_datp, datp, "e", e_out
  mod_datp, datp, "x", x_out
  mod_datp, datp, "y", y_out
  mod_datp, datp, "x_tit", "Energy [meV]"
  mod_datp, datp, "y_tit", "$2\theta$]"
  give_datp, datp

  return, w_out

END