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