how to solve trigonometric equations

Here I have 5 equations:

eq1= DE*cos(alpha)-AC*cos(beta)-CE*cos(sigma)==-AD;
eq2 = DE*sin(alpha)-AC*sin(beta)-CE*sin(sigma)==0;
eq3 = CE*cos(alpha)+CB*cos(beta)-FG*cos(theta)-GE*sin(lambda)-BF*cos(delta)==0;
eq4= CE*sin(alpha)-CB*cos(beta)-FG*sin(theta)+GE*cos(lambda)-BF*sin(delta)==0;
eq5 = sigma - lambda == 1.815;

I've already specified other parameters as follows:

DE = 35.39;
AC = [35.51;40.51];
BF = [20.82;24.52];
CE = 4.31;
GE = 15.29;
FG = 3.5;
AD = 4.75;
CB = 6.5;

Where AC and BF are lengths that are variable with upper and lower bounds.

The only known angle will be beta which is also an input angle with range from 0 deg to 80 deg. I've tried using fsolve but this function can only deal with inputs that are constant not 'range'. Are there any methods to obtain the rest of five unknown angles as results of also a range of angle on the plot?

Something similar to this:

[Daniel Frisch]

Well, just solve the equation system several times: - AC=35.51; BF=20.82; - AC=35.51; BF=24.52; - AC=40.51; BF=20.82; - AC=40.51; BF=24.52; Then find the minimum and maximum of the results and you have your range.

[David Goodmanson]

Hello Jenny, The first two equations have only two unknown angles, so alpha and sigma can be solved independently of the other variables. However, it appears that in the range of 0 to 3pi/2, only some values of beta allow a solution for alpha and sigma. Some other values of beta do not satisfy a certain triangle inequality and do not provide a solution for alpha and sigma. Is that expected?

[Jenny Wang]

>David Goodmanson That's why I am trying to implement those constants to the code first and verifying the mechanism design. Based on journal paper, all the angles if solved, the results are supposed to be continuous within the specified range. The results of alpha, lambda + certain angle, and theta plotted verses beta should be something similar to this:

So far what I did was using fsolve to solve all the functions together but I can only use scatterplot not regular plot because the results aren't continuous...

    F =  [DE*cos(alpha)-AC*cos(beta)-CE*cos(sigma)+AD;
    DE*sin(alpha)-AC*sin(beta)-CE*sin(sigma);
    CE*cos(alpha)+CB*cos(beta)-FG*cos(theta)-GE*sin(lambda)-BF*cos(delta);
    CE*sin(alpha)-CB*cos(beta)-FG*sin(theta)+GE*cos(lambda)-BF*sin(delta);
    sigma - lambda-1.815];

[David Goodmanson]

Yes it looks like the results can be pretty strange. What is being varied to give the different plot lines in the graphs? Also, you have a CB*cos(beta) term in both eqn. 3 and eqn. 4. Should one of those be a sin(beta) term?

[Jenny Wang]

>David Goodmanson

Right...even bigger mistake I made...the equations should be as follows

F =  [DE*cos(alpha)-AC*cos(beta)-CE*cos(sigma)+AD;
  DE*sin(alpha)-AC*sin(beta)-CE*sin(sigma);
  CE*cos(alpha)+CB*cos(beta)-FG*cos(theta)-GE*sin(lambda)-BF*cos(delta);
  CE*sin(alpha)-CB*sin(beta)-FG*sin(theta)+GE*cos(lambda)-BF*sin(delta);
  sigma - lambda-1.815];

In the figures, the horizontal axis will be beta. And the varied factors that give different plot lines are the variable AC and BF, which thus gives the upper and lower bound of the plot.

Answer 1

[John D'Errico]

You have insufficient information to solve those equations.

There are 5 equations. Unknowns? There are 6 unknowns: {alpha, beta, sigma, lambda, delta, theta}.

So at best there will be infinitely many solutions, an entire manifold of solutions. Worse, two of the coefficients are themselves unknowns, as all you are willing to say is they live in intervals.

Anyone who tells you to set AC and BF to their limits, then solve the problem is simply wrong, because you can never solve the equations for a solution even if AC and BF are fixed and known. And no, interval arithmetic will also fail you here, again because no solution exists.

You can never solve this system of equations as it is. This is not a problem that can be solved EVER using fsolve, nor will solve or vpasolve succeed. Sorry.

Were you to derive a 6th independent equation, independent of the other 5, then you could start to try solving the problem.

[David Goodmanson]

From the original question:

"The only known angle will be beta which is also an input angle with range from 0 to 3pi/2"

[Jenny Wang]

>David Goodmanson Sorry, I made a mistake, beta should ranges from 0 deg to 80 deg. Don't know what I was thinking when typing the question :P

>John D'Errico Although there are six unknowns in the equations, one of them, beta, will be known, which ranges from 0 deg to 80 deg. If this range is known, will these equations possibly be solved?

[José-Luis]

There are still infinitely many points between 0 and 80.

[Jenny Wang]

>José-Luis

Yes. So as you can see in the figures above, those lines are supposed to be connected dots which are the solved results each times in my understanding...

[John D'Errico]

Yes. There are an infinite number of points. An infinite manifold of points. In fact, there will be more than just a simple surface. There will be a complete volume of points. Beta is NOT known. It ranges over a range that you will allow. But so are two of your coefficients. So you have three parameters that are allowed to range over a THREE dimensionalvolume. Fix those three to ANY number in their ranges, and you will POSSIBLY get some solution. But those three parameters mean that there is a 3-dimensional volume of solutions. That volume will live in the 5 dimensional space of the 5 unknown parameters. It is not as simple as the plot that you show.

My point is that for ANY combination of those three values, you will get 5 other parameters. So you essentially need to plot alpha as a function of (beta,AC,BF). Then plot sigma(beta,AC,BF), theta(beta,AC,BF), etc. But each of these plots is a FOUR dimensional plot.

So, you will need to loop over the possible values of AC, BF, beta. Now you can solve for the other 5 parameters. For ANY value of beta, compute the min and max of alpha. Do this for all AC and BF in their ranges. That gives you two curves, each as a function of beta. It does not give you the complete understanding of that volume. But it will give you the 2-d plot you are looking for.

[Jenny Wang]

>John D'Errico

Then I'm not really sure how they graphed those figures in this paper.(I'll be glad if you can open the file and checked if you're available)Compared to the equations shown in this paper, the 5th one was the one I derived based on my design. Probably like you suggested, they only showed the 2D image of the plot.

Answer 2

[John BG]

Hi Jenny

this is John BG <mailto:[e-mail] [e-mail]>

1.

The equations in page 409, 3rd page, of the document you have supplied, just above Figure 4, are not complete.

The authors supply the main equations, but there are missing boundary equations needed to narrow down results, for instance the rigid (constant) angle µ.

The key comment is

.

' .. a MATLAB model is created ..'

.

The equations (2) to (5) seem to be used in the model, but fur such unknown model to supply graphs 5 to 7, it has to take more equations into account that (2) to (5).

Probably a MATLAB SIMSCAPE model has been built, or with Autodesk, and Figure 4

.

% 002

.

shows precisely a couple views of the simulated hardware of model, which prodces Figures 5 to 7.

If one attempts to work just with supplied equations only, as already commented by different readers, the direct assault doesn't work.

2.

One can then try combining equations, like (2) and (3)

with AC spring extended (ABC=AC=40.51) and BF spring compressed (BF=20.82)

close all; clear all; clc

    DE = 35.39;
    % AC = [35.51;40.51];  % spring compressed 35.51 spring extended 40.51
    AC=40.51;
    AB=35.51;
    % BF = [20.82;24.52];   % again, spring
    BF=20.82;
    CE = 4.31;
    GE = 15.29;
    FG = 3.5;
    AD = 4.75;
    CB = 6.5;

    % AB=35;DE=35;CE=5;AC=40;AD=5; 

    da=1/10*pi/180*2;a=[0:da:pi/2]; 
    b=a;
    [A,B]=meshgrid(a,b);
    Y=acos(-1/CE*(AC*cos(B)-AD-DE*cos(A)))-asin(1/CE*(DE*sin(A)-AC*sin(B)));
    Yimag=imag(Y);
    hs=surf(A,B,Yimag)
    hs.EdgeColor='none'

    campos([.8 .8 5])
.

% 101

.

Now with AC spring compressed (ABC=AB=AC=35.51) and BF spring extended (BF=24.52)

.

close all; clear all; clc

DE = 35.39;
% AC = [35.51;40.51];  % spring compressed 35.51  spring extended 40.51
AC=35.51;
AB=35.51;
% BF = [20.82;24.52];   % again, spring
BF=24.52;
CE = 4.31;
GE = 15.29;
FG = 3.5;
AD = 4.75;
CB = 6.5;

% AB=35;DE=35;CE=5;AC=40;AD=5; 

da=1/10*pi/180*2;a=[0:da:pi/2]; 
b=a;
[A,B]=meshgrid(a,b);
Y=acos(-1/CE*(AC*cos(B)-AD-DE*cos(A)))-asin(1/CE*(DE*sin(A)-AC*sin(B)));
Yimag=imag(Y);
hs=surf(A,B,Yimag)
hs.EdgeColor='none'

campos([.8 .8 5])
.

% 102

.

One notes that combining (2) and (3), although there's a clear narrow diagonal strip of null imaginary values of the resulting angle sigma (just considering levers ABC and DE, joined with CE), such angles are possible, but the complete model has a far more narrower imag(sigma)=0 gap between alpha and beta, really close but not null, can be appreciated around mid range graph Figure 5.

.

% 003

.

3.

So, if in order to generate the graphs you want to build the model with Simscape https://uk.mathworks.com/products/simscape.html

a start point could from scratch, with an empty Simscape model.

Or open already existing Simscape examples like

https://uk.mathworks.com/help/physmod/sm/examples/double-pendulum-in-simulink-and-simscape-multibody.html

4.

There are people who have attempted to model similar moving linkage systems direrctly with MATLAB only, like

https://uk.mathworks.com/matlabcentral/fileexchange/24237-watt-linkage-animation?s_tid=srchtitle

attached function wattLinkage.m

attempting to modify these kind of examples can be difficult because there are no blocks to play with, all the equations have to be available.

5.

Another example, 4 linkages, with similar approach

https://uk.mathworks.com/matlabcentral/fileexchange/32294-four-bar-linkage--all-configurations-?s_tid=srchtitle

This example, attached script Four_bar_linkage.m has the same configuration (horizontal) of the section ABC and DE.

6.

Yet another of interest

https://uk.mathworks.com/matlabcentral/fileexchange/8690-linkage-mechanism--mechanical-engineering-?s_tid=srchtitle

Attached Linkage_mechanism.zip

the test lines

   xxx = mnism(3,5,6,6,8)
   xxx = animism(3,5,6,7,8,1,10)

show again a similar configuration to a partial of the sought model.

.

% 103

.

7.

Jed's Linkage solver helps build the equations

https://uk.mathworks.com/matlabcentral/fileexchange/27814-linkage-solver?s_tid=srchtitle

Linkage_solver.m

8.

Another option is Do It Yourself:

The linkage diagram to reproduce, animate and then read the angles of interest, to obtain the 3 graphs, is the following

.

% 001

.

The 1st section, related to Figure 5 graph, with 4 sample points

clear all;close all;clc;
% format long;
% stepper motor angle = input angle = beta = b: angle between beam ABC and vertical knucle plane
%
% angles that tell whether good grip can be achieved:
% Proximal Phalanx angle = alpha = a
% angle MP = beta in diagram = b : metacarpophalangeal or 1st phalanx. Spelling error in text, the word metacarphopahalangeal 

doesn't exist, page 1
% PIP = sigma in diagram = s: proximal phalanx or 2nd phaanx
% DIP = delta in diagram = d: distal phalanx or 3rd phalanx

% alpha, a: angle 

AB=35;DE=35;CE=5;AC=40;AD=5;  % beams, test values, slightly simplified from supplied figures

dv=.5;v1=[0:dv:180];circ1=[sind(v1);cosd(v1)];  % r=1 centre [0 0] reference circle points [x1 x2 x3 ..; y1 y2 y3 ..]
v2=[0:dv:360];circ2=[sind(v2);cosd(v2)];

cntrABC=[0;50];cntrDE=[0;50+AD];  % static pivot points on vertical knucle plane

% closest point that point E can get to knucle plane, with cosines trigonometric law, checked it matches calculating beta swapping 

beams not shown
% S1/sin(a1)=S2/sin(a2)=S3/sin(a3) 
% A1^2=S2^2+S3^2-2*S2*S3*cos(a1)
% translated to the diagram: S2=S3=DE S1=CE a1=b0
a0=acosd(-1/(2*DE^2)*(CE^2-2*DE^2))
cntrCE1=DE*[sind(a0);cosd(a0)]+cntrDE            %    [4.9872 89.6429]

cntrCE2=[27.77;76.31];cntrCE3=[34.78;58.96];cntrCE4=[25.6;31.13];cntrCE5=[0;20]; % sample points along circle DE radius

circAB=AB*circ1+cntrABC;
circAC=AC*circ1+cntrABC;
circDE=DE*circ1+cntrDE;

circCE1=CE*circ2+cntrCE1;  % sample circles along circle BE radius
circCE2=CE*circ2+cntrCE2;
circCE3=CE*circ2+cntrCE3;
circCE4=CE*circ2+cntrCE4;
circCE5=CE*circ2+cntrCE5;

hf=figure;ax=gca

plot(zeros(1,2),[0 120],'k--','LineWidth',1)  % knucle or vertical plane
hold all
plot(circAB(1,:),circAB(2,:),'Color',[.1 .1 .8],'LineWidth',1)  % beam ABC circle
plot(circAC(1,:),circAC(2,:),'Color',[.1 .1 .8],'LineWidth',1)

plot(circDE(1,:),circDE(2,:),'Color',[.8 .1 .1],'LineWidth',1)    % beam DE circle

plot([cntrDE(1) cntrCE1(1)],[cntrDE(2) cntrCE1(2)],'Color',[.8 .5 .5],'LineWidth',4) % beam DE, sample position
plot([cntrDE(1) cntrCE2(1)],[cntrDE(2) cntrCE2(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot([cntrDE(1) cntrCE3(1)],[cntrDE(2) cntrCE3(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot([cntrDE(1) cntrCE4(1)],[cntrDE(2) cntrCE4(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot([cntrDE(1) cntrCE5(1)],[cntrDE(2) cntrCE5(2)],'Color',[.8 .5 .5],'LineWidth',4)

plot(circCE1(1,:),circCE1(2,:),'Color',[.8 .5 .5],'LineWidth',1)  % sample circles
plot(circCE2(1,:),circCE2(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circCE3(1,:),circCE3(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circCE4(1,:),circCE4(2,:),'Color',[.8 .5 .5],'LineWidth',1)
plot(circCE5(1,:),circCE5(2,:),'Color',[.8 .5 .5],'LineWidth',1)

plot(cntrABC(1),cntrABC(2),'bo')  % beam ABC circle centre
plot(cntrDE(1),cntrDE(2),'ro')      % beam DE circle centre
text(cntrABC(1)-5,cntrABC(2),'A')  % text point A
text(cntrDE(1)-5,cntrDE(2),'D')      % text point D

plot(cntrCE1(1),cntrCE1(2),'ro')  % sample circles CE centres
plot(cntrCE2(1),cntrCE2(2),'ro')
plot(cntrCE3(1),cntrCE3(2),'ro')
plot(cntrCE4(1),cntrCE4(2),'ro')
plot(cntrCE5(1),cntrCE5(2),'ro')

% the spring pushes out, expands AC towards AD, not AD backwards to AC

% 1st sample position beta=0 then alpha is not zero
[x_CE1_AC,y_CE1_AC]=intersections(circCE1(1,:),circCE1(2,:),circAC(1,:),circAC(2,:));  % when motor input angle b=0, a is already a 

few degree ahead, the Proximal Phalanx angle in 1st graph, Y axis in 1st graph, effect of the spring pushing. a passive device like a 

spring substantially limits the grimping capacity, thus the top weight that can be gripped.
[xmin nxmin]=min(x_CE1_AC);                                                 % choose intersection point closer to knucle plane
x_CE1_AC=x_CE1_AC(nxmin);y_CE1_AC=y_CE1_AC(nxmin);
plot(x_CE1_AC,y_CE1_AC,'bd')                                                  % joint CE AC
plot(cntrCE1(1),cntrCE1(2),'go')                                                     % join CE DE
plot([x_CE1_AC cntrCE1(1)],[y_CE1_AC cntrCE1(2)],'Color',[.8 .5 .5],'LineWidth',4)  % beam CE
plot([cntrABC(1) x_CE1_AC],[cntrABC(2) y_CE1_AC],'Color',[.5 .5 .8],'LineWidth',4) % beam ABC, 1st sample position
% when beta=0 0, alpha=MP angle is not zero.
b1=0
a1=a0 % angle when all up

% 2nd sample position
[x_CE2_AC,y_CE2_AC]=intersections(circCE2(1,:),circCE2(2,:),circAC(1,:),circAC(2,:));  
[xmin nxmin]=min(x_CE2_AC); % choose intersection point closer to knucle plane
x_CE2_AC=x_CE2_AC(nxmin);y_CE2_AC=y_CE2_AC(nxmin);
plot(x_CE2_AC,y_CE2_AC,'bd')     % joint CE AC
plot(cntrCE2(1),cntrCE2(2),'go')      % joint CE DE
plot([x_CE2_AC cntrCE2(1)],[y_CE2_AC cntrCE2(2)],'Color',[.8 .5 .5],'LineWidth',4)  % beam CE
plot([cntrABC(1) x_CE2_AC],[cntrABC(2) y_CE2_AC],'Color',[.5 .5 .8],'LineWidth',4) % beam ABC, 1st sample position
b2=atand((x_CE2_AC-cntrABC(1))/(y_CE2_AC-cntrABC(2)))
a2=atand((cntrCE2(1)-cntrDE(1))/(cntrCE2(2)-cntrDE(2)))

b2-(a2+a0)  % just shy of half degree different alpha-beta, when beta ~ 45° matching green blue graph lines in 1st graph page 4
%  =
%    -0.4198

% 3rd sample position
[x_CE3_AC,y_CE3_AC]=intersections(circCE3(1,:),circCE3(2,:),circAC(1,:),circAC(2,:));  
[xmin nxmin]=min(x_CE3_AC); % choose intersection point closer to knucle plane
x_CE3_AC=x_CE3_AC(nxmin);y_CE3_AC=y_CE3_AC(nxmin);
plot(x_CE3_AC,y_CE3_AC,'bd')     % joint CE AC
plot(cntrCE3(1),cntrCE3(2),'go')      % joint CE DE
plot([x_CE3_AC cntrCE3(1)],[y_CE3_AC cntrCE3(2)],'Color',[.8 .5 .5],'LineWidth',4)  % beam CE
plot([cntrABC(1) x_CE3_AC],[cntrABC(2) y_CE3_AC],'Color',[.5 .5 .8],'LineWidth',4) % beam ABC, 1st sample position
b3=atand((x_CE3_AC-cntrABC(1))/(y_CE3_AC-cntrABC(2)))
a3=atand((cntrCE3(1)-cntrDE(1))/(cntrCE3(2)-cntrDE(2)))
b3-(a3+a0) 

hf.Position=[200 200 1000 600]
ax.DataAspectRatio=[1 1 1]
ax.XLim=[-10 150];
ax.YLim=[0 100];
ax.PlotBoxAspectRatio=[1 1 1]

.

% 010

.

gripping with beta>90° not plotted in document, it may either damage some fragile objects or the spring may reach full run and not pushing so no grip, too small items may drop

.

% MIP
% EG=15 but because in linkage diagram, for visual convenience EG length
% similar to DE length, I am choosing EG=DE
EG=DE
% BF=[20 24]  % again BF has another spring, up to 24
BF1=40
%FG=4   % short beam, again from linkage diagram, despite FG value < CE, in diagram look same size so FG=CE
FG=CE
% rigid angle u not supplied, but from linkage diagram, counting pixels,
% triangle EQG complementary angle to u is  ~ 17° thus u=107°
u=109

% [x_CE1_AC,y_CE1_AC] rotate CW around [cntrCE1(1),cntrCE1(2)]
G01=[cosd(-u) -sind(-u);sind(-u) cosd(-u)]*[x_CE1_AC-cntrCE1(1);y_CE1_AC-cntrCE1(2)]+[cntrCE1(1);cntrCE1(2)]
plot(G01(1),G01(2),'go')
m1=(G01(2)-cntrCE1(2))/(G01(1)-cntrCE1(1))
G1=[cntrCE1(1)+EG*cosd(atand(m1));cntrCE1(2)+EG*sind(atand(m1))]
plot([cntrCE1(1) G1(1)],[cntrCE1(2) G1(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(G1(1),G1(2),'ro')             % 1st sample point G acquired
% BF also has another spring pushing forward, but there's no spring
% compression point like B on ABC, assuming BF fully extended

% [x_CE2_AC,y_CE2_AC] rotate CW around [cntrCE2(1),cntrCE2(2)]
G02=[cosd(-u) -sind(-u);sind(-u) cosd(-u)]*[x_CE2_AC-cntrCE2(1);y_CE2_AC-cntrCE2(2)]+[cntrCE2(1);cntrCE2(2)]
plot(G02(1),G02(2),'go')
m1=(G02(2)-cntrCE2(2))/(G02(1)-cntrCE2(1))
G2=[cntrCE2(1)+EG*cosd(atand(m1));cntrCE2(2)+EG*sind(atand(m1))]
plot([cntrCE2(1) G2(1)],[cntrCE2(2) G2(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(G2(1),G2(2),'ro')  

% [x_CE3_AC,y_CE3_AC] rotate CW around [cntrCE3(1),cntrCE3(2)]
G03=[cosd(-u) -sind(-u);sind(-u) cosd(-u)]*[x_CE3_AC-cntrCE3(1);y_CE3_AC-cntrCE3(2)]+[cntrCE3(1);cntrCE3(2)]
plot(G03(1),G03(2),'go')
m1=(G03(2)-cntrCE3(2))/(G03(1)-cntrCE3(1))
G3=[cntrCE3(1)+EG*cosd(atand(m1));cntrCE3(2)+EG*sind(atand(m1))]
plot([cntrCE3(1) G3(1)],[cntrCE3(2) G3(2)],'Color',[.8 .5 .5],'LineWidth',4)
plot(G3(1),G3(2),'ro')  

ax.YLim=[0 140];

.

% 011

.

% calculate 1st set of sigmas s, mCE: slope section between points C and E
% reviewing the diagram:  if m_angle<0 s=90°+abs(m_angle); if m_angle>0 s=90-m_angle+180

a1CE=atand((cntrCE1(2)-y_CE1_AC)/(cntrCE1(1)-x_CE1_AC))
if a1CE<0 && abs(a1CE)<90 s1=90+abs(a1CE); end
if a1CE<0 && abs(a1CE)==90 s1=90; end
if a1CE>0 s1=180+(90-abs(a1CE)); end

a2CE=atand((cntrCE2(2)-y_CE2_AC)/(cntrCE2(1)-x_CE2_AC))
if a2CE<0 && abs(a2CE)<90 s2=90+abs(a2CE); end
if a2CE<0 && abs(a2CE)==90 s2=90; end
if a2CE>0 s2=180+(90-abs(a2CE)); end

a3CE=atand((cntrCE3(2)-y_CE3_AC)/(cntrCE3(1)-x_CE3_AC))
if a3CE<0 && abs(a3CE)<90 s3=90+abs(a3CE); end
if a3CE<0 && abs(a3CE)==90 s3=90; end
if a3CE>0 s3=180+(90-abs(a3CE)); end
s1
s2
s3

s1 =
   94.0961
s2 =
  153.5055
s3 =
  218.0812

% now going for 1st set of F points

    circFG1=FG*circ2+G1;     % small beam circles
    circFG2=FG*circ2+G2; 
    circFG3=FG*circ2+G3; 
    plot(circFG1(1,:),circFG1(2,:),'Color',[.8 .5 .5],'LineWidth',1) 
    plot(circFG2(1,:),circFG2(2,:),'Color',[.8 .5 .5],'LineWidth',1) 
    plot(circFG3(1,:),circFG3(2,:),'Color',[.8 .5 .5],'LineWidth',1) 

    v3=[0:dv:45];v4=[45:dv:130];v5=[135:dv:170];  % avoiding overcrowding schematic, now just arcs
    circ3=[sind(v3);cosd(v3)];   % support arcs
    circ4=[sind(v4);cosd(v4)];
    circ5=[sind(v5);cosd(v5)];

    B1=AB*[sind(b1);cosd(b1)]+cntrABC   % sample B points
    plot(B1(1),B1(2),'b*')
    B2=AB*[sind(b2);cosd(b2)]+cntrABC
    plot(B2(1),B2(2),'b*')
    B3=AB*[sind(b3);cosd(b3)]+cntrABC
    plot(B3(1),B3(2),'b*')

    circBF1=BF1*circ3+B1             % support circles
    circBF2=BF1*circ4+B2
    circBF3=BF1*circ5+B3
    plot(circBF1(1,:),circBF1(2,:),'Color',[.7 .7 .7])
    plot(circBF2(1,:),circBF2(2,:),'Color',[.7 .7 .7])
    plot(circBF3(1,:),circBF3(2,:),'Color',[.7 .7 .7])

    [x_FG1_BF,y_FG1_BF]=intersections(circBF1(1,:),circBF1(2,:),circFG1(1,:),circFG1(2,:));  
    [ymax nymax]=max(y_FG1_BF);                                                 % choose intersection point with higher y
    x_FG1_BF=x_FG1_BF(nymax);y_FG1_BF=y_FG1_BF(nymax);        % F point acquired
    plot(x_FG1_BF,y_FG1_BF,'bd')                                                  % joint FG BF
    plot([x_FG1_BF B1(1)],[y_FG1_BF B1(2)],'Color',[.5 .5 .8],'LineWidth',4)  % beam BF

.

% 012

.

I have exhausted the time allocated to your question without feed back from you yet.

The above lines attached test_123.m

Jenny would you please comment how would you like to proceed in order to obtain graphs 5 to 7?

Do you take it from here, or would you like to see how to complete the points of interest, assess angles, and compare to graph?

Therefore

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>