I modeled the forward kinematics of this machine:
I don't care about the Z movements (the ups and downs of the pen), just the XY. As you can see there are 2 motors aligned in the same axes driving the arms and there are 3 joints that define the constraints of the system. I decided to model the machine using vectors placing the motors on 0,0:
The vectors A,B,C,D represents the arms of the machine, and I know their module and their angles since in the forward kinematics thoose are the inputs. In order to find the position of the pen, the vector U, I can say:
$$ \overrightarrow{A} + \overrightarrow{B} + \overrightarrow{D} = \overrightarrow{U} \\ \\ \widehat{D}= -\widehat{A} $$
The second equation comes from a limitation of the system, the arm D will always be parallel to the arm A, so in vector terms, the unitary vector of D is equal to the unitary vector of A, the negative is to indicate the different direction.
Putting both together: $$ \overrightarrow{A} + \overrightarrow{B} - \lVert D \rVert \widehat{A} = \overrightarrow{X} $$
And from this, I get 2 equations that define the position of the pen X,Y:
$$ X=\lVert A \rVert*cos( \alpha )+\lVert B \rVert*cos( \beta )-\lVert D \rVert*cos( \alpha ) \\ Y=\lVert A \rVert*sin( \alpha )+\lVert B \rVert*sin( \beta )-\lVert D \rVert*sin( \alpha ) $$
Where Alpha is the angle of the arm A and Beta the angle of the arm B.
To me, it looks correct, by using the power of imagination if you lock one arm in place and rotate the other the pen should describe a circular motion, which is what I got using the equations in a spreadsheet:
Now, let's input the know values, the module of vector A is 60, and the module of vector B is 60.5, that give us 2 beautiful equations:
$$ X=60.5*\cos( \beta )-60*\cos( \alpha ) \\ Y=60.5*\sin( \beta )-60*\sin( \alpha ) $$
Now as you may know, what is really useful here is to know the angles Alpha and Beta given a combination of X and Y, or in other words, the reverse kinematics. Well, I just need to put Alpha and Beta in terms of X and Y from the 2 equations we got, simple right?. It turns out that the more I went into the equation the crazier it got, so much that I had to make another question on this page to know if I was doing something wrong. It cannot be this complicated.
Looking online I found a guy that did exactly what I'm trying but unfortunately, he did no share the equations. What he shares tho was the reasoning for finding the inverse kinematics using circles instead of vectors like I did to find the mathematical model (see the Kinematics paragraph from the last link) but to be honest, I did not understand what he did, maybe you can shine a light here.
By the way, I'm building the thing, your help will give life to a robot!
I'll take $A, B, D$ to be the lengths of the vectors of concern to simplify notation. Original equations are
$$\begin{aligned} X &= (A-D)\cos\alpha + B \cos \beta\\ Y &= (A-D)\sin\alpha + B \sin \beta \end{aligned}$$
Rearrange to solve for $\alpha,$
$$\begin{aligned} X - B \cos \beta &= (A-D)\cos\alpha\\ Y - B \sin \beta &= (A-D)\sin\alpha \end{aligned}$$
Assuming $\beta, \alpha$ are a solution to these equations, they must satisfy the nonlinear constraint
$$ (X - B \cos \beta)^2 + (Y - B \sin \beta)^2 = (A-D)^2 $$
which is computed by the squared sum of the two equations. Observe that geometrically this is the circle of radius $A-D$ centred on the point $(B\cos\beta, B\sin\beta).$
Solutions to the above constraint can be found by expanding and collecting the terms in $\beta.$ Start by expanding the constraint,
$$ X^2 - 2X B \cos\beta + B^2\cos^2\beta + Y^2 - 2YB\sin\beta + B^2\sin^2\beta= (A-D)^2 $$
Then collecting our knowns
$$ - 2X B \cos\beta - 2YB\sin\beta = (A-D)^2 - X^2 - Y^2 - B^2 $$
Dividing by common constants
$$ X \cos\beta + Y\sin\beta = \frac{X^2 + Y^2 + B^2 -(A-D)^2}{2B} $$
Now we are getting somewhere. This is, again, a very geometric equation. This is saying that the vector $(X,Y)$ has a certain dot product with respect to the vector $(\cos\beta,\sin\beta).$ I'll take the algebraic approach again, however. Let us define a helpful constant,
$$ F = \frac{X^2 + Y^2 + B^2 - (A-D)^2}{2B}, $$
which simplifies our equation to
$$ X \cos \beta + Y \sin \beta = F $$
It is known that we can reduce the left-hand side equation to the expression $R \cos(\beta + \theta)$ where $\theta, R$ are functions of the coefficients. In particular, you can verify that the reduce the equation to
$$ \sqrt{X^2 + Y^2}\cos\left( \beta + \arctan\left(-Y, X\right)\right) = F$$
You are now off to the races. $\beta$ can be solved explicitly
$$ \beta = \arccos\frac{F}{\sqrt{X^2 + Y^2}} - \arctan\left(-Y, X\right) $$
Once you have $\beta,$ $\alpha$ may be produced by our first few equations. In particular you will find
$$ \alpha = \arctan\left(Y-B\cos\beta, X-B\cos\beta\right) $$
Even if I made some mistakes, hopefully this equips you with some tools to help solve the problem! Cheers!