Commit 1f3e5053 authored by Didier Richard's avatar Didier Richard
Browse files

Added IN5 macros. Plugin IN5 ok from today.

parent 30e869d9
......@@ -1122,12 +1122,12 @@ pro P_AFTER_REALIZE_DID ,sepben,sepdon,sepdid
i =0
if sepdon gt 0L then begin & widget_control,bad_id=i,sepdon ,get_value=j
if i eq 0 then begin & wset,j & erase,50 & did_etc.spDon=j
if sys_dep('VERSION') ge '6.4' then widget_control,sepdon,Tooltip='Physical Memory used'
if sys_dep('VERSION') ge '6.4' then widget_control,sepdon,Tooltip='Physical Memory: 1Gb/tick'
endif & endif
if sepdid gt 0L then begin & widget_control,bad_id=i,sepdid ,get_value=j
if i eq 0 then begin & wset,j & erase,50 & did_etc.spDid=j
if sys_dep('VERSION') ge '6.4' then widget_control,sepdid,Tooltip='Physical Memory used'
if sys_dep('VERSION') ge '6.4' then widget_control,sepdid,Tooltip='Physical Memory: 1Gb/tick'
endif & endif
if (did_wd gt 0) or (!D.name eq 'Z') then loadct , tcol
......@@ -6474,15 +6474,18 @@ IF n_tags(alter) gt 0 then begin
if (Other gt 0.) and (Other+IdlMem lt did_etc.spRamT) then begin ; We have potentialy good info AND run into physical memory
;.spDid is the top BarId , .spDon is the bottom BarId , .spSiz is bytarr(600,4)
szB = (size(did_etc.spSiz))[1] & Bcol =[ 27 , 118 , 63] ; total size in pixel, col= green grass
szO = round(szB * Other / TotMem -1) >0 & Ocol =[113 , 0 , 0] ; Other size in pixel, col= red dark
szI = round(szB * IdlMem / TotMem -1) >0 & Icol =[186 , 61 , 0] ; Idl size in pixel, col= orange
did_etc.spSiz[*,*,0] = Bcol[0] & did_etc.spSiz[*,*,1] = Bcol[1] & did_etc.spSiz[*,*,2] = Bcol[2]
did_etc.spSiz[0:szO,*,0]= Ocol[0] & did_etc.spSiz[0:szO,*,1]= Ocol[1] & did_etc.spSiz[0:szO,*,2]= Ocol[2]
szO = round(szB * Other / TotMem -1) >0 & Ocol =[113 , 0 , 0] ; Other size in pixel, col= dark red
szI = round(szB * IdlMem / TotMem -1) >0 & Icol =[186 , 61 , 0] ; Idl size in pixel, col= red
did_etc.spSiz[*,*,0] = Bcol[0] & did_etc.spSiz[*,*,1] = Bcol[1] & did_etc.spSiz[*,*,2] = Bcol[2] ;Total green
did_etc.spSiz[0:szO,*,0]= Ocol[0] & did_etc.spSiz[0:szO,*,1]= Ocol[1] & did_etc.spSiz[0:szO,*,2]= Ocol[2] ;Other dark red
I0 = szO+1
I1 =(szO+szI) > I0
I2 = I1+1 & J=[1,2]
did_etc.spSiz[I0:I1,*,0]= Icol[0] & did_etc.spSiz[I0:I1,*,1]= Icol[1] & did_etc.spSiz[I0:I1,*,2]= Icol[2]
did_etc.spSiz[I2 ,J,0]= Icol[0] & did_etc.spSiz[I2 ,J,1]= Icol[1] & did_etc.spSiz[I2 ,J,2]= Icol[2]
I1 =(szO+szI+1) < (szB-1)
did_etc.spSiz[I0:I1,*,0]= Icol[0] & did_etc.spSiz[I0:I1,*,1]= Icol[1] & did_etc.spSiz[I0:I1,*,2]= Icol[2] ;Lamp red
I2 = [I1-3,I1-2,I1-1,I1] & J=[0,3]
did_etc.spSiz[I2 ,J,0]= Bcol[0] & did_etc.spSiz[I2 ,J,1]= Bcol[1] & did_etc.spSiz[I2 ,J,2]= Bcol[2] ;Arrow
I2 = [I1-1,I1] & J=[2]
did_etc.spSiz[I2 ,J,0]= Bcol[0] & did_etc.spSiz[I2 ,J,1]= Bcol[1] & did_etc.spSiz[I2 ,J,2]= Bcol[2] ;Arrow
did_etc.spSiz[I2 ,1,0]= 255 ;Arrow end red
;1Gb/tick
nticks = TotMem / 1000. & btick = round( szb / nticks) > 1
iticks =((findgen(nticks)+1) * btick) & iticks = iticks[where(iticks lt szb) >0] & J=[0,1]
......@@ -7999,6 +8002,7 @@ pro W_store ,W=wii ,ALL=all
lamp_sys=W_lamp_sys()
if keyword_set(ALL) then FOR i=1,lamp_sys do A_store, W=(i) $
else A_store, W=wii
SnapShotWks_event,0, alter={doit:'Memory'}
end
pro W_restore ,W=wii ,ALL=all
......@@ -8010,7 +8014,7 @@ pro W_restore ,W=wii ,ALL=all
lamp_sys=W_lamp_sys()
if keyword_set(ALL) then FOR i=1,lamp_sys do A_restore, W=(i) $
else A_restore, W=wii
SnapShotWks_event,0, alter={doit:'Memory'}
if n_elements(wii) ne 1 then wii=1 & to_don_history, wii,-1,' ',/nojournal
end
......@@ -8041,6 +8045,7 @@ if keyword_set(ALL) then begin & FOR i=1,lamp_sys do A_clear, W=(i)
n=0 & o=0 & p=0 & q=0 & r=0 & s=0 & t=0 & u=0 & v=0 & w=0 & x=0 & y=0 & z=0
logo,1
endif else A_clear, W=wii
SnapShotWks_event,0, alter={doit:'Memory'}
end
pro A_exchange ,W=wii
......
......@@ -4,7 +4,7 @@ pro Show_Version, wid, version, bil
;** return lamp version , keep some pseudo colors for lamp windows at start
;-
@lamp.cbk
version='14 Aug 2015 for idl -> 8.5'
version='09 Oct 2015 for idl -> 8.5'
if wid gt 0 then begin
bil=widget_base( wid,/row)
......@@ -436,7 +436,7 @@ common for_users, a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z
lamp_ben[17]=lamp_mac2
;did
lamp_did =lonarr(10)
if lamp_siz ge 900 then $
if lamp_siz ge 800 then $
sepdid =widget_draw( lamp_tmp2,xsize=600,ysize=4) else sepdid=0
lamp_did[0]=widget_base( lamp_tmp2,/frame,resource_name='did')
......
FUNCTION straight_2d, win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, texture=txtr, integrate=integ, Weights=Iweight
FUNCTION straight_2d, win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, texture=txtr, integrate=integ, Weights=Iweight, verbose=verbose
;+
;******* ***********
;**
......@@ -89,9 +89,9 @@ 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
print, 'straightening ' ,h_pixels,' , ',XX[0],' ->',XX[win_s1],regul
;print,'number of scans=',strtrim(n_acc_sc,2),' ,straightening ' ,h_pixels,' , ',XX(0),' ->',XX(win_s1),regul
wait ,.01
if keyword_set(verbose) then print, 'straightening ' ,h_pixels,' , ',XX[0],' ->',XX[win_s1],regul
if keyword_set(verbose) then wait ,.01
;print,'number of scans=',strtrim(n_acc_sc,2),' ,straightening ' ,h_pixels,' , ',XX[0],' ->',XX[win_s1],regul
;Calculation of the rays for each 2theta (needed for textures)
if keyword_set(txtr) then begin
......@@ -241,7 +241,7 @@ if not keyword_set(integ) then begin
Eout =temporary(Eout) * totIDM
if Eok then datp.e=Eout
print, 'straightening done ',round((systime(1)-st1)*100)/100.,'s'
if keyword_set(verbose) then print, 'straightening done ',round((systime(1)-st1)*100)/100.,'s'
if pasp then if (regu or Eok) then give_datp,datp
return, Wout
......@@ -250,8 +250,8 @@ endif
v=0 & e=0 & vv=0 & ee=0 & sk=0 & ik=0 & ik1=0 & PA=0 & PA1=0 ;FREEs
print, 'integrating ',round(systime(1)-st1),'s'
wait ,.01
if keyword_set(verbose) then print, 'integrating ',round(systime(1)-st1),'s'
if keyword_set(verbose) then wait ,.01
;sum over the height to output 1d array
media = median(Iweight)>1. & Iweight=media/(Iweight>2.)
......
; Absolute cross section:
; ------------------------
; We were hoping to have the data in fully normalised units (i.e. mb sr^-1 meV^-1 f.u.^-1)
; f.u. = formula unit
; We usually don't obtain the data in absolute units on IN5. If one use the in5_vnorm and in5_t2e
; together (both have to be used in this case ) one can obtain the absolute cross-section by reference to the vanadium.
; One has then to modify the in5_vnorm in the template1.prox :
;
; Absolute cross section calculation:
;
; w3 = in5_vnorm(w1,w2,tv=0.5,Ns=7.22305,ts=0.5, /verbose)
;
; Absolute cross section:
; Ns is sample number density (in 10^22 at/cm^-3)
; tv and ts are total volume [cm^3] seen by the beam that might be
; approximated by the effective sample thickness [cm] if the sample are
; flat slab and/or the the beam cross section is not known.
;
;
; tv and ts are vanadium and sample volume or surface or thickness in the beam. One can obtain the volume from the
; mass but it is made as such because it depends also on the quantity of sample in the beam (and it is more easy to handle
; slab samples ).
;
;
; template1.prox:
; ------------------
; ---------------------------------------------------------------------
;
; @template1
;
; ---------------------------------------------------------------------
; Sample
; ---------------------------------------------------------------------
w1 = rdrun(c)
w1 = in5_remove_spectra(w1,badSpectra=s, /verbose)
; ---------------------------------------------------------------------
; vanadium
; ----------------------------------------------------------
w2 = rdrun(e)
w2 = in5_remove_spectra(w2,badSpectra=s, /verbose)
; ---------------------------------------------------------------------
; 3) Normalize by the elastic peak of the vanadium
; (vnorm correct the vana from DWF first.)
;
; LaFePO sample: I think we want in terms of 'formula unit' of LaFePO so it is
; Ns = 1.49 [10^22 f.u./cm^3]. (Ns = 5.98 for atoms/cm^3)
; ts = 1.77 [cm^3]
;
; LaZnPO sample:
; Ns = 1.37 [10^22 f.u./cm^3]. f.u. = formula unit
; (Ns = 5.98 for atoms/cm^3)
; ts = 1.37 [cm^3]
;
; vanadium: tv = 1.34 (Nv = 7.223 mbut integrated in the code)
; ---------------------------------------------------------------------
; LaZnPO
; w3 = in5_vnorm(w1,w2, tv=1.34 ,Ns=1.37,ts=1.37, /verbose)
; LaFePO
w3 = in5_vnorm(w1,w2, tv=1.34 ,Ns=1.49,ts=1.77, /verbose)
; ---------------------------------------------------------------------
; 4) Correction for detector wavelength-efficiency +
; remove flat bkgd.
; ---------------------------------------------------------------------
w4 = in5_corrtof(w3, /deteff, /psd, /verbose)
; ---------------------------------------------------------------------
; 5) Conversion to S(2*theta, omega)
; ---------------------------------------------------------------------
w5 = in5_t2e(w4,/verbose)
; ---------------------------------------------------------------------
; 7) Force constant energy bins for convolution in the fit
; ---------------------------------------------------------------------
w6 = in5_energy_rebin(w5,de=d,Emin=-1.0,/force, /verbose)
; ---------------------------------------------------------------------
; Store S(2*theta, omega) into dpe format
; ---------------------------------------------------------------------
in5_outputspe, w6, file='spe/'+strmid(c,4,strlen(c)), dZCells=241, /phx, /verbose
;+
; NAME:
; CMREPLICATE
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
; Replicates an array or scalar into a larger array, as REPLICATE does.
;
; CALLING SEQUENCE:
; ARRAY = CMREPLICATE(VALUE, DIMS)
;
; DESCRIPTION:
;
; The CMREPLICATE function constructs an array, which is filled with
; the specified VALUE template. CMREPLICATE is very similar to the
; built-in IDL function REPLICATE. However there are two
; differences:
;
; * the VALUE can be either scalar or an ARRAY.
;
; * the dimensions are specified as a single vector rather than
; individual function arguments.
;
; For example, if VALUE is a 2x2 array, and DIMS is [3,4], then the
; resulting array will be 2x2x3x4.
;
; INPUTS:
;
; VALUE - a scalar or array template of any type, to be replicated.
; NOTE: These two calls do not produce the same result:
; ARRAY = CMREPLICATE( 1, DIMS)
; ARRAY = CMREPLICATE([1], DIMS)
; In the first case the output dimensions will be DIMS and
; in the second case the output dimensions will be 1xDIMS
; (except for structures). That is, a vector of length 1 is
; considered to be different from a scalar.
;
; DIMS - Dimensions of output array (which are combined with the
; dimensions of the input VALUE template). If DIMS is not
; specified then VALUE is returned unchanged.
;
; RETURNS:
; The resulting replicated array.
;
; EXAMPLE:
; x = [0,1,2]
; help, cmreplicate(x, [2,2])
; <Expression> INT = Array[3, 2, 2]
; Explanation: The 3-vector x is replicated 2x2 times.
;
; x = 5L
; help, cmreplicate(x, [2,2])
; <Expression> LONG = Array[2, 2]
; Explanation: The scalar x is replicated 2x2 times.
;
; SEE ALSO:
;
; REPLICATE
;
; MODIFICATION HISTORY:
; Written, CM, 11 Feb 2000
; Fixed case when ARRAY is array of structs, CM, 23 Feb 2000
; Apparently IDL 5.3 can't return from execute(). Fixed, CM, 24 Feb
; 2000
; Corrected small typos in documentation, CM, 22 Jun 2000
; Removed EXECUTE() call by using feature of REBIN() new in IDL 5.6,
; (thanks to Dick Jackson) CM, 24 Apr 2009
; Remove some legacy code no longer needed after above change
; (RETVAL variable no longer defined; thanks to A. van Engelen),
; CM, 08 Jul 2009
;
;-
; Copyright (C) 2000, 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function cmreplicate, array, dims
if n_params() EQ 0 then begin
message, 'RARRAY = CMREPLICATE(ARRAY, DIMS)', /info
return, 0L
endif
on_error, 2
if n_elements(dims) EQ 0 then return, array
if n_elements(array) EQ 0 then $
message, 'ERROR: ARRAY must have at least one element'
;; Construct new dimensions, being careful about scalars
sz = size(array)
type = sz[sz[0]+1]
if sz[0] EQ 0 then return, make_array(value=array, dimension=dims)
onedims = [sz[1:sz[0]], dims*0+1] ;; For REFORM, to extend # of dims.
newdims = [sz[1:sz[0]], dims ] ;; For REBIN, to enlarge # of dims.
nnewdims = n_elements(newdims)
if nnewdims GT 8 then $
message, 'ERROR: resulting array would have too many dimensions.'
if type NE 7 AND type NE 8 AND type NE 10 AND type NE 11 then begin
;; Handle numeric types
;; Passing REBIN(...,ARRAY) is a feature introduced in
;; IDL 5.6, and will be incompatible with previous versions.
array1 = rebin(reform([array], onedims), newdims)
return, array1
endif else begin
;; Handle strings, structures, pointers, and objects separately
;; Handle structures, which are never scalars
if type EQ 8 AND sz[0] EQ 1 AND n_elements(array) EQ 1 then $
return, reform(make_array(value=array, dimension=dims), dims, /over)
nold = n_elements(array)
nadd = 1L
for i = 0L, n_elements(dims)-1 do nadd = nadd * round(dims[i])
if type EQ 8 then $
array1 = make_array(value=array[0], dimension=[nold,nadd]) $
else $
array1 = make_array(type=type, dimension=[nold,nadd], /nozero)
array1 = reform(array1, [nold,nadd], /overwrite)
array2 = reform([array], n_elements(array))
;; Efficient copying, done by powers of two
array1[0,0] = temporary(array2)
stride = 1L ;; stride increase by a factor of two in each iteration
i = 1L & nleft = nadd - 1
while nleft GT stride do begin
array1[0,i] = array1[*,0:stride-1] ;; Note sneaky IDL optimization
i = i + stride & nleft = nleft - stride
stride = stride * 2
endwhile
if nleft GT 0 then array1[0,i] = array1[*,0:nleft-1]
return, reform(array1, newdims, /overwrite)
endelse
end
function convolve, image, psf, FT_PSF=psf_FT, FT_IMAGE=imFT, NO_FT=noft, $
CORRELATE=correlate, AUTO_CORRELATION=auto
;+
; NAME:
; CONVOLVE
; PURPOSE:
; Convolution of an image with a Point Spread Function (PSF)
; CATEGORY:
; Image Processing, Math
; EXPLANATION:
; The default is to compute the convolution using a product of
; Fourier transforms (for speed).
;
; CALLING SEQUENCE:
;
; psf = psf_Gaussian( NPIXEL=21, FWHM=5, /NORMALIZE )
; imconv = convolve( image1, psf, FT_PSF = psf_FT )
; or:
; correl = convolve( image1, image2, /CORREL )
; or:
; correl = convolve( image, /AUTO )
;
; INPUTS:
; image = 2-D array (matrix) to be convolved with psf
; psf = the Point Spread Function, (size < or = to size of image).
;
; OPTIONAL INPUT KEYWORDS:
;
; FT_PSF = passes out/in the Fourier transform of the PSF,
; (so that it can be re-used the next time function is called).
; FT_IMAGE = passes out/in the Fourier transform of image.
;
; /CORRELATE uses the conjugate of the Fourier transform of PSF,
; to compute the cross-correlation of image and PSF,
; (equivalent to IDL function convol() with NO rotation of PSF)
;
; /AUTO_CORR computes the auto-correlation function of image using FFT.
;
; /NO_FT overrides the use of FFT, using IDL function convol() instead.
; (then PSF is rotated by 180 degrees to give same result)
; METHOD:
; When using FFT, PSF is centered & expanded to size of image.
; HISTORY:
; written, Frank Varosi, NASA/GSFC 1992.
; Converted to IDL V5.0 W. Landsman September 1997
;-
sp = size( psf_FT ) & sif = size( imFT )
sim = size( image ) & sc = sim/2 & npix = N_elements( image )
if (sim[0] NE 2) OR keyword_set( noft ) then begin
if keyword_set( auto ) then begin
message,"auto-correlation only for images with FFT",/INF
return, image
endif else if keyword_set( correlate ) then $
return, convol( image, psf ) $
else return, convol( image, rotate( psf, 2 ) )
endif
if (sif[0] NE 2) OR (sif[sif[0]+1] NE 6) OR $
(sif[1] NE sim[1]) OR (sif[2] NE sim[2]) then imFT = FFT( image,-1 )
if keyword_set( auto ) then $
return, shift( npix*float( FFT( imFT*conj( imFT ),1 ) ), sc[1],sc[2] )
if (sp[0] NE 2) OR (sp[sp[0]+1] NE 6) OR $
(sp[1] NE sim[1]) OR (sp[2] NE sim[2]) then begin
sp = size( psf )
if (sp[0] NE 2) then begin
message,"must supply PSF matrix (2nd arg.)",/INFO
return, image
endif
Loc = ( sc - sp/2 ) > 0 ;center PSF in new array,
s = (sp/2 - sc) > 0 ;handle all cases: smaller or bigger
L = (s + sim-1) < (sp-1)
psf_FT = complexarr( sim[1], sim[2] )
psf_FT[ Loc[1], Loc[2] ] = psf[ s[1]:L[1], s[2]:L[2] ]
psf_FT = FFT( psf_FT, -1, /OVERWRITE )
endif
if keyword_set( correlate ) then $
conv = npix * float( FFT( imFT * conj( psf_FT ), 1 ) ) $
else conv = npix * float( FFT( imFT * psf_FT, 1 ) )
sc = sc + (sim MOD 2) ;shift correction for odd size images.
return, shift( conv, sc[1], sc[2] )
end
FUNCTION debyescherrer, win, interpol=inter, distance=e_d, datp=datp, $
verbose=verbose
; ----------------------------------------------------------------------------------
;
;+
;
; NAME:
;
;
; debyescherrer
;
;
; USAGE:
;
; w2 = debyescherrer(w1, distance=400.0)
;
; ARGUMENTS:
;
; distance: detector distance in [cm]
; /interpol: set keyword if interpolation needed
; /verbose: set keyword if verbosity required
;
;
;
; win is a 2D 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 DS cones should be straightened out
; For each detector height calculate a vector of 2theta
; angles and interpolate onto equatorial plane
; using INTERPOL
;
; h is the physical height of the multidetector
; r is the physical radius of the multidetector
; set them
;
;
; VERSION HISTORY:
;
; Johnson MR. 2005?
; DRD, 2006 improve speed
; DRD, 2008/09 improve for IN5 data
; JOR, 2008/09 try to improve (?) change the name, improve speed
; by changing the last loop by a *
;-
;
; ----------------------------------------------------------------------------------
;
;
; The call is w6 = debyescherrer(w1, distance=400.0, /verbose)
;
; 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 DS cones should be straightened out
; For each detector height calculate a vector of 2theta
; angles and interpolate onto equatorial plane
;
; h is the physical height of the multidetector
; r is the physical radius of the multidetector
; set them
; ----------------------------------------------------------------------------------
FORWARD_FUNCTION interpol
st1=systime(1)
if (!version.release lt '5.3') then inter=1 ;***VALUE_LOCATE function***
; ----------------------------------------------------------------------------------
; h= 30cm, defined in reading D2B and carried forward in Y
; h= 8cm, defined in reading SALSA and carried forward in Y
; h= 40cm, defined in reading D19 and carried forward in Y
; h=300cm, defined in reading IN5 and carried forward in Y
; ----------------------------------------------------------------------------------
r =130.0 ;D2B
if n_elements(e_d) eq 1 then $
r =e_d $;for SALSA 100.0 , D19 76.3 , IN5 400.
else $
print,'DEBYESCHERRER: Warning: distance not given -> set to ',strtrim(r,2),' cm!'
if (size(win))[0] eq 3 then kz=(size(win))[3]-1 else kz=0
; ----------------------------------------------------------------------------------
; get variables associated with win
; ----------------------------------------------------------------------------------
pasp=1 & if n_elements(datp) eq 0 then take_datp,datp else pasp=0
; ----------------------------------------------------------------------------------
; 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.) & indPhiEqZero = wzero[0]
; ----------------------------------------------------------------------------------
;wzero=lindgen((size(datp.x))(1))
; ----------------------------------------------------------------------------------
rcos2t = r*cos(XX[wzero]*!pi/180.0)
rsin2t2 = (r*sin(XX[wzero]*!pi/180.0))^2
RAD2DEG = 180.0/!pi
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
; ----------------------------------------------------------------------------------
wout = win*0.
win_s1 =(size(win))[1] -1
rc_s1 = n_elements(rcos2t)-1
; ----------------------------------------------------------------------------------
; Find how many times data has been accumulated in input workspace
; ----------------------------------------------------------------------------------
n_acc_sc=1 & Eok=0
if n_elements(datp.e) gt 1 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,'DEBYESCHERRER: number of scans = ',strtrim(n_acc_sc,2)
if keyword_set(verbose) then $
print,'DEBYESCHERRER: straightening ' ,strtrim(h_pixels,2),' pixels * ',strtrim(rc_s1+1,2),' angles: ',strtrim(XX[0],2),' ->',strtrim(XX[win_s1],2),' deg.'
; ----------------------------------------------------------------------------------
; For each detector height calculate the 2theta vector
; ----------------------------------------------------------------------------------
for j=0,h_pixels-1 do begin
zq = sqrt(rsin2t2 + y2[j]) ; height -h/2 to +h/2
Alpha = atan(zq,rcos2t)*RAD2DEG ; calculate 2theta values for the equatorial plane
if Alpha[0] lt Alpha[rc_s1] then $
xdx = where((XX ge Alpha[0]) and (XX le Alpha[rc_s1])) $ ;find corresponding indices Alpha->XX
else $
xdx = where((XX le Alpha[0]) and (XX ge Alpha[rc_s1])) ;(ascending or descending)
if xdx[0] ge 0 then begin
nxe = n_elements(xdx)-1
phi = XX[xdx]