Commit 5d64d110 authored by Miguel Angel Gonzalez's avatar Miguel Angel Gonzalez

Changes done by Thomas Hansen

parent f70ba93f
Function mergd20,wn,nosum=nosum,nointerpolation=noint,threshold=thresh,$
sdev=sdev,noprint=noprint,bad=bad,mult=bad_mult,bin=bin,$
norm=norm,datp=datp,min_count=min_count
norm=norm,datp=datp,min_count=min_count,xr=xr
;+
; Merge twotheta scan of n steps (PSD cell width must be 0.1 deg!)
; W1 = mergd20 (W2,/nosum)
......@@ -16,8 +16,8 @@ compile_opt idl2
COMMON C_LAMP_INFO
bad_mult=INTARR(1600)
IF NOT KEYWORD_SET(bin) THEN bin=.05 ELSE bin=bin/2.
bad_mult=INTARR(3200)
;IF NOT KEYWORD_SET(bin) THEN bin=.05 ELSE bin=bin/2.
IF NOT KEYWORD_SET(thresh) THEN thresh=3.
IF N_ELEMENTS(wn) GT 1 THEN BEGIN
n=two
......@@ -32,6 +32,13 @@ IF N_ELEMENTS(n) EQ 1 THEN BEGIN
x=d.x
scans=N_ELEMENTS(d.x[0,*])
cells=N_ELEMENTS(d.x[*,0])
bad_mult=INTARR(cells);
IF NOT KEYWORD_SET(bin) THEN BEGIN
bin=ABS(median(x[0:N_elements(x[*,0])-2,*]-x[1:*,*])/2.)
print,"bin:",bin
ENDIF ELSE bin=bin/2.
cellsperdegree=ROUND(1./ABS(median(x[0:N_elements(x[*,0])-2,*]-x[1:*,*])/2.))/2
print,"Cells per degree:",cellsperdegree
IF N_ELEMENTS(w[0,*]) GT 2 THEN BEGIN
;############# Normalisation #############
norm=TOTAL(d.n[0,0,*])/N_ELEMENTS(d.n[0,0,*])
......@@ -42,31 +49,38 @@ IF N_ELEMENTS(n) EQ 1 THEN BEGIN
n[*,*,i]=n[*,*,i]*norm/d.n[0,0,i]
ENDFOR
;############# 2Theta ####################
delta=(10.*d.x[0,*]-FLOOR(10.*d.x[0,*]))/10.
delta=(cellsperdegree*d.x[0,*]-FLOOR(cellsperdegree*d.x[0,*]))/cellsperdegree
test=ABS(median(delta)-total(delta)/N_ELEMENTS(d.x[0,*]))
delta=(10.*(d.x[0,*]-.05)-FLOOR(10.*(d.x[0,*]-.05)))/10.
IF test LT ABS(median(delta+.05)-total(delta+.05)/N_ELEMENTS(d.x[0,*])) THEN BEGIN
delta=(10.*d.x[0,*]-FLOOR(10.*d.x[0,*]))/10.
delta=(cellsperdegree*(d.x[0,*]-0.5/cellsperdegree)-FLOOR(cellsperdegree*(d.x[0,*]-.5/cellsperdegree)))/cellsperdegree
IF test LT ABS(median(delta+.5/cellsperdegree)-total(delta+.5/cellsperdegree)/N_ELEMENTS(d.x[0,*])) THEN BEGIN
delta=(cellsperdegree*d.x[0,*]-FLOOR(cellsperdegree*d.x[0,*]))/cellsperdegree
delta=MEDIAN(delta)
x=FLOOR(10.*d.x)/10.+delta
x=FLOOR(cellsperdegree*d.x)/(1.*cellsperdegree)+delta
print,'+ ','Delta: ',delta
ENDIF ELSE BEGIN
delta=MEDIAN(delta+.05)
delta=((10.*delta)-FLOOR(10.*delta))/10.
x=FLOOR(10.*(d.x-.05))/10.+delta
delta=MEDIAN(delta+.5/cellsperdegree)
delta=((cellsperdegree*delta)-FLOOR(cellsperdegree*delta))/cellsperdegree
x=FLOOR(cellsperdegree*(d.x-.5/cellsperdegree))/(1.*cellsperdegree)+delta
print,'- ','Delta: ',delta
ENDELSE
FOR i=0,N_ELEMENTS(w[0,*])-1 DO BEGIN
IF ABS(d.x[0,i]-x[0,i]) GE bin THEN BEGIN
IF NOT KEYWORD_SET(noint) THEN BEGIN
PRINT,'Step',STRCOMPRESS(i),' outside limit, interpolated: ' ,d.x[0,i],x[0,i]
PRINT,'Step',STRCOMPRESS(i),' outside limit, interpolated: ' ,d.x[0,i],x[0,i],ABS(d.x[0,i]-x[0,i]),bin
w[*,i]=INTERPOL(w[*,i],d.x[*,i],x[*,i])
d.e[*,i]=INTERPOL(d.e[*,i],d.x[*,i],x[*,i])
d.x[*,i]=x[*,i]
ENDIF ELSE PRINT,'Step',STRCOMPRESS(i),' outside limit, rejected: ' ,d.x[0,i],x[0,i]
ENDIF ;ELSE PRINT,'Step',STRCOMPRESS(i),' inside limit: ' ,d.x(0,i),x(0,i)
ENDFOR
xrange=[FLOOR(10.*min(d.x))/10.+delta,CEIL(10.*max(d.x))/10.+delta]
xrange=[FLOOR(cellsperdegree*min(d.x))/cellsperdegree+delta,CEIL(cellsperdegree*max(d.x))/cellsperdegree+delta]
xrange=round(xrange*20.)/20.
;help,xr,xrange,cellsperdegree
if keyword_set(xr) then begin
xrange=xr
cellsperdegree=20
endif
;help,xr,xrange
x=0
e=0
y=0
......@@ -76,14 +90,24 @@ IF N_ELEMENTS(n) EQ 1 THEN BEGIN
IF KEYWORD_SET(sdev) THEN min_count=min_count>2 ELSE min_count=min_count>1
ENDELSE
IF NOT KEYWORD_SET(nosum) THEN BEGIN
FOR i=xrange[0],xrange[1],0.1 DO BEGIN
;print,xrange[0],xrange[1],1./cellsperdegree
stepwidth=1./cellsperdegree
;FOR i=xrange[0]*,xrange[1],1./cellsperdegree DO BEGIN
;FOR j=xrange[0]*20,xrange[1]*20 DO BEGIN
FOR j=xrange[0]*cellsperdegree,xrange[1]*cellsperdegree DO BEGIN
;i=j/20.
i=j/cellsperdegree
; print,i
; FOR i=xrange[0],xrange[1],.05 DO BEGIN
index=WHERE(d.x GE i-bin AND d.x LT i+bin,count)
;print,i,bin,count
IF count GT 0 THEN BEGIN
tmp_w = w[index]
tmp_m = MEDIAN(tmp_w)
tmp_e = d.e[index]
tmp_i= WHERE(ABS(tmp_w-tmp_m) LE thresh*tmp_e,count)
tmp_n= WHERE(ABS(tmp_w-tmp_m) GT thresh*tmp_e,count_n)
;print,i,bin,count,tmp_m,median(tmp_e)
IF count EQ 1 THEN BEGIN
tmp=min(tmp_w[WHERE(tmp_w GE tmp_m)]-tmp_m,tmp_pos)
tmp=max(tmp_w[WHERE(tmp_w LE tmp_m)]-tmp_m,tmp_neg)
......@@ -108,6 +132,7 @@ IF N_ELEMENTS(n) EQ 1 THEN BEGIN
IF tmp_0 GT 0 THEN BEGIN
mean=TOTAL(w[index])/count
x=[x,i]
;print,i
y=[y,mean]
IF count GT 1 THEN BEGIN
sdev=SQRT(TOTAL((w[index]-mean)^2))/((count-1)>1)
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment