straight_2d.pro 13.3 KB
Newer Older
1
2
3
FUNCTION straight_2d, win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, $
                           texture=txtr, integrate=integ, Weights=Iweight, $
                           verbose=verbose
Didier Richard's avatar
Didier Richard committed
4
5
6
7
;+
;******* ***********
;**
;** User callable
8
;** The call is w6 = straight_2d(w1 [,distance=samp2detect_Cm] [,/texture])
Didier Richard's avatar
Didier Richard committed
9
;**
10
11
12
13
14
15
16
17
;** 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.
Didier Richard's avatar
Didier Richard committed
18
19
20
21
;**
;** !!! WORKS BETTER WITHOUT INTERPOL (& with regular angles).
;** !!! Do not use the same input/output workspace.
;**
22
23
24
25
26
27
;** 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.
Didier Richard's avatar
Didier Richard committed
28
;**
29
;** Dec 2008: Correct vertical position for textured diffractions (use keyword /TEXTURE)
Didier Richard's avatar
Didier Richard committed
30
;**
31
32
;** Dec 2009: !!! Replace -1 (no data) by -1E-9
;** Jan 2009: (texture) Replace arrows*fY/fX  by  arrows*fX/fY
Didier Richard's avatar
Didier Richard committed
33
;**
34
35
;** Nov 2011: Join straight_2d and straight_1d, adding an integrate keyword
;**           Compute a weighting vector for integration (returned in keyword)
Didier Richard's avatar
Didier Richard committed
36
;**
37
;** Nov 2013: Loop in C code implementation
Didier Richard's avatar
Didier Richard committed
38
39
;-

40
41
FORWARD_FUNCTION interpol, regular
st1 = systime(1)
Didier Richard's avatar
Didier Richard committed
42
43
44
45
46
47

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

48
49
50
51
52
53
; 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
54
55
56
57
   inst = ''
   tags = tag_names(datp)
   idx = where(tags eq 'OTHER_TIT', found)
   if found then inst = strlowcase((strsplit(datp.other_tit, /extract))[0])
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
   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
75
lamp_journal, 'Starting straight_2d with r = ' + string(r)
76
77
78
79
80

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
81
82
83
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) 
84
85
86
87

; Function works only with regular 2thetha !!!
regul = "" 
regu = 0
Didier Richard's avatar
Didier Richard committed
88
if not keyword_set(inter) then begin
89
90
91
92
93
94
95
96
97
98
   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
Didier Richard's avatar
Didier Richard committed
99

100
101
; 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
Didier Richard's avatar
Didier Richard committed
102

103
104
105
106
; Don't use negative angles !!!
wzero = where(XX ge 0.)
w0 = wzero[0]
nx = n_elements(XX)
Didier Richard's avatar
Didier Richard committed
107

108
109
110
rcos2t  = r*cos(XX[wzero]*!pi/180.0)
rsin2t2 = (r*sin(XX[wzero]*!pi/180.0))^2
y2      = datp.y^2
Didier Richard's avatar
Didier Richard committed
111

112
113
114
115
; 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
Didier Richard's avatar
Didier Richard committed
116

117
; Define output array and initialise it
Didier Richard's avatar
Didier Richard committed
118
119
Wout    = win*0.
Eout    = 0
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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) )
Didier Richard's avatar
Didier Richard committed
135
endif
136
137
138
139
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)
140
; Need to use float(!const.rtod) to avoid error (would need to explore this issue)
141
if keyword_set(txtr) then begin
142
143
144
145
  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
146
147
148
149
150
151
152
153
154
155
156
  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
Didier Richard's avatar
Didier Richard committed
157
158

   if keyword_set(txtr) then begin
159
160
161
162
163
164
165
166
      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

167
168
   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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253

   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
Didier Richard's avatar
Didier Richard committed
254
            endelse
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308

         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,*,*])
Didier Richard's avatar
Didier Richard committed
309
310

if not keyword_set(integ) then begin
311
312
313
314
315
316
   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
Didier Richard's avatar
Didier Richard committed
317
318
319
320
endif

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

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

Didier Richard's avatar
Didier Richard committed
324
;sum over the height to output 1d array
325
326
media   = median(Iweight)>1. 
Iweight = media/(Iweight>2.)
Didier Richard's avatar
Didier Richard committed
327

328
329
Wout_1d = reform(total(wout,2))     
WOUT=0 ;FREEs
Didier Richard's avatar
Didier Richard committed
330

331
for z=0,kz do Wout_1d[*,z] *= Iweight
Didier Richard's avatar
Didier Richard committed
332

333
334
totIDM  = totIN  / total(Wout_1d) ;Keep same count
Wout_1d = temporary(Wout_1d) * totIDM
Didier Richard's avatar
Didier Richard committed
335

336
337
338
339
340
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
Didier Richard's avatar
Didier Richard committed
341

342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
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' 
Didier Richard's avatar
Didier Richard committed
357

358
if pasp then give_datp, datp
Didier Richard's avatar
Didier Richard committed
359
360
361
362

return, wout_1d

end