Merge pull request #54 from exorad/corrections

Corrections

gcm_toolkit/core/writer.py

@@ -128,7 +128,7 @@ def writer_setup(typ):
         Writer.file_name = None
 
 
-def write_status(tag, message):
+def write_status(tag, message, end='\n'):
     """
     Write a status output that is auto layouted according to
     _writer set up. Careful: an ERROR status will call a
@@ -141,6 +141,8 @@ def write_status(tag, message):
         'STAT', 'INFO', 'E-INFO', 'WARN', 'ERROR'
     message : str
         The message to be printed.
+    end : str, optional
+        end of string, default='\n'
     """
     writer = Writer()
     # define colors and layout according to tags
@@ -166,7 +168,7 @@ def write_status(tag, message):
     if Writer.on:
         # write message to terminal
         if Writer.file_name is None:
-            print(color + liner + Writer.ENDC)
+            print(color + liner + Writer.ENDC, end=end)
 
         # write message to file
         else:

gcm_toolkit/utils/clouds.py

@@ -150,6 +150,7 @@ def eff_func(volume_fraction, wavelength, use_bruggemann=False):
     for wav in range(len_w):
         for hig in range(len_h):
             if use_bruggemann:
+                wrt.write_status('INFO', 'Bruggemann progress: ' + str(round((hig + wav*len_h)/len_h/len_w*100, 1)) + '  %', end='\r')
                 # initial guess using linear approximation
                 eff_0 = np.zeros((2,))
                 eff_0[0] = sum(work[:, hig, wav, 0] * work[:, hig, wav, 1])
@@ -329,6 +330,7 @@ def patch_cloud_mix_opa(self, clouds=None):
         self.line_struc_kappas[:, :, :-1, :] = lsk
         self.line_struc_kappas[:, :, -1, :] = np.ones_like(self.line_struc_kappas[:, :, -1, :]) \
                                               * clouds[0][np.newaxis, :, :]
+        self.continuum_opa_scat_emis += clouds[1]
 
 
 def patch_delete_clouds(self, clouds=None):

gcm_toolkit/utils/interface.py

@@ -344,6 +344,8 @@ def __init__(self, tools, prt):
         super().__init__(tools)
         self.prt = prt
 
+        wrt.write_status("STAT", "Created Interface class for gcm_toolkit and petitRADTRANS")
+
     def set_data(self, time, tag=None, regrid_lowres=False,
                  terminator_avg=False, lat_points=1, lon_resolution=10):
         """
@@ -365,6 +367,17 @@ def set_data(self, time, tag=None, regrid_lowres=False,
             longitudinal opening angle for terminator avareging
         """
 
+        wrt.write_status("STAT", "Selected data set for petitRADTRANS interface")
+        wrt.write_status("INFO", "time: " + str(time))
+        if tag is not None:
+            wrt.write_status("INFO", "tag: " + tag)
+        else:
+            wrt.write_status("INFO", "No tag given")
+
+        if terminator_avg:
+            wrt.write_status("INFO", "Data set ready for transmission spectrum calculation")
+
+
         self._set_data_common(time, tag=tag, regrid_lowres=regrid_lowres,
                               terminator_avg=terminator_avg, lat_points=lat_points,
                               lon_resolution=lon_resolution)
@@ -578,7 +591,8 @@ def calc_transit_spectrum(
             pressure_0=None,
             mass_frac=None,
             clouds=None,
-            use_bruggemann=False
+            use_bruggemann=False,
+            asymmetric=False,
     ):
         """
         Calculate the transit spectrum. This function avarages T-p profiles in
@@ -615,6 +629,8 @@ def calc_transit_spectrum(
         use_bruggemann: bool, optional
             If this flag is set to true, bruggemann mixing is used. This is
             more accurate but takes much longer to calculate.
+        asymmetric : bool, optional
+            If true, the spectra will be calculated for the morning and evening limb speratly.
 
         Returns
         -------
@@ -643,6 +659,22 @@ def calc_transit_spectrum(
         if pressure_0 is None:
             pressure_0 = np.max(self.prt.press)*1e-6
 
+        # output
+        wrt.write_status("STAT", "Calculate transmission spectra")
+        wrt.write_status("INFO", "gravity: " + str(gravity))
+        wrt.write_status("INFO", "rplanet: " + str(rplanet))
+        if clouds is not None:
+            if isinstance(clouds, bool):
+                wrt.write_status("INFO", "Using GCM cloud structure")
+            else:
+                wrt.write_status("INFO", "Using given cloud structure")
+            if use_bruggemann:
+                wrt.write_status("INFO", "Bruggemann mixing used")
+            else:
+                wrt.write_status("INFO", "LLL mixing used")
+
+
+
         # check if clouds are wished
         do_clouds = False
         if clouds is not None:
@@ -660,23 +692,38 @@ def calc_transit_spectrum(
 
         # get profile for each latitude and each terminator
         spectra_list = []
+        output_list = []
+        del_lat = 180/len(self.dsi[c['lat']].values)
         for i, lat in enumerate(self.dsi[c['lat']].values):
             # add morning terminator spectra
             spectra_list.append(self._get_1_transit_spectra([i, lat], [0, -90], mass_frac, gravity,
                                                             mmw, rplanet, pressure_0,
                                                             do_clouds, clouds, use_bruggemann))
+            output_list.append([(-90, (i+0.5)*del_lat - 90), spectra_list[-1]])
 
             # add evening terminator spectra
             spectra_list.append(self._get_1_transit_spectra([i, lat], [1, 90], mass_frac, gravity,
                                                             mmw, rplanet, pressure_0,
                                                             do_clouds, clouds, use_bruggemann))
-
-        # avarage over all profiles (area avaraged)
-        spectra = (np.asarray(spectra_list))**2/len(np.asarray(spectra_list))
-        spectra = np.sqrt(np.sum(spectra, axis=0))
+            output_list.append([(90, (i+0.5)*del_lat - 90), spectra_list[-1]])
 
         # calcualte wavelengths in micron
         wavelengths = 29979245800.0/self.prt.freq/1e-4
+
+        # calcualte spectra either for each limb seperatly or together
+        if not asymmetric:
+            spectra = (np.asarray(spectra_list))**2/len(np.asarray(spectra_list))
+            spectra = np.sqrt(np.sum(spectra, axis=0))
+        else:
+            spectra = np.zeros((2, len(wavelengths)))
+            for pos, spec in output_list:
+                if pos[0] == -90:
+                    spectra[0, :] += (np.asarray(spec))**2/len(np.asarray(spectra_list))*2
+                else:
+                    spectra[1, :] += (np.asarray(spec))**2/len(np.asarray(spectra_list))*2
+
+            spectra = np.sqrt(spectra)
+
 
         # return the final spectra
         return wavelengths, spectra
@@ -794,7 +841,7 @@ def _get_1_transit_spectra(self, lat, lon, mass_frac, gravity, mmw, rplanet,
             if key not in abus:
                 # if a species is missing, check for partial matches
                 for k_abus in abus:
-                    if k_abus in key:
+                    if k_abus+'_' in key:
                         abus[key] = abus[k_abus]
                         break
                 else: