FUNCTION straight_2d, win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, $
                           texture=txtr, integrate=integ, Weights=Iweight, $
                           verbose=verbose
;+
;******* ***********
;**
;** User callable
;** The call is w6 = straight_2d(w1 [,distance=samp2detect_Cm] [,/texture])
;**
;** Win  is a 2D (3D) diffractogram, i.e. I(2theta,h)
;**      where h is the physical height of the multidetector.
;**
;** Wout is a 2D diffractogram, i.e. I(2theta,h)
;** in which the Debye-Scherrer cones have been straightened out.
;**
;** For each detector height, it calculates a vector of 2theta
;** angles and interpolates onto equatorial plane.
;**
;** !!! WORKS BETTER WITHOUT INTERPOL (& with regular angles).
;** !!! Do not use the same input/output workspace.
;**
;** Set the right values for:
;**    h = physical height of the multidetector --> in Y
;**    r = physical radius of the multidetector --> in call
;**
;** Oct 2008: Care with 2thetha non regular if inter is not set,
;**           deltaX is added to better control the regular function.
;**
;** Dec 2008: Correct vertical position for textured diffractions (use keyword /TEXTURE)
;**
;** Dec 2009: !!! Replace -1 (no data) by -1E-9
;** Jan 2009: (texture) Replace arrows*fY/fX  by  arrows*fX/fY
;**
;** Nov 2011: Join straight_2d and straight_1d, adding an integrate keyword
;**           Compute a weighting vector for integration (returned in keyword)
;**
;** Nov 2013: Loop in C code implementation
;-

FORWARD_FUNCTION interpol, regular
st1 = systime(1)

lamp_loop_so, SFdll ;We may have loops in C code

if (!version.release lt '5.3') then inter=1 $ ;***unknown VALUE_LOCATE function***
else if keyword_set(txtr)      then inter=0

; Get variables associated with win
pasp = 1
if n_elements(datp) eq 0 then take_datp, datp else pasp = 0

; Detector distance
if n_elements(e_d) eq 1 then r = e_d else begin
   inst = ''
   tags = tag_names(datp)
   idx = where(tags eq 'OTHER_TIT', found)
   if found then inst = strlowcase((strsplit(datp.other_tit, /extract))[0])
   case inst of  
      'd2b':   r = 130.0
      'salsa': r = 100.0
      'd19':   r = 76.3
      'in5':   r = 400.0
      else: begin
               lamp_journal, '------------------------------------------------------------------'
               lamp_journal, 'ERROR !!!'
               lamp_journal, 'Instrument ' + inst + ' is not in case list.'
               lamp_journal, 'Sample to detector distance not known!'
               lamp_journal, 'Use distance keyword in call to overcome this issue.'
               lamp_journal, 'Returning empty workspace!'
               lamp_journal, '------------------------------------------------------------------'
               return, 0
            end
   endcase 
endelse
lamp_journal, 'Starting straight_2d with r = ' + string(r)

if (size(win))[0] eq 3 then kz = (size(win))[3]-1 else kz = 0


; Deal with the Hexapod case in Salsa, where the point information is stored as a string chain
tags = tag_names(datp)
idx = where(tags eq 'Z', found)
if found then if typename(datp.z) eq 'STRING' then mod_datp, datp, 'z', fix(indgen(n_elements(datp.z)) + 1) 

; Function works only with regular 2thetha !!!
regul = "" 
regu = 0
if not keyword_set(inter) then begin
   Sxyz = regular(win, dY=0, dX=deltaX, datp=datp, /GetNewSize) 
   regul = " regular"
   if Sxyz[0] ne (size(win))[1] then begin
      win = regular(win, dY=0, dX=deltaX, datp=datp) 
      regul = " non-regular"
      datp.x_tit = datp.x_tit + ' (regular)' 
      regu=1
      if (pasp) then give_datp, datp, /second  
   endif  
endif

; If x is bidim then take the middle. Not the best correction, to be improved!!!
if (size(datp.x))[0] eq 2 then XX = datp.x[*,(((size(datp.x))[2])-1)/2] else XX = datp.X

; Don't use negative angles !!!
wzero = where(XX ge 0.)
w0 = wzero[0]
nx = n_elements(XX)

rcos2t  = r*cos(XX[wzero]*!pi/180.0)
rsin2t2 = (r*sin(XX[wzero]*!pi/180.0))^2
y2      = datp.y^2

; 2D array, number of pixels vertically defined here (21)
h_pixels     = (size(datp.y))[1]
glob_low_idx = fltarr(h_pixels)
glob_hi_idx  = glob_low_idx

; Define output array and initialise it
Wout    = win*0.
Eout    = 0
win_s1  = (size(win))[1] - 1
rc_s1   = n_elements(rcos2t) - 1
Iweight = XX*0.                   ;Keep weights for integration
vv      = fltarr(kz+1) & ee=vv    ;temporaries

; Find how many times data have been accumulated in input workspace
n_acc_sc = 1
Eok = 0
if n_elements(datp.e) eq n_elements(win) then begin
  Eok = 1
  Eout = Wout
  idx = where(datp.e[*,h_pixels/2,0] gt 0.)
  if idx[0] ge 0 then n_acc_sc = ceil( total(win[idx,h_pixels/2,0] / $
    datp.e[idx,h_pixels/2]^2,0) / $
    n_elements(idx) )
endif
if keyword_set(verbose) then print, 'Straightening ', h_pixels, ' , ', XX[0], ' ->', XX[win_s1], regul
if keyword_set(verbose) then wait, 0.01

; Calculation of the rays for each 2-theta (needed for textures)
; Need to use float(!const.rtod) to avoid error (would need to explore this issue)
if keyword_set(txtr) then begin
  zq     = sqrt(rsin2t2 + y2[0])                 ; height -h/2 to +h/2 (j=0)
  xh     = atan(zq, rcos2t) * float(!const.rtod) ; calculate 2theta values for the equatorial plane xh
  XHcm   = sin(xh/float(!const.rtod)) * r
  XXcm   = sin(XX[wzero]/float(!const.rtod)) * r
  fY     = abs(datp.y[1] - datp.y[0])        ; delta Y
  fX     = abs((XX[nx-1]-XX[0]) / nx)        ; delta X
  corde  = abs(2.0*datp.y[0])                ; corde of the arcs
  arrows = (abs((XHcm-XXcm)*fX/fY))>(1e-8)   ; arrows of the arcs normalized with XYfct-> *fX/fY
  rays   = (4.*arrows^2+corde^2)/(8.*arrows) ; rays of the arcs
  ;arcs   = rays*2.*asin(corde/(2.*rays))     ; length of the arcs unused
endif
jof = rcos2t * 0.

; For each detector height, calculate the 2-theta vector
for j = 0L, h_pixels-1 do begin

   if keyword_set(txtr) then begin
      jof[*] = ((rays*asin(datp.y[j]/rays))/fY +(h_pixels-1)/2.)>0.<(h_pixels-1.) ;New J index (for texture)
      jof1   = floor(jof) 
      jof2   = ceil(jof) 
      ;vertical weighting
      PAj2   = jof-jof1 
      PAj1   = 1.0 - PAj2       
   endif else jof[*] = j

   zq = sqrt(rsin2t2 +  y2[j])              ; height -h/2 to +h/2
   xh = atan(zq,rcos2t)* float(!const.rtod) ; calculate 2theta values for the equatorial plane xh

   if xh[0] lt xh[rc_s1] then xdx = where((XX ge xh[0]) and (XX le xh[rc_s1])) $ ;find corresponding indices xh->XX
                         else xdx = where((XX le xh[0]) and (XX ge xh[rc_s1]))   ;(ascending or descending)

   if xdx[0] ge 0 then begin
    
      nxe = n_elements(xdx)-1
      xeq = XX[xdx]
      ;keep first and last 2theta index (straight_1d)
      glob_low_idx[j] = xdx[0] 
      glob_hi_idx[j]  = xdx[nxe] 

      ;INTERPOLATE ONTO EQUATORIAL PLANE
      if keyword_set(inter) then begin

         for KK = 0L,kz do begin
            v = win[wzero,j,KK]
            ;Will prevent from interpol bounds artifacts
            v[0,*] = v[1,*] 
            v[rc_s1,*] = v[rc_s1-1,*] 
            wout[xdx,j,KK] = INTERPOL(v,xh,xeq)>0.
            if Eok then Eout[xdx,j,KK] = INTERPOL(datp.e[wzero,j,KK],xh,xeq)>0.  
         endfor

      endif else begin

         sk = VALUE_LOCATE(xeq, xh)  
         ind = where((sk ge 0 ) and (sk le nxe-1))
         sk = sk[ind] 
         ik = sk+xdx[0] 
         ik1 = ik+1
         ;horizontal  weighting
         PA = ((xh[ind]-xeq[sk])/(xeq[sk+1]-xeq[sk]))>0.<1.  
         PA1 = 1.0 - PA      
         v = win[wzero,j,*] 
         if Eok then e = datp.e[wzero,j,*]

         if keyword_set(txtr) then begin
          
            ;horizontal & vertical weighting for texture
            PA11 = PA1 * PAj1[ind] 
            PA01 = PA * PAj1[ind] 
            PA12 = PA1 * PAj2[ind] 
            PA02=PA*PAj2[ind]   
            if SFdll gt ' '  then begin        ;Try the C code if available
               dims_wout = SIZE(wout, /DIMENSIONS)
               dims_v_e  = SIZE(v,    /DIMENSIONS) ; same size for e and v
               if n_elements(dims_wout) eq 2 then begin
                  dims_wout=[dims_wout,1]
                  dims_v_e =[dims_v_e ,1,1] 
               endif
               if not Eok then e = 0
               catch, stat 
               if stat ne 0 then bid=0 else begin
                  bid = CALL_EXTERNAL(SFdll, 'STR_textu', $
                         n_elements(sk), ind, Eok, jof1, jof2, ik, ik1, $
                         dims_wout[0], dims_wout[1], dims_wout[2], wout, Eout, $
                         dims_v_e[0], dims_v_e[1], dims_v_e[2], v, e, $
                         PA01, PA02, PA11, PA12, $
                         VALUE=[ 1,0,1,0,0,0,0,  1,1,1,0,0,  1,1,1,0,0,  0,0,0,0], $
                         /I_VALUE, /CDECL)
               endelse
               catch, /cancel
               if (bid ne 1) then begin
                  lamp_loop_so, /reset
                  SFdll = '' 
               endif
            endif else begin ;Use the IDL code
               for k = 0L, n_elements(sk)-1 do begin 
                  kk = ind[k]
                  jj1 = jof1[kk] 
                  jj2 = jof2[kk] 
                  vv[*] = v[kk,*,*] 
                  if Eok then ee[*] = e[kk,*,*]
                  wout[ik[ k],jj1,*] += (vv*PA11[k])
                  wout[ik1[k],jj1,*] += (vv*PA01[k])
                  wout[ik[ k],jj2,*] += (vv*PA12[k])
                  wout[ik1[k],jj2,*] += (vv*PA02[k])
                  if Eok then begin
                     Eout[ik[ k],jj1,*] += (ee*PA11[k]) ;Increase errors
                     Eout[ik1[k],jj1,*] += (ee*PA01[k])
                     Eout[ik[ k],jj2,*] += (ee*PA12[k])
                     Eout[ik1[k],jj2,*] += (ee*PA02[k])
                  endif
               endfor
            endelse

         endif else begin ;no texture case: horizontal weighting only

            if SFdll gt ' '  then begin ;Try the C code if there.
               dims_wout = SIZE(wout, /DIMENSIONS)
               dims_v_e  = SIZE(v,    /DIMENSIONS) ; same size for e and v
               if n_elements(dims_wout) eq 2 then begin
                  dims_wout = [dims_wout,1]
                  dims_v_e  = [dims_v_e ,1,1] 
               endif
               if not Eok then e = 0
               catch, stat 
               if stat ne 0 then bid=0 else begin
                  bid = CALL_EXTERNAL(SFdll, 'STR_proj', $
                         n_elements(sk), ind, ik, ik1, j, Eok, $
                         dims_wout[0], dims_wout[1], dims_wout[2], wout, Eout, Iweight, $
                         dims_v_e[0], dims_v_e[1], dims_v_e[2], v, e, PA, PA1, $
                         VALUE=[ 1,0,0,0,1,1,  1,1,1,0,0,0,  1,1,1,0,0,0,0], $
                         /I_VALUE, /CDECL)
               endelse
               if (bid ne 1) then begin
                  lamp_loop_so, /reset 
                  SFdll = '' 
               endif
            endif else begin ;Use the IDL code
               for k = 0L, n_elements(sk)-1 do begin 
                  kk = ind[k]
                  vv[*] = v[kk,*,*] 
                  if Eok then ee[*] = e[kk,*,*]
                  wout[ik[ k], j ,*] += (vv*PA1[k])
                  wout[ik1[k], j ,*] += (vv*PA[ k])
                  if Eok then begin
                     Eout[ik[ k], j ,*] += (ee*PA1[k])   ;Increase errors
                     Eout[ik1[k], j ,*] += (ee*PA[ k])
                  endif
                  Iweight[ik[ k]] +=  PA1[k]       ;Calculate weights for integration
                  Iweight[ik1[k]] +=  PA[ k]
               endfor
            endelse

         endelse ;end if keyword_set(txtr) else block
         
      endelse ;end if keyword_set(inter) else block
      
      if not keyword_set(integ) then begin
         if xdx[0]   gt 0      then wout[0:xdx[0]-1       ,[floor(jof[0])    ,ceil(jof[0])    ],*] = -1E-9 ;No data
         if xdx[nxe] lt win_s1 then wout[xdx[nxe]+1:win_s1,[floor(jof[rc_s1]),ceil(jof[rc_s1])],*] = -1E-9
      endif
      
   endif ;end if xdx[0]>0 block
   
endfor ;end loop over h_pixels

totIN = total(win[wzero,*,*])

if not keyword_set(integ) then begin
   totIDM = totIN / total(Wout) ;Keep same count
   Wout   = temporary(Wout) * totIDM
   Eout   = temporary(Eout) * totIDM
   if Eok then datp.e = Eout
   if pasp then if (regu or Eok) then give_datp, datp
   return, Wout
endif

v=0 & e=0 & vv=0 & ee=0 & sk=0 & ik=0 & ik1=0 & PA=0 & PA1=0 ;FREEs

if keyword_set(verbose) then print, 'Total time = ', round((systime(1)-st1)*100)/100., ' s'
if keyword_set(verbose) then wait, 0.01

;sum over the height to output 1d array
media   = median(Iweight)>1. 
Iweight = media/(Iweight>2.)

Wout_1d = reform(total(wout,2))     
WOUT=0 ;FREEs

for z=0,kz do Wout_1d[*,z] *= Iweight

totIDM  = totIN  / total(Wout_1d) ;Keep same count
Wout_1d = temporary(Wout_1d) * totIDM

if Eok then begin
   Eout_1d = reform(sqrt(total(Eout^2,2)))  
   EOUT=0 ;FREEs
   for z=0,kz do Eout_1d[*,z] *= (Iweight * totIDM)
endif

wout_1d = wout_1d[wzero,*]
if Eok then Eout_1d = Eout_1d[wzero,*] else Eout_1d = sqrt(wout_1d/n_acc_sc)
mod_datp, datp,'e', Eout_1d
mod_datp, datp,'x', XX[wzero]

szz = size(datp.z) 
if szz[0] eq 2 then mod_datp, datp, 'y', reform(datp.z,szz[1]*szz[2]) $
               else mod_datp, datp, 'y', datp.z
syy=size(datp.y) 
if syy[0] eq 1 then if datp.y[0] eq datp.y[syy[1]-1] then mod_datp, datp, 'y', findgen(syy[1])+1.
if (size(wout_1d))[0] eq 2 then  begin
   mod_datp, datp, 'y', datp.z 
   datp.y_tit = datp.z_tit 
endif
if (size(wout_1d))[0] eq 1 then datp.y_tit='Integrations' 

if pasp then give_datp, datp

return, wout_1d

end