diff --git a/include/libcloudph++/lgrngn/opts_init.hpp b/include/libcloudph++/lgrngn/opts_init.hpp
index f3df5974..619fed87 100644
--- a/include/libcloudph++/lgrngn/opts_init.hpp
+++ b/include/libcloudph++/lgrngn/opts_init.hpp
@@ -30,7 +30,7 @@ namespace libcloudphxx
       // uses shared_ptr to make opts_init copyable
       typedef std::unordered_map<
         real_t,                // kappa
-        std::shared_ptr<unary_function<real_t>> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
+        std::shared_ptr<unary_function<double>> // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
       > dry_distros_t;
       dry_distros_t dry_distros;
 
@@ -163,7 +163,7 @@ namespace libcloudphxx
       typedef std::unordered_map<
         real_t,                // kappa
         std::tuple<
-          std::shared_ptr<unary_function<real_t>>, // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
+          std::shared_ptr<unary_function<double>>, // n(ln(rd)) @ STP; alternatively it's n(ln(rd)) independent of rhod if aerosol_independent_of_rhod=true
           std::pair<real_t, real_t>, // kappa range of CCN considered to belong to this distribution, ranges of different members of the map need to be exclusive (TODO: add a check of this)
           std::pair<real_t, real_t>  // range of altitudes at which this relaxation acts
         >
diff --git a/src/detail/eval_and_oper.hpp b/src/detail/eval_and_oper.hpp
index 50d38077..00d15059 100644
--- a/src/detail/eval_and_oper.hpp
+++ b/src/detail/eval_and_oper.hpp
@@ -10,12 +10,12 @@ namespace libcloudphxx
       template <typename real_t>
       struct eval_and_mul
       {   
-        const common::unary_function<real_t> &fun;
+        const common::unary_function<double> &fun;
         const real_t &mul;
 
         // ctor
         eval_and_mul(
-          const common::unary_function<real_t> &fun, 
+          const common::unary_function<double> &fun, 
           const real_t &mul
         )
           : fun(fun), mul(mul)
@@ -32,12 +32,12 @@ namespace libcloudphxx
       template <typename real_t>
       struct eval_and_add
       {   
-        const common::unary_function<real_t> &fun;
+        const common::unary_function<double> &fun;
         const real_t &add;
 
         // ctor
         eval_and_add(
-          const common::unary_function<real_t> &fun, 
+          const common::unary_function<double> &fun, 
           const real_t &add
         )
           : fun(fun), add(add)
diff --git a/src/impl/particles_impl.ipp b/src/impl/particles_impl.ipp
index b07d3332..78b8abfc 100644
--- a/src/impl/particles_impl.ipp
+++ b/src/impl/particles_impl.ipp
@@ -422,9 +422,9 @@ namespace libcloudphxx
 
       void sanity_checks();
       void init_SD_with_distros();
-      void init_SD_with_distros_sd_conc(const common::unary_function<real_t> &, const real_t &);
-      void init_SD_with_distros_tail(const common::unary_function<real_t> &, const real_t);
-      void init_SD_with_distros_const_multi(const common::unary_function<real_t> &);
+      void init_SD_with_distros_sd_conc(const common::unary_function<double> &, const real_t &);
+      void init_SD_with_distros_tail(const common::unary_function<double> &, const real_t);
+      void init_SD_with_distros_const_multi(const common::unary_function<double> &);
       void init_SD_with_distros_finalize(const real_t &, const bool unravel_ijk = true);
       void init_SD_with_sizes();
       void init_sanity_check(
@@ -436,23 +436,23 @@ namespace libcloudphxx
 
       void init_dry_sd_conc();
       void init_dry_const_multi(
-        const common::unary_function<real_t> &n_of_lnrd
+        const common::unary_function<double> &n_of_lnrd
       );
       void init_dry_dry_sizes(real_t);
 
       void init_n_sd_conc(
-        const common::unary_function<real_t> &n_of_lnrd
+        const common::unary_function<double> &n_of_lnrd
       );
       void init_n_const_multi(const thrust_size_t &);
       void init_n_dry_sizes(const real_t &conc, const thrust_size_t &sd_conc);
 
       void dist_analysis_sd_conc(
-        const common::unary_function<real_t> &n_of_lnrd,
+        const common::unary_function<double> &n_of_lnrd,
         const n_t sd_conc,
         const real_t dt = 1.
       );
       void dist_analysis_const_multi(
-        const common::unary_function<real_t> &n_of_lnrd 
+        const common::unary_function<double> &n_of_lnrd 
       );
       void reserve_hskpng_npart();
       void init_ijk();
@@ -460,8 +460,8 @@ namespace libcloudphxx
       void init_kappa(const real_t &);
       void init_incloud_time();
       void init_count_num_sd_conc(const real_t & = 1);
-      void init_count_num_const_multi(const common::unary_function<real_t> &);
-      void init_count_num_const_multi(const common::unary_function<real_t> &, const thrust_size_t &);
+      void init_count_num_const_multi(const common::unary_function<double> &);
+      void init_count_num_const_multi(const common::unary_function<double> &, const thrust_size_t &);
       void init_count_num_dry_sizes(const std::pair<real_t, int> &);
       void init_count_num_hlpr(const real_t &, const thrust_size_t &);
       void init_count_num_src(const thrust_size_t &);
diff --git a/src/impl/particles_impl_dist_analysis.ipp b/src/impl/particles_impl_dist_analysis.ipp
index ac63c974..6b8a991f 100644
--- a/src/impl/particles_impl_dist_analysis.ipp
+++ b/src/impl/particles_impl_dist_analysis.ipp
@@ -15,7 +15,7 @@ namespace libcloudphxx
     // init
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::dist_analysis_sd_conc(
-      const common::unary_function<real_t> &n_of_lnrd_stp,
+      const common::unary_function<double> &n_of_lnrd_stp,
       const impl::n_t sd_conc,
       const real_t dt
     )
@@ -78,7 +78,7 @@ namespace libcloudphxx
 
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::dist_analysis_const_multi(
-      const common::unary_function<real_t> &n_of_lnrd_stp
+      const common::unary_function<double> &n_of_lnrd_stp
     )
     {
       if(opts_init.rd_min >= 0 && opts_init.rd_max >= 0) // user-defined bin edges
diff --git a/src/impl/particles_impl_init_SD_with_distros_const_multi.ipp b/src/impl/particles_impl_init_SD_with_distros_const_multi.ipp
index 31b332ab..8f269fd9 100644
--- a/src/impl/particles_impl_init_SD_with_distros_const_multi.ipp
+++ b/src/impl/particles_impl_init_SD_with_distros_const_multi.ipp
@@ -11,7 +11,7 @@ namespace libcloudphxx
   namespace lgrngn
   {
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::init_SD_with_distros_const_multi(const common::unary_function<real_t> &fun)
+    void particles_t<real_t, device>::impl::init_SD_with_distros_const_multi(const common::unary_function<double> &fun)
     {
       // analyze the distribution, TODO: just did it
       dist_analysis_const_multi(fun);
diff --git a/src/impl/particles_impl_init_SD_with_distros_sd_conc.ipp b/src/impl/particles_impl_init_SD_with_distros_sd_conc.ipp
index fbdf45bd..c41d179e 100644
--- a/src/impl/particles_impl_init_SD_with_distros_sd_conc.ipp
+++ b/src/impl/particles_impl_init_SD_with_distros_sd_conc.ipp
@@ -11,7 +11,7 @@ namespace libcloudphxx
   namespace lgrngn
   {
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::init_SD_with_distros_sd_conc(const common::unary_function<real_t> &fun, const real_t &tot_lnrd_rng)
+    void particles_t<real_t, device>::impl::init_SD_with_distros_sd_conc(const common::unary_function<double> &fun, const real_t &tot_lnrd_rng)
     {
       // analyze the distribution, TODO: just did it in init_SD_with_distros
       dist_analysis_sd_conc(
diff --git a/src/impl/particles_impl_init_SD_with_distros_tail.ipp b/src/impl/particles_impl_init_SD_with_distros_tail.ipp
index 7d35d727..c2a2ae8b 100644
--- a/src/impl/particles_impl_init_SD_with_distros_tail.ipp
+++ b/src/impl/particles_impl_init_SD_with_distros_tail.ipp
@@ -11,7 +11,7 @@ namespace libcloudphxx
   namespace lgrngn
   {
     template <typename real_t, backend_t device>
-    void particles_t<real_t, device>::impl::init_SD_with_distros_tail(const common::unary_function<real_t> &fun, const real_t log_rd_min_init)
+    void particles_t<real_t, device>::impl::init_SD_with_distros_tail(const common::unary_function<double> &fun, const real_t log_rd_min_init)
     {
       dist_analysis_const_multi(fun);
  
diff --git a/src/impl/particles_impl_init_count_num.ipp b/src/impl/particles_impl_init_count_num.ipp
index ad3dcb41..401a375e 100644
--- a/src/impl/particles_impl_init_count_num.ipp
+++ b/src/impl/particles_impl_init_count_num.ipp
@@ -15,10 +15,10 @@ namespace libcloudphxx
       // calculate numerical integral with the trapezoidal rule
       // TODO: use thrust
       template<typename real_t>
-      real_t integrate(const common::unary_function<real_t> &fun, const real_t &min, const real_t &max, const real_t &bin_size)
+      real_t integrate(const common::unary_function<double> &fun, const real_t &min, const real_t &max, const real_t &bin_size)
       {
         const int n = (max - min) / bin_size; //no of bins
-        real_t integral = (fun(min) + fun(max)) / 2.;
+        double integral = (fun(min) + fun(max)) / 2.;
 
         for(int i=1; i<n; ++i)
           integral += fun(min + i * bin_size);
@@ -86,7 +86,7 @@ namespace libcloudphxx
 
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::init_count_num_const_multi(
-      const common::unary_function<real_t> &n_of_lnrd_stp
+      const common::unary_function<double> &n_of_lnrd_stp
     )
     {
       const real_t integral = detail::integrate(n_of_lnrd_stp, log_rd_min, log_rd_max, config.bin_precision);
@@ -95,7 +95,7 @@ namespace libcloudphxx
 
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::init_count_num_const_multi(
-      const common::unary_function<real_t> &n_of_lnrd_stp,
+      const common::unary_function<double> &n_of_lnrd_stp,
       const thrust_size_t &const_multi
     )
     {
diff --git a/src/impl/particles_impl_init_dry_const_multi.ipp b/src/impl/particles_impl_init_dry_const_multi.ipp
index 28f86cf4..c6534215 100644
--- a/src/impl/particles_impl_init_dry_const_multi.ipp
+++ b/src/impl/particles_impl_init_dry_const_multi.ipp
@@ -24,7 +24,7 @@ namespace libcloudphxx
     {
       // calculate cumulative distribution function
       template<typename real_t, typename vec_t>
-      void calc_CDF(const common::unary_function<real_t> &fun, const real_t &min, const real_t &max, const real_t &bin_size, vec_t &vec)
+      void calc_CDF(const common::unary_function<double> &fun, const real_t &min, const real_t &max, const real_t &bin_size, vec_t &vec)
       {
         const thrust_size_t n = (max - min) / bin_size + 1; //no of points at which cdf will be calculated
         vec.resize(n);
@@ -46,7 +46,7 @@ namespace libcloudphxx
 
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::init_dry_const_multi(
-      const common::unary_function<real_t> &n_of_lnrd_stp 
+      const common::unary_function<double> &n_of_lnrd_stp 
     )
     {
       // calculate cumulative distribution function
diff --git a/src/impl/particles_impl_init_n.ipp b/src/impl/particles_impl_init_n.ipp
index 50c4880a..df079e5e 100644
--- a/src/impl/particles_impl_init_n.ipp
+++ b/src/impl/particles_impl_init_n.ipp
@@ -22,14 +22,14 @@ namespace libcloudphxx
     namespace detail
     {
       template <typename real_t>
-      struct eval_and_multiply
+      struct eval_log_of_cubic_root_and_multiply
       {   
-        const common::unary_function<real_t> &fun;
+        const common::unary_function<double> &fun;
         const real_t &mul;
 
         // ctor
-        eval_and_multiply(
-          const common::unary_function<real_t> &fun, 
+        eval_log_of_cubic_root_and_multiply(
+          const common::unary_function<double> &fun, 
           const real_t &mul
         )
           : fun(fun), mul(mul)
@@ -46,7 +46,7 @@ namespace libcloudphxx
     // init
     template <typename real_t, backend_t device>
     void particles_t<real_t, device>::impl::init_n_sd_conc(
-      const common::unary_function<real_t> &n_of_lnrd_stp 
+      const common::unary_function<double> &n_of_lnrd_stp 
     )
     {
       // temporary space on the host 
@@ -72,7 +72,7 @@ namespace libcloudphxx
       thrust::transform(
         tmp_real.begin(), tmp_real.end(), // input 
         tmp_real.begin(),                 // output
-        detail::eval_and_multiply<real_t>(n_of_lnrd_stp, multiplier)
+        detail::eval_log_of_cubic_root_and_multiply<real_t>(n_of_lnrd_stp, multiplier)
       );
 
       {