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