function r = solve_cos(a, b, k, h)
%%
% r = solve_cos(a, b, k, h)
% 求解方程:-k*r + h = a * cos(b*r), r >= 0 的第一个正数解
%
% 输入:
% a 和 b 是余弦参数, a > 0, b > 0
% k 是直线斜率的绝对值,k > 0
% h 是直线的截距, h > 0
%
% 返回值:
% r 是直线与余弦函数的第一个交点
% 精度有k决定,k越大精度越高,k越小精度越低
%%
r0 = (h - a) / k; %直线z = -k*r + h与z = a的交点横坐标
T = 2 * pi / b; %余弦函数的周期
n = floor(r0 / T); %交点位于第n个周期中
r1 = n * T; %第n个周期的起点(波峰)横坐标
r2 = r1 + T / 2; %第n个周期的中点(波谷)横坐标
r3 = r1 + T; %第n个周期的终点(波峰)横坐标
if k < a*b %存在多解情况
r4 = (asin(k / a / b) + 2 * pi * n) / b;
%与直线平行与余弦函数相切的直线与余弦函数的切点横坐标
z4 = a * cos(b * r4); %切点纵坐标
h4 = k * r4 + z4; %切线截距
if h4 > h %直线交于余弦左臂,有多个交点
r = Iteration(a, b, k, h, r0);
%迭代法求解
elseif h4 < h %直线交于余弦右臂
r = Dichotomy(a, b, k, h, r2, r3);
%二分法求解
else %直线与余弦相切
r = r0;
end
else %只有唯一解
if k * r2 - h > a %直线交于余弦左臂
r = Dichotomy(a, b, k, h, r2, r3);
%二分法求解
elseif k * r2 - h < a %直线交于余弦右臂
r = Iteration(a, b, k, h, r0);
%迭代法求解
else %直线交于余弦谷底
r = r2;
end
end