Fixed units

Fixed the Jy/pixel conversion.
Added a check kwargs key
Save option to fits file for the maps.
made a tutorial
updated the test we_observe

maria/__init__.py

@@ -97,8 +97,6 @@ def get_pointing(pointing_name, **kwargs):
 def get_site(site_name, **kwargs):
     return Site(**get_site_config(site_name, **kwargs))
 
-
-
 class TOD:
 
     def __init__(self):
@@ -139,6 +137,30 @@ def to_hdf(self, filename):
     def plot(self):
         pass
 
+class KeyNotFoundError(Exception):
+    def __init__(self, invalid_keys):
+        super().__init__(f"The key \'{invalid_keys}\' is not in the database. ")
+
+def check_nested_keys(keys_found, data, keys):
+
+    for key in data.keys():
+        for i in range(len(keys)):
+            if keys[i] in data[key].keys():
+                keys_found[i] = True
+
+def check_json_file_for_key(keys_found, file_path, *keys_to_check):
+    with open(file_path, 'r') as json_file:
+        data = json.load(json_file)       
+        return check_nested_keys(keys_found, data, keys_to_check)
+
+def test_multiple_json_files(files_to_test, *keys_to_find):   
+    keys_found = np.zeros(len(keys_to_find)).astype(bool)
+
+    for file_path in files_to_test:
+        check_json_file_for_key(keys_found, file_path, *keys_to_find)
+
+    if np.sum(keys_found) != len(keys_found):
+        raise KeyNotFoundError(np.array(keys_to_find)[~keys_found])
 
 class Simulation(base.BaseSimulation):
     """
@@ -162,15 +184,21 @@ def __init__(self,
         if isinstance(site, str):
             site = get_site(site, **kwargs)
 
+        json_pattern = os.path.join(here+'/configs', '*.json')
+        json_files    = glob.glob(json_pattern)
+        test_multiple_json_files(json_files, *tuple(kwargs.keys()))
+
         super().__init__(array, pointing, site)
 
+        self.atm_model = atm_model
         if atm_model in ["linear_angular", "LA"]:
             self.atm_sim = atmosphere.LinearAngularSimulation(array, pointing, site, **kwargs)
         elif atm_model in ["kolmogorov_taylor", "KT"]:
             self.atm_sim = atmosphere.KolmogorovTaylorSimulation(array, pointing, site, **kwargs)
         else:
             self.atm_sim = None
 
+        self.map_file = map_file
         if map_file is not None:
             self.map_sim = sky.MapSimulation(array, pointing, site, map_file, **kwargs)
         else:
@@ -192,6 +220,7 @@ def run(self):
         tod.ra   = self.pointing.ra
         tod.dec  = self.pointing.dec
 
+        # number of bands are lost here
         tod.data = np.zeros((self.array.n_det, self.pointing.n_time))
 
         if self.atm_sim is not None:
@@ -207,3 +236,46 @@ def run(self):
                         'altitude': self.site.altitude}
 
         return tod
+
+    def save_maps(self, mapper):
+
+            self.map_sim.map_data["header"]['comment'] = 'Made Synthetic observations via maria code'
+            self.map_sim.map_data["header"]['comment'] = 'Changed resolution and size of the output map'
+            self.map_sim.map_data["header"]['CDELT1']  = np.rad2deg(mapper.resolution)
+            self.map_sim.map_data["header"]['CDELT2']  = np.rad2deg(mapper.resolution)
+            self.map_sim.map_data["header"]['CRPIX1']  = mapper.maps[list(mapper.maps.keys())[0]].shape[0]
+            self.map_sim.map_data["header"]['CRPIX2']  = mapper.maps[list(mapper.maps.keys())[0]].shape[1]
+
+            self.map_sim.map_data["header"]['comment'] = 'Changed pointing location of the output map'
+            self.map_sim.map_data["header"]['CRVAL1']  = self.pointing.ra.mean()
+            self.map_sim.map_data["header"]['CRVAL2']  = self.pointing.dec.mean()
+
+            self.map_sim.map_data["header"]['comment'] = 'Changed spectral position of the output map'
+            self.map_sim.map_data["header"]['CTYPE3']  = 'FREQ    '
+            self.map_sim.map_data["header"]['CUNIT3']  = 'Hz      '
+            self.map_sim.map_data["header"]['CRPIX3']  = 1.000000000000E+00
+
+            self.map_sim.map_data["header"]['BTYPE']   = 'Intensity'
+            if self.map_sim.map_data["inbright"] == 'Jy/pixel': 
+                self.map_sim.map_data["header"]['BUNIT']   = 'Jy/beam '   
+            else: 
+                self.map_sim.map_data["header"]['BUNIT']   = 'Kelvin RJ'   
+
+            for i in range(len(self.array.ubands)):
+
+                self.map_sim.map_data["header"]['CRVAL3']  = self.array.detectors[i][0]
+                self.map_sim.map_data["header"]['CDELT3']  = self.array.detectors[i][1]
+
+                save_map = mapper.maps[list(mapper.maps.keys())[i]] 
+
+                if self.map_sim.map_data["inbright"] == 'Jy/pixel': 
+                    save_map *= utils.KbrightToJyPix(self.map_data["header"]['CRVAL3'], 
+                                                     self.map_data["header"]['CDELT1'], 
+                                                     self.map_data["header"]['CDELT2']
+                                                    )
+                fits.writeto( filename = here+'/outputs/atm_model_'+str(self.atm_model)+'_file_'+ str(self.map_file.split('/')[-1].split('.fits')[0])+'.fits', 
+                                  data = save_map, 
+                                header = self.map_sim.map_data["header"],
+                             overwrite = True 
+                            )
+

maria/cmb_sky.py

@@ -0,0 +1,57 @@
+
+# class CMBSimulation(base.BaseSimulation):
+#     """
+#     This simulates scanning over celestial sources.
+#     """
+#     def __init__(self, array, pointing, site, map_file, add_cmb=False, **kwargs):
+#         super().__init__(array, pointing, site)
+
+#         pass    
+
+#     def _get_CMBPS(self):
+
+#         import camb
+
+#         pars = camb.CAMBparams()
+#         pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.06, omk=0, tau=0.06)
+#         pars.InitPower.set_params(As=2e-9, ns=0.965, r=0)
+
+#         # correct mode would l=129600 for 5"
+#         pars.set_for_lmax(5000, lens_potential_accuracy=0) 
+
+#         results = camb.get_results(pars)
+#         powers = results.get_cmb_power_spectra(pars, CMB_unit="K")["total"][:, 0]
+
+#         self.CMB_PS = np.empty((len(np.unique(self.array.band)), len(powers)))
+#         for i in range(len(np.unique(self.array.band))):
+#             self.CMB_PS[i] = powers
+
+
+#     def _cmb_imager(self, bandnumber=0):
+
+#         import pymaster as nmt
+
+#         nx, ny = self.map_data[0].shape
+#         Lx = nx * np.deg2rad(self.sky_data["incell"])
+#         Ly = ny * np.deg2rad(self.sky_data["incell"])
+
+#         self.CMB_map = nmt.synfast_flat(
+#             nx,
+#             ny,
+#             Lx,
+#             Ly,
+#             np.array([self.CMB_PS[bandnumber]]),
+#             [0],
+#             beam=None,
+#             seed=self.pointing.seed,
+#         )[0]
+
+#         self.CMB_map += utils.Tcmb
+#         self.CMB_map *= utils.KcmbToKbright(np.unique(self.array.band_center)[bandnumber])
+
+#     #self._cmb_imager(i)
+#         #         cmb_data = sp.interpolate.RegularGridInterpolator(
+#         #                     (map_x, map_y), self.CMB_map, bounds_error=False, fill_value=0
+#         #                     )((lam_x, lam_y))
+#         #         self.noise    = self.lam.temperature + cmb_data
+#         #         self.combined = self.map_data + self.lam.temperature + cmb_data

maria/configs/pointings.json

@@ -9,19 +9,22 @@
         "pointing_frame": "az_el", 
         "pointing_units": "degrees", 
         "scan_period": 120, 
-        "sample_rate": 20,  
+        "sample_rate": 50,  
         "seed": 42
     },
 
     "DAISY_2deg": {
         "description": "",
         "start_time": "2022-07-01T08:00:00",
         "scan_pattern": "daisy",  
+        "integration_time": 60,
         "scan_options": {"k": 3.1416},
+        "pointing_center": [4, 10],
         "pointing_throws": [2, 2],
         "pointing_units": "degrees", 
+        "pointing_frame": "ra_dec",
         "scan_period": 60,  
-        "sample_rate": 20, 
+        "sample_rate": 50, 
         "seed": 42
     },
 
@@ -36,7 +39,7 @@
         "pointing_frame": "ra_dec",
         "pointing_units": "degrees", 
         "scan_period": 60,  
-        "sample_rate": 20, 
+        "sample_rate": 50, 
         "seed": 42
     },
 
@@ -51,7 +54,7 @@
         "pointing_frame": "az_el",
         "pointing_units": "degrees", 
         "scan_period": 60,  
-        "sample_rate": 20, 
+        "sample_rate": 50, 
         "seed": 42
     }
 }

maria/configs/sky.json

@@ -0,0 +1,9 @@
+{
+    "Input": {
+        "map_center": [[0.0,0.0]],
+        "map_res": 0.5,
+        "map_inbright": 1.0,
+        "map_units": "KRJ"
+    }
+
+}

maria/mappers.py

@@ -64,8 +64,6 @@ def add_tods(self, tods):
         for tod in np.atleast_1d(tods):
             self.tods.append(self.expand_tod(tod))
 
-
-
 
 class RawBinMapper(BaseMapper):
 
@@ -112,16 +110,4 @@ def run(self):
                 self.map_sums[band] += map_sum
                 self.map_cnts[band] += map_cnt
 
-        self.maps[band] = self.map_sums[band] / self.map_cnts[band]
-
-
-
-
-
-
-
-
-
-
-
-
+        self.maps[band] = self.map_sums[band] / self.map_cnts[band]