dW1_fit = input('W1 slope: ');
dW2_fit = input('W2 slope: ');
dO1_fit = input('O1 slope: ');
dO2_fit = input('O2 slope: ');
        theta1_meas =  input('theta_1: ');
        theta2_meas =  input('theta_2: ');
        theta1 = 54.736*pi/180;
        theta2 = -54.736*pi/180;
                theta1_meas = theta1_meas*pi/180; % An uncalibrated measure of the angle of slope1
        theta2_meas = theta2_meas*pi/180; % An uncalibrated measure of the angle of slope2
        
        s1_meas = tan(theta1_meas); % the uncalibrated slope1
        s2_meas = tan(theta2_meas); % the uncalibrated slope2

        gamma = ((s1_meas -s2_meas) + sqrt((s1_meas-s2_meas)^2 - 4*s1_meas*s2_meas*tan(theta1-theta2)^2))/(2*s1_meas*s2_meas*tan(theta1-theta2));

        disp(['The z-piezo correction factor gamma = ',num2str(gamma)]);

        tilt_angle = atan(gamma*s1_meas)-theta1;
        % This is the tilt angle of the surface normal with respect to the Z axis.
        % If the sample was perfectly flat, the tilt angle = 0.
        disp(['The sample tilt angle = ',num2str(tilt_angle*180/pi)]);

        theta1_true = atan(tan(theta1_meas)*gamma);
        theta2_true = atan(tan(theta2_meas)*gamma);
        % These are the true angles of the slopes (corrected)
        disp(['The measured theta1 = ',num2str(theta1_meas*180/pi)]);
        disp(['The true theta1 = ',num2str(theta1_true*180/pi)]);
        disp(['The measured theta2 = ',num2str(theta2_meas*180/pi)]);
        disp(['The true theta2 = ',num2str(theta2_true*180/pi)]);

        p = dW2_fit/dW1_fit;
        q = (dO1_fit-dO2_fit)/dW2_fit;
        
        for mu1 = 0.05:0.001:0.6
                     
            k = p*mu1/(cos(theta1_true)^2 - (mu1^2)*sin(theta1_true)^2);
            % k is the kappa from the 1996 paper
            % To solve for mu2 in terms of mu1, use the quadratic formula.
            a2 = k*sin(theta2_true)^2;
            b2 = 1;
            c2 = (-k*sin(2*theta2_true)^2)/(4*sin(theta2_true)^2);
            
            mu2 = (-b2 + sqrt(b2^2 - 4*a2*c2))/(2*a2);
            % This equation is in terms of mu1 (symbolic variable)                    
    
            gg = 0.5*((1/mu1 + mu1)*sin(2*theta1_true)*1/p - (1/mu2 + mu2)*sin(2*theta2_true)) - q;
            % This equation has mu1 and mu2 in it, so mu2 needs to be substituted with
            % mu2(mu1) from above. Also, all the other parameters (p, q,
            % theta1_true,...) have to be substituted as well.
            
            if gg < 0
                break
            end
        end
              
        mu1_found = mu1;
        disp(['The friction coefficient on slope 1 = ',num2str(mu1_found)]);
        
        mu2_found = round(mu2*1000)/1000;
        disp(['The friction coefficient on slope 2 = ',num2str(mu2_found)]);
        
        S = mu1_found/(dW1_fit*(cos(theta1_true)^2-(mu1_found^2)*sin(theta1_true)^2));
        disp(['The S (alpha/beta) conversion factor = ',num2str(S)]);       

        disp(sprintf(['\nTo convert your lateral signal T[V] to nN, here''s the equation:\n T[V]*S*k_normal[nN/nm]*defl_sens[nm/V] = T[nN]']));