diff --git a/README.md b/README.md
index 845e1a0..68a5678 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,7 @@
 # Introduction
 
 `ecx` is a Fortran library providing formulas for electrochemistry.
+C API allows usage from C, or can be used as a basis for other wrappers. Python wrapper allows easy usage from Python.
 
 To use `ecx` within your `fpm <https://github.com/fortran-lang/fpm>`_ project,
 add the following to your `fpm.toml` file:
diff --git a/src/capi.f90 b/src/capi.f90
index dc95524..37a2e3c 100644
--- a/src/capi.f90
+++ b/src/capi.f90
@@ -1,5 +1,7 @@
 module capi
     !! Main module for ECX library: C API.
     use capi__version
+    use capi__eis
+    use capi__kinetics
 
 end module
diff --git a/src/capi_eis.f90 b/src/capi_eis.f90
new file mode 100644
index 0000000..2634d6d
--- /dev/null
+++ b/src/capi_eis.f90
@@ -0,0 +1,39 @@
+module capi__eis
+    use ecx__eis
+    implicit none
+    
+subroutine capi_z(p, w, zout, e, k, n, errstat, errmsg)bind(C, name="ecx_eis_z")
+    !! Compute the complex impedance for the given element.
+    implicit none
+    integer(c_size_t), intent(in), value :: n
+        !! Size of w
+    integer(c_size_t), intent(in), value :: k
+        !! Size of p
+    character(len=1,kind=c_char), intent(in), value :: e
+        !! Electrochemical element: R, C, L, Q, O, T, G
+    real(c_double), intent(in) :: p(k)
+        !! Parameters.
+    real(c_double), intent(in) :: w(n)
+        !! Angular frequencies in rad.s-1
+    complex(c_double_complex), intent(out) :: zout(n)
+        !! Complex impedance in Ohms.
+    integer(c_int), intent(out) :: errstat
+        !! Error status
+    type(c_ptr), intent(out) :: errmsg
+        !! errmsg Error message
+    
+    character(len=:), pointer :: fptr
+
+    call z(p, w, zout, e, errstat, fptr)
+
+    if(allocated(errmsg_c))then
+        deallocate(errmsg_c)
+    endif
+    allocate(character(len=len(fptr)+1) :: errmsg_c)
+
+    errmsg_c = fptr // c_null_char
+    errmsg = c_loc(errmsg_c)
+
+end subroutine
+
+end module
diff --git a/src/ecx_kinetics_capi.f90 b/src/capi_kinetics.f90
similarity index 86%
rename from src/ecx_kinetics_capi.f90
rename to src/capi_kinetics.f90
index cc29ffc..fdde17c 100644
--- a/src/ecx_kinetics_capi.f90
+++ b/src/capi_kinetics.f90
@@ -1,5 +1,6 @@
-module ecx__kinetics_capi
-    use iso_fortran_env
+module capi__kinetics    
+    !! Kinetics: C API.
+    use stdlib_kinds, only, dp, int32
     use iso_c_binding
     use ecx__kinetics
     implicit none
@@ -7,7 +8,7 @@ module ecx__kinetics_capi
 
 contains
 
-pure function ecx_kinetics_capi_nernst(E0, z, aox, vox, nox, ared, vred, nred, T)result(E)bind(C)
+pure function capi_nernst(E0, z, aox, vox, nox, ared, vred, nred, T)result(E)bind(C)
     !! Compute the Nernst electrochemical potential in V.
     implicit none
     real(c_double), intent(in), value :: E0
@@ -33,7 +34,7 @@ pure function ecx_kinetics_capi_nernst(E0, z, aox, vox, nox, ared, vred, nred, T
 
 end function
 
-pure subroutine ecx_kinetics_capi_sbv(U, OCV, j0, aa, ac, za, zc, A, T, I, n)bind(c)
+pure subroutine capi_sbv(U, OCV, j0, aa, ac, za, zc, A, T, I, n)bind(c)
     !! Compute Butler Volmer equation without mass transport.
     ! arguments
     integer(c_size_t), intent(in), value :: n
@@ -63,7 +64,7 @@ pure subroutine ecx_kinetics_capi_sbv(U, OCV, j0, aa, ac, za, zc, A, T, I, n)bin
 
 end subroutine
 
-pure subroutine ecx_kinetics_capi_bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A, T, I, n)bind(c)
+pure subroutine capi_bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A, T, I, n)bind(c)
     !! Compute Butler Volmer equation without mass transport.
     ! arguments
     integer(c_size_t), intent(in), value :: n
@@ -74,9 +75,9 @@ pure subroutine ecx_kinetics_capi_bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A,
         !! Potential in volts.
     real(c_double), intent(in), value :: j0
         !! Exchange current density in A.cm-2
-    real(real64), intent(in), value :: jdla
+    real(4), intent(in), value :: jdla
         !! Anodic diffusion limiting current density in A.cm-2.
-    real(real64), intent(in), value :: jdlc
+    real(4), intent(in), value :: jdlc
         !! Cathodic diffusion limiting current density in A.cm-2.
     real(c_double), intent(in), value :: aa
         !! Anodic transfert coefficient.
@@ -93,7 +94,7 @@ pure subroutine ecx_kinetics_capi_bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A,
     real(c_double), intent(out) :: I(n)
         !! Current in A.
 
-    I = ecx_kinetics_bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A, T)
+    I = bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A, T)
 end subroutine
 
-end module
\ No newline at end of file
+end module
diff --git a/src/capi_version.f90 b/src/capi_version.f90
index 77beba2..8fbc504 100644
--- a/src/capi_version.f90
+++ b/src/capi_version.f90
@@ -1,5 +1,5 @@
 module capi__version
-    !! version
+    !! version: C API.
     use iso_c_binding, only: c_loc, c_ptr, c_null_char
     use ecx__version, only: get_version
     implicit none
diff --git a/src/ecx_core.f90 b/src/ecx_core.f90
index 469c8f6..a622bea 100644
--- a/src/ecx_core.f90
+++ b/src/ecx_core.f90
@@ -1,21 +1,18 @@
-!> @file
-!! @brief Core module.
-
-!> @brief Core module.
 module ecx__core
-    use iso_fortran_env
+    !! Core.
     use stdlib_kinds, only: dp, int32
     use iso_c_binding, only: c_ptr, c_loc, c_double, c_size_t
     use ieee_arithmetic
-    use codata, only: PLANCK_CONSTANT_IN_EV_HZ, &
+    use codata, only: PI_dp &,
+                        PLANCK_CONSTANT_IN_EV_HZ, &
                       SPEED_OF_LIGHT_IN_VACUUM, &
                       BOLTZMANN_CONSTANT_IN_EV_K
     use stdlib_math, only: linspace, logspace
     implicit none
     private
     
-    real(real64), parameter :: PI = 4.0d0*datan(1.0d0) !! PI
-    real(real64), parameter :: T_K=273.15d0 !! 0°C in Kelvin.
+    real(dp), parameter :: PI = 4.0d0*datan(1.0d0) !! PI
+    real(dp), parameter :: T_K=273.15d0 !! 0°C in Kelvin.
     real(dp), parameter :: kB_eV = BOLTZMANN_CONSTANT_IN_EV_K%value
     real(dp), parameter :: h_eV = PLANCK_CONSTANT_IN_EV_HZ%value 
     real(dp), parameter :: c = SPEED_OF_LIGHT_IN_VACUUM%value
@@ -35,13 +32,13 @@ module ecx__core
 pure elemental function roundn(x, n)result(r)
     !! Round x to n digits.
     implicit none
-    real(real64), intent(in) :: x
+    real(dp), intent(in) :: x
         !! Number to be rounded.
     integer(int32), intent(in) :: n
         !! Number of digits.s
-    real(real64) :: r
+    real(dp) :: r
         !! Rounded number
-    real(real64) :: fac
+    real(dp) :: fac
 
     fac = 10**n
     r = nint(x*fac, kind=kind(x)) / fac
@@ -50,18 +47,18 @@ pure elemental function roundn(x, n)result(r)
 function assertEqual(x1, x2, n)result(r)
     !! Assert if two numbers are equal.
     implicit none
-    real(real64), intent(in) :: x1
+    real(dp), intent(in) :: x1
         !! First number to be compared.
-    real(real64), intent(in) :: x2
+    real(dp), intent(in) :: x2
         !! Second number to be compared.
     integer(int32), intent(in) :: n
         !! Number of digits.
     logical :: r
         !! Comparison result.
 
-    real(real64) :: fac
-    real(real64) :: ix1
-    real(real64) :: ix2
+    real(dp) :: fac
+    real(dp) :: ix1
+    real(dp) :: ix2
     
     if(ieee_is_nan(x1) .or. ieee_is_nan(x2))then
         r = .false.
@@ -75,14 +72,14 @@ function assertEqual(x1, x2, n)result(r)
 
 pure subroutine ecx_core_linspace(start, end, x)
     !! Linear spaced 1d-array.
-    real(real64), intent(in) :: start
+    real(dp), intent(in) :: start
         !! Starting value.
-    real(real64), intent(in) :: end
+    real(dp), intent(in) :: end
         !! Ending value (included).
     integer(int32), intent(out) :: x(:)
         !! 1d-array where to put the linear spaced values.
 
-    real(real64) :: dx
+    real(dp) :: dx
     integer(int32) :: n, i
 
     n = size(x)
@@ -96,9 +93,9 @@ pure subroutine ecx_core_linspace(start, end, x)
 
 pure subroutine ecx_core_logspace(start, end, x)
     !! Log spaced 1d-array.
-    real(real64), intent(in) :: start
+    real(dp), intent(in) :: start
         !! Starting value.
-    real(real64), intent(in) :: end
+    real(dp), intent(in) :: end
         !! Ending value (included).
     integer(int32), intent(out) :: x(:)
         !! 1d-array where to put the log spaced values.
@@ -109,9 +106,9 @@ pure subroutine ecx_core_logspace(start, end, x)
 pure elemental function nm2eV(lambda)result(E)
     !! Convert wavelength to energy
     implicit none
-    real(real64), intent(in) :: lambda
+    real(dp), intent(in) :: lambda
         !! Wavelength in nm.
-    real(real64) :: E
+    real(dp) :: E
         !! Energy in eV.
 
     E = h_eV * c / (lambda*1.0d-9)
@@ -133,9 +130,9 @@ pure subroutine capi_nm2eV(lambda, E, n)bind(C, name="ecx_core_nm2eV")
 pure elemental function eV2nm(E)result(lambda)
     !! Convert wavelength to energy
     implicit none
-    real(real64), intent(in) :: E
+    real(dp), intent(in) :: E
         !! Energy in eV.
-    real(real64) :: lambda
+    real(dp) :: lambda
         !! Wavelength in nm.
 
     lambda = h_eV * c / (E * 1.0d-9)
@@ -144,9 +141,9 @@ pure elemental function eV2nm(E)result(lambda)
 pure elemental function deg2rad(theta)result(phase)
     !! Converts degrees to rad.
     implicit none
-    real(real64), intent(in) :: theta
+    real(dp), intent(in) :: theta
         !! Angle in degrees.
-    real(real64) :: phase
+    real(dp) :: phase
         !! Angle in rad.
     phase = theta * PI / 180.0d0
 end 
@@ -154,9 +151,9 @@ pure elemental function deg2rad(theta)result(phase)
 pure elemental function rad2deg(phase)result(theta)
     !! Converts degrees to rad.
     implicit none
-    real(real64), intent(in) :: phase
+    real(dp), intent(in) :: phase
         !! Angle in rad.
-    real(real64) :: theta
+    real(dp) :: theta
         !! Angle in degrees.
     theta = phase * 180.0d0 / PI
 end 
@@ -164,9 +161,9 @@ pure elemental function rad2deg(phase)result(theta)
 pure elemental function kTe(T)result(r)
     !! Compute the thermal voltage.
     implicit none
-    real(real64), intent(in) :: T
+    real(dp), intent(in) :: T
         !! Temperature in °C.
-    real(real64) :: r
+    real(dp) :: r
         !! Thermal voltage in V.
     
     r = (T+T_K) * kB_eV
diff --git a/src/ecx_eis.f90 b/src/ecx_eis.f90
index 2cc454f..20bbb4e 100644
--- a/src/ecx_eis.f90
+++ b/src/ecx_eis.f90
@@ -1,6 +1,6 @@
 module ecx__eis
-    !! EIS module
-    use iso_fortran_env
+    !! EIS module.
+    use stdlib_kinds, only: dp, int32
     use iso_c_binding, only: c_double, c_int, c_double_complex, c_size_t, c_char, c_loc, c_ptr, c_null_char
     use ieee_arithmetic, only: ieee_quiet_nan, ieee_value
     use ecx__core
diff --git a/src/ecx_kinetics.f90 b/src/ecx_kinetics.f90
index 94500a2..352b5e3 100644
--- a/src/ecx_kinetics.f90
+++ b/src/ecx_kinetics.f90
@@ -1,34 +1,34 @@
 module ecx__kinetics
     !! Module for computing kinetics using the Butler-Volmer equations.
-    use iso_fortran_env
+    use stdlib_kinds only: dp, int32
     use ecx__core
     implicit none
     private
 
-    public :: ecx_kinetics_nernst
-    public :: ecx_kinetics_sbv, ecx_kinetics_bv
+    public :: nernst
+    public :: sbv, bv
 
 contains
 
-pure function ecx_kinetics_nernst(E0, z, aox, vox, ared, vred, T)result(E)
+pure function nernst(E0, z, aox, vox, ared, vred, T)result(E)
     !! Compute the Nernst electrochemical potential in V.
     implicit none
-    real(real64), intent(in) :: E0
+    real(dp), intent(in) :: E0
         !! Standard electrochemical potential in V.
     integer(int32), intent(in) :: z
         !! Number of exchanged electrons.
-    real(real64), intent(in) :: aox(:)
+    real(dp), intent(in) :: aox(:)
         !! Activities of the oxidants.
-    real(real64), intent(in) :: vox(:)
+    real(dp), intent(in) :: vox(:)
         !! Coefficients for the oxidants.
-    real(real64), intent(in) :: ared(:)
+    real(dp), intent(in) :: ared(:)
         !! Activities of the reductants
-    real(real64), intent(in) :: vred(:)
+    real(dp), intent(in) :: vred(:)
         !! Coefficients for the reductants.
-    real(real64), intent(in) :: T
+    real(dp), intent(in) :: T
         !! Temperature in °C.
 
-    real(real64) :: E, ox, red, kTe_
+    real(dp) :: E, ox, red, kTe_
     
     kTe_ = kTe(T)
     ox = product(aox**vox)
@@ -38,67 +38,67 @@ pure function ecx_kinetics_nernst(E0, z, aox, vox, ared, vred, T)result(E)
 
 end function
 
-pure elemental function ecx_kinetics_sbv(U, OCV, j0, aa, ac, za, zc, A, T)result(I)
+pure elemental function sbv(U, OCV, j0, aa, ac, za, zc, A, T)result(I)
     !! Compute Butler Volmer equation without mass transport.
-    real(real64), intent(in) :: OCV
+    real(dp), intent(in) :: OCV
         !! Open Circuit Voltage in V.
-    real(real64), intent(in) :: U
+    real(dp), intent(in) :: U
         !! Electrochemical potential in V.
-    real(real64), intent(in) :: j0
+    real(dp), intent(in) :: j0
         !! Exchange current density in A.cm-2.
-    real(real64), intent(in) :: aa
+    real(dp), intent(in) :: aa
         !! Anodic transfer coefficient.
-    real(real64), intent(in) :: ac
+    real(dp), intent(in) :: ac
         !! Cathodic transfer coefficient.
-    real(real64), intent(in) :: za
+    real(dp), intent(in) :: za
         !! Number of exchnaged electrons in the anodic branch.
-    real(real64), intent(in) :: zc
+    real(dp), intent(in) :: zc
         !! Number of exchnaged electrons in the cathodic branch.
-    real(real64), intent(in) :: A
+    real(dp), intent(in) :: A
         !! Area in cm2.
-    real(real64), intent(in) :: T
+    real(dp), intent(in) :: T
         !! Temperature in °C.
 
-    real(real64) :: I
+    real(dp) :: I
 
-    real(real64) :: kTe_
+    real(dp) :: kTe_
     
     kTe_ = kTe(T)
 
     I = A * j0 * (exp(aa * za * (U - OCV) / kTe_) - exp(-ac * zc * (U - OCV) / kTe_));
 end function
 
-pure elemental function ecx_kinetics_bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A, T)result(I)
+pure elemental function bv(U, OCV, j0, jdla, jdlc, aa, ac, za, zc, A, T)result(I)
     !! Compute Butler Volmer equation with mass transport.
     implicit none
-    real(real64), intent(in) :: OCV
+    real(dp), intent(in) :: OCV
         !! Open Circuit Voltage in V.
-    real(real64), intent(in) :: U
+    real(dp), intent(in) :: U
         !! Electrochemical potential in V.
-    real(real64), intent(in) :: j0
+    real(dp), intent(in) :: j0
         !! Exchange current density in A.cm-2.
-    real(real64), intent(in) :: jdla
+    real(dp), intent(in) :: jdla
         !! Anodic diffusion limiting current density in A.cm-2.
-    real(real64), intent(in) :: jdlc
+    real(dp), intent(in) :: jdlc
         !! Cathodic diffusion limiting current density in A.cm-2.
-    real(real64), intent(in) :: aa
+    real(dp), intent(in) :: aa
         !! Anodic transfer coefficient.
-    real(real64), intent(in) :: ac
+    real(dp), intent(in) :: ac
         !! Cathodic transfer coefficient.
-    real(real64), intent(in) :: za
+    real(dp), intent(in) :: za
         !! Number of exchnaged electrons in the anodic branch.
-    real(real64), intent(in) :: zc
+    real(dp), intent(in) :: zc
         !! Number of exchnaged electrons in the cathodic branch.
-    real(real64), intent(in) :: A
+    real(dp), intent(in) :: A
         !! Area in cm2.
-    real(real64), intent(in) :: T
+    real(dp), intent(in) :: T
         !! Temperature in °C.
 
-    real(real64) :: I, kTe_, num, denom
+    real(dp) :: I, kTe_, num, denom
 
     kTe_ = kTe(T)
     
-    num = ecx_kinetics_sbv(U, OCV, j0, aa, ac, za, zc, 1.0d0, T)
+    num = sbv(U, OCV, j0, aa, ac, za, zc, 1.0d0, T)
     denom = 1 + j0 / jdla * exp(aa * za * (U - OCV) / kTe_) + j0 / jdlc * exp(-ac * zc * (U - OCV) / kTe_);
 
     I = A * num / denom;
diff --git a/src/ecx_pec.f90 b/src/ecx_pec.f90
index 3493ffb..e5770dd 100644
--- a/src/ecx_pec.f90
+++ b/src/ecx_pec.f90
@@ -1,9 +1,7 @@
-!> @file
-!! @brief PEC Module.
-
-!> @brief PEC Module.
 module ecx__pec
+    !! PEC.
     use iso_fortran_env
+    use stdlib_kinds, only: dp
     use ecx__core
     implicit none
     private
@@ -13,18 +11,18 @@ module ecx__pec
 contains
 
 
-!> @brief Compute the not scaled absorbance coefficient.
-!! @param[in] hv Light energy in eV.
-!! @param[in] Eg Bandgap in eV.
-!! @param[in] n Exponent for direct (1/2) or indirect transition (2)
-!! @return Absorbance coefficient in eV.
 pure elemental function alpha(hv, Eg, n)result(res)
+    !> Compute the not scaled absorbance coefficient.
     implicit none
     real(real64), intent(in) :: hv
+        !! Light energy in eV.
     real(real64), intent(in) :: Eg
+        !! Bandgap in eV.
     real(real64), intent(in) :: n
+        !! Exponent for direct (1/2) or indirect transition (2)
 
     real(real64) :: res
+        !! Absorbance coefficient in eV.
 
     real(real64) :: d
     d = hv - Eg
@@ -36,21 +34,21 @@ pure elemental function alpha(hv, Eg, n)result(res)
     endif
 end function
 
-!> @brief Compute the complex photocurrent
-!! @param[in] hv Light energy in eV.
-!! @param[in] K Scaling factor for absorbance in .
-!! @param[in] Eg Bandgap in eV.
-!! @param[in] theta Phase in degrees.
-!! @param[in] n Transition type: n=1/2 for direct transition and n=2 for indirect transition
-!! @return iph Complex photocurrent.
 pure elemental function iph(hv, K, Eg, theta, n)result(res)
+    !! Compute the complex photocurrent
     implicit none
     real(real64), intent(in) :: hv
+        !! Light energy in eV.
     real(real64), intent(in) :: K
+        !! Scaling factor for absorbance in .
     real(real64), intent(in) :: Eg
+        !! Bandgap in eV.
     real(real64), intent(in) :: theta
+        !! Phase in degrees.
     real(real64), intent(in) :: n 
+        !! Transition type: n=1/2 for direct transition and n=2 for indirect transition
     complex(real64) :: res
+        !! Complex photocurrent.
 
     real(real64) :: re, im, mod, phase