% Version 1.0 released by David Romps on September 12, 2017.
% Version 1.1 vectorized lcl.R, released on May 24, 2021.
% 
% When using this code, please cite:
% 
% @article{16lcl,
%   Title   = {Exact expression for the lifting condensation level},
%   Author  = {David M. Romps},
%   Journal = {Journal of the Atmospheric Sciences},
%   Year    = {2017},
%   Month   = dec,
%   Number  = {12},
%   Pages   = {3891--3900},
%   Volume  = {74}
% }
%
% This lcl function returns the height of the lifting condensation level
% (LCL) in meters.  The inputs are:
% - p in Pascals
% - T in Kelvins
% - Exactly one of rh, rhl, and rhs (dimensionless, from 0 to 1):
%    * The value of rh is interpreted to be the relative humidity with
%      respect to liquid water if T >= 273.15 K and with respect to ice if
%      T < 273.15 K. 
%    * The value of rhl is interpreted to be the relative humidity with
%      respect to liquid water
%    * The value of rhs is interpreted to be the relative humidity with
%      respect to ice
% - return_ldl is an optional logical flag.  If true, the lifting deposition
%   level (LDL) is returned instead of the LCL. 
% - return_min_lcl_ldl is an optional logical flag.  If true, the minimum of the
%   LCL and LDL is returned.
%
% Since MATLAB does not allow optional input variables or default values
% for input variables, this function has been written to only allow the 
% input of rhl.  If you like unsolicited advice, do not use MATLAB.  R
% and Python are free and better.

if (abs(lcl(1e5,300,.5,false,false)/( 1433.844139279)-1) < 1e-10) && ...
   (abs(lcl(1e5,200,.5,false,false)/( 542.8017712435)-1) < 1e-10) && ...
   (abs(lcl(1e5,300,.5,true ,false)/( 1639.249726127)-1) < 1e-10) && ...
   (abs(lcl(1e5,200,.5,true ,false)/(-8.609834216556)-1) < 1e-10);
   disp('Success')
else
   disp('Failure')
end
