Commit 52dfba69 authored by Gonzalez, Miguel's avatar Gonzalez, Miguel
Browse files

Changes to integrate and straight routines to handle better the needs for...

Changes to integrate and straight routines to handle better the needs for Salsa, including the hexapod scan (where datp.z is a string array)
parent 393fc5f4
function str_integrate, w_in,r
function str_integrate, w_in, r
;+
;******* *************
;** User callable
;** function to read a workspace containing a number of 2D data sets
;** Integrate the 2D data to 1D
;** Typical call w2=str_integrate(w1)
;** Typical call w2 = str_integrate(w1)
;-
take_datp,datp
;find dimensions/size of w_in
nx =(size(datp.x))[1]
nspec =(size(datp.z))[1] ; & help,datp,/struct & print, datp.z
w_out =fltarr(nx,nspec) ; array for integrated data
take_datp, datp
;integrate data
for ispec=0,nspec-1 do begin
w_out[*,ispec]=total(w_in[*,*,ispec],2)
; Find dimensions/size of w_in
nx = (size(datp.x))[1]
nspec = (size(datp.z))[1] ; & help,datp,/struct & print, datp.z
w_out = fltarr(nx, nspec) ; array for integrated data
; Integrate data
for ispec = 0, nspec-1 do begin
w_out[*,ispec] = total(w_in[*,*,ispec], 2)
endfor
;return 2D workspace
mod_datp,datp,'y',datp.z
datp.y_tit='scan'
give_datp,datp
; Deal with the Hexapod case in Salsa, where the point information is stored as a string chain
datp.y_tit = 'scan'
if typename(datp.z) eq 'STRING' then mod_datp, datp, 'y', fix(indgen(nspec) + 1) $
else mod_datp, datp, 'y', datp.z
; Return 2D workspace
give_datp, datp
return, w_out
end
FUNCTION straight_1d, win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, texture=txtr, integrate=integ, Weights=Iweight
FUNCTION straight_1d, 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_1d(w1)
;** The call is w6 = straight_1d(w1)
;**
;** win is a 2D (3D) diffractogram i.e. I(2theta,h)
;** where h is the physical height of the multidetector
;** wout is a 1D diffractogram i.e. sum over h of 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
;** Win is a 2D (3D) diffractogram, i.e. I(2theta,h)
;** where h is the physical height of the multidetector.
;**
;** h is the physical height of the multidetector
;** r is the physical radius of the multidetector
;** set them
;** Wout is a 1D diffractogram, i.e. sum over h of I(2theta,h)
;** in which the Debye-Scherrer cones have been straightened out.
;**
;** Nov 2011 call straight_2d adding the integrate keyword
;** For each detector height, it calculates a vector of 2theta
;** angles and "interpolates" onto equatorial plane
;**
;** Set the right values for:
;** h = physical height of the multidetector --> in Y
;** r = physical radius of the multidetector --> in call
;**
;** Nov 2011: Call straight_2d adding the integrate keyword
;-
return, straight_2d(win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, /integrate, Weights=Iweight)
return, straight_2d(win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, /integrate, Weights=Iweight, verbose=verbose)
end
FUNCTION straight_2d, win, interpol=inter, distance=e_d, datp=datp, dX=deltaX, texture=txtr, integrate=integ, Weights=Iweight, verbose=verbose
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])
;** 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 DS cones should be straightened out
;** For each detector height calculate a vector of 2theta
;** angles and interpolate onto equatorial plane
;** 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.
;**
;** h is the physical height of the multidetector
;** r is the physical radius of the multidetector
;** set them
;** 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.
;**
;** 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 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
;**
;** 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 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
;** Nov 2013: Loop in C code implementation
;-
FORWARD_FUNCTION interpol,regular
st1=systime(1)
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
;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.
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
;works only with regular 2thetha !!!
regul="" & regu=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 = 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 for ' + inst + ' 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
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
; Regular will fail if an axis value is not numeric (e.g. Hexapod scan in Salsa)
if typename(datp.x) eq 'STRING' or typename(datp.y) eq 'STRING' or typename(datp.z) eq 'STRING' then begin
lamp_journal, 'Non numeric axis (hexapod scan?)'
lamp_journal, 'Convert to numeric axis before calling straight_2d,'
lamp_journal, 'e.g. z1 = indgen(n_elements(z1)).'
lamp_journal, 'Returning input workspace!'
return, win
endif
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
; 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)
; 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
todeg = 180.0/!pi
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
; 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
;if (pasp) then give_w,0 ;(clear output workspace so we gain memory)
; Define output array and initialise it
Wout = win*0.
Eout = 0
win_s1 =(size(win))[1] -1
rc_s1 = n_elements(rcos2t)-1
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 has 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))
; 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 ,.01
;print,'number of scans=',strtrim(n_acc_sc,2),' ,straightening ' ,h_pixels,' , ',XX[0],' ->',XX[win_s1],regul
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 2theta (needed for textures)
if keyword_set(txtr) then begin
; Calculation of the rays for each 2-theta (needed for textures)
if keyword_set(txtr) then begin
zq = sqrt(rsin2t2 + y2[0]) ; height -h/2 to +h/2 (j=0)
xh = atan(zq,rcos2t)*todeg ; calculate 2theta values for the equatorial plane xh
XHcm = sin(xh/todeg)*r & XXcm=sin(XX[wzero]/todeg)*r ; we work in centimeters
fY = abs(datp.y[1]-datp.y[0]) & fX=abs((XX[nx-1]-XX[0])/nx); delta Y,X
corde = abs(datp.y[0]*2.) ; corde of the arcs
;arrows =(abs((XHcm-XXcm)*fY/fX))>(1e-8) ; arrows of the arcs normalized with XYfct-> *fY/fX (bad!)
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)) ; lenght of the arcs unused
endif
jof=rcos2t*0.
xh = atan(zq, rcos2t) * !const.rtod ; calculate 2theta values for the equatorial plane xh
XHcm = sin(xh/!const.rtod) * r
XXcm = sin(XX[wzero]/!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 2theta vector
for j=0L,h_pixels-1 do begin
; 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) & PAj2=(jof-jof1) & PAj1=1.-PAj2 ;vertical weighting
endif else jof[*]=j
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)* !const.rtod ; calculate 2theta values for the equatorial plane xh
zq =sqrt(rsin2t2 + y2[j]) ; height -h/2 to +h/2
xh =atan(zq,rcos2t)*todeg ; 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 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
if xdx[0] ge 0 then begin & nxe=n_elements(xdx)-1
xeq=XX[xdx]
glob_low_idx[j]=xdx[0] & glob_hi_idx[ j]=xdx[nxe] ;keep first and last 2theta index (straight_1d)
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]
v[0,*]=v[1,*] & v[rc_s1,*]=v[rc_s1-1,*] ;Will prevent from interpol bounds artifacts
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. ;Ok I know!
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
;sk =value_locate(xh, xeq)>0L<(rc_s1-1) & sk1=sk+1 ;(xh,v) (xeq,sk,r)
;PA =(xeq-xh(sk))/(xh(sk1)-xh(sk)) & PA1=1.-PA
;v =win(wzero,j,*)
;v(0,*,*)=v(1,*,*) & v(rc_s1,*,*)=v(rc_s1-1,*,*) ;Will prevent from interpol bounds artifacts
;if Eok then e=datp.e(wzero,j,*)
;for k=0L,nxe do begin
; wout(xdx(k),j,*)=(v(sk(k),*,*)*PA1(k) +v(sk1(k),*,*)*PA(k))>0.
; if Eok then Eout(xdx(k),j,*)=(e(sk(k),*,*)*PA1(k) +e(sk1(k),*,*)*PA(k))>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
PA =((xh[ind]-xeq[sk])/(xeq[sk+1]-xeq[sk]))>0.<1. & PA1=1.-PA ;*** horizontal weighting
v =win[wzero,j,*] & if Eok then e=datp.e[wzero,j,*]
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
PA11=PA1*PAj1[ind] & PA01=PA*PAj1[ind] & PA12=PA1*PAj2[ind] & PA02=PA*PAj2[ind] ;*** hori & vert weighting for texture
if SFdll gt ' ' then begin ;Try the C code if there.
;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 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)
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
catch, /cancel
if (bid ne 1) then begin
lamp_loop_so, /reset
SFdll = ''
endif
if SFdll le ' ' then 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,*,*]
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])
......@@ -187,31 +253,36 @@ for j=0L,h_pixels-1 do begin
Eout[ik1[k],jj2,*] += (ee*PA02[k])
endif
endfor
endelse
endif
endif else begin ;*** horizontal weighting only
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 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)
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
if (bid ne 1) then begin
lamp_loop_so, /reset
SFdll = ''
endif
if SFdll le ' ' then 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,*,*]
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
......@@ -221,65 +292,73 @@ for j=0L,h_pixels-1 do begin
Iweight[ik[ k]] += PA1[k] ;Calculate weights for integration
Iweight[ik1[k]] += PA[ k]
endfor
endelse
endif
endelse ;end if keyword_set(txtr) else block
endelse ;end if keyword_set(inter) else block
endelse ;keyword_set(txtr)
endelse
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
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
endfor
totIN =total(win[wzero,*,*])
endif ;end if xdx[0]>0 block
if not keyword_set(integ) then begin
totIDM=totIN / total(Wout) ;Keep same count
Wout =temporary(Wout) * totIDM
Eout =temporary(Eout) * totIDM
endfor ;end loop over h_pixels
if Eok then datp.e=Eout
if keyword_set(verbose) then print, 'straightening done ',round((systime(1)-st1)*100)/100.,'s'
totIN = total(win[wzero,*,*])
if pasp then if (regu or Eok) then give_datp,datp
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, 'integrating ',round(systime(1)-st1),'s'
if keyword_set(verbose) then wait ,.01
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.)
media = median(Iweight)>1. & Iweight=media/(Iweight>2.)
Wout_1d = reform(total(wout,2))
WOUT=0 ;FREEs
Wout_1d= reform( total(wout,2)) & WOUT=0 ;FREEs
for z=0,kz do Wout_1d[*,z] *= Iweight
for z=0,kz do Wout_1d[*,z] *= Iweight
totIDM =totIN / total(Wout_1d) ;Keep same count
Wout_1d=temporary(Wout_1d)* totIDM
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
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]
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