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