-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathRun_MEDEA.py
More file actions
57 lines (47 loc) · 1.77 KB
/
Run_MEDEA.py
File metadata and controls
57 lines (47 loc) · 1.77 KB
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
"""
File: medea/Run_MEDEA.py
Author: Joshua J. Hibbard
Date: January 2024
Description: Example script opening Cryo-coefficients and the Cryo-basis
for a particular experiment and giving them
to MEDEA for an analytical horizontal dipole suspended above
a PEC. In this case we use the flat horizon Cryo-basis,
and the corresponding Cryo-coefficients for this basis
and the analytical dipole beam.
"""
import numpy as np
from medea import BeamEmulator
import h5py
import os
input_file_path = os.getenv('MEDEA')
cryo_coefficients_path = input_file_path + \
'/input/cryo_coeff_flat_horizon_horizontal_dipole_PEC.hdf5'
cryo_basis_path = input_file_path + \
'/input/cryo_basis_flat_horizon_nside32.hdf5'
with h5py.File(cryo_basis_path, 'r') as cfile:
cryo_basis = cfile['Basis']['Transform_to_map'][()]
unmasked_indices = cfile['Basis']['nonmasked_indices'][()]
frequencies = np.linspace(50,99,50)
hyper_parameter_array = np.round(np.linspace(1,3,21),3) #in meters between 1 and 3
nside = 32
cryo_coefficients = []
with h5py.File(cryo_coefficients_path, 'r') as cofile:
for par in hyper_parameter_array:
par_by_frequency = []
for frequency in frequencies:
par_by_frequency.append(cofile[str(par)]['coefficients']\
['Freq_'+str(int(frequency))][()])
cryo_coefficients.append(par_by_frequency)
cryo_coefficients = np.array(cryo_coefficients)
print(cryo_coefficients.shape)
medea_emulator = BeamEmulator(frequencies, nside, cryo_basis, \
cryo_coefficients, hyper_parameter_array,\
unmasked_indices, interpolation_order=3, use_gpr=False)
print(medea_emulator([1.234]).shape)
try:
import healpy as hp
import matplotlib.pyplot as plot
hp.orthview(medea_emulator([1.234])[0,:])
plot.show()
except:
raise ImportError('Install healpy to view interpolated beam maps!')