Add option to boats.get_model_spectrum() to interpolate over stellar grid.

[jeroenbouwman]

I would like to request an additional feature to the boats.get_model_spectrum() function. Currently it returns the closest stellar model in the grid to the specified stellar parameters, It would be great if an interpolated model could be returned rater than the closets models. Or perhaps better to have the option to either have the closets or interpolated model returned. This would then be similar as done for the calculation of the limbdarkening coefficients.

[gmorello]

Hi Jeroen, I created a new branch "jb_boats" to test this new functionality. You can add the following line to any boats_example file: star_database_interpolation !nearest seq_linear If not added, the current default is nearest. seq_linear is same scheme adopted for the limb-darkening coefficients. It worked on one example. Let me know if it is ok, and/or if you note any bugs.

[jeroenbouwman]

Hi Giuseppe, how does this work with the get_model_spectrum() function? If I look at the code this still returns the nearest neighbor file in the model grid, or am I mistaken? I run this part of the boats functionality in the following way:

`

    params = [InputParameter['Tstar'] * u.K,
              InputParameter['logg'],
              InputParameter['star_metallicity']]

    model_wavelengths, model_fluxes = \
        boats.get_model_spectrum(InputParameter['stellar_models_grids'],
                                 params=params)

`

will this now return a interpolated stellar model?

[jeroenbouwman]

I see you added the new functionality to the process_configuration_eclipse (_transit) function. Would it be possible to incorporate this into the get_model_spectrum() rather than in the higher level process_configuration_eclipse() function?

[gmorello]

Hi Jeroen, I tried incorporating it in get_model_spectrum, but it was causing much more conflicts to resolve. I think the quickest option to use it without the configuration file is to copy into your own function the procedure implemented in process_configuration_XX, under the elif star_database_interpolation=='seq_linear':

[gmorello]

Let me know if something is not clear, and I can fix it. I may try againg to modify the get_model_spectrum directly, but it would require more work.

[jeroenbouwman]

Hi, that is indeed often a problem. I would have thought that adding an optional keyword star_database_interpolation to the get_model_spectrum() function would solve this as you have the same if else construct in process_configuration_XX as in get_model_spectrum what the model grids are concerned. I can certainly use the solution you created directly in my code. I could also edit the get_model_spectrum function and propose a solution for backwards compatibility which you can then test.

[jeroenbouwman]

Based on your solution, I propose the following substitution:

def get_model_spectrum(models_grid, params=None, file_to_read=None,
                       nearest_warning=None, star_database_interpolation='nearest'):
    """
    This function returns the model spectrum from user file, blackbody calculation or built-in dataset

    :param str models_grid: the choice of stellar_models_grid or planet_models_grid
    :argument quantity array params: default is None
    :argument str file_to_read: user file to read, default is None
    :argument bool nearest_warning: default is None
    :argument str star_database_interpolation: default = 'nearest'
    :return: the model wavelengths (in Angstrom) and the corresponding flux (in erg/(cm^2 s A))
    :rtype: quantity array, quantity array
    """
    if models_grid=='Userfile':
        try:
            model_file = pickle.load(open(file_to_read, 'rb'))
            model_wavelengths = model_file['wavelengths'].to(u.Angstrom)
            if model_file['fluxes'].unit == u.dimensionless_unscaled:
                model_fluxes = model_file['fluxes']
            else:
                model_fluxes = \
                    model_file['fluxes'].to(u.erg / (u.cm**2 * u.second * u.Angstrom),
                                            equivalencies=u.spectral_density(model_wavelengths))
        except:
            print('ERROR: Something went wrong when reading the file', file_to_read)
            exit()            
    elif models_grid=='Blackbody':
        blackbody_temperature = params[0]
        model_wavelengths = get_waves_fromR(1.0, 2000000.0) * u.Angstrom
        model_fluxes = blackbody_lambda(model_wavelengths, blackbody_temperature)
        model_fluxes *= np.pi * u.sr
    else:
        # Reading file names and stellar parameters from the selected database
        [star_files_grid, star_params_grid] = get_stellar_grid_parameters(models_grid) 
        star_effective_temperature =  params[0]
        star_log_gravity =  params[1]
        star_metallicity =  params[2]
        if  star_database_interpolation=='nearest':
            neighbour_index = \
                get_nearest_file_index(star_effective_temperature.value,
                                       star_log_gravity, star_metallicity,
                                       star_params_grid, models_grid)
            if nearest_warning:
                print('WARNING: Adopting nearest model in the', models_grid,
                      'grid: Teff=', star_params_grid[neighbour_index,0]*u.Kelvin,
                      ', logg=', star_params_grid[neighbour_index,1],
                      ', [M/H]=', star_params_grid[neighbour_index,2])
            neighbour_model = \
                databases[models_grid].get_file_content(dbx_file=star_files_grid[neighbour_index])
            model_wavelengths = neighbour_model['wavelengths']
            model_fluxes = neighbour_model['fluxes']
        elif star_database_interpolation=='seq_linear':
            target_name = 'this target'
            neighbour_files_indices = \
                get_neighbour_files_indices(target_name,
                                            star_effective_temperature.value,
                                            star_log_gravity, star_metallicity,
                                            star_params_grid, models_grid)
            neighbour_star_params = np.array([])
            neighbour_star_model_wavelengths = np.array([])
            neighbour_star_model_fluxes = np.array([])
            for i in neighbour_files_indices:
                neighbour_star_params = \
                    my_vstack(neighbour_star_params, star_params_grid[i,:])
                neighbour_model = \
                    databases[models_grid].get_file_content(dbx_file=star_files_grid[i])
                neigh_star_model_wavelengths = neighbour_model['wavelengths']
                neigh_star_model_fluxes = neighbour_model['fluxes']
                neighbour_star_model_wavelengths = \
                    my_vstack(neighbour_star_model_wavelengths,
                              neigh_star_model_wavelengths )
                neighbour_star_model_fluxes = \
                    my_vstack(neighbour_star_model_fluxes,
                              neigh_star_model_fluxes)
            [model_wavelengths, model_fluxes] = \
                interp_model_spectra(star_effective_temperature.value,
                                     star_log_gravity, star_metallicity,
                                     neighbour_star_params,
                                     neighbour_star_model_wavelengths,
                                     neighbour_star_model_fluxes)
    return model_wavelengths, model_fluxes

[gmorello]

Thanks Jeroen, your solution worked and it is now implemented in the default branch. I will wait some time before uploading to Pypi.