-
Notifications
You must be signed in to change notification settings - Fork 2
/
OctofitterInterferometry.jl
294 lines (256 loc) · 10.8 KB
/
OctofitterInterferometry.jl
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
module OctofitterInterferometry
using Octofitter
using PlanetOrbits
using Tables, TypedTables
using FITSIO, OIFITS
# Querying a FITSIO file by HDU name works, but if there are multiple columns,
# it automatically selects the *last* matching HDU.
# Unforuantely GRAVITY stores their data with two sets of named HDUs per file.
# The first set are from the data fiber, and the last set are from the fringe
# tracker.
# We therefore use this function to find the *first* matching HDU number
# given an HDU name.
# Note: this is certainly not threadsafe, as the underlying library is not either.
function get_first_named_hdu(fits::FITS, hduname::AbstractString)
for i in 1:length(fits)
FITSIO.fits_movabs_hdu(fits.fitsfile, i)
hduname′ = something(FITSIO.fits_try_read_extname(fits.fitsfile), "")
if hduname == hduname′
return fits[i]
end
end
throw(KeyError("FITS HDU with name \"$hduname\" not found."))
end
const required_cols = (:epoch, :u, :v, :cps_data, :dcps, :vis2_data, :dvis2, :index_cps1, :index_cps2, :index_cps3, :band, :use_vis2)
struct InterferometryLikelihood{TTable<:Table} <: Octofitter.AbstractLikelihood
table::TTable
function InterferometryLikelihood(observations...)
input_table = Table(observations...)
if :filename ∈ Tables.columnnames(input_table)
rows = map(eachrow(input_table)) do row
row = only(row)
FITS(row.filename, "r") do f
wavs = read(OIDataBlock, f["OI_WAVELENGTH",10])
vis2s = read(OIDataBlock, f["OI_VIS2",10])
cps = read(OIDataBlock, f["OI_T3",10])
#read data
eff_wave = wavs.eff_wave
vis2 = vis2s.vis2data
vis2_err = vis2s.vis2err
ut = vis2s.ucoord
vt = vis2s.vcoord
vis2_index = vis2s.sta_index
cp = cps.t3phi
cp_err = cps.t3phierr
cp_index = cps.sta_index
#convert u,v to units of wavelength
u = ut ./ eff_wave' # Units of inverse wavelength
v = vt ./ eff_wave' # Units of inverse wavelength
# These say what baseline (cp1) should be added to (cp2) and then subtract (cp3)
# to get a closure phase in our modelling.
cp_inds1, cp_inds2, cp_inds3 = cp_indices(;vis2_index, cp_index)
return (;
row.epoch,
row.band,
row.use_vis2,
u,
v,
eff_wave=eff_wave,
cps_data=transpose(cp),
dcps=transpose(cp_err),
vis2_data=transpose(vis2),
dvis2=transpose(vis2_err),
index_cps1=cp_inds1,
index_cps2=cp_inds2,
index_cps3=cp_inds3,
)
end
end
table = Table(rows)
else
table = input_table
end
if !issubset(required_cols, Tables.columnnames(table))
error("Expected columns $vis_cols")
end
return new{typeof(table)}(table)
end
end
InterferometryLikelihood(observations::NamedTuple...) = InterferometryLikelihood(observations)
export InterferometryLikelihood
"""
Visibliitiy modelling likelihood for point sources.
"""
function Octofitter.ln_like(vis::InterferometryLikelihood, θ_system, orbits, num_epochs::Val{L}=Val(length(vis.table))) where {L}
T = typeof(θ_system.M)
ll = zero(T)
# Access the data here:
epochs = vis.table.epoch
band = vis.table.band
u = vis.table.u
v = vis.table.v
cps_data = vis.table.cps_data
dcps = vis.table.dcps
vis2_data = vis.table.vis2_data
dvis2 = vis.table.dvis2
i_cps1 = vis.table.index_cps1
i_cps2 = vis.table.index_cps2
i_cps3 = vis.table.index_cps3
use_vis2 = vis.table.use_vis2
# Loop through planets
for j in eachindex(epochs)
epoch = epochs[j]
this_band = band[j]
cvis = zeros(complex(T), size(u[j]))
norm_factor = 0.0 #to normalize complex visibilities
for i in eachindex(orbits)
# All parameters relevant to this planet
θ_planet = θ_system.planets[i]
# Get model contrast parameter in this band (band provided as a symbol, e.g. :L along with data in table row.)
contrast = getproperty(θ_planet, this_band) #secondary/primary
#contrast = contrast[i]
# orbit object pre-created from above parameters (shared between all likelihood functions)
orbit = orbits[i]
sol = orbitsolve(orbit, epoch)
Δra = raoff(sol) # in mas
Δdec = decoff(sol) # in mas
#add complex visibilities from all planets at a single epoch
cvis_bin!(cvis; Δdec, Δra, contrast, u=u[j], v=v[j])
norm_factor += contrast
end
cvis .+= 1.0 #add contribution from the primary primary
cvis .*= 1.0 / (1.0 + norm_factor)
#compute closure phase
cps = closurephase(vis=cvis, index_cps1=i_cps1[j], index_cps2=i_cps2[j], index_cps3=i_cps3[j])
lnlike_v2 = 0.0
if use_vis2[j]
#compute squared visibilities
vis2 = abs.(cvis) .^ 2
const_v2 = -sum(log.(2π*(dvis2[j] .^ 2))) / 2
#calculate vis2 ln likelihood
lnlike_v2 = lnlike_v2 .+ -0.5 * sum((vis2_data[j] .- vis2) .^ 2 ./ dvis2[j] .^ 2) .+ const_v2
end
#calculate cp ln likelihood
# TODO: move mask into object creation and add warning there
mask = dcps[j] .> 0
const_cp = -sum(log.(2π*(dcps[j][mask] .^ 2))) / 2
lnlike_cp = -0.5 * sum((cps_data[j][mask] .- cps[mask]) .^ 2 ./ dcps[j][mask] .^ 2) .+ const_cp
if any(!isfinite, const_cp)
@show any(dcps[j] .< 0)
@show any(!isfinite, dcps[j])
@show any(2π*(dcps[j] .^ 2) .<= 0)
end
# Accumulate into likelihood
ll = ll .+ lnlike_cp .+ lnlike_v2
end
return ll
end
function cvis_bin!(cvis; Δdec, Δra, contrast, u, v)
#u,v: baselines [wavelengths]
#Δdec: dec offset [mas]
#Δra: ra offset [mas]
#contrast: secondary/primary contrast
##################################
#returns complex visibilities of a point source at position Δdec,Δra
l2 = contrast
# phase-factor
cvis .+= l2 .* (
cos.(-2π*(u * Δra + v * Δdec) * π / (180 * 3600 * 1000)) .+
sin.(-2π*(u * Δra + v * Δdec) * π / (180 * 3600 * 1000))im
)
return cvis
end
function closurephase(; vis::AbstractArray, index_cps1::AbstractArray, index_cps2::AbstractArray, index_cps3::AbstractArray)
#vis: complex visibilities
#i_cps1,i_cps2,i_cps3: closure triangle indices
##################################
#returns closure phases [degrees]
realt = real.(vis)
imagt = imag.(vis)
visphi = rad2deg.(atan.(imagt, realt)) #convert to degrees
visphi = mod.(visphi .+ 10980.0, 360.0) .- 180.0 #phase in range (-180,180]
# TODO: replace with rem2pi(, RoundNearest)
cp = @views visphi[index_cps1, :] .+ visphi[index_cps2, :] .- visphi[index_cps3, :]
return cp
end
"""
Extracts indices for calculating closure phases from visibility and closure phase station indices
"""
function cp_indices(; vis2_index::Matrix{<:Int64}, cp_index::Matrix{<:Int64})
i_cps1 = zeros(Int64, size(cp_index)[2])
i_cps2 = zeros(Int64, size(cp_index)[2])
i_cps3 = zeros(Int64, size(cp_index)[2])
nh = maximum(vis2_index) #number of stations making up your interferometer
nb = Int64(nh * (nh - 1) / 2) #number of baselines
ncp = Int64(nh * (nh - 1) * (nh - 2) / 6) #total number of closure phases
for i in range(1, size(cp_index)[2]), j in range(1, size(vis2_index)[2])
if cp_index[1, i] == vis2_index[1, j] && cp_index[2, i] == vis2_index[2, j]
if floor((j - 1) / nb) == floor((i - 1) / ncp)
i_cps1[i] = j
end
end
if cp_index[2, i] == vis2_index[1, j] && cp_index[3, i] == vis2_index[2, j]
if floor((j - 1) / nb) == floor((i - 1) / ncp)
i_cps2[i] = j
end
end
if cp_index[1, i] == vis2_index[1, j] && cp_index[3, i] == vis2_index[2, j]
if floor((j - 1) / nb) == floor((i - 1) / ncp)
i_cps3[i] = j
end
end
end
return i_cps1, i_cps2, i_cps3
end
# Generate new observations for a system of possibly multiple planets
function Octofitter.generate_from_params(like::InterferometryLikelihood, θ_system, orbits::Vector{<:AbstractOrbit})
# # Get epochs, uncertainties, and planet masses from observations and parameters
epochs = like.table.epoch
bands = like.table.band
u = like.table.u
v = like.table.v
cps_data = like.table.cps_data
dcps = like.table.dcps
vis2_data = like.table.vis2_data
dvis2 = like.table.dvis2
i_cps1 = like.table.index_cps1
i_cps2 = like.table.index_cps2
i_cps3 = like.table.index_cps3
use_vis2 = like.table.use_vis2
cp_all = Any[]
vis2_all = Any[]
for j in eachindex(epochs)
band = bands[j]
epoch = epochs[j]
cvis = zeros(Complex{Float64}, length(u[j]))
norm_factor = 0.0 # to normalize complex visibilities
for i in eachindex(orbits)
# All parameters relevant to this planet
θ_planet = θ_system.planets[i]
# orbit object pre-created from above parameters (shared between all likelihood functions)
orbit = orbits[i]
contrast = getproperty(θ_planet, band) #secondary/primary
sol = orbitsolve(orbit, epoch)
Δra = raoff(sol) # in mas
Δdec = decoff(sol) # in mas
# Accumulate into cvis
cvis_bin!(cvis; Δdec, Δra, contrast, u=u[j], v=v[j])
norm_factor += contrast
end
cvis .+= 1.0 #add contribution from the primary
cvis ./= 1.0 + norm_factor
# compute closure phase
cp = closurephase(vis=cvis, index_cps1=i_cps1[j], index_cps2=i_cps2[j], index_cps3=i_cps3[j])
#compute squared visibilities
vis2 = abs.(cvis) .^ 2
cp_all = append!(cp_all, [cp])
vis2_all = append!(vis2_all, [vis2])
end
new_vis_like_table = Table(epoch=epochs, u=u, v=v, cps_data=cp_all, dcps=dcps,
vis2_data=vis2_all, dvis2=dvis2, index_cps1=i_cps1,
index_cps2=i_cps2, index_cps3=i_cps3, band=bands, use_vis2=use_vis2)
# return with same number of rows: band, epoch
# position(s) of point sources according to orbits, θ_system
return InterferometryLikelihood(new_vis_like_table)
end
end