-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstix_sim_dirty_map_from_vs.pro
146 lines (121 loc) · 5.99 KB
/
stix_sim_dirty_map_from_vs.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
145
146
FUNCTION stix_sim_dirty_map_from_vs, visibilities, weights=weights, npix=npix, pixsize=pixsize, xc=xc, yc=yc, xvalues=xvalues, yvalues=yvalues
;+
; :project:
; STIX imaging system simple simulation
;
; :description:
; This function takes 30 complexe visibilities and create a dirty map
;
; :inputs:
; visibilities: array of 30 complex numbers (visibilities), typically calculated with the function stix_sim_calculate_visibilities
;
; :keywords:
; weights: in, default is 30 identical values: array of 30 values containing the weights to apply to each visibility (collimator)
; npix: in, default is 64: number of pixels in each dimension of the image
; pixsize: in, default is 4: pixel size in arcsec
; xc: in, default is 0: center of the map in arcsec
; yc: in, default is 0: idem
; xvalues, out, : array containing the x values for each pixel
; yvalues, out, : array containing the y values for each pixel
;
; :history:
; 2018/07/30, Sophie Musset (University of Minnesota), initial release
;-
;-------------------------------------------------------------
; define default values for weight, several possibilities
;-------------------------------------------------------------
; DEFAULT, weights, replicate(1./n_elements(visibilities), n_elements(visibilities))
; DEFAULT, weights, exp( alog10([replicate(100.,3), replicate(10.*9./10.,3), replicate(10.*8./10.,3), $
; replicate(10.*7./10.,3), replicate(10.*6./10.,3), replicate(10.*5./10.,3), $
; replicate(10.*4./10.,3), replicate(10.*3./10.,3), replicate(10.*2./10.,3), $
; replicate(10.*1./10.,3)]) )
pitches = [replicate(0.0380,3),replicate(0.0543,3),replicate(0.0777,3),$
replicate(0.1112,3),replicate(0.1590,3),replicate(0.2275,3),$
replicate(0.3254,3),replicate(0.4655,3),replicate(0.6659,3),$
replicate(0.9526,3)]
DEFAULT, weights, 1/pitches
;-------------------------------------------------------------
; define other default values and restore useful variables
;-------------------------------------------------------------
DEFAULT, npix, 400
DEFAULT, pixsize, 1
DEFAULT, xc, 0
DEFAULT, yc, 0
restore, 'C:\Users\SMusset\Documents\GitHub\STIX\stix_visibility_uv.sav'
; this restore the u and v variables
facto=2*!pi
;-------------------------------------------------------------
; calculate the x and y values of the image
;-------------------------------------------------------------
xvalues = indgen(npix)*pixsize - (npix*pixsize/2 - xc)
yvalues = indgen(npix)*pixsize - (npix*pixsize/2 - yc)
;-------------------------------------------------------------
; calculate the amplitude and phase of each visibitily
;-------------------------------------------------------------
real = real_part(visibilities)
imgn = imaginary(visibilities)
amplitude = sqrt( real^2 + imgn^2 )
phase = atan(visibilities, /phase)
;print, phase
;phase = dblarr(30)
;FOR k=0,29 DO BEGIN
; phasee = stix_sim_complex_phase(visibilities[k])
; phase[k] = phasee
;ENDFOR
;-------------------------------------------------------------
; construct the dirty image
;-------------------------------------------------------------
; initialisation
dmap = dblarr(npix,npix)
; calculation of the value in each pixel
FOR i=0, npix-1 DO BEGIN
FOR j=0, npix-1 DO BEGIN
dmap[i,j] = total(weights*amplitude*cos(facto*(u*xvalues[i]+v*yvalues[j])+phase))/total(weights)
ENDFOR
ENDFOR
;-------------------------------------------------------------
; do intermediate maps (for info plot)
;-------------------------------------------------------------
dmap1 = dblarr(npix,npix)
FOR i=0, npix-1 DO BEGIN
FOR j=0, npix-1 DO BEGIN
dmap1[i,j] = total(weights[0]*amplitude[0]*cos(facto*(u[0]*xvalues[i]+v[0]*yvalues[j])+phase[0]))/total(weights[0])
ENDFOR
ENDFOR
dmap2 = dblarr(npix,npix)
FOR i=0, npix-1 DO BEGIN
FOR j=0, npix-1 DO BEGIN
dmap2[i,j] = total(weights[0:1]*amplitude[0:1]*cos(facto*(u[0:1]*xvalues[i]+v[0:1]*yvalues[j])+phase[0:1]))/total(weights[0:1])
ENDFOR
ENDFOR
dmap4 = dblarr(npix,npix)
FOR i=0, npix-1 DO BEGIN
FOR j=0, npix-1 DO BEGIN
dmap4[i,j] = total(weights[0:3]*amplitude[0:3]*cos(facto*(u[0:3]*xvalues[i]+v[0:3]*yvalues[j])+phase[0:3]))/total(weights[0:3])
ENDFOR
ENDFOR
dmap8 = dblarr(npix,npix)
FOR i=0, npix-1 DO BEGIN
FOR j=0, npix-1 DO BEGIN
dmap8[i,j] = total(weights[0:7]*amplitude[0:7]*cos(facto*(u[0:7]*xvalues[i]+v[0:7]*yvalues[j])+phase[0:7]))/total(weights[0:7])
ENDFOR
ENDFOR
dmap16 = dblarr(npix,npix)
FOR i=0, npix-1 DO BEGIN
FOR j=0, npix-1 DO BEGIN
dmap16[i,j] = total(weights[0:15]*amplitude[0:15]*cos(facto*(u[0:15]*xvalues[i]+v[0:15]*yvalues[j])+phase[0:15]))/total(weights[0:15])
ENDFOR
ENDFOR
;-------------------------------------------------------------
; for info, plot images with different number of visibilities
;-------------------------------------------------------------
i=image(dmap1, xvalues, yvalues, margin=0.1, title='1 vis', layout=[5,1,1], axis_style=1, rgb_table=5, dimensions = [2100,500])
i=image(dmap2, xvalues, yvalues, margin=0.1, title='2 vis', layout=[5,1,2], axis_style=1, rgb_table=5, /current)
i=image(dmap4, xvalues, yvalues, margin=0.1, title='4 vis', layout=[5,1,3], axis_style=1, rgb_table=5, /current)
i=image(dmap8, xvalues, yvalues, margin=0.1, title='8 vis', layout=[5,1,4], axis_style=1, rgb_table=5, /current)
i=image(dmap16, xvalues, yvalues, margin=0.1, title='16 vis', layout=[5,1,5], axis_style=1, rgb_table=5, /current)
;-------------------------------------------------------------
; return the dirty image
;-------------------------------------------------------------
RETURN, dmap
END