Removed calls to GSL macros and functions.

src/filter.c

@@ -77,7 +77,7 @@ int main(
     /* Triangle... */
     else if (type == 1) {
       ff[fn] = 1.0 - fabs(nu - center) / fwhm;
-      ff[fn] = GSL_MAX(ff[fn], 0.0);
+      ff[fn] = MAX(ff[fn], 0.0);
     }
 
     /* Gaussian... */

src/invert.c

@@ -139,7 +139,7 @@ int main(
     atm.np = n = 0;
     for (i = 0; i < NMAX; i++) {
       ndata[i] = 0;
-      x[i] = y[i] = GSL_NAN;
+      x[i] = y[i] = NAN;
     }
 
     /* Loop over lines... */
@@ -162,10 +162,10 @@ int main(
 
 	/* Get maxima... */
 	if (data == 1) {
-	  x[i] = (gsl_finite(x[i]) ? GSL_MAX(x[i], obs_sim) : obs_sim);
-	  y[i] = (gsl_finite(y[i]) ? GSL_MAX(y[i], obs_meas) : obs_meas);
+	  x[i] = (isfinite(x[i]) ? MAX(x[i], obs_sim) : obs_sim);
+	  y[i] = (isfinite(y[i]) ? MAX(y[i], obs_meas) : obs_meas);
 	  y_err[i] = obs_err;
-	  if (gsl_finite(x[i]) && gsl_finite(y[i]))
+	  if (isfinite(x[i]) && isfinite(y[i]))
 	    ndata[i] = 1;
 	}
 
@@ -224,15 +224,15 @@ int main(
 	if (ndata[i] > 0) {
 	  x[i] /= ndata[i];
 	  y[i] /= ndata[i];
-	  y_err[i] = sqrt(GSL_MAX(y_err[i] / ndata[i] - POW2(y[i]), 0.0))
+	  y_err[i] = sqrt(MAX(y_err[i] / ndata[i] - POW2(y[i]), 0.0))
 	    / sqrt(ndata[i]);	/* standard error! */
 	}
 
     /* Filter data... */
     n = 0;
     for (i = 0; i < NMAX; i++)
-      if (ndata[i] > 0 && gsl_finite(x[i]) && gsl_finite(y[i])
-	  && gsl_finite(y_err[i])) {
+      if (ndata[i] > 0 && isfinite(x[i]) && isfinite(y[i])
+	  && isfinite(y_err[i])) {
 	x2[n] = x[i];
 	y2[n] = y[i];
 	y2_err[n] = y_err[i];

src/jurassic.c

@@ -3009,7 +3009,7 @@ void copy_obs(
   if (init)
     for (int id = 0; id < ctl->nd; id++)
       for (int ir = 0; ir < obs_dest->nr; ir++)
-	if (gsl_finite(obs_dest->rad[id][ir])) {
+	if (isfinite(obs_dest->rad[id][ir])) {
 	  obs_dest->rad[id][ir] = 0;
 	  obs_dest->tau[id][ir] = 0;
 	}
@@ -3044,7 +3044,7 @@ void formod(
   /* Save observation mask... */
   for (int id = 0; id < ctl->nd; id++)
     for (int ir = 0; ir < obs->nr; ir++)
-      mask[id * NR + ir] = !gsl_finite(obs->rad[id][ir]);
+      mask[id * NR + ir] = !isfinite(obs->rad[id][ir]);
 
   /* Hydrostatic equilibrium... */
   hydrostatic(ctl, atm);
@@ -3071,7 +3071,7 @@ void formod(
   for (int id = 0; id < ctl->nd; id++)
     for (int ir = 0; ir < obs->nr; ir++)
       if (mask[id * NR + ir])
-	obs->rad[id][ir] = GSL_NAN;
+	obs->rad[id][ir] = NAN;
 
   /* Free... */
   free(mask);
@@ -3158,8 +3158,8 @@ void formod_fov(
 
     /* Get radiance and transmittance profiles... */
     int nz = 0;
-    for (int ir2 = GSL_MAX(ir - NFOV, 0);
-	 ir2 < GSL_MIN(ir + 1 + NFOV, obs->nr); ir2++)
+    for (int ir2 = MAX(ir - NFOV, 0);
+	 ir2 < MIN(ir + 1 + NFOV, obs->nr); ir2++)
       if (obs->time[ir2] == obs->time[ir]) {
 	z[nz] = obs2->vpz[ir2];
 	for (int id = 0; id < ctl->nd; id++) {
@@ -3654,7 +3654,7 @@ void init_srcfunc(
     /* Get minimum grid spacing... */
     double dnu = 1.0;
     for (int i = 1; i < n; i++)
-      dnu = GSL_MIN(dnu, nu[i] - nu[i - 1]);
+      dnu = MIN(dnu, nu[i] - nu[i - 1]);
 
     /* Compute source function table... */
 #pragma omp parallel for default(none) shared(ctl,tbl,id,nu,f,n,dnu)
@@ -3780,7 +3780,7 @@ void intpol_tbl_cga(
 		       tbl->p[id][ig][ipr + 1], eps11, los->cgp[ip][ig]);
 
 	  /* Check emssivity range... */
-	  eps00 = GSL_MAX(GSL_MIN(eps00, 1), 0);
+	  eps00 = MAX(MIN(eps00, 1), 0);
 
 	  /* Determine segment emissivity... */
 	  eps = 1 - (1 - eps00) / tau_path[id][ig];
@@ -3875,7 +3875,7 @@ void intpol_tbl_ega(
 		      tbl->p[id][ig][ipr + 1], eps11, los->p[ip]);
 
 	  /* Check emssivity range... */
-	  eps00 = GSL_MAX(GSL_MIN(eps00, 1), 0);
+	  eps00 = MAX(MIN(eps00, 1), 0);
 
 	  /* Determine segment emissivity... */
 	  eps = 1 - (1 - eps00) / tau_path[id][ig];
@@ -4051,11 +4051,11 @@ void kernel(
     /* Set perturbation size... */
     double h;
     if (iqa[j] == IDXP)
-      h = GSL_MAX(fabs(0.01 * gsl_vector_get(x0, j)), 1e-7);
+      h = MAX(fabs(0.01 * gsl_vector_get(x0, j)), 1e-7);
     else if (iqa[j] == IDXT)
       h = 1.0;
     else if (iqa[j] >= IDXQ(0) && iqa[j] < IDXQ(ctl->ng))
-      h = GSL_MAX(fabs(0.01 * gsl_vector_get(x0, j)), 1e-15);
+      h = MAX(fabs(0.01 * gsl_vector_get(x0, j)), 1e-15);
     else if (iqa[j] >= IDXK(0) && iqa[j] < IDXK(ctl->nw))
       h = 1e-4;
     else if (iqa[j] == IDXCLZ || iqa[j] == IDXCLDZ)
@@ -4189,7 +4189,7 @@ size_t obs2y(
   /* Determine measurement vector... */
   for (int ir = 0; ir < obs->nr; ir++)
     for (int id = 0; id < ctl->nd; id++)
-      if (gsl_finite(obs->rad[id][ir])) {
+      if (isfinite(obs->rad[id][ir])) {
 	if (y != NULL)
 	  gsl_vector_set(y, m, obs->rad[id][ir]);
 	if (ida != NULL)
@@ -4237,12 +4237,12 @@ void raytrace(
   /* Get altitude range of atmospheric data... */
   gsl_stats_minmax(&zmin, &zmax, atm->z, 1, (size_t) atm->np);
   if (ctl->nsf > 0) {
-    zmin = GSL_MAX(atm->sfz, zmin);
+    zmin = MAX(atm->sfz, zmin);
     if (atm->sfp > 0) {
       int ip = locate_irr(atm->p, atm->np, atm->sfp);
       double zip = LIN(log(atm->p[ip]), atm->z[ip],
 		       log(atm->p[ip + 1]), atm->z[ip + 1], log(atm->sfp));
-      zmin = GSL_MAX(zip, zmin);
+      zmin = MAX(zip, zmin);
     }
   }
 
@@ -4297,7 +4297,7 @@ void raytrace(
 	xh[i] = x[i] / norm;
       double cosa = fabs(DOTP(ex0, xh));
       if (cosa != 0)
-	ds = GSL_MIN(ctl->rayds, ctl->raydz / cosa);
+	ds = MIN(ctl->rayds, ctl->raydz / cosa);
     }
 
     /* Determine geolocation... */
@@ -5996,7 +5996,7 @@ void y2obs(
   /* Decompose measurement vector... */
   for (int ir = 0; ir < obs->nr; ir++)
     for (int id = 0; id < ctl->nd; id++)
-      if (gsl_finite(obs->rad[id][ir])) {
+      if (isfinite(obs->rad[id][ir])) {
 	obs->rad[id][ir] = gsl_vector_get(y, m);
 	m++;
       }

src/jurassic.h

@@ -107,6 +107,14 @@
       ERRMSG("Error while writing!");					\
   }
 
+/*! @brief Macro to determine the maximum of two values. */
+#define MAX(a,b)				\
+  (((a)>(b))?(a):(b))
+
+/*! Macro to determine the minimum of two values. */
+#define MIN(a,b)				\
+  (((a)<(b))?(a):(b))
+
 /*! Compute linear interpolation. */
 #define LIN(x0, y0, x1, y1, x)			\
   ((y0)+((y1)-(y0))/((x1)-(x0))*((x)-(x0)))

src/raytrace.c

@@ -98,7 +98,7 @@ int main(
       for (int ig = 0; ig < ctl.ng; ig++)
 	atm2.q[ig][ip] = los.q[ip][ig];
       for (int iw = 0; iw < ctl.nw; iw++)
-	atm2.k[iw][ip] = GSL_NAN;
+	atm2.k[iw][ip] = NAN;
     }
 
     /* Save data... */

src/retrieval.c

@@ -608,22 +608,22 @@ void optimal_estimation(
 
       /* Check atmospheric state... */
       for (ip = 0; ip < atm_i->np; ip++) {
-	atm_i->p[ip] = GSL_MIN(GSL_MAX(atm_i->p[ip], 5e-7), 5e4);
-	atm_i->t[ip] = GSL_MIN(GSL_MAX(atm_i->t[ip], 100), 400);
+	atm_i->p[ip] = MIN(MAX(atm_i->p[ip], 5e-7), 5e4);
+	atm_i->t[ip] = MIN(MAX(atm_i->t[ip], 100), 400);
 	for (ig = 0; ig < ctl->ng; ig++)
-	  atm_i->q[ig][ip] = GSL_MIN(GSL_MAX(atm_i->q[ig][ip], 0), 1);
+	  atm_i->q[ig][ip] = MIN(MAX(atm_i->q[ig][ip], 0), 1);
 	for (iw = 0; iw < ctl->nw; iw++)
-	  atm_i->k[iw][ip] = GSL_MAX(atm_i->k[iw][ip], 0);
+	  atm_i->k[iw][ip] = MAX(atm_i->k[iw][ip], 0);
       }
-      atm_i->clz = GSL_MAX(atm_i->clz, 0);
-      atm_i->cldz = GSL_MAX(atm_i->cldz, 0.1);
+      atm_i->clz = MAX(atm_i->clz, 0);
+      atm_i->cldz = MAX(atm_i->cldz, 0.1);
       for (icl = 0; icl < ctl->ncl; icl++)
-	atm_i->clk[icl] = GSL_MAX(atm_i->clk[icl], 0);
-      atm_i->sfz = GSL_MAX(atm_i->sfz, 0);
-      atm_i->sfp = GSL_MAX(atm_i->sfp, 0);
-      atm_i->sft = GSL_MIN(GSL_MAX(atm_i->sft, 100), 400);
+	atm_i->clk[icl] = MAX(atm_i->clk[icl], 0);
+      atm_i->sfz = MAX(atm_i->sfz, 0);
+      atm_i->sfp = MAX(atm_i->sfp, 0);
+      atm_i->sft = MIN(MAX(atm_i->sft, 100), 400);
       for (isf = 0; isf < ctl->nsf; isf++)
-	atm_i->sfeps[isf] = GSL_MIN(GSL_MAX(atm_i->sfeps[isf], 0), 1);
+	atm_i->sfeps[isf] = MIN(MAX(atm_i->sfeps[isf], 0), 1);
 
       /* Forward calculation... */
       formod(ctl, atm_i, obs_i);
@@ -953,7 +953,7 @@ void set_cov_meas(
   for (int ir = 0; ir < obs_err.nr; ir++)
     for (int id = 0; id < ctl->nd; id++)
       obs_err.rad[id][ir]
-	= (gsl_finite(obs->rad[id][ir]) ? ret->err_noise[id] : GSL_NAN);
+	= (isfinite(obs->rad[id][ir]) ? ret->err_noise[id] : NAN);
   obs2y(ctl, &obs_err, sig_noise, NULL, NULL);
 
   /* Forward model error (always considered in retrieval fit)... */