How to solve a symbolic angle for an equation without complex numbers

This program takes a diameter of a nitinol coil and calculates the distance the spring would be pulled apart at a force (F_A). The problem is that using the syms function with these complicated formulas gives me complex numbers when I do vpa(solve(function==value, angle_to_be_solved_for). Also, I can't use just a solve() to get a value I can put into a tand() function to get the other values afterwards. To summarize, I need a way to solve for the angle in degrees without getting a complex number. The code in question is bolded and everything else feeds in values.

clear;
clc;

Input Values
%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio

%Wire and Coil Diameter
D=input('The Mean diameter of the Spring in mm;  ');
D=D/1000;

d_net=0;
length_f=0.121;
%Iterating final length to find a diameter for force
while length_f >= 0.01 & d_net <= 0.07857
length_f=length_f-0.001;
dl=105;
    %dl=input('The Displacement of the Spring in mm;  ');
dl=dl/1000;

angle_f=90;
g=zeros(40);
while (g<=0.08) & (angle_f>=0)
angle_f=angle_f-0.1; %iterated final angle
pitch_f=tand(angle_f)*pi*D; %Coil pitch
 %Length of unsprung spring set to forearm
length_i=length_f-dl; %Final length = initial + displacement
n=length_f/pitch_f; %Number of loops in coil
angle_i=atand(length_i/(pi*n*D));

%Shear Strain and Detwinned Martensite Fraction
d=150*(10^-6); %The diameter of the Wire (microns converted to m)

%Calculate martensite percent from range of shear strain

g=d./D.*cosd(angle_i).^2.*(sind(angle_f-(2*angle_i))+...
    sind(angle_i))./(cosd(angle_f-(2*angle_i)).^2.*...
    (cosd(angle_f-(2*angle_i)).^2.+...
    (sind(angle_f-(2*angle_i)).^2)./(1+v)));
end 

angle_range=angle_i:0.1:angle_f;
g_range=d./D.*cosd(angle_i).^2.*(sind(angle_range-(2*angle_i))+...
    sind(angle_i))./(cosd(angle_range-(2*angle_i)).^2.*...
    (cosd(angle_range-(2*angle_i)).^2.+...
    (sind(angle_range-(2*angle_i)).^2)./(1+v)));

%plot(angle_range,g_range)


%Defining critical shear stresses and variable shear stress
gcr_s=0.01;
gcr_f=0.12;

%Detwinned martensite fraction
ESg=0.5.*cos(pi.*(g_range-gcr_f)./(gcr_s-gcr_f))+0.5;
eigth=round(length(ESg)*0.01/g);
ESg(1:eigth)=0;
%plot(g_range,ESg)

%Final angles of austentite and martensite phase finish (spring range)
syms angle_A_f

%Displacement of the two metal phases for 100g hystersis force
F_A=G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
    ^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f);
F_M=G_M*d^4/(8*D^3*n)*(cosd(angle_i)^3)/(cosd(angle_f)^2*(cosd(angle_f)...
    ^2+sind(angle_f)^2/(1+v)))*dl-pi*d^3/(8*D)*G_M*ESg*0.06; %residual g
angle_A_f=(solve(F_A==.981, angle_A_f));       %MARKED
d_A=n*pi*D*tand(angle_A_f);                             %MARKED
d_net=dl-d_A;
end

disp(d_net)

Answer 1

[Alan Stevens]

You could find the angle with fzero:

%Input Values
%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio

%Wire and Coil Diameter
%D=input('The Mean diameter of the Spring in mm;  ');
D = 5;
D=D/1000;

d_net=0;
length_f=0.121;
%Iterating final length to find a diameter for force
while (length_f >= 0.01) & (d_net <= 0.07857)
    length_f=length_f-0.001;
    dl=105;
        %dl=input('The Displacement of the Spring in mm;  ');
    dl=dl/1000;

    angle_f=90;
    g=zeros(40);
    while (g<=0.08) & (angle_f>=0)
        angle_f=angle_f-0.1; %iterated final angle
        pitch_f=tand(angle_f)*pi*D; %Coil pitch
         %Length of unsprung spring set to forearm
        length_i=length_f-dl; %Final length = initial + displacement
        n=length_f/pitch_f; %Number of loops in coil
        angle_i=atand(length_i/(pi*n*D));

        %Shear Strain and Detwinned Martensite Fraction
        d=150*(10^-6); %The diameter of the Wire (microns converted to m)

        %Calculate martensite percent from range of shear strain

        g=d./D.*cosd(angle_i).^2.*(sind(angle_f-(2*angle_i))+...
            sind(angle_i))./(cosd(angle_f-(2*angle_i)).^2.*...
            (cosd(angle_f-(2*angle_i)).^2.+...
            (sind(angle_f-(2*angle_i)).^2)./(1+v)));
    end 

    angle_range=angle_i:0.1:angle_f;
    g_range=d./D.*cosd(angle_i).^2.*(sind(angle_range-(2*angle_i))+...
        sind(angle_i))./(cosd(angle_range-(2*angle_i)).^2.*...
        (cosd(angle_range-(2*angle_i)).^2.+...
        (sind(angle_range-(2*angle_i)).^2)./(1+v)));

    %plot(angle_range,g_range)


    %Defining critical shear stresses and variable shear stress
    gcr_s=0.01;
    gcr_f=0.12;

    %Detwinned martensite fraction
    ESg=0.5.*cos(pi.*(g_range-gcr_f)./(gcr_s-gcr_f))+0.5;
    eigth=round(length(ESg)*0.01/g);
    ESg(1:eigth)=0;
    %plot(g_range,ESg)

    %Final angles of austentite and martensite phase finish (spring range)
    angle_A_f0 = 30;
    angle_A_f = fzero(@(angle_A_f) findangle(angle_A_f,G_A,d,D,angle_i,v,n),angle_A_f0);
    %Displacement of the two metal phases for 100g hystersis force
    %F_A=G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
    %   ^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f);
    %F_M=G_M*d^4/(8*D^3*n)*(cosd(angle_i)^3)/(cosd(angle_f)^2*(cosd(angle_f)...
    %    ^2+sind(angle_f)^2/(1+v)))*dl-pi*d^3/(8*D)*G_M*ESg*0.06; %residual g
    %angle_A_f=(solve(F_A==.981, angle_A_f));       %MARKED
    d_A=n*pi*D*tand(angle_A_f);                             %MARKED
    d_net=dl-d_A;
end
disp(angle_A_f)
disp(d_net)

function Z = findangle(angle_A_f,G_A,d,D,angle_i,v,n)

         Z = G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
        ^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f) - 0.981;
end

[Vincent Pipitone]

It looks like the angle can be found that way, by setting the function equal to zero and solving for it. I was able to find the austentite angle by this method and it seems to make sense, thank you.

Answer 2

[Walter Roberson]

You can use vpasolve() with range [-inf inf]

Caution: if you copy this code back, fix the entry for D back to use input()... the testing facility I used does not accept input() so I used rand()

%Shear_max=560*(10^6); %yield stress from data sheet
g_max_A=172.369*(10^6); %max cyclic shear stress austentite MPa (25ksi)
g_max_M=68.9476*(10^6); %max cyclic shear stress martensite MPa (10ksi)
G_M=10.8*(10^9); %Shear Modulus of the cold Material
G_A=28.8*(10^9); %Shear Modulus of the hot Material
v=0.33; %Poisson's ratio

%Wire and Coil Diameter
%D=input('The Mean diameter of the Spring in mm;  ');

D = rand() * 10;      %<==== FOR ONLINE TESTING

D=D/1000;

d_net=0;
length_f=0.121;
%Iterating final length to find a diameter for force
while length_f >= 0.01 & d_net <= 0.07857
length_f=length_f-0.001;
dl=105;
    %dl=input('The Displacement of the Spring in mm;  ');
dl=dl/1000;

angle_f=90;
g=zeros(40);
while (g<=0.08) & (angle_f>=0)
angle_f=angle_f-0.1; %iterated final angle
pitch_f=tand(angle_f)*pi*D; %Coil pitch
 %Length of unsprung spring set to forearm
length_i=length_f-dl; %Final length = initial + displacement
n=length_f/pitch_f; %Number of loops in coil
angle_i=atand(length_i/(pi*n*D));

%Shear Strain and Detwinned Martensite Fraction
d=150*(10^-6); %The diameter of the Wire (microns converted to m)

%Calculate martensite percent from range of shear strain

g=d./D.*cosd(angle_i).^2.*(sind(angle_f-(2*angle_i))+...
    sind(angle_i))./(cosd(angle_f-(2*angle_i)).^2.*...
    (cosd(angle_f-(2*angle_i)).^2.+...
    (sind(angle_f-(2*angle_i)).^2)./(1+v)));
end 

angle_range=angle_i:0.1:angle_f;
g_range=d./D.*cosd(angle_i).^2.*(sind(angle_range-(2*angle_i))+...
    sind(angle_i))./(cosd(angle_range-(2*angle_i)).^2.*...
    (cosd(angle_range-(2*angle_i)).^2.+...
    (sind(angle_range-(2*angle_i)).^2)./(1+v)));

%plot(angle_range,g_range)


%Defining critical shear stresses and variable shear stress
gcr_s=0.01;
gcr_f=0.12;

%Detwinned martensite fraction
ESg=0.5.*cos(pi.*(g_range-gcr_f)./(gcr_s-gcr_f))+0.5;
eigth=round(length(ESg)*0.01/g);
ESg(1:eigth)=0;
%plot(g_range,ESg)

%Final angles of austentite and martensite phase finish (spring range)
syms angle_A_f

%Displacement of the two metal phases for 100g hystersis force
F_A=G_A*d^4/(8*D^3*n)*cosd(angle_i)^3/(cosd(angle_A_f)^2*(cosd(angle_A_f)...
    ^2+sind(angle_A_f)^2/(1+v)))*n*pi*D*tand(angle_A_f);
F_M=G_M*d^4/(8*D^3*n)*(cosd(angle_i)^3)/(cosd(angle_f)^2*(cosd(angle_f)...
    ^2+sind(angle_f)^2/(1+v)))*dl-pi*d^3/(8*D)*G_M*ESg*0.06; %residual g
angle_A_f=(vpasolve(F_A==.981, angle_A_f, [-inf inf]));       % <====== HERE
d_A=n*pi*D*tand(angle_A_f);
d_net=dl-d_A;
end

disp(vpa(d_net))

91.943581578108842623752027152327

[Vincent Pipitone]

Your answer helped me with the next version of the program. Should I make a new thread about it?

[Walter Roberson]

Sure