function ci = cbrt(r)
%cube root calculator
%   - - - AUGI - - -
%   2007-12-26

%   Cube Root algorithm
%   This M-file will compute the "simple" cube root of any number. Since
%   MATLAB uses a combination of logs and antilogs to compute [non-integer]
%   powers of numbers, this algorithm will do it. This algorithm will
%   computer the cuberoot by analysing the angle of the original number in
%   the complex plane, and find the first cube root closest to that angle
%   as the solution. If the solution and the number have the same angle in
%   the complex plane, then that is the root choosen.


if (imag(r)==0)
    %   simple compute of negative and positive cube roots
    ci = sign(r)*(abs(r))^(1/3);
else
    %   RECTANGULAR TO POLAR
    len = abs(r);
    phi = angle(r);
    %  angle must be from 0 to 2pi
    if (phi < 0);
        phi = (2*pi) + phi;
    end

    l = len^(1/3);% radius of the cuberoot
    a1 = phi/3;% first angle of the cuberoot
    a2 = a1 + (2/3)*pi;%    second angle of cuberoot
    a3 = a1 + (4/3)*pi;%    third angle cuberoot
    %   ! angles larger than 2pi
    %   a1 cannot be larger than 2pi, phi is always between 0 and 2pi
    %   a2 cannot be larger than 2pi, a1 will always be less than (2/3)pi
    %   a3 cannot be larger than 2pi, a1 will always be less than (2/3)pi
    
    %   determines which angle is closest to phi
    if (abs(phi-a1)>abs(phi-a2))
        a1 = a2;
    end
    if (abs(phi-a1)>abs(phi-a3))
        a1 = a3;
    end
    ci = l*(cos(a1) + sqrt(-1)*sin(a1));
end

function xi = cardano(a, b, c, d)
%zeros soultion cubic 3rd third
%   - - - AUGI - - -
%   2007-12-25

%   Cardano Cubic Equation Solver
%   This M-file will solve the roots of a cubic function given the function
%   is of the form: f(x) = a*x^3 + b*x^2 + c*x + d
%IMPORTANT: this program is not tested, nor is expected to function with
%   arrays, lists, or matrices. it expects numbers.
%WARNING: This mathematical algorithm has proved to be inaccurate in
%   certain cases. The fault has been attributed to the method MATLAB uses
%   to compute non-integer powers by using logs and antilogs. For any n,
%   three distinct cube-roots exist, like for any n, there are two distinct
%   square-roots.


%   BEGIN CARDANO ALGORITHM
tol = 10^-6;%   PROGRAM TOLERANCE, DO NOT CHANGE!
j = sqrt(-1);

q = ( (3*a*c) - (b^2)) / (9 * (a^2) );
r = ( (9*a*b*c) - (27*(a^2)*d) - 2*(b^3)) / (54*(a^3));

s = cbrt(r + sqrt(q^3 + r^2));
t = cbrt(r - sqrt(q^3 + r^2));

%   SOLUTIONS
x1 = s + t - (b / (3*a));
x2 = -(1/2)*(s+t) - (b / (3*a));
x3 = -(1/2)*(s+t) - (b / (3*a));
if (abs(s-t) > tol);
    x2 = x2 + (sqrt(3)/2)*(s-t)*j;
    x3 = x3 - (sqrt(3)/2)*(s-t)*j;
end

xi = sort([x1 x2 x3]);

%   END CARDANO ALGORITHM

%strain_analysis.m
%read write excel xls xlsx

%   - - - AUGI - - -
%   2007-12-25
%   - - - CLASS - - - / - - - PROFESSOR - - -
%   Strain and Deflection Lab analysis

%   This utility will analyze data from an Excel spreadsheet that is beyond
%   the capabilities of MS Excel. The data in question is the transverse
%   deflection of a bar when the length to the measured transverse
%   deflection is unknown.

%   This M-file utilizes the ability to read and write to Excel Documents.


clc;%   clear command window
disp('-------------------- BEGIN strain_analysis --------------------')

d = [0;0];% variable for clamp-gage length, list

%   SET DOCUMENT TO BE ANALYZED, INCLUDE PATH
xlsx = 'D:\PATH\file.xlsx';
disp('Document:');disp(xlsx);disp(' ')


disp('BEGIN COLLECTING CONSTANTS & DATA FROM DOCUMENT')
l = xlsread(xlsx, 'Sheet1', 'A1');% Length [weight to gage] (mm)
E = xlsread(xlsx, 'Sheet1', 'A4');% Modulus of Elasticity (N/mm^2)
I = xlsread(xlsx, 'Sheet1', 'A6');% Area Moment of Inertia (mm^4)
P = xlsread(xlsx, 'Sheet1', 'E12:E23');%    Load (N), list
v = xlsread(xlsx, 'Sheet1', 'D12:D23');%    [Transverse] Deflection (mm), list
disp('----- END COLLECTING CONSTANTS & DATA FROM DOCUMENT');disp(' ')


disp('BEGIN DATA ANALYSIS')
%   (1)   v = (-P / 6EI) x^2 (3(l+x) - x)
%   0 = (2)x^3 + 3lx^2 + 0x + 6EIv/P
%   0 = a3 x^3 + a2 x^2 + a1 x^1 + a0 x^0
a3 = 2;
a2 = 3*l;
a1 = 0;
for n=1:12;
    a0 = 6*E*I*(v(n,1)) / (P(n,1));
    xi = cardano(a3,a2,a1,a0);
    %   WARNING: As the Cardano alrogithm has proven unreliable with
    %   certain values, strain_analysis has tested the Cardano
    %   Algorithm results. Each positive result of clamp-gage length has
    %   been shown to match equation (1), line 37
    for t=1:3;% remove impossible values
        if (imag(xi(t)) ~= 0);
            xi(t) = -1;
        end%    removes all complex solutions
        if (real(xi(t)) > 100);
            xi(t) = -1;
        end%    removes values longer than clamp-gage length, 100(mm)
        if (real(xi(t)) < 0);
            xi(t) = -1;
        end%  removes negative lengths
    end
    xi = sort(xi);
    d(n,1) = xi(3);
end
disp('----- END DATA ANALYSIS');disp(' ')


xlswrite(xlsx, d, 'Sheet1', 'I12:I23');
disp('Values written to:');disp(xlsx)


disp('-------------------- END strain_analysis --------------------')