diff --git a/gcm_toolkit/core/writer.py b/gcm_toolkit/core/writer.py index c4edf61..32491b7 100644 --- a/gcm_toolkit/core/writer.py +++ b/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: diff --git a/gcm_toolkit/utils/clouds.py b/gcm_toolkit/utils/clouds.py index 227ccf6..71328ce 100644 --- a/gcm_toolkit/utils/clouds.py +++ b/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): diff --git a/gcm_toolkit/utils/interface.py b/gcm_toolkit/utils/interface.py index 875e5cc..f75452b 100644 --- a/gcm_toolkit/utils/interface.py +++ b/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: