htgop

commit 69d3ddf421892235c04c476ac25f4dc53fe5e6f4
parent ac5a4071d4ab1e6bbf07ef5144d7ca4875708980
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Mon, 30 Jul 2018 16:55:55 +0200

Implement htgop_layer_fetch_sw_ka

Diffstat:
M
src/htgop.c
 | 
70
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
M
src/htgop.h
 | 
27
+++++++++++++++++++++++++++

2 files changed, 94 insertions(+), 3 deletions(-)
diff --git a/src/htgop.c b/src/htgop.c
@@ -94,9 +94,9 @@ trapezoidal_integration
 }
 
 /* The following 3 functions are used to fit the CIE Xbar, Ybar and Zbar curved
- * has defined by the 1931 standard observers. These analytical fits are
- * propsed by C. Wyman, P. P. Sloan & P. Shirley in "Simple Analytic
- * Approximations to the CIE XYZ Color Matching Functions" - JCGT 2013. */
+ * has defined by the 1931 standard. These analytical fits are propsed by C.
+ * Wyman, P. P. Sloan & P. Shirley in "Simple Analytic Approximations to the
+ * CIE XYZ Color Matching Functions" - JCGT 2013. */
 static FINLINE double
 fit_x_bar_1931(const double lambda)
 {
@@ -471,6 +471,7 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
       log_err(htgop,
         "%s:%lu: The levels must be sorted in strict ascending order "
         "wrt their height.\n", rdr.name, rdr.iline);
+      res = RES_BAD_ARG;
       goto error;
     }
   }
@@ -482,6 +483,12 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
 
   /* Parse the length of the tabulation */
   CALL(cstr_to_ulong(read_line(&rdr), &tab_len));
+  if(!tab_len) {
+    log_err(htgop,
+      "%s:%lu: The tabulation length cannot be null.\n", rdr.name, rdr.iline);
+    res = RES_BAD_ARG;
+    goto error;
+  }
 
   /* Allocate the tabulated xH2O of each layer */
   FOR_EACH(ilay, 0, nlays) {
@@ -498,6 +505,7 @@ load_stream(struct htgop* htgop, FILE* stream, const char* stream_name)
         log_err(htgop,
           "%s:%lu: The tabulated x H2O must be sorted in strict ascending "
           "order.\n", rdr.name, rdr.iline);
+        res = RES_BAD_ARG;
         goto error;
       }
     }
@@ -993,6 +1001,62 @@ htgop_position_to_layer_id
   return RES_OK;
 }
 
+res_T
+htgop_layer_fetch_sw_ka
+  (const struct htgop_layer* layer,
+   const size_t ispecint,
+   const size_t iquad_point,
+   const double x_h2o, /* H2O mol / mixture mol */
+   double* out_ka)
+{
+  struct htgop_layer_sw_spectral_interval specint;
+  struct htgop_layer_sw_spectral_interval_tab tab;
+  double* find;
+  double x1, x2, k1, k2;
+  double ka;
+  res_T res = RES_OK;
+
+  if(!layer || !out_ka) return RES_BAD_ARG;
+
+  find = search_lower_bound(&x_h2o, layer->x_h2o_tab, layer->tab_length,
+    sizeof(double), cmp_dbl);
+  if(!find) return RES_BAD_ARG;
+
+  res = htgop_layer_get_sw_spectral_interval(layer, ispecint, &specint);
+  if(res != RES_OK) return res;
+  res = htgop_layer_sw_spectral_interval_get_tab(&specint, iquad_point, &tab);
+  if(res != RES_OK) return res;
+
+  ASSERT(layer->tab_length);
+  ASSERT(tab.tab_length == layer->tab_length);
+
+  if(find == layer->x_h2o_tab) {
+    ASSERT(layer->x_h2o_tab[0] >= x_h2o);
+    x2 = layer->x_h2o_tab[0];
+    k2 = tab.ka_tab[0];
+    ka = k2 / x2*x_h2o;
+  } else {
+    const size_t i = (size_t)(find - layer->x_h2o_tab);
+    ASSERT(i < layer->tab_length);
+    ASSERT(layer->x_h2o_tab[i-0] >= x_h2o);
+    ASSERT(layer->x_h2o_tab[i-1] <  x_h2o);
+    x1 = layer->x_h2o_tab[i-1];
+    x2 = layer->x_h2o_tab[i-0];
+    k1 = tab.ka_tab[i-1];
+    k2 = tab.ka_tab[i-0];
+
+    if(!k1 || !k2) {
+      ka = k1 + (k2-k1)/(x2-x1) * (x_h2o-x1);
+    } else {
+      const double alpha = (log(k2) - log(k1))/(log(x1) - log(x2));
+      const double beta = log(k1) - alpha * log(x1);
+      ka = exp(alpha * log(x_h2o) + beta);
+    }
+  }
+  *out_ka = ka;
+  return RES_OK;
+}
+
 /*******************************************************************************
  * Local functions
  ******************************************************************************/
diff --git a/src/htgop.h b/src/htgop.h
@@ -239,6 +239,33 @@ htgop_position_to_layer_id
    const double pos,
    size_t* ilayer);
 
+/* TODO comment the interpolation procedure */
+HTGOP_API res_T
+htgop_layer_fetch_lw_ka
+  (const struct htgop_layer* layer,
+   const size_t ispecint,
+   const size_t iquad_point,
+   const double x_h2o, /* H2O mol / mixture mol */
+   double* ka);
+
+/* TODO comment the interpolation procedure */
+HTGOP_API res_T
+htgop_layer_fetch_sw_ka
+  (const struct htgop_layer* layer,
+   const size_t ispecint,
+   const size_t iquad_point,
+   const double x_h2o, /* H2O mol / mixture mol */
+   double* ka); 
+
+/* TODO comment the interpolation procedure */
+HTGOP_API res_T
+htgop_layer_fetch_sw_ks
+  (const struct htgop_layer* layer,
+   const size_t ispecint,
+   const size_t iquad_point,
+   const double x_h2o, /* H2O mol / mixture mol */
+   double* ks);
+
 END_DECLS
 
 #endif /* HTGOP */