/***********************************************************
* Smithsonian Astrophysical Observatory
* Submillimeter Receiver Laboratory
* am
*
* spectra.c                       S. Paine rev. 2023 June 23
*
* Functions for computing various output spectra including
* those derived from the optical depth and radiance spectra.
************************************************************/

#include <math.h>
#include <float.h>
#include <stdlib.h>

#include "abscoeff.h"
#include "am_types.h"
#include "column.h"
#include "errlog.h"
#include "math_const.h"
#include "model.h"
#include "output.h"
#include "phys_const.h"
#include "planck.h"
#include "specfunc.h"
#include "transform.h"
#include "units.h"


/***********************************************************
* int compare_spectral_subgrid_ranges(model_t *model, model_t *lmodel)
*
* Purpose:
*   Checks whether the spectral subgrid ranges being
*   computed have changed from the last model computation
*   to the current one.  This might happen when computing
*   IF spectra if the LO frequency has changed.
*
* Arguments:
*   model_t *model  - pointer to model data structure
*   model_t *lmodel - pointer to a model data structure
*       containing scalar data from a prior computation.
*
* Return:
*   1 if the subgrid range has changed or if lmodel == NULL
*   0 if there have been no changes
************************************************************/

int compare_spectral_subgrid_ranges(model_t *model, model_t *lmodel)
{
    return
        (lmodel == NULL) ||
        (model->isub[0].min != lmodel->isub[0].min) ||
        (model->isub[0].max != lmodel->isub[0].max) ||
        (model->isub[1].min != lmodel->isub[1].min) ||
        (model->isub[1].max != lmodel->isub[1].max);
}   /* compare_spectral_subgrid_ranges() */


/***********************************************************
* void compute_k_out(model_t *model)
*
* Purpose:
*   Copies the first spectral absorption coefficient from
*   the model to the array model.k_out.  A warning is logged
*   if the model contains zero or more than one absorption
*   coefficient.  Appropriate units are set according to
*   whether the absorption coefficient is a molecular one or
*   a binary one.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_k_out(model_t *model)
{
    int lnum, cnum, knum;
    int subgrid;
    int abscoeff_count = 0;

    for (lnum = 0; lnum < model->nlayers; ++lnum) {
        layer_t *layer = model->layer[lnum];
        for (cnum = 0; cnum < layer->ncols; ++cnum) {
            column_t *column = layer->column[cnum];
            for (knum = 0; knum < column->n_abscoeffs; ++knum) {
                abscoeff_t *abscoeff = column->abscoeff[knum];
                if (abscoeff->k == NULL)
                    continue;
                if (abscoeff_count == 0) {
                    output[OUTPUT_K].default_unitnum =
                        k_type[abscoeff->k_typenum].unitnum;
                    /*
                     * If no output unit was set by the user, assign
                     * it here.  Otherwise, if a user unit was
                     * assigned, correct with a warning if
                     * needed.
                     */
                    if (output[OUTPUT_K].unitnum == UNIT_AUTO) {
                        output[OUTPUT_K].unitnum =
                            k_type[abscoeff->k_typenum].unitnum;
                    } else if (
                            unit_tab[output[OUTPUT_K].unitnum].group !=
                            unit_tab[output[OUTPUT_K].default_unitnum].group) {
                        errlog(200,0);
                        output[OUTPUT_K].unitnum =
                            output[OUTPUT_K].default_unitnum;
                    }
                    output[OUTPUT_K].k_typenum = abscoeff->k_typenum;
                    for (subgrid = 0; subgrid <= 1; ++subgrid) {
                        gridsize_t i;
                        gridsize_t imin = model->isub[subgrid].min;
                        gridsize_t imax = model->isub[subgrid].max;
                        for (i = imin; i <= imax; ++i)
                            model->k_out[i] = abscoeff->k[i];
                    }
                }
                ++abscoeff_count;
            }
        }
    }
    /*
     * If the model contains no spectral absorption coefficients,
     * set k_out to zero and log a warning.
     */
    if (abscoeff_count == 0) {
        for (subgrid = 0; subgrid <= 1; ++subgrid) {
            gridsize_t i;
            gridsize_t imin = model->isub[subgrid].min;
            gridsize_t imax = model->isub[subgrid].max;
            for (i = imin; i <= imax; ++i)
                model->k_out[i] = 0.0;
        }
        errlog(198, 0);
    }
    /*
     * If the model contains more than one spectral absorption
     * coefficent, log a warning
     */
    if (abscoeff_count > 1)
        errlog(199, abscoeff_count);
    return;
}   /* compute_k_out() */


/***********************************************************
* int compute_delay_spectrum(model_t *model)
*
* Purpose:
*   Computes the excess delay spectrum, which is divided
*   into two contributions.  Optical transitions contribute
*   an essentially frequency-independent refractivity, which
*   produces a delay proportional to column density only.
*   Lower-lying (mainly rotational) transitions contribute a
*   dispersive delay which depends on all the factors which
*   influence these spectra, including P, T, mixing ratios
*   and column densities.
*
*   In this function, a first pass through the model
*   computes the total optical delay.  In am this is
*   associated with particular h2o and dry_air column types.
*   Then, the dispersive delay associated with all types of
*   absorption which have been included in the model is
*   computed by Hilbert transformation of the opacity
*   spectrum.  Both contributions are then added together
*   to produce the total spectral delay L.
*
* Arguments:
*   model_t *model - pointer to model structure
*
* Return:
*   0 if successful, 1 otherwise
************************************************************/

int compute_delay_spectrum(model_t *model)
{
    gridsize_t i, if0, imax, jr, ji, np2;
    double Lopt;
    double *z;
    int lnum, cnum;

    if (model->L == NULL) {
        errlog(71, 0);
        return 1;
    }
    /*
     * Compute the optical delay.  This is the constant part of the
     * delay spectrum arising from optical transitions.
     */
    Lopt = 0.0;
    for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) {
        for (cnum = 0; cnum < model->layer[lnum]->ncols; ++cnum) {
            Lopt += column_optical_delay(model, lnum, cnum);
        }
    }
    for (lnum = model->path_mid; lnum >= model->path_end; --lnum) {
        for (cnum = 0; cnum < model->layer[lnum]->ncols; ++cnum) {
            Lopt += column_optical_delay(model, lnum, cnum);
        }
    }
    if0 = (gridsize_t)(0.5 + model->f[0] / model->df);
    np2 = 2 * model->nLpad;
    /*
     * Allocate temporary space for a complex array z[0..2*nLpad-1]
     */
    if ((z = (double*)malloc(np2 * sizeof(double))) == NULL) {
        errlog(68, 0);
        return 1;
    }
    /*
     * The imaginary part ni of the refractive index n = nr + i * ni,
     * integrated over the propagation path, is
     *
     *   integral(ni * ds) = c * tau / (4 * pi * f).
     *
     * Put this into the real part of z, antisymmetrically about f = 0.
     * Zero fill the rest of the array.
     */
    for (i = 0; i < if0; ++i) {
        jr = i << 1;
        ji = jr + 1;
        z[jr] = 0.0;
        z[ji] = 0.0;
    }
    jr = if0 << 1;
    ji = jr + 1;
    z[jr] = (model->f[0] < F_EPSILON) ?
        0.0 : (model->tau[0] * C_ON_FOURPI / model->f[0]);
    z[ji] = 0.0;
    for (i = 1; i < model->ngrid; ++i) {
        jr = (i + if0) << 1;
        ji = jr + 1;
        z[jr] = model->tau[i] * C_ON_FOURPI / model->f[i];
        z[ji] = 0.0;
    }
    imax = model->nLpad - if0 - model->ngrid;
    for (i = if0 + model->ngrid; i <= imax; ++i) {
        jr = i << 1;
        ji = jr + 1;
        z[jr] = 0.0;
        z[ji] = 0.0;
    }
    for (i = model->ngrid - 1; i > 0; --i) {
        jr = (model->nLpad - if0 - i) << 1;
        ji = jr + 1;
        z[jr] = -model->tau[i] * C_ON_FOURPI / model->f[i];
        z[ji] = 0.0;
    }
    for (i = model->nLpad - if0; i < model->nLpad; ++i) {
        jr = i << 1;
        ji = jr + 1;
        z[jr] = 0.0;
        z[ji] = 0.0;
    }
    /*
     * Compute the delay spectrum L = integral[(nr-1) * dz], which is
     * a Hilbert transform of the integrated imaginary part.
     */
    hilbert(z, (unsigned long)model->nLpad);
    /*
     * Combine the delay spectrum and the optical delay.
     */
    for (i = 0; i < model->ngrid; ++i) {
        jr = (i + if0) << 1;
        ji = jr + 1;
        model->L[i] = z[ji] + Lopt;
    }
    free(z);
    return 0;
}   /* compute_delay_spectrum() */


/***********************************************************
* void compute_free_space_loss(model_t *model)
*
* Purpose:
*   Computes the natural log of the free space path loss.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_free_space_loss(model_t *model)
{
    if (model->path_begin == 0 || model->path_end == 0) {
        /*
         * If one path endpoint is at infinity, the free space
         * path loss is not meaningful.  Log an error and write
         * zeros into tau_fsl.
         */
        int subgrid;
        errlog(194, 0);
        for (subgrid = 0; subgrid <= 1; ++subgrid) {
            gridsize_t i;
            gridsize_t imin = model->isub[subgrid].min;
            gridsize_t imax = model->isub[subgrid].max;
            for (i = imin; i <= imax; ++i)
                model->tau_fsl[i] = 0.0;
        }
    } else {
        int subgrid;
        double d = source_to_obs_path_distance(model);
        for (subgrid = 0; subgrid <= 1; ++subgrid) {
            gridsize_t i;
            gridsize_t imin = model->isub[subgrid].min;
            gridsize_t imax = model->isub[subgrid].max;
            for (i = imin; i <= imax; ++i) {
                double f   = model->f[i];
                double fsl = FOURPI_ON_C * f * d;
                fsl *= fsl;
                if (fsl < 1.0) {
                    /*
                     * (4 pi d) < lambda.  Clamp fsl at 1
                     * (i.e. tau_fsl = 0), and log a warning
                     * unless f == 0.0.
                     */
                    model->tau_fsl[i] = 0.0;
                    if (f != 0.0)
                        errlog(195, 0);
                } else {
                    model->tau_fsl[i] = log(fsl);
                }
            }
        }
    }
    return;
}   /* compute_free_space_loss() */


/***********************************************************
* int compute_if_spectra(model_t *model)
*
* Purpose:
*   Computes IF spectra from spectra computed on the model
*   frequency grid.
*
* Arguments:
*   model_t *model - pointer to model structure
*
* Return:
*   0 if successful, 1 otherwise
************************************************************/

int compute_if_spectra(model_t *model)
{
    gridsize_t i;

    if (!model->ifmode)
        return 0;
    for (i = 0; i < OUTPUT_END_OF_TABLE; ++i) {
        if ((output[i].flags & OUTPUT_ACTIVE) && (i != OUTPUT_FREQUENCY)) {
            double *s_mod;  /* model spectrum, on model frequency grid */
            double *s_if;   /* if spectrum, on if frequency grid       */
            gridsize_t j;
            s_mod = output[i].spectrum;
            s_if = model->fif;  /* use fif array for scratch space */
            if (model->ifmode & IFMODE_DSB) {
                double su, sl, wtu, wtl;
                gridsize_t ju, jl;
                wtu = model->dsb_utol_ratio / (1.0 + model->dsb_utol_ratio);
                wtl = 1.0 - wtu;
                for (j = 0,
                    ju = model->iusb.min, jl = model->ilsb.max;
                    j < model->nif;
                    ++j, ++ju, --jl) {
                    su = s_mod[ju]
                        + model->fif_delta * (s_mod[ju+1] - s_mod[ju]);
                    sl = s_mod[jl]
                        + model->fif_delta * (s_mod[jl+1] - s_mod[jl]);
                    s_if[j] = wtu * su + wtl * sl;
                }
            } else if (model->ifmode & IFMODE_USB) {
                gridsize_t ju;
                for (j = 0, ju = model->iusb.min; j < model->nif; ++j, ++ju) {
                    s_if[j] = s_mod[ju]
                        + model->fif_delta * (s_mod[ju+1] - s_mod[ju]);
                }
            } else if (model->ifmode & IFMODE_LSB) {
                gridsize_t jl;
                for (j = 0, jl = model->ilsb.max; j < model->nif; ++j, --jl) {
                    s_if[j] = s_mod[jl]
                        + model->fif_delta * (s_mod[jl+1] - s_mod[jl]);
                }
            } else {
                errlog(86, __LINE__);
                return 1;
            }
            /* copy IF spectrum array back to model spectrum array */
            for (j = 0; j < model->nif; ++j)
                s_mod[j] = s_if[j];
        }
    }
    /* fill in the fif array */
    for (i = 0; i < model->nif; ++i)
        model->fif[i] = model->fif_0 + i * model->df;
    return 0;
}   /* compute_if_spectra() */


/***********************************************************
* void compute_radiance_difference_spectrum(model_t *model, model_t *lmodel)
*
* Purpose:
*   Computes the radiance difference spectrum, defined as
*  
*    I_diff = I - I_ref
*  
*   where
*  
*    I_ref = B(T_ref)
*  
*   I_diff models the output of spectrometers which
*   effectively measure the difference in radiance between
*   the scene and a reference blackbody load at T = T_ref.
*   Examples are instruments that chop between reference
*   load and scene, and rapid-scan Fourier spectrometers
*   that have one input port of the interferometer
*   terminated by a reference load.
*
* Arguments:
*   model_t *model - pointer to model structure
*   model_t *lmodel - pointer to a model data structure
*       containing scalar data from a prior computation.
************************************************************/

void compute_radiance_difference_spectrum(model_t *model, model_t *lmodel)
{
    int subgrid;

    /*
     * Compute I_ref if needed.
     */
    if (lmodel == NULL ||
        fabs(model->Tref - lmodel->Tref) > (DBL_EPSILON * model->Tref) ||
        compare_spectral_subgrid_ranges(model, lmodel)) {
        for (subgrid = 0; subgrid <= 1; ++subgrid) {
            gridsize_t imin = model->isub[subgrid].min;
            gridsize_t imax = model->isub[subgrid].max;
            B_Planck(
                &model->I_ref[imin],
                &model->f[imin],
                model->df,
                imax - imin + 1,
                model->Tref);
        }
    }
    for (subgrid = 0; subgrid <= 1; ++subgrid) {
        gridsize_t i;
        gridsize_t imin = model->isub[subgrid].min;
        gridsize_t imax = model->isub[subgrid].max;
        for (i = imin; i <= imax; ++i)
            model->I_diff[i] = model->I[i] - model->I_ref[i];
    }
    return;
}   /* compute_radiance_difference_spectrum() */


/***********************************************************
* void compute_spectral_Tsys(model_t *model)
*
* Purpose:
*   Computes a system temperature spectrum.  This is a
*   shifted and scaled version of the Rayleigh-Jeans
*   temperature spectrum, defined as
*
*     Tsys[i] = (Trj[i] + Trx) * rx_gain_factor
*
*   where rx_gain_factor is intended to be used for modeling
*   fluctuations in receiver gain when fitting series of
*   spectra.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_spectral_Tsys(model_t *model)
{
    int subgrid;
    for (subgrid = 0; subgrid <= 1; ++subgrid) {
        gridsize_t i;
        gridsize_t imin = model->isub[subgrid].min;
        gridsize_t imax = model->isub[subgrid].max;
        for (i = imin; i <= imax; ++i)
            model->Tsys[i] =
                (model->Trj[i] + model->Trx) * model->rx_gain_factor;
    }
    return;
}   /* compute_spectral_Tsys() */


/***********************************************************
* void compute_spectral_Y_factor(model_t *model)
*
* Purpose:
*   Compute the spectral Y-factor, defined as
*
*     Y[i] = {(Trj[i] + Trx) / [Trj(Tref) + Trx]} * rx_gain_factor
*
*   where Trj(Tref) is a Rayleigh-Jeans brightness
*   temperature computed for a physical reference
*   temperature Tref.  Note that if Y is convolved with an
*   ILS, we are making the implicit assumption here that the
*   frequency dependence of Trj(Tref) is negligible across
*   the width of the ILS.  The receiver noise temperature
*   Trx is assumed to be already expressed in terms of an
*   equivalent Rayleigh-Jeans source brightness temperature.
*
*   As in the case of Tsys, rx_gain_factor is intended to be
*   used for modeling fluctuations in receiver gain when
*   fitting series of spectra.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_spectral_Y_factor(model_t *model)
{
    /*
     * If the denominator in the Y-factor is too small, log an error
     * and set the spectrum to zero.
     */
    if ((model->Tref + model->Trx) < DBL_EPSILON) {
        int subgrid;
        errlog(110, 0);
        for (subgrid = 0; subgrid <= 1; ++subgrid) {
            gridsize_t i;
            gridsize_t imin = model->isub[subgrid].min;
            gridsize_t imax = model->isub[subgrid].max;
            for (i = imin; i <= imax; ++i)
                model->Y[i] = 0.0;
        }
    } else {
        int subgrid;
        for (subgrid = 0; subgrid <= 1; ++subgrid) {
            gridsize_t i;
            gridsize_t imin = model->isub[subgrid].min;
            gridsize_t imax = model->isub[subgrid].max;
            for (i = imin; i <= imax; ++i)
                model->Y[i] = model->rx_gain_factor *
                    (model->Trj[i] + model->Trx) /
                    (T_Rayleigh_Jeans(model->f[i], model->Tref) + model->Trx);
        }
    }
    return;
}   /* compute_spectral_Y_factor() */


/***********************************************************
* void compute_Tb(model_t *model)
*
* Purpose:
*   Compute Planck brightness temperature spectrum from the
*   spectral radiance.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_Tb(model_t *model)
{
    int subgrid;

    for (subgrid = 0; subgrid <= 1; ++subgrid) {
        gridsize_t i;
        gridsize_t imin = model->isub[subgrid].min;
        gridsize_t imax = model->isub[subgrid].max;
        /*
         * The first point is handled specially, in case the frequency
         * grid starts at zero.  lim Tb, f->0 = T0, so set Tb(f=0) = T0.
         * This prevents ringing if Tb gets convolved with an ILS.
         */
        if (imin == 0) {
            model->Tb[0] = model->f[0] < F_EPSILON ?
                model->T0 : T_Planck(model->f[0], model->I[0]);
            ++imin;
        }
        #ifdef _OPENMP 
        #pragma omp parallel for schedule(guided, 128)\
            if (model->isub[0].max - model->isub[0].min >= 16384)
        #endif
        for (i = imin; i <= imax; ++i)
            model->Tb[i] = T_Planck(model->f[i], model->I[i]);
    }
    return;
}   /* compute_Tb() */


/***********************************************************
* void compute_transmittance(model_t *model)
*
* Purpose:
*   Compute spectral transmittance from opacity.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_transmittance(model_t *model)
{
    int subgrid;

    for (subgrid = 0; subgrid <= 1; ++subgrid) {
        gridsize_t i;
        gridsize_t imin = model->isub[subgrid].min;
        gridsize_t imax = model->isub[subgrid].max;
        #ifdef _OPENMP
        #pragma omp parallel for schedule(guided, 128)\
            if (imax - imin >= 16384)
        #endif
        for (i = imin; i <=imax; ++i)
            model->tx[i] = am_exp(-model->tau[i]);
    }
    return;
}   /* compute_transmittance() */


/***********************************************************
* void compute_Trj(model_t *model)
*
* Purpose:
*   Compute Rayleigh-Jeans brightness temperature spectrum
*   from the spectral radiance.
*
* Arguments:
*   model_t *model - pointer to model structure
************************************************************/

void compute_Trj(model_t *model)
{
    int subgrid;

    for (subgrid = 0; subgrid <= 1; ++subgrid) {
        gridsize_t i;
        gridsize_t imin = model->isub[subgrid].min;
        gridsize_t imax = model->isub[subgrid].max;
        /*
         * The first point is handled specially, in case the frequency
         * grid starts at zero.  lim Trj, f->0 = T0, so set Tb(f=0) = T0.
         * This prevents ringing if Trj gets convolved with an ILS.
         */
        if (imin == 0) {
            model->Trj[0] = model->f[0] < F_EPSILON ?
                model->T0 : model->I[0] / (TWOKB_ON_CSQUARED * model->f2[0]);
            ++imin;
        }
        /*
         * This loop was found not to gain significantly from being
         * parallelized.
         */
        for (i = imin; i <= imax; ++i)
            model->Trj[i] = model->I[i] / (TWOKB_ON_CSQUARED * model->f2[i]);
    }
    return;
}   /* compute_Trj() */


/***********************************************************
* static int set_IF_spectrum_subgrid_ranges(model_t *model)
*
* Purpose:
*
*   If IF spectra are being computed, this function sets the
*   ranges of grid points used to compute LSB and/or USB
*   spectra.
*
*   If no IF range has been specified, the IF grid range is
*   set to the maximum range supported by the model grid.
*   Otherwise, it is set to the intersection of the
*   specified IF grid and the range supported by the model
*   grid.
*
*   The model frequency grid is always aligned to f = 0,
*   whereas the IF frequency grid will be aligned to
*   fif = 0.  That is, the origin of the IF frequency grid
*   is flo, with the same grid interval df as the model
*   grid:
*
*                                           f = flo
*                ->|       |<- df             |
*   f   <- |       |       |       |       |       |       | ->
*   fif <-    |       |       |       |       |       |      ->
*                ->|  |<- df * fif_delta      |
*                                           fif = 0
*
*   Here, positive (USB) and negative (LSB) IF frequencies
*   are shown.  However, the IF frequency grid is always
*   output as positive.
*
*   Besides the grid range, this function sets three
*   additional parameters in the model structure:
*
*     model->fif_delta - normalized offset between a point
*       on the IF grid, and the model grid point at or below
*       the IF grid point.
*     model->fif_0 -  the lowest IF frequency
*     model->nif - the number of IF frequency grid points.
*
* Arguments:
*   model_t *model - pointer to model structure
*
* Return:
*   0 if successful, 1 otherwise
************************************************************/

int set_IF_spectrum_subgrid_ranges(model_t *model)
{
    double lsb_fmin = 0.0, lsb_fmax = 0.0;
    double usb_fmin = 0.0, usb_fmax = 0.0;
    int fif_defined = (model->fif_max >= 0.0);

    /*
     * At least two grid points are needed for IF spectra
     */
    if (model->ngrid < 2) {
        errlog(88, 0);
        return 1;
    }
    /*
     * LSB frequency range
     */
    if (model->ifmode & (IFMODE_LSB | IFMODE_DSB)) {
        lsb_fmax = model->flo - model->f[0];
        if (lsb_fmax < 0.0) {
            errlog(85, 0);
            return 1;
        }
        if (fif_defined && (lsb_fmax > model->fif_max * (1.0 + DBL_EPSILON)))
            lsb_fmax = model->fif_max * (1.0 + DBL_EPSILON);
        lsb_fmin = model->flo - model->f[model->ngrid - 1];
        if (lsb_fmin < 0.0) 
            lsb_fmin = 0.0;
        if (fif_defined && (lsb_fmin < model->fif_min)) {
            if (model->fif_min > lsb_fmax) {
                errlog(85, 0);
                return 1;
            }
            lsb_fmin = model->fif_min;
        }
    }
    /*
     * USB frequency range.
     */
    if (model->ifmode & (IFMODE_USB | IFMODE_DSB)) {
        usb_fmax = model->f[model->ngrid - 1] - model->flo;
        if (usb_fmax < 0.0) {
            errlog(85, 0);
            return 1;
        }
        if (fif_defined && (usb_fmax > model->fif_max * (1.0 + DBL_EPSILON)))
            usb_fmax = model->fif_max * (1.0 + DBL_EPSILON);
        usb_fmin = model->f[0] - model->flo;
        if (usb_fmin < 0.0)
            usb_fmin = 0.0;
        if (fif_defined && (usb_fmin < model->fif_min)) {
            if (model->fif_min > usb_fmax) {
                errlog(85, 0);
                return 1;
            }
            usb_fmin = model->fif_min;
        }
    }
    /*
     * For DSB, find the overlap between the LSB and USB ranges.
     */
    if (model->ifmode & (IFMODE_DSB)) {
        if (lsb_fmax < usb_fmax)
            usb_fmax = lsb_fmax;
        else
            lsb_fmax = usb_fmax;
        if (lsb_fmin > usb_fmin)
            usb_fmin = lsb_fmin;
        else
            lsb_fmin = usb_fmin;
    }
    /*
     * Model grid indices covering the LSB and USB spectra.  The
     * ranges ilsb and iusb cover the model grid points just below
     * each corresponding IF grid point.  Note that, in principle,
     * there could be an index rounding inconsistency, in the DSB
     * case, between LSB and USB, but this would occur near
     * fif_delta = 0.0, so the resulting interpolation error would be
     * negligible.
     */
    if (model->ifmode & (IFMODE_LSB | IFMODE_DSB)) {
        double f;
        model->fif_0 =
            model->df * ceil((lsb_fmin * (1.0 - DBL_EPSILON)) / model->df);
        model->nif = (gridsize_t)ceil(((lsb_fmax - model->fif_0) *
            (1.0 - DBL_EPSILON)) / model->df);
        f = model->flo - model->fif_0;
        model->ilsb.max = (gridsize_t)floor(((f - model->f[0]) *
            (1.0 + DBL_EPSILON)) / model->df);
        model->fif_delta = (f - model->f[model->ilsb.max]) / model->df;
        model->ilsb.min = model->ilsb.max - (model->nif - 1);
    }
    if (model->ifmode & (IFMODE_USB | IFMODE_DSB)) {
        double f;
        model->fif_0 = model->df *
            ceil((usb_fmin * (1.0 - DBL_EPSILON)) / model->df);
        model->nif = (gridsize_t)ceil(((usb_fmax - model->fif_0) *
            (1.0 - DBL_EPSILON)) / model->df);
        f = model->flo + model->fif_0;
        model->iusb.min = (gridsize_t)floor(((f - model->f[0]) *
            (1.0 + DBL_EPSILON)) / model->df);
        model->fif_delta = (f - model->f[model->iusb.min]) / model->df;
        model->iusb.max = model->iusb.min + (model->nif - 1);
    }
    return 0;
}   /* set_IF_spectrum_subgrid_ranges() */


/***********************************************************
* int set_spectral_subgrid_ranges(model_t *model, model_t *lmodel)
*
* Purpose:
*
*   This function sets ranges for two subgrids of the model
*   frequency grid.  For normal, LSB, and USB spectra, one
*   of these subgrid ranges is empty.  For DSB spectra, if
*   the USB and LSB ranges do not overlap, the two subgrid
*   ranges correspond to USB and LSB.  If they do overlap,
*   they are consolidated into a single subgrid.
*
*   Note that the IF sideband ranges computed in
*   set_IF_spectrum_subgrid_ranges() are the range of model
*   grid indices for the model grid points just at or below
*   each corresponding IF frequency point, so the ranges set
*   here include an extra point at the top for
*   interpolation.
*
*   These abbreviated sub ranges are used to accelerate
*   computations at the layer and column levels, which
*   typically occur repetitively during fits or Jacobian
*   computations.  (This optimization is most important for
*   DSB spectra.)
*
*   If an ILS is being applied, or if a delay spectrum is
*   being computed, the model will always be computed over
*   the entire frequency grid; in these cases any difference
*   between the model and IF grid ranges is assumed to be
*   spectral padding desired by the user.
*
* Arguments:
*   model_t *model - pointer to model structure
*
* Return:
*   0 if successful, 1 otherwise
************************************************************/

int set_spectral_subgrid_ranges(model_t *model, model_t *lmodel)
{
    /*
     * For IF spectra, initialize or adjust subgrid ranges as needed.
     */
    if (model->ifmode) {
        if ((lmodel == NULL) ||
            (fabs(model->flo - lmodel->flo)
                > DBL_EPSILON * model->flo) ||
            (fabs(model->fif_min - lmodel->fif_min)
                > DBL_EPSILON * model->fif_min) ||
            (fabs(model->fif_max - lmodel->fif_max)
                > DBL_EPSILON * model->fif_max)) {
            if (set_IF_spectrum_subgrid_ranges(model))
                return 1;
        }
    }
    if ((model->ifmode == 0) || (model->ils != NULL) || (model->L != NULL)) {
        model->isub[0].min = 0;
        model->isub[0].max = model->ngrid - 1;
        model->isub[1].min = 0;
        model->isub[1].max = -1;
    } else if (model->ifmode & IFMODE_LSB) {
        model->isub[0].min = model->ilsb.min;
        model->isub[0].max = model->ilsb.max + 1;
        model->isub[1].min = 0;
        model->isub[1].max = -1;
    } else if (model->ifmode & IFMODE_USB) {
        model->isub[0].min = model->iusb.min;
        model->isub[0].max = model->iusb.max + 1;
        model->isub[1].min = 0;
        model->isub[1].max = -1;
    } else if (model->ifmode & IFMODE_DSB) {
        if (model->ilsb.max + 1 >= model->iusb.min) {
            model->isub[0].min = model->ilsb.min;
            model->isub[0].max = model->iusb.max + 1;
            model->isub[1].min = 0;
            model->isub[1].max = -1;
        } else {
            model->isub[0].min = model->ilsb.min;
            model->isub[0].max = model->ilsb.max + 1;
            model->isub[1].min = model->iusb.min;
            model->isub[1].max = model->iusb.max + 1;
        }
    } else {
        errlog(86, __LINE__);
        return 1;
    }
    return 0;
}   /* set_spectral_subgrid_ranges() */