-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgconv.pro
144 lines (133 loc) · 4.31 KB
/
gconv.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
pro gconv,x,y,fwhm,nx,ny,ppr=ppr,nsigma=nsigma,inter=inter,original=original
;+
; Gaussian convolution
;
; IN: x -fltarr independent variable
; y -fltarr array of data to smooth
; fwhm -float FWHM of the Gaussian kernel (pixels)
;
; OUT: nx -fltarr new sampled x-axis
; ny -fltarr smoothed data
;
; KEYWORDS: ppr -float pixels per resolution element
; (default is ppr=3.)
;
; nsigma -float how far the convolution goes
;
; inter if set, the x-axis is interpolated
; to force a constant step; if the step
; already constant, 'inter' will have no
; effect
;
; original if set, the original sampling is kept
; for the output arrays, unless inter is
; also set and an interpolation is
; performed, in which case, the output
; x-axis will just have the same
; number of points as the input x, but
; resampled to constant step.
;
; EXAMPLES: 1) Smooth an input spectrum with a Gaussian of FWHM=2.0 AA
; assuming the spectrum is given on a uniform (linearly sampled)
; wavelength scale with a step of 0.1 AA.
;
; IDL> fwhm=2.0/0.1
; IDL> gconv,w,f,fwhm,w2,f2
;
; 2) Smooth the same spectrum with a Gaussian of FWHM=50. km/s
;
; IDL> step=(max(alog(w))-min(alog(w)))/n_elements(w)
; IDL> fwhm=50./299792.458/step
; IDL> gconv,alog(w),f,fwhm,w2,f2,/inter
;
; Because of the keyword inter, the spectrum will be internally
; interpolated to a step in ln (w) of
; step=(max(alog(w))-min(alog(w)))/n_elements(w)
; Then the 50 km/s (or V/c = 0.000166782) is equivalent to
; 50./299792.458/step pixels.
;
; C. Allende Prieto, Univ. of Texas, March 2006
; "" , MSSL/UCL, July 2009 -- changed to double, avoid psf_gaussian
; , MSSL/UCL, August 2009 -- modified to avoid rate=0
; , IAC, January 2010 -- changed indgen -> indgen(,/long)
; , IAC, January 2011 -- bug fixed (checking linear or log x scale)
; , IAC, November 2012 -- changed ppr type from integer to float
;-
if (N_params() lt 4) then begin
message,'use -- gconv,x,y,fwhm,nx,ny[,ppr=ppr,nsigma=nsigma',/cont
message,' inter=inter,original=original]',/cont
return
endif
;initialize
tol=1d-3
nx=0.0d0
ny=0.0d0
if keyword_set(inter) then begin
nel=n_elements(x)
xinput=x
yinput=y
minx=min(x)
step=(max(x)-minx)/nel
x=dindgen(nel)*step+minx
y=interpol(yinput,xinput,x)
endif else begin
;checking that the input x is evenly sampled
nel=n_elements(x) ; nel input pixels
step=abs(x-shift(x,1))
med=median(step[1:nel-1])
if max(abs(step[1:nel-1]-med)) gt tol then begin
;not linearly spaced, maybe logarithmically
lx=alog10(x)
step=abs(lx-shift(lx,1))
med=median(step[1:nel-1])
if max(abs(step[1:nel-1]-med)) gt tol then begin
message,'not log spaced, cannot handle it without interpolation',/cont
return
endif
endif
endelse
;set default params
if not keyword_set(ppr) then ppr=3.0 ;sample output function with ppr
;points per FWHM
if not keyword_set(nsigma) then nsigma=3.0 ;how far the convolution goes
;Gaussian kernel ; units are pixels
sigma=fwhm/2.d0/sqrt(-2.d0*alog(0.5d0)) ; equivalent sigma for kernel
grange=floor(nsigma*sigma) ; range for kernel (-range:range)
;psf=psf_gaussian(npix=2*grange+1,fwhm=fwhm,/normalize,ndimen=1)
psf=1.d0/sqrt(2.d0*!dpi)/sigma*exp(-((dindgen(2*grange+1)-grange)/sigma)^2/2.d0)
psf=psf/total(psf)
;select output x-axis
if keyword_set(original) then begin
;keep input sampling on output arrays
sampling=indgen(nel,/long)
endif else begin
;rate=floor(fwhm)/ppr ; old version, whih ppr being integer
rate=floor(fwhm/ppr) ; 1/rate pixels will be kept
;make sure that rate is at least 1
if rate lt 1 then rate=1
sampling=uniq(indgen(nel,/long)/rate)-rate/2
endelse
;trim edges
wtrim=where(sampling ge grange[0] and sampling lt nel-grange[0])
if max(wtrim) le 0 then begin
message,'The Gaussian is too wide for the input range',/cont
return
endif
sampling=sampling[wtrim]
nx=x[sampling] ; new x-axis
nel2=n_elements(nx) ; nel2 output pixels
ny=dblarr(nel2)
;convolve
for i=0l,nel2-1 do begin
ny[i]=total(psf*y[sampling[i]-grange:sampling[i]+grange])
;print,'y=',y[sampling[i]-grange:sampling[i]+grange]
;print,'psf=',psf
;print,'ny[i]=',ny[i]
;stop
endfor
;return the original values of x and y if resampled
if keyword_set(inter) then begin
x=xinput
y=yinput
endif
end