/*********************************************************** * Smithsonian Astrophysical Observatory * Submillimeter Receiver Laboratory * am * * model.c S. Paine rev. 2024 July 16 * * Atmospheric model and radiative transfer computations. ************************************************************/ #include #include #include "am_sysdep.h" #include "am_types.h" #include "column.h" #include "errlog.h" #include "ils.h" #include "interp.h" #include "layer.h" #include "math_const.h" #include "model.h" #include "output.h" #include "phys_const.h" #include "rt.h" #include "spectra.h" #include "units.h" /* * The following two constants control iterative solution * for the tangent height of paths below the astronomical * horizon. Convergence is deemed to have occurred when the * absolute value of the fractional change in ztan from one * iteration to the next is less than ZTAN_TOL. * ZTAN_ITER_MAX is the maximum number of iterations * allowed, after which convergence is deemed to have * failed. */ #define ZTAN_TOL 1e-10 #define ZTAN_ITER_MAX 100 /* * If the sine of the zenith angle computed from the * refractive invariant at the base of a layer exceeds * unity, then we're in a duct on the layer below. However, * especially near a tangent point, there may be a small * numerical overshoot. If sin(za) exceeds unity by less * than SIN_ZA_DUCT_TOL, it is simply clamped to unity. * * This tolerance can be somewhat smaller than ZTAN_TOL, * above, since ZTAN_TOL enters the refractive invariant via * the refractivity n-1. */ #define SIN_ZA_DUCT_TOL 1e-12 /* * For thin layers, TRIG_TOL sets the point at which path * trigonometry computations cannot be done to useful * accuracy at double precision, and a thin layer limit * is used instead. */ #define TRIG_TOL 1e-12 /* * Constants for set_hydrostatic_pressure_limits() * * The product of TOA_EPSILON and the base pressure of the * highest-altitude model definition layer is considered the * "top of atmosphere" pressure level. A user-defined * source, observing, or tangent level placed at a lower * pressure will generate a warning, but will be * accommodated. A user-defined level set at a height that * would correspond to a lower pressure (such as the orbital * height zobs of a satellite) is placed at P = 0. When * such a level is inserted into the model, it is sorted * into height order with any other such levels. * * The product of P_OVERSHOOT_LIMIT and the base pressure of * the lowest model definition layer is the maximum allowed * pressure for any extrapolated level. */ #ifndef TOA_EPSILON #define TOA_EPSILON 0.1353352832366127 #endif #ifndef P_OVERSHOOT_LIMIT #define P_OVERSHOOT_LIMIT 2.718281828459045 #endif static int compute_airmass_and_refraction(model_t*); static int compute_refractive_invariant(model_t*); static int copy_model_scalars(model_t*, model_t*); static int find_n_from_z(double*, const model_t*, const double); static int find_P_from_z(double*, const model_t*, const double); static int find_z_from_P(double*, const model_t*, const double); static void set_hydrostatic_pressure_limits(model_t*); static int set_user_defined_levels(model_t*, model_t*); static int set_propagation_path_indices(model_t*); static int solve_for_tan_level(model_t*, model_t*); /*********************************************************** * static int compute_airmass_and_refraction(model_t *model) * * Purpose: * Compute airmass and refraction for the propagation path. * * An excellent and highly readable general reference on * this topic is the review: * * A. T. Young 2006, "Understanding astronomical * refraction." The Observatory 126:82-115. * * For spherical ray tracing, we use the approximations * from: * * S. N. Kivalov 2007, "Improved ray tracing air mass * numbers model." Applied Optics 46:7091-7098. * * The Kivalov formulation has the very useful property * that the path segment in a layer is computed using the * mean of the radii and zenith angles at the upper and * lower layer boundaries. This enables a midpoint * integration of the path as a function of height, * avoiding singularities at the tangent level. * * For plane-parallel geometry, a similar approximation is * used for the airmass, which is computed as the secant of * the mean of the zenith angles at the upper and lower * layer boundaries. * * Arguments: * model_t *model - pointer to model data structure * * Return: * 1 if there were errors, 0 otherwise. ************************************************************/ static int compute_airmass_and_refraction(model_t *model) { if (model->geometry & GEOMETRY_PLANE_PARALLEL) { double sin_za, za_hi; int lnum; sin_za = model->p; if (sin_za > 1.0) { /* duct */ errlog(172, 0); return 1; } za_hi = asin(sin_za); for (lnum = 0; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; double za_lo, za_bar; sin_za = model->p / layer->nbase; if (sin_za > 1.0 + SIN_ZA_DUCT_TOL) { /* duct */ errlog(172, lnum); return 1; } sin_za = sin_za > 1.0 ? 1.0 : sin_za; za_lo = layer->za_base = asin(sin_za); za_bar = 0.5 * (za_hi + za_lo); za_hi = za_lo; /* * An instrument or plume layer on a hydrostatic * model is assigned a fixed airmass factor of unity. */ if ((model->PTmode & PTMODE_HYDROSTATIC) && (layer->h > 0.0)) layer->m = 1.0; else layer->m = 1.0 / cos(za_bar); } } else if (model->geometry & (GEOMETRY_SPHERICAL | GEOMETRY_LIMB)) { double z_hi = 0.0; /* inits to suppress compiler warning */ double r_hi = 0.0; double za_hi = 0.0; int lnum; /* * Here, scan down through the layer stack to find the * highest layer with its base at finite height and upper * boundary at infinity. The airmass of this layer is * asymptotically equal to unity. The central angle dphi * is asymptotically equal to the zenith angle at the * layer base; in the total refraction computation these * will cancel, and there will be no net refractive * bending in the top layer of spherical models. */ for (lnum = 0; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; if (layer->Pbase > 0.0 || layer->type != LAYER_TYPE_DEF) { double sin_za; z_hi = layer->zbase; r_hi = model->R0 + z_hi; sin_za = model->p / (layer->nbase * r_hi); if (sin_za > 1.0 + SIN_ZA_DUCT_TOL) { /* duct */ errlog(172, lnum); return 1; } sin_za = sin_za > 1.0 ? 1.0 : sin_za; za_hi = asin(sin_za); layer->m = 1.0; layer->dphi = za_hi; layer->za_base = za_hi; break; } /* * Dummy layer at infinity. */ layer->m = 1.0; layer->dphi = 0.0; layer->za_base = 0.0; } /* * Log a warning if the entire path lies within a single * layer with upper bound at infinity (P == 0), since the * user might not have expected this. */ if (lnum == model->path_mid) errlog(169, 0); /* * Now work down through the remaining layers. */ for (++lnum; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; double z_lo, dz; double r_lo; double sin_za; double za_lo, za_bar; /* * A layer with h > 0 in a hydrostatic model * (spherical models are always hydrostatic) is an * embedded non-hydrostatic instrument or plume * layer. It is treated as a zero-thickness layer * having no effect on the path geometry. The * airmass factor is set to unity. */ if (layer->h > 0.0) { layer->m = 1.0; layer->dphi = 0.0; layer->za_base = model->layer[lnum - 1]->za_base; continue; } z_lo = layer->zbase; dz = z_hi - z_lo; r_lo = model->R0 + z_lo; sin_za = model->p / (layer->nbase * r_lo); if (sin_za > 1.0 + SIN_ZA_DUCT_TOL) { /* duct */ errlog(172, lnum); return 1; } sin_za = sin_za > 1.0 ? 1.0 : sin_za; za_lo = asin(sin_za); za_bar = 0.5 * (za_lo + za_hi); /* * Thin layer limits are used when the computation * cannot be done at double precision. The airmass * is undefined on a zero-thickness layer, and is * simply set to zero. */ if ((PI_ON_TWO - za_bar < TRIG_TOL) || (dz < r_lo * TRIG_TOL)) { layer->m = 0.0; layer->dphi = 0.0; layer->za_base = model->layer[lnum - 1]->za_base; continue; } /* * For a normal layer, compute dphi and m using * Kivalov Eq.(14) and Eq.(16). */ layer->dphi = 2.0 * atan(dz * tan(za_bar) / (r_lo + r_hi)); layer->m = sqrt(r_lo * r_lo + r_hi * r_hi - 2.0 * r_lo * r_hi * cos(layer->dphi)) / dz; layer->za_base = za_lo; z_hi = z_lo; r_hi = r_lo; za_hi = za_lo; } } else { return 1; /* model geometry flags were invalid */ } return 0; } /* compute_airmass_and_refraction() */ /*********************************************************** * int compute_model(model_t *model, model_t *lmodel_in) * * Purpose: * This is the top-level function for setting up the * atmospheric model and running radiative transfer * computations. * * If lmodel_in != NULL, and lmodel_in->ngrid is nonzero, * then this is an update of a prior model computation, and * *lmodel_in holds the values of scalar model variables at * the time of that prior computation. In this case, * expensive spectral computations are only redone if they * depend on changed scalar variables. (Scalar variables * are always recomputed, because these computations take * relatively negligible time.) * * After all computations are finished, scalar variables in * *lmodel_in are updated with the corresponding values * from *model to save the most recently computed model * state. * * Arguments: * model_t *model - pointer to model data structure * model_t *lmodel_in - pointer to a model data structure * holding scalar data from a prior * computation, or NULL for a single- * pass computation. * * Return: * 1 if there were errors, 0 otherwise. ************************************************************/ int compute_model(model_t *model, model_t *lmodel_in) { model_t *lmodel; double tmark = 0.0; if (model == NULL) { errlog(76, 0); return 1; } /* * *lmodel_in is ignored for scalar variable dependency * checking if it does not yet reflect a completed model * computation (i.e., if ngrid == 0 for non-NULL lmodel_in). * This is done by setting lmodel to NULL. Otherwise, * lmodel is just a copy of lmodel_in. */ lmodel = (lmodel_in != NULL && lmodel_in->ngrid == 0) ? NULL : lmodel_in; /* * Start the model run performance timer. */ if (model->log_runtimes) model->runtime = tmark = am_timer(0.0); /* * Set up the atmospheric model. This step may insert or * reorder model layers to support interpolated observing, * source, or tangent levels. Because of this, we operate on * lmodel_in also, even though it is not being used here for * dependency checking. */ if (setup_atmospheric_model(model, lmodel_in)) return 1; /* * Log the model setup timing. Start the layer performance * timers now that all layers are in place. */ if (model->log_runtimes) { tmark += (model->am_runtime = am_timer(tmark)); reset_layer_performance_timers(model); } /* * In some cases (IF spectra, particularly), we don't need to * compute all spectral channels. Here, the sub-grids to be * computed are defined (or updated). */ if (set_spectral_subgrid_ranges(model, lmodel)) return 1; /* * Compute the optical depth (opacity) spectrum for the * entire propagation path. */ if (compute_opacity_spectrum(model, lmodel)) return 1; if (model->log_runtimes) tmark += (model->od_runtime = am_timer(tmark)); /* * Do the radiative transfer computation along the * propagation path to compute the spectral radiance. */ compute_spectral_radiance(model, lmodel); if (model->log_runtimes) tmark += (model->rt_runtime = am_timer(tmark)); /* * Compute other spectra derived from radiance and opacity as * needed. */ if (output[OUTPUT_TRANSMITTANCE].flags & OUTPUT_ACTIVE) compute_transmittance(model); if (output[OUTPUT_TB_PLANCK].flags & OUTPUT_ACTIVE) compute_Tb(model); if (output[OUTPUT_TB_RAYLEIGH_JEANS].flags & OUTPUT_ACTIVE || output[OUTPUT_Y].flags & OUTPUT_ACTIVE || output[OUTPUT_TSYS].flags & OUTPUT_ACTIVE) compute_Trj(model); if (output[OUTPUT_TSYS].flags & OUTPUT_ACTIVE) compute_spectral_Tsys(model); if (output[OUTPUT_Y].flags & OUTPUT_ACTIVE) compute_spectral_Y_factor(model); if (output[OUTPUT_RADIANCE_DIFF].flags & OUTPUT_ACTIVE) compute_radiance_difference_spectrum(model, lmodel); if (output[OUTPUT_DELAY].flags & OUTPUT_ACTIVE) { if (compute_delay_spectrum(model)) return 1; } /* * Compute other output spectra as needed. */ if (output[OUTPUT_FREE_SPACE_LOSS].flags & OUTPUT_ACTIVE) compute_free_space_loss(model); if (output[OUTPUT_K].flags & OUTPUT_ACTIVE) compute_k_out(model); /* * Apply modifications to spectra that model the response of * spectrometers or receivers. */ if (model->ils != NULL) { if (apply_instrumental_line_shape(model, lmodel)) return 1; } if (model->ifmode) { compute_if_spectra(model); } if (model->log_runtimes) tmark += (model->spec_runtime = am_timer(tmark)); /* * Finished. Log the total model compute time and update * *lmodel_in. */ if (model->log_runtimes) model->runtime = am_timer(model->runtime); if (lmodel_in != NULL) return copy_model_scalars(model, lmodel_in); else return 0; } /* compute_model() */ /*********************************************************** * static int compute_refractive_invariant(model_t *model) * * Purpose: * Computes the refractive invariant, which is n * sin(za) * in plane-parallel mode, and n * r * sin(za) in spherical * mode. This function assumes that layer heights have * already been computed. * * Arguments: * model_t *model - pointer to model structure * * Return: * 0 if successful, 1 otherwise ************************************************************/ static int compute_refractive_invariant(model_t *model) { int lnum; /* * Find the index number of the observing level for * plane-parallel or spherical models, or of the tangent * level for limb models. The default observing or tangent * level if not otherwise defined is the base of the lowest * defined model layer. */ if (model->geometry & GEOMETRY_LIMB) lnum = get_tan_lnum(model); else lnum = get_obs_lnum(model); /* * Compute the refractive invariant. */ if (model->geometry & GEOMETRY_PLANE_PARALLEL) { model->p = model->layer[lnum]->nbase * sin(model->za); } else if (model->geometry & GEOMETRY_SPHERICAL) { double r = model->R0 + model->layer[lnum]->zbase; double sin_za = sin(model->za); /* * Here, clean up numerical rounding error for the * common nadir case, za == PI. */ if (fabs(sin_za) < DBL_EPSILON) sin_za = 0.0; model->p = model->layer[lnum]->nbase * r * sin_za; } else if (model->geometry & GEOMETRY_LIMB) { double r = model->R0 + model->layer[lnum]->zbase; model->p = model->layer[lnum]->nbase * r; } else { return 1; } return 0; } /* compute_refractive_invariant() */ /*********************************************************** * static int copy_model_scalars(model_t *model, model_t *lmodel) * * Purpose: * Copies scalar variables from *model to *lmodel, where * *lmodel is a model structure recording the state of * a prior model computation. The variables that need to * be copied are those that are compared from one model * computation to the next when determining whether * expensive spectral computations need to be redone, and * those that are involved in insertion and removal of * interpolated layers. * * Arguments: * model_t *model, *lmodel - pointers to model structures * * Return: * 0 if successful, 1 otherwise ************************************************************/ static int copy_model_scalars(model_t *model, model_t *lmodel) { int lnum, cnum, knum; if (model->nlayers != lmodel->nlayers) { errlog(35, 0); return 1; } lmodel->flo = model->flo; lmodel->fif_min = model->fif_min; lmodel->fif_max = model->fif_max; lmodel->ils_fwhm = model->ils_fwhm; lmodel->ils_fif = model->ils_fif; lmodel->dsb_utol_ratio = model->dsb_utol_ratio; lmodel->T0 = model->T0; lmodel->Tref = model->Tref; lmodel->sec_za = model->sec_za; lmodel->isub[0] = model->isub[0]; lmodel->isub[1] = model->isub[1]; lmodel->ifmode = model->ifmode; lmodel->ilsmode = model->ilsmode; lmodel->ils_typenum = model->ils_typenum; /* * Layer variables are copied only from layers that lie * within the range covered by the propagation path indices. * Outside this range they are not guaranteed to be updated. */ for (lnum = model->path_min; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; layer_t *llayer = lmodel->layer[lnum]; if (layer->ncols != llayer->ncols) { errlog(35, 0); return 1; } llayer->P = layer->P; llayer->T = layer->T; llayer->dP = layer->dP; llayer->dP_def = layer->dP_def; llayer->Pbase = layer->Pbase; llayer->Tbase = layer->Tbase; llayer->h = layer->h; llayer->m = layer->m; llayer->type = layer->type; for (cnum = 0; cnum < layer->ncols; ++cnum) { column_t *column = layer->column[cnum]; column_t *lcolumn = llayer->column[cnum]; if (column->n_abscoeffs != lcolumn->n_abscoeffs) { errlog(35, 0); return 1; } lcolumn->N = column->N; lcolumn->N_def = column->N_def; lcolumn->N_scaled = column->N_scaled; lcolumn->xvmr = column->xvmr; lcolumn->vmr_scaled = column->vmr_scaled; lcolumn->RH = column->RH; for (knum = 0; knum < column->n_abscoeffs; ++knum) { abscoeff_t *abscoeff = column->abscoeff[knum]; abscoeff_t *labscoeff = lcolumn->abscoeff[knum]; labscoeff->vmr_selfbroad = abscoeff->vmr_selfbroad; } } } /* * ngrid is initialized to 0, but is always non-zero in a * valid model. It is copied here to indicate that valid * model data have been copied to *lmodel. */ lmodel->ngrid = model->ngrid; return 0; } /* copy_model_scalars() */ /*********************************************************** * int find_n_from_z( * double *n, const model_t *model, const double z) * * Purpose: * Given a height z, finds the refractive index at that * height by interpolating in density (P/T) between model * definition levels. If z lies outside the model * definition levels, n-1 is scaled by P/T. * * This interpolation scheme ensures stability of path * geometry computations. It avoids inadvertent creation * of ducts caused by composition changes above thin * interpolated layers, and maintains continuity of the * refractive index profile when solving for the tangent * level of paths below the astronomical horizon. * * This function assumes that layer base heights are up to * date on model definition layers, and that model pressure * limits have been set. * * Arguments: * double *n - pointer to receive refractive index * model_t *model - pointer to model data structure * const double z - height * * Return: * 0 on success, 1 on failure ************************************************************/ int find_n_from_z( double *n, const model_t *model, const double z) { int lnum_lo, lnum_hi; layer_t *layer_lo = NULL; layer_t *layer_hi = NULL; double P; if (model->Pmax == 0.0) { /* * There are no model definition layers at non-zero * pressure. */ *n = 1.0; return 0; } /* * find_P_from_z() clamps P at a maximum extrapolated value * of model->Pmax. Below that level, we assume constant P * for the purpose of computing n. This situation arises * when the tangent point lies below the planetary surface. * This is normal for nadir sounding paths, and will not * cause trouble since in such cases the source level is on * the near side of the tangent point and no tangent level * will be inserted. For this reason, we ignore the error * return unless P == 0.0. (Note that a non-error return * with P == 0.0 is also possible.) */ if (find_P_from_z(&P, model, z) && P == 0.0) return 1; /* * Scan for two model definition layers with nonzero base * pressures that have zbase values that bracket z. Here, * lnum_hi and lnum_lo refer to indices of model definition * levels above and below the interpolated level, * respectively. */ lnum_hi = -1; for (lnum_lo = 0; lnum_lo < model->nlayers; ++lnum_lo) { layer_lo = model->layer[lnum_lo]; if (layer_lo->type != LAYER_TYPE_DEF) { /* skip if not a model definition layer */ continue; } else if (layer_lo->Pbase == 0.0) { /* skip zero-pressure model definition layer */ continue; } else if (layer_lo->dP == 0.0 && layer_lo->h >= 0.0) { /* skip plume layer */ continue; } else if (layer_lo->zbase > z) { /* update upper z boundary */ lnum_hi = lnum_lo; layer_hi = model->layer[lnum_hi]; continue; } else { /* found lower z boundary */ break; } } if (lnum_hi < 0 && lnum_lo == model->nlayers) { /* * No model definition layers above or below z with * nonzero base pressure. Not a normal case. */ *n = 1.0; } else if (lnum_hi < 0) { /* * z is in the highest model definition layer with * nonzero base pressure. Extrapolate isothermally from * the base of this layer upwards. This is equivalant to * scaling n-1 by pressure. */ *n = layer_lo->nbase - 1.0; *n *= P / layer_lo->Pbase; *n += 1.0; } else if (lnum_hi < lnum_lo && lnum_lo < model->nlayers) { /* * z is bracketed by two model definition levels. * Compute n by interpolating in P/T. */ double T; if (model->PTmode & PTMODE_TBASE) { T = log_x_interp( layer_lo->Pbase, layer_lo->Tbase, layer_hi->Pbase, layer_hi->Tbase, P); } else { T = layer_lo->T; } *n = lin_interp( layer_lo->Pbase / layer_lo->Tbase, layer_lo->nbase, layer_hi->Pbase / layer_hi->Tbase, layer_hi->nbase, P / T); } else { /* * z is below the lowest model definition level. Assume * constant molecular refractivity and scale n-1 by P/T. */ double gamma = layer_base_lapse_rate(model, lnum_hi); double T = layer_hi->Tbase + gamma * (layer_hi->zbase - z); *n = layer_hi->nbase - 1.0; *n *= (P / layer_hi->Pbase) * (layer_hi->Tbase / T); *n += 1.0; } return 0; } /* find_n_from_z() */ /*********************************************************** * static int find_P_from_z( * double *P, const model_t *model, const double z) * * Purpose: * Given a height z, finds the corresponding pressure level * in a hydrostatic model, by interpolating log(P) linearly * in z. This function assumes that layer base heights are * up to date on model definition layers, and that model * pressure limits have been set. * * If z is not bracketed by model definition levels, P is * extrapolated isothermally at the top of the atmosphere, * and extrapolated according to the extrapolation mode * setting at the base of the atmosphere. * * If the pressure corresponding to z is lower than the * nominal top-of-atmosphere pressure model->Ptoa, then the * height z is regarded as lying outside the atmosphere, * and this function returns P = 0. * * If the pressure corresponding to z is higher than the * base pressure extrapolation limit model->Pmax, P is * clamped at model->Pmax. * * Arguments: * double *P - pointer to receive computed pressure * model_t *model - pointer to model data structure * const double z - height * * Return: * 0 on success, 1 on failure or if P is clamped at Pmax ************************************************************/ static int find_P_from_z( double *P, const model_t *model, const double z) { int lnum_lo, lnum_hi; layer_t *layer_lo = NULL; layer_t *layer_hi = NULL; if (model->Pmax == 0.0) { /* * There are no model definition layers at non-zero * pressure. */ *P = 0.0; return 0; } /* * Scan for two model definition layers with nonzero base * pressures that have zbase values that bracket z. Here, * lnum_hi and lnum_low refer to indices of model definition * levels above and below the interpolated level, * respectively. */ lnum_hi = -1; for (lnum_lo = 0; lnum_lo < model->nlayers; ++lnum_lo) { layer_lo = model->layer[lnum_lo]; if (layer_lo->type != LAYER_TYPE_DEF) { /* skip if not a model definition layer */ continue; } else if (layer_lo->Pbase == 0.0) { /* skip zero-pressure model definition layer */ continue; } else if (layer_lo->dP == 0.0 && layer_lo->h >= 0.0) { /* skip plume layer */ continue; } else if (layer_lo->zbase > z) { /* update upper z boundary */ lnum_hi = lnum_lo; layer_hi = model->layer[lnum_hi]; continue; } else { /* found lower z boundary */ break; } } if (lnum_hi < 0 && lnum_lo == model->nlayers) { /* * No model definition layers above or below z with * nonzero base pressure. Not a normal case. */ *P = 0.0; } else if (lnum_hi < 0) { /* * z is in the highest model definition layer with * nonzero base pressure. Extrapolate isothermally from * the base of this layer upwards. If the pressure is * below the top-of-atmosphere extrapolation limit, set * it to zero. */ double R, T, x; T = model->PTmode & PTMODE_TBASE ? layer_lo->Tbase : layer_lo->T; if (T == 0.0) { errlog(188, 0); *P = 0.0; return 1; } R = PCONST_KB / (layer_lo->Mair * PCONST_AMU); R *= unit_tab[UNIT_M].factor * unit_tab[UNIT_M].factor; x = (layer_lo->zbase - z) * layer_lo->gbase / (R * T); *P = layer_lo->Pbase * exp(x); if (*P < model->Ptoa) *P = 0.0; } else if (lnum_hi < lnum_lo && lnum_lo < model->nlayers) { /* * z is bracketed by two model definition levels. * Compute P by linearly interpolating log(P) in z. */ *P = log_y_interp( layer_lo->zbase, layer_lo->Pbase, layer_hi->zbase, layer_hi->Pbase, z); } else { /* * z is below the lowest model definition level. Get the * lapse rate gamma corresponding to the base extension * mode setting and extrapolate. */ double gamma, dz, Tbar, R, x; gamma = layer_base_lapse_rate(model, lnum_hi); dz = layer_hi->zbase - z; Tbar = layer_hi->Tbase + gamma * 0.5 * dz; if (Tbar <= 0.0) { errlog(188, 0); *P = 0.0; return 1; } R = PCONST_KB / (layer_hi->Mair * PCONST_AMU); R *= unit_tab[UNIT_M].factor * unit_tab[UNIT_M].factor; x = (layer_hi->zbase - z) * layer_hi->gbase / (R * Tbar); *P = layer_hi->Pbase * exp(x); if (*P > model->Pmax) { *P = model->Pmax; return 1; } } return 0; } /* find_P_from_z() */ /*********************************************************** * static int find_z_from_P( * double *z, const model_t *model, const double P) * * Purpose: * Given a pressure level P, finds the corresponding height * z by interpolating linearly in log(P) between model * definition levels. This function assumes that base * level heights have been computed on all model definition * layers, and model atmosphere pressure limits have been * set. * * If P is not bracketed by model definition levels, z is * found by isothermal extrapolation. * * Return: * 0 on success, 1 on failure ************************************************************/ static int find_z_from_P( double *z, const model_t *model, const double P) { int lnum_lo, lnum_hi; layer_t *layer_lo = NULL; layer_t *layer_hi = NULL; if (P <= 0.0) { /* * Fail for non-positive P */ errlog(163, 0); return 1; } else if (P < model->Ptoa) { /* * If P is more than two scale heights above the * highest model definition level, log a warning. */ errlog(164, 0); } else if (P > model->Pmax) { /* * Fail if P would extrapolate more than one scale * height below the lowest model definition level. * This also catches the case when there are no model * definition levels at non-zero pressure. */ errlog(171, 0); *z = 0.0; return 1; } /* * Scan for two model definition layers having Pbase * values that bracket P. */ lnum_hi = -1; for (lnum_lo = 0; lnum_lo < model->nlayers; ++lnum_lo) { layer_lo = model->layer[lnum_lo]; if (layer_lo->type != LAYER_TYPE_DEF) { /* skip if not a model definition layer */ continue; } else if (layer_lo->Pbase < P) { /* update upper P boundary */ lnum_hi = lnum_lo; layer_hi = model->layer[lnum_hi]; } else if (layer_lo->Pbase >= P) { /* found lower P boundary */ break; } } if (lnum_hi < 0 && lnum_lo == model->nlayers) { /* * No model definition layers above or below P. This * case should already have been trapped by the tests * above. */ errlog(205, __LINE__); return 1; } else if (lnum_hi < 0) { /* * P is in the highest model definition layer in the * atmosphere. Extrapolate isothermally from the base of * this layer upwards to find z. */ double R, T; T = model->PTmode & PTMODE_TBASE ? layer_lo->Tbase : layer_lo->T; R = PCONST_KB / (layer_lo->Mair * PCONST_AMU); R *= unit_tab[UNIT_M].factor * unit_tab[UNIT_M].factor; *z = layer_lo->zbase - (R * T / layer_lo->gbase) * log(P / layer_lo->Pbase); } else if (lnum_hi < lnum_lo && lnum_lo < model->nlayers) { /* * P is bracketed by two model definition levels. * Compute z by linearly interpolating z in log(P). */ *z = log_x_interp( layer_lo->Pbase, layer_lo->zbase, layer_hi->Pbase, layer_hi->zbase, P); } else { /* * P is below the lowest model definition level. Get the * lapse rate gamma corresponding to the base extension * mode setting and extrapolate. */ double gamma, R, x, dz; gamma = layer_base_lapse_rate(model, lnum_hi); R = PCONST_KB / (layer_hi->Mair * PCONST_AMU); R *= unit_tab[UNIT_M].factor * unit_tab[UNIT_M].factor; x = log(P / layer_hi->Pbase); dz = x * R * layer_hi->Tbase; dz /= layer_hi->gbase - 0.5 * R * gamma * x; *z = layer_hi->zbase - dz; } return 0; } /* find_z_from_P() */ /*********************************************************** * static void set_hydrostatic_pressure_limits(model_t *model) * * Purpose: * Set pressure limits for hydrostatic models. These are * the top of atmosphere pressure Ptoa, defined as the * product of TOA_EPSILON and the base pressure of the * highest model definition layer; and the pressure * extrapolation limit Pmax, defined as the product of * P_OVERSHOOT_LIMIT and the base pressure of the lowest * model layer. See the definitions of these constants * near the top of this file for details. * * Arguments: * model_t *model - pointer to model data structure ************************************************************/ static void set_hydrostatic_pressure_limits(model_t *model) { int lnum; /* * If this is not a hydrostatic model, do nothing. */ if (!(model->PTmode & PTMODE_HYDROSTATIC)) return; /* * Scan for the highest model definition layer having a * non-zero base pressure. */ model->Ptoa = 0.0; for (lnum = 0; lnum < model->nlayers; ++lnum) { layer_t *layer = model->layer[lnum]; if (layer->type == LAYER_TYPE_DEF && layer->Pbase > 0.0) { model->Ptoa = layer->Pbase; break; } } model->Ptoa *= TOA_EPSILON; /* * Scan for the lowest model definition layer having a * non-zero base pressure. */ model->Pmax = 0.0; for (lnum = model->nlayers - 1; lnum >= 0; --lnum) { layer_t *layer = model->layer[lnum]; if (layer->type == LAYER_TYPE_DEF && layer->Pbase > 0.0) { model->Pmax = layer->Pbase; break; } } model->Pmax *= P_OVERSHOOT_LIMIT; return; } /* set_hydrostatic_pressure_limits() */ /*********************************************************** * static int set_user_defined_levels(model_t *model, model_t *lmodel) * * Purpose: * Inserts any user-defined interpolated source, observing, * and tangent levels into the model, inserting the * additional layers needed to split the associated parent * model definition layers. If any of these layers were * carried over from a prior computation, they will be * adjusted as needed. This adjustment will include * deletion and re-insertion if the level has moved to a * different parent layer, changing the layer ordering. * * 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: * 0 on success, 1 on error. ************************************************************/ static int set_user_defined_levels(model_t *model, model_t *lmodel) { set_hydrostatic_pressure_limits(model); if (model->geometry & GEOMETRY_POBS_USER_DEFINED) { double zobs, nobs; if ( find_z_from_P(&zobs, model, model->Pobs) || find_n_from_z(&nobs, model, zobs) || insert_interpolated_level( model, lmodel, LAYER_TYPE_OBS, zobs, model->Pobs, nobs)) { errlog(190, 0); return 1; } } if (model->geometry & GEOMETRY_ZOBS_USER_DEFINED) { double Pobs, nobs; if (find_P_from_z(&Pobs, model, model->zobs)) { errlog(171, 0); errlog(190, 0); return 1; } if ( find_n_from_z(&nobs, model, model->zobs) || insert_interpolated_level( model, lmodel, LAYER_TYPE_OBS, model->zobs, Pobs, nobs)) { errlog(190, 0); return 1; } } if (model->geometry & GEOMETRY_PSOURCE_USER_DEFINED) { double zsource, nsource; if ( find_z_from_P(&zsource, model, model->Psource) || find_n_from_z(&nsource, model, zsource) || insert_interpolated_level( model, lmodel, LAYER_TYPE_SOURCE, zsource, model->Psource, nsource)) { errlog(191, 0); return 1; } } if (model->geometry & GEOMETRY_ZSOURCE_USER_DEFINED) { double Psource, nsource; if (find_P_from_z(&Psource, model, model->zsource)) { errlog(171, 0); errlog(191, 0); return 1; } if ( find_n_from_z(&nsource, model, model->zsource) || insert_interpolated_level( model, lmodel, LAYER_TYPE_SOURCE, model->zsource, Psource, nsource)) { errlog(191, 0); return 1; } } if (model->geometry & GEOMETRY_PTAN_USER_DEFINED) { double ztan, ntan; if ( find_z_from_P(&ztan, model, model->Ptan) || find_n_from_z(&ntan, model, ztan) || insert_interpolated_level( model, lmodel, LAYER_TYPE_TAN, ztan, model->Ptan, ntan)) { errlog(192, 0); return 1; } } if (model->geometry & GEOMETRY_ZTAN_USER_DEFINED) { double Ptan, ntan; if (find_P_from_z(&Ptan, model, model->ztan)) { errlog(171, 0); errlog(192, 0); return 1; } if ( find_n_from_z(&ntan, model, model->ztan) || insert_interpolated_level( model, lmodel, LAYER_TYPE_TAN, model->ztan, Ptan, ntan)) { errlog(192, 0); return 1; } } return 0; } /* set_user_defined_levels() */ /*********************************************************** * static int set_propagation_path_indices(model_t *model) * * Purpose: * Sets up the layer index ranges for the propagation path. * * The propagation path is divided across the tangent point * into upward and downward segments. Both of these * segments occur for limb paths and paths below the * astronomical horizon through the tangent point. * * Paths above the horizon (which includes all plane- * parallel paths), and paths below the horizon from a * source point on the near side of the tangent point, have * a single upward or downward segment. In these cases, * the unused segment is assigned reversed index ranges * relative to the propagation direction so it will always * be skipped. * * An error can occur if the propagation path does not * intersect a user-defined source level. * * Arguments: * model_t *model - pointer to model data structure * * Return: * 0 on success, 1 on error. ************************************************************/ static int set_propagation_path_indices(model_t *model) { int obs_lnum = get_obs_lnum(model); int source_lnum = get_source_lnum(model); int tan_lnum = get_tan_lnum(model); int base_lnum = get_lnum_by_type(model, LAYER_TYPE_DEF); /* * Moderate extrapolation of observing, source, or tangent * levels below the base level of the lowest defined model is * allowed. Nevertheless, log a warning to alert the user, * unless an extrapolation mode has been explicitly specified * in an "extend" statement. */ if (!(model->PTmode & PTMODE_EXTEND_USER_DEFINED)) { if (obs_lnum > base_lnum) errlog(165, 0); if (source_lnum > base_lnum) errlog(177, 0); if (tan_lnum > base_lnum) errlog(166, 0); } if (model->geometry & GEOMETRY_LIMB) { /* * Limb path. The tangent point is either user-defined, * or defaults to the bottom model layer. */ model->path_begin = 0; model->path_mid = tan_lnum; model->path_end = 0; } else if (fabs(model->za) <= PI_ON_TWO) { /* * Path on or above the astronomical horizon. Fail if * the source level is below the astronomical horizon. */ if (source_lnum > obs_lnum) { errlog(178, 0); return 1; } else { model->path_begin = source_lnum + 1; model->path_mid = obs_lnum; model->path_end = obs_lnum + 1; } } else if (!(model->geometry & GEOMETRY_SOURCE_NEAR)) { /* * Path through tangent point below the astronomical * horizon to a source on the far side of the tangent * point. Fail if the source level is below the tangent * level. */ if (source_lnum > tan_lnum) { errlog(179, 0); return 1; } else { model->path_begin = source_lnum + 1; model->path_mid = tan_lnum; model->path_end = obs_lnum + 1; } } else { /* * Path below the astronomical horizon to a source on the * near side of the tangent point. Fail if the source * level isn't between the tangent and observing levels. * * A tangent level in this case could be below the * planetary surface, and will not have been inserted in * the model, so check heights instead of layer indices. */ double zobs = model->layer[obs_lnum]->zbase; double zsource = model->layer[source_lnum]->zbase; if (zsource > zobs || zsource < model->ztan) { errlog(180, 0); return 1; } else { model->path_begin = source_lnum + 1; model->path_mid = source_lnum; model->path_end = obs_lnum + 1; } } if (model->geometry & GEOMETRY_REVERSE) { /* * For reverse geometry, swap path_begin and path_end. * This has the effect of reversing both the order and * the direction of the path. */ int temp = model->path_end; model->path_end = model->path_begin; model->path_begin = temp; } /* * Set the minimum layer index for the path. */ if (model->path_begin < model->path_end) model->path_min = model->path_begin; else model->path_min = model->path_end; return 0; } /* set_propagation_path_indices() */ /*********************************************************** * int setup_atmospheric_model(model_t *model, model_t *lmodel) * * Purpose: * Sets up the atmospheric model and path geometry. * * 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: * 0 on success, 1 on error. ************************************************************/ int setup_atmospheric_model(model_t *model, model_t *lmodel) { /* * If there are no model layers, set up path ranges that * cover no layers and return. */ if (model->nlayers <= 0) { model->path_begin = 0; model->path_mid = -1; model->path_end = 0; return 0; } /* * Set the observer zenith angle za and sec(za). */ if (model->geometry & GEOMETRY_SEC_ZA_USER_DEFINED) { model->za = acos(1.0 / model->sec_za); } else { if (model->geometry & GEOMETRY_PLANE_PARALLEL) { /* * Plane-parallel models are restricted to zenith * angles less than 90 degrees. */ if ((PI_ON_TWO - fabs(model->za)) < FLT_EPSILON) { errlog(161, 0); return 1; } } model->sec_za = 1.0 / cos(model->za); } /* * The atmospheric model is anchored to the layers as defined * in the model configuration file (the "model definition * layers"). In hydrostatic models, additional interpolated * or extrapolated levels may subsequently be added for * user-defined source, observing, or tangent levels. These * require the addition of new model layers created by * splitting a model definition layer (a "parent layer") into * parts, or by extrapolating a parent layer below the layer * stack. * * The first pass through the layer stack is carried out on * model definition layers only. Any added layers carried * over from a prior computation are ignored, and parent * layers are restored to their unmodified state. This pass * computes scalar parameters including mixing ratios, * geometric heights, and base level refractive indicies, * which are inherited by or interpolated on user-defined * levels and their associated layers. */ restore_def_layer_explicit_col_densities(model); set_def_layer_pressures(model); set_def_layer_temperatures(model); set_def_layer_vmrs_and_apply_Nscale_factors(model); set_def_layer_heights_and_col_densities(model); set_def_layer_refractive_indices(model); /* * Add user-defined interpolated source, observing, or * tangent levels and associated layers. Existing * interpolated layers carried over from a prior computation * are deleted and reinserted if the layer ordering will * change, otherwise layers are adjusted in place. * * Because these operations affect the number of layers or * their ordering, they are carried out on both *model and * *lmodel, but level positions and other variables are * updated on *model only. On initial layer insertion, or * in the case of a layer deletion and re-insertion, the * default initialization of parameters on *lmodel will * trigger (re)computation of all spectra on that layer. */ if (set_user_defined_levels(model, lmodel)) return 1; /* * Now that we know the observer level height and zenith * angle (or the tangent level height in limb mode) we have * what we need to compute the refractive invariant. */ if (compute_refractive_invariant(model)) return 1; /* * If we're in spherical geometry mode, and looking below the * astronomical horizon, iteratively solve for and insert the * tangent level. */ if ((model->geometry & GEOMETRY_SPHERICAL) && (fabs(model->za) > PI_ON_TWO)) { if (solve_for_tan_level(model, lmodel)) return 1; } /* * With all interpolated levels now in place, set and adjust * remaining variables on interpolated and parent layers. */ set_interpolated_layer_variables(model); compute_adjusted_selfbroad_vmrs(model); /* * Compute the propagation path geometry. */ if (set_propagation_path_indices(model)) return 1; if (compute_airmass_and_refraction(model)) return 1; return 0; } /* setup_atmospheric_model() */ /*********************************************************** * static int solve_for_tan_level( * model_t *model, model_t *lmodel) * * Purpose: * In spherical geometry mode, for paths below the * astronomical horizon, this function solves for the * tangent height and inserts a tangent level. Because the * refractive index n varies with altitude, this has to be * done iteratively. * * Arguments: * model_t *model - pointer to model data structure * * Return: * 0 on success, 1 on error or failure to converge. ************************************************************/ static int solve_for_tan_level(model_t *model, model_t *lmodel) { double ntan, ztan; int iter, lnum; /* * If a tangent layer already exists from a prior * computation, start the iteration at the base of that * layer. Otherwise, start at the base of the lowest model * layer. */ if ((lnum = get_lnum_by_type(model, LAYER_TYPE_TAN)) < 0) lnum = model->nlayers - 1; ztan = model->layer[lnum]->zbase; for (iter = 0; iter < ZTAN_ITER_MAX; ++iter) { double n, ztan_last; if (find_n_from_z(&n, model, ztan)) { errlog(168, 0); return 1; } ztan_last = ztan; ztan = (model->p / n) - model->R0; if (fabs((ztan - ztan_last) / ztan) < ZTAN_TOL) break; } /* * The iteration may fail to converge to a solution due to * ducting, extrapolation beyond defined layers, or for some * other reason. If so, log an error and give up. */ if (iter == ZTAN_ITER_MAX) { errlog(168, 0); return 1; } /* * If the tangent level pressure is higher than the * extrapolation limit, and the source is on the far side of * the tangent point from the observer, log an error and * fail. */ if (find_P_from_z(&model->Ptan, model, ztan) && !(model->geometry & GEOMETRY_SOURCE_NEAR)) { errlog(196, 0); return 1; } /* * If the tangent level height is below the planetary surface * and the source is on the far side of the tangent point * from the observer, log a warning, but continue * nonetheless. */ if (ztan < 0.0 && !(model->geometry & GEOMETRY_SOURCE_NEAR)) errlog(167, 0); model->ztan = ztan; if (find_n_from_z(&ntan, model, ztan)) return 1; /* * Insert a tangent level if the source level is on the far * side of the tangent point. */ if (!(model->geometry & GEOMETRY_SOURCE_NEAR) && insert_interpolated_level( model, lmodel, LAYER_TYPE_TAN, model->ztan, model->Ptan, ntan)) { errlog(182, 0); return 1; } return 0; } /* solve_for_tan_level() */ /*********************************************************** * double source_to_obs_geometric_distance(model_t *model) * * Purpose: * Computes the geometric (i.e. unrefracted) distance * between source and observer. If either path endpoint * lies at infinity, the return value is DBL_MAX. * * Arguments: * model_t *model - pointer to model data structure * * Return: * geometric distance ************************************************************/ double source_to_obs_geometric_distance(model_t *model) { if (model->path_begin == 0 || model->path_end == 0) return DBL_MAX; if (model->geometry & GEOMETRY_PLANE_PARALLEL) { double x = 0.0; double z = 0.0; int lnum; for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; double dz; if (layer->h > 0.0) dz = layer->h; else dz = model->layer[lnum - 1]->zbase - layer->zbase; x += dz * tan(layer->za_base); z += dz; } for (lnum = model->path_mid; lnum >= model->path_end; --lnum) { layer_t *layer = model->layer[lnum]; double dz; if (layer->h > 0.0) dz = layer->h; else dz = model->layer[lnum - 1]->zbase - layer->zbase; x += dz * tan(layer->za_base); z += dz; } return sqrt(x * x + z * z); } else { double phi = 0.0; double Robs, Rsource; int lnum; for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) phi += model->layer[lnum]->dphi; for (lnum = model->path_mid ; lnum >= model->path_end; --lnum) phi += model->layer[lnum]->dphi; lnum = get_obs_lnum(model); Robs = model->R0 + model->layer[lnum]->zbase; lnum = get_source_lnum(model); Rsource = model->R0 + model->layer[lnum]->zbase; return sqrt(Robs * Robs + Rsource * Rsource - 2.0 * Robs * Rsource * cos(phi)); } } /* source_to_obs_geometric_distance() */ /*********************************************************** * double source_to_obs_path_distance(model_t *model) * * Purpose: * Computes the total refracted path distance from source * to observer. If either path endpoint lies at infinity, * the return value is DBL_MAX. * * Arguments: * model_t *model - pointer to model data structure * * Return: * path distance ************************************************************/ double source_to_obs_path_distance(model_t *model) { double d = 0.0; int lnum; if (model->path_begin == 0 || model->path_end == 0) return DBL_MAX; for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; if (layer->h > 0.0) d += layer->h; else d += layer->m * (model->layer[lnum - 1]->zbase - layer->zbase); } for (lnum = model->path_mid; lnum >= model->path_end; --lnum) { layer_t *layer = model->layer[lnum]; if (layer->h > 0.0) d += layer->h; else d += layer->m * (model->layer[lnum - 1]->zbase - layer->zbase); } return d; } /* source_to_obs_path_distance() */ /*********************************************************** * double source_to_obs_projected_distance(model_t *model) * * Purpose: * In plane-parallel models, computes the length of the * horizontal plane projection of the path from source to * observer. * * In spherical models, computes the projected great circle * map distance along the planetary surface between source * and observer nadir points. * * Arguments: * model_t *model - pointer to model data structure * * Return: * Projected map distance ************************************************************/ double source_to_obs_projected_distance(model_t *model) { if (model->geometry & GEOMETRY_PLANE_PARALLEL) { double x = 0.0; int lnum; /* * If either path endpoint is at infinity in a * plane-parallel model, the projected path length * is infinite unless the zenith angle is zero. */ if (model->path_begin == 0 || model->path_end == 0) return fabs(model->za) < DBL_MIN ? 0.0 : DBL_MAX; for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; double dz; if (layer->h > 0.0) dz = layer->h; else dz = model->layer[lnum - 1]->zbase - layer->zbase; x += dz * tan(layer->za_base); } for (lnum = model->path_mid; lnum >= model->path_end; --lnum) { layer_t *layer = model->layer[lnum]; double dz; if (layer->h > 0.0) dz = layer->h; else dz = model->layer[lnum - 1]->zbase - layer->zbase; x += dz * tan(layer->za_base); } return x; } else { double phi = 0.0; int lnum; for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) phi += model->layer[lnum]->dphi; for (lnum = model->path_mid ; lnum >= model->path_end; --lnum) phi += model->layer[lnum]->dphi; return phi * model->R0; } } /* source_to_obs_projected_distance() */ /*********************************************************** * double total_airmass(model_t *model) * * Purpose: * Computes the total relative airmass by column density, * irrespective of species, for the propagation path. * * Airmass is normalized to the sum of the zenith column * densities of all layers in the path from source to * observer. If the path traverses a tangent point, layers * that are traversed twice are only counted once in the * zenith column density sum. * * Arguments: * model_t *model - pointer to model data structure * * Return: * Total nominal airmass, or 0.0 if undefined. ************************************************************/ double total_airmass(model_t *model) { double N = 0.0; double N_los = 0.0; int lnum; int path_begin = model->path_begin; int path_mid = model->path_mid; int path_end = model->path_end; int path_min = path_begin < path_end ? path_begin : path_end; /* * Zenith total column density. */ for (lnum = path_mid; lnum >= path_min; --lnum) { layer_t *layer = model->layer[lnum]; int cnum; for (cnum = 0; cnum < layer->ncols; ++cnum) { column_t *column = layer->column[cnum]; int ctype = column->col_typenum; if (!(col_type[ctype].flags & COL_PARAMETRIC)) N += column->N_scaled; } } /* * Line-of-sight total column density. */ for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) { layer_t *layer = model->layer[lnum]; int cnum; for (cnum = 0; cnum < layer->ncols; ++cnum) { column_t *column = layer->column[cnum]; int ctype = column->col_typenum; if (!(col_type[ctype].flags & COL_PARAMETRIC)) N_los += column->N_scaled * layer->m; } } for (lnum = model->path_mid; lnum >= model->path_end; --lnum) { layer_t *layer = model->layer[lnum]; int cnum; for (cnum = 0; cnum < layer->ncols; ++cnum) { column_t *column = layer->column[cnum]; int ctype = column->col_typenum; if (!(col_type[ctype].flags & COL_PARAMETRIC)) N_los += column->N_scaled * layer->m; } } /* * The airmass is the ratio of line-of-sight to zenith column * density. */ return N > 0.0 ? N_los / N : 0.0; } /* total_airmass() */ /*********************************************************** * double total_refraction(model_t *model) * * Purpose: * Computes the total refraction through the model * propagation path. For spherical paths, the sign of the * refractive bending is considered positive when the line * of sight is concave towards the planetary center, * irrespective of propagation direction. Similarly, for * plane-parallel geometry, the sign is positive for * paths for which the zenith angle is smaller at lower * altitude, irrespective of propagation direction. * * * Arguments: * model_t *model - pointer to model data structure * * Return: * Total refraction in radians. ************************************************************/ double total_refraction(model_t *model) { double R = 0.0; int lnum; if (model->geometry & GEOMETRY_PLANE_PARALLEL) { /* * In the plane parallel case, the total refraction is * just the difference in zenith angles between the two * ends of a path segment. */ if (model->path_begin <= model->path_mid) { R += model->layer[model->path_begin]->za_base; R -= model->layer[model->path_mid]->za_base; } if (model->path_mid >= model->path_end) { R += model->layer[model->path_end]->za_base; R -= model->layer[model->path_mid]->za_base; } } else { /* * In the spherical case, the refractive bending is zero * for a layer with an upper boundary pressure of zero, * so all such layers are skipped. (One way to think of * this is that, formally, such a layer has a pressure * gradient asymptotically equal to zero, and hence no * refractive index gradient.) */ for (lnum = model->path_begin; lnum <= model->path_mid; ++lnum) { if (lnum != 0 && model->layer[lnum - 1]->Pbase > 0.0) R += model->layer[lnum - 1]->za_base - model->layer[lnum]->za_base + model->layer[lnum]->dphi; } for (lnum = model->path_mid; lnum >= model->path_end; --lnum) { if (lnum != 0 && model->layer[lnum - 1]->Pbase > 0.0) R += model->layer[lnum - 1]->za_base - model->layer[lnum]->za_base + model->layer[lnum]->dphi; } } return R; } /* total_refraction() */