-
Notifications
You must be signed in to change notification settings - Fork 0
/
lstdrv_example.pro~
90 lines (64 loc) · 2.02 KB
/
lstdrv_example.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
; synthetic sources
x_src = [10, 12]
y_src = [10, 10]
f_src = [1, 2]
;x_src = [10]
;y_src = [10]
;f_src = [1]
n_src = n_elements(x_src)
; map
naxis1 = 20
naxis2 = 20
make_2d, findgen(naxis1), findgen(naxis2), x_pix, y_pix
map = fltarr(naxis1, naxis2)
n_pix = naxis1*naxis2
; sources in map
map[x_src, y_src] = f_src
; PRF
paxis1 = 5
paxis2 = 5
fwhm = 2
prf=psf_gaussian(npixel=[paxis1,paxis2],fwhm=fwhm)
; Convolution
map = convol(map, prf)
; check maximum pixel is where we expect it to be
tmp = max(map, ij)
whereimage, size(map), ij, i, j
print, i, j
; matrix of pixel <-->src positions
LSTDRV_MATRIX, X_PIX, Y_PIX, X_SRC, Y_SRC, DX, DY, DR
; getting PRF values at dx, dy offsets
; dx and dy are defined relative to the centre of the PRF image
; we need them in IDL pixel coordinates
px = dx+paxis1/2
py = dy+paxis2/2
good = where(px GE 0 AND px LT paxis1 AND py GE 0 AND py LT paxis2, ngood, comp=bad, ncomp=nbad)
PRF_matrix = fltarr(n_pix, n_src)
IF ngood GT 0 THEN PRF_matrix[good] = prf[px[good], py[good]]
IF nbad GT 0 THEN PRF_matrix[bad] = 0
; check my making up an image using the matrix approach
map2 = PRF_matrix#f_src
map2 = reform(map2, naxis1, naxis2)
tmp = max(map2, ij)
whereimage, size(map2), ij, i, j
good = where(map GT 0)
good2 = where(map2 GT 0)
; quick check
icplot, map-map2
;----------------------------------------------------------------------
; writing out data
;----------------------------------------------------------------------
name = 'eg1'
get_lun, unit
openw, unit, 'eg1.data'
printf, unit, '# '+string(n_pix)+' elelement DATA vector'
data = reform(map, n_pix)
printf, unit, data
free_lun, unit
get_lun, unit
openw, unit, 'eg1.matrix'
printf, unit, '# '+string(n_pix)+' X'+string(n_src)+' elelement N_DATA X N_SRC response matrix. Formatted as a sparse matrix i,j, M(i,j) with 0=<i<N_DATA, 0=<j<N_SRC and only if M(i,j) non zero'
make_2d, indgen(n_pix), indgen(n_src), i, j
good = where(PRF_matrix NE 0, ngood)
FOR k=0, ngood-1 DO printf, unit, i[good[k]], j[good[k]], PRF_matrix[good[k]]
free_lun, unit