function [F] = F_func(x,d,p,H_plus)

W_a4 = x(2);
W_b4 = x(3);
K_a1 = p(9);
K_a2 = p(10);
K_w =  p(11);

        
F = W_b4*(K_a1/H_plus + (2*K_a1*K_a2)/H_plus^2)/(1+K_a1/H_plus+K_a1*K_a2/H_plus^2)...
    +W_a4+K_w/H_plus - H_plus;

end