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
|
/* Given a line with a continuum, subtract the continuum.
To do this, first a straight line is fit to the lower
and upper 10% (or at least 5 values each) of the vector,
this fit subtracted and a second fit then done to the
values "close" (within 1 sigma) to this initial fit and
again subtracted.
*/
function subcontinuum, line {
low = round(nelements(line) / 10)
if (low < 5) low = 5
high = nelements(line) - low
xvalues = yvalues = doublearray(low * 2)
for x=1, low {
xvalues[x] = x
yvalues[x] = line[x]
xvalues[2*low - x + 1] = nelements(line) - x + 1
yvalues[2*low - x + 1] = line[nelements(line) - x + 1]
}
fit = polyfitxy(xvalues, yvalues, 1)
initialguess = line - fit[1] - fit[2] * [1:nelements(line)]
deviation = stddev(initialguess)
xvalues = yvalues = doublearray(nelements(line))
yvalues = line
counter = 0
for x=1, nelements(line) {
if (abs(initialguess[x]) < deviation) {
counter++
xvalues[counter] = x
yvalues[counter] = line[x]
}
}
xvalues = xvalues[1:counter]
yvalues = yvalues[1:counter]
fit = polyfitxy(xvalues, yvalues, 1)
subcontinuum = line - fit[1] - fit[2] * [1:nelements(line)]
}
function moments, x, y {
moments = [0, 0]
indices = where(y > 0)
if (nelements(indices) > 0) {
xx = x[indices]
yy = y[indices]
sumY = total(yy)
sumXY = total(xx * yy)
moments[1] = sumXY / sumY
moments[2] = sqrt(total(yy*(xx-moments[1])^2)/total(yy))
}
}
function clipsigma, line, deviation {
clipsigma = line
clipsigma[where(abs(line) < deviation)] = 0
}
userDialog "myvelmap",["cube","wavelength center","fwhm"],["fits","real","real"],["","2.13",".001"]
userDialog "sin",["argument","Do this in radians","Do this in degrees"],["fits","/rad","/deg"],["","1","0"]
function myvelmap, cube, center, fwhm {
myvelmap = doublearray(naxis1(cube), naxis2(cube), 8)
centralwave = getfitskey(cube, "CRVAL3")
centralpix = getfitskey(cube, "CRPIX3")
dispersion = getfitskey(cube, "CDELT3")
zc = round((center - centralwave) / dispersion + centralpix)
z1 = zc - round(fwhm / dispersion * 5)
z2 = zc + round(fwhm / dispersion * 5)
if (z1 < 1) z1 = 1
if (z2 > naxis3(cube)) z2 = naxis3(cube)
xvalues = cube[1,1,z1:z2]
centralwave = getfitskey(xvalues, "CRVAL1")
centralpix = getfitskey(xvalues, "CRPIX1")
dispersion = getfitskey(xvalues, "CDELT1")
xvalues = centralwave + ([1:nelements(xvalues)] - centralpix) * dispersion
errors = xvalues * 0 + 1
estimate = doublearray(4)
for x=1, naxis1(cube) {
for y=1, naxis2(cube) {
line = cube[x, y, z1:z2]
line = subcontinuum(line)
//line[where(abs(line) < stddev(line))] = 0
sss = stddev(line)
line = clipsigma(line, sss)
//estimate[1] = 0
//estimate[2] = max(line)
//estimate[3] = centralwave + (xmax(line) - centralpix) * dispersion
//estimate[4] = fwhm
//myvelmap[x,y,*] = gaussfit(xvalues, line, errors, estimate)
myvelmap[x,y,3:4] = moments(xvalues, line)
}
print x
}
}
function mysin, arg {
mysin = sin( arg ) * 2
}
function mymysin, arg {
mymysin = mysin(arg * 2)
}
|
