<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introdcution
This trial describes how to create edge and screw dislocations in iron BCC strating with one unitcell containing two atoms  
Background
The elastic solution for displacement field of dislocations is provided in the paper Dislocation Displacement Fields in Anisotropic Media.
Theoritical
The paper mentioned in backgroud subsection deals with only one dislocation. Here we describe how to extend the solution to periodic array of dislocations. Since we are dealing with linear elasticity we can superpose (sum up) the displacement field of all the individual dislocations. Looking at the Eqs. (2-8)  of abovementioned reference this boils done to finding a closed form soloution for 
$$\sum_{m=-\infty}^{\infty} \log\left(z-ma \right).$$
Where $z= x+yi$ and $a$ is a real number, equivakent to $\mathbf{H}_{00}$ that defines the periodicity of dislocations on x direction. 
Let us simplify the problem a bit further. Since this is the component displacement field we can add or subtract constant term so for each $\log\left(z-ma \right)$ we subtract a factor of $log\left(a \right)$, leading to
$$\sum_{m=-\infty}^{\infty} \log\left(\frac{z}{a}-m \right).$$
Lets change $z/a$ to $z$ and when we arrive the solution we will change ot back
$$\sum_{m=-\infty}^{\infty} \log\left(z-m \right).$$
Objective is to find a closed form solution for
$$f\left(z\right)=\sum_{m=-\infty}^{\infty} \log\left(z-m \right).$$
First note that
$$
f'\left(z\right)=\frac{1}{z}+\sum_{m=1}^{\infty}\frac{1}{z-m}+\frac{1}{z+m},
$$
and also
$$
\frac{1}{z\mp m}=\mp \frac{1}{m}\sum_{n=0}^{\infty}
\left(\pm \frac{z}{m}\right)^n.
$$
This leads to
$$
\frac{1}{z-m}+\frac{1}{z+m}=-\frac{2}{z}\sum_{n=1}^{\infty}\left(\frac{z}{m}\right)^{2n},
$$
and subsequently
$$
f'\left(z\right)=\frac{1}{z}-\frac{2}{z}\sum_{n=1}^{\infty}\left(z\right)^{2n}\sum_{m=1}^{\infty}m^{-2n},
$$
$$
=\frac{1}{z}-\frac{2}{z}\sum_{n=1}^{\infty}\left(z\right)^{2n}\zeta\left(2n\right).
$$
Where $\zeta$ is Riemann zeta function. Since $\zeta\left(0\right)=-1/2$, it simplifies to
Step1:  Block the output of all cores except for one
Step2:  Define an md.export_cfg object
md.export_cfg has a call method that we can use to create quick snapshots of our simulation box
Step3:  Screw dislocation
Step4:  Create a $\langle110\rangle\times\langle112\rangle\times\frac{1}{2}\langle111\rangle$ cell
create a $\langle110\rangle\times\langle112\rangle\times\langle111\rangle$ cell
Since mapp4py.md.atoms.cell_chenge() only accepts integer values start by creating a $\langle110\rangle\times\langle112\rangle\times\langle111\rangle$ cell
Step5:  Remove half of the atoms and readjust the position of remaining
Now one needs to cut the cell in half in $[111]$ direction. We can achive this in three steps
Step6:  Readjust the postions
Step7:  Replicating the unit cell
Step8:  putting it all together
Step9:  Edge dislocation
Step10:  putting it all together
Step11:  Putting it all together
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
import mapp4py
from mapp4py import md
from lib.elasticity import rot, cubic, resize, displace, HirthEdge, HirthScrew
Explanation: Introdcution
This trial describes how to create edge and screw dislocations in iron BCC strating with one unitcell containing two atoms  
Background
The elastic solution for displacement field of dislocations is provided in the paper Dislocation Displacement Fields in Anisotropic Media.
Theoritical
The paper mentioned in backgroud subsection deals with only one dislocation. Here we describe how to extend the solution to periodic array of dislocations. Since we are dealing with linear elasticity we can superpose (sum up) the displacement field of all the individual dislocations. Looking at the Eqs. (2-8)  of abovementioned reference this boils done to finding a closed form soloution for 
$$\sum_{m=-\infty}^{\infty} \log\left(z-ma \right).$$
Where $z= x+yi$ and $a$ is a real number, equivakent to $\mathbf{H}_{00}$ that defines the periodicity of dislocations on x direction. 
Let us simplify the problem a bit further. Since this is the component displacement field we can add or subtract constant term so for each $\log\left(z-ma \right)$ we subtract a factor of $log\left(a \right)$, leading to
$$\sum_{m=-\infty}^{\infty} \log\left(\frac{z}{a}-m \right).$$
Lets change $z/a$ to $z$ and when we arrive the solution we will change ot back
$$\sum_{m=-\infty}^{\infty} \log\left(z-m \right).$$
Objective is to find a closed form solution for
$$f\left(z\right)=\sum_{m=-\infty}^{\infty} \log\left(z-m \right).$$
First note that
$$
f'\left(z\right)=\frac{1}{z}+\sum_{m=1}^{\infty}\frac{1}{z-m}+\frac{1}{z+m},
$$
and also
$$
\frac{1}{z\mp m}=\mp \frac{1}{m}\sum_{n=0}^{\infty}
\left(\pm \frac{z}{m}\right)^n.
$$
This leads to
$$
\frac{1}{z-m}+\frac{1}{z+m}=-\frac{2}{z}\sum_{n=1}^{\infty}\left(\frac{z}{m}\right)^{2n},
$$
and subsequently
$$
f'\left(z\right)=\frac{1}{z}-\frac{2}{z}\sum_{n=1}^{\infty}\left(z\right)^{2n}\sum_{m=1}^{\infty}m^{-2n},
$$
$$
=\frac{1}{z}-\frac{2}{z}\sum_{n=1}^{\infty}\left(z\right)^{2n}\zeta\left(2n\right).
$$
Where $\zeta$ is Riemann zeta function. Since $\zeta\left(0\right)=-1/2$, it simplifies to:
$$
f'\left(z\right)=-\frac{2}{z}\sum_{n=0}^{\infty}\left(z\right)^{2n}\zeta\left(2n\right)
$$
Note that
$$
-\frac{\pi z\cot\left(\pi z\right)}{2}=\sum_{n=0}^{\infty}z^{2n} \zeta\left(2n\right)
$$
I have no idea how I figured this out but it is true. Therefore,
$$
f'\left(z\right)=\pi\cot\left(\pi z\right).
$$
At this point one can naively assume that the problem is solved (like I did) and the answer is something like:
$$
f\left(z\right)=\log\left[\sin\left(\pi z\right)\right]+C,
$$
Where $C$ is a constant. However, after checking this against numerical vlaues you will see that this is completely wrong. 
The issue here is that startegy was wrong at the very begining. The sum of the displacelment of infinte dislocations will not converge since we have infinite discountinuity in displacement field. In other words they do not cancel each other they feed each other.
But there is still a way to salvage this. Luckily, displacement is relative quantity and we are dealing with crystals. We can easily add a discontinuity in form an integer number burger vectors to a displacement field and nothing will be affected. 
So here is the trick: We will focus only on the displacement field of one unit cell dislocation (number 0). At each iteration we add two dislocation to its left and right. 
At $n$th iterations we add a discontinuity of the form
$$
-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] \pi i
$$
and a constant of the form:
$$
-2\log n.
$$
In other words and we need to evaluate: 
$$
\lim_{m\to\infty}\sum_{n=-m}^{m}
\biggl{
\log\left(z-n\right)
-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] \pi i 
-2\log\left(n \right)
\biggr} + \pi,
$$
which simplifies to 
$$
\lim_{m\to\infty}\sum_{n=-m}^{m}\log\left(z-n\right)
-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] m \pi i 
-2\log\left(\frac{m!!}{\sqrt{\pi}} \right)
$$
Note that we added an extra $\pi$ to displacement field for aesthetic reasons. After a lot of manipulations and tricks (meaning I dont't remember how I got here) we arrive at the following relation:
$$
\lim_{m\to\infty}\sum_{n=-m}^{m}\log\left(z-n\right)
-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] m \pi i 
-2\log\left(\frac{m!!}{\sqrt{\pi}} \right)=\log\left[\sin\left(\pi z\right)\right]
$$
However, this is only valid when 
$$-1/2 \le\mathrm{Re}\left(z\right)\lt 1/2.$$ 
If one exceeds this domain the answer is:
$$
\boxed{
\log\left[\sin\left(\pi z\right)\right]-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right]\left \lceil{\mathrm{Re}\left(\frac{z}{2}\right)}-\frac{3}{4}\right \rceil 2 \pi i
}
$$
Where $\lceil . \rceil$ is the cieling function. Of course there is probably a nicer form. Feel free to derive it
Final formulation
To account for peridicity of dislocations in $x$ direction, the expression $\log\left(z\right)$ in Eqs(2-7) of the paper, it should be replaced by:
$$\lim_{m\to\infty}\sum_{n=-m}^{m}\log\left(z-na\right)
-\mathrm{Sign}\left[\mathrm{Im}\left(z\right)\right] m \pi i 
-2\log\left(\frac{m\,\,!!}{\sqrt{\pi}} \right),$$
which has the closed form:
$$
\boxed{
\log\left[\sin\left(\pi\frac{z}{a}\right)\right]-\mathrm{Sign}\left[\mathrm{Im}\left(\frac{z}{a}\right)\right]\left \lceil{\mathrm{Re}\left(\frac{z}{2a}\right)}-\frac{3}{4}\right \rceil 2 \pi i.
}
$$
Preperation
Import packages
End of explanation
from mapp4py import mpi
if mpi().rank!=0:
    with open(os.devnull, 'w') as f:
        sys.stdout = f;
Explanation: Block the output of all cores except for one
End of explanation
xprt = md.export_cfg("");
Explanation: Define an md.export_cfg object
md.export_cfg has a call method that we can use to create quick snapshots of our simulation box
End of explanation
sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
nlyrs_fxd=2
a=sim.H[0][0];
b_norm=0.5*a*np.sqrt(3.0);
b=np.array([1.0,1.0,1.0])
s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
Explanation: Screw dislocation
End of explanation
sim.cell_change([[1,-1,0],[1,1,-2],[1,1,1]])
Explanation: Create a $\langle110\rangle\times\langle112\rangle\times\frac{1}{2}\langle111\rangle$ cell
create a $\langle110\rangle\times\langle112\rangle\times\langle111\rangle$ cell
Since mapp4py.md.atoms.cell_chenge() only accepts integer values start by creating a $\langle110\rangle\times\langle112\rangle\times\langle111\rangle$ cell
End of explanation
H=np.array(sim.H);
def _(x):
    if x[2] > 0.5*H[2, 2] - 1.0e-8:
        return False;
    else:
        x[2]*=2.0;
sim.do(_);
_ = np.full((3,3), 0.0)
_[2, 2] = - 0.5
sim.strain(_)
Explanation: Remove half of the atoms and readjust the position of remaining
Now one needs to cut the cell in half in $[111]$ direction. We can achive this in three steps:
Remove the atoms that are above located above $\frac{1}{2}[111]$
Double the position of the remiaing atoms in the said direction
Shrink the box affinly to half on that direction
End of explanation
displace(sim,np.array([sim.H[0][0]/6.0,sim.H[1][1]/6.0,0.0]))
Explanation: Readjust the postions
End of explanation
max_natms=100000
H=np.array(sim.H);
n_per_area=sim.natms/(H[0,0] * H[1,1]);
_ =np.sqrt(max_natms/n_per_area);
N0 = np.array([
    np.around(_ / sim.H[0][0]),
    np.around(_ / sim.H[1][1]), 
    1], dtype=np.int32)
sim *= N0;
H = np.array(sim.H);
H_new = np.array(sim.H);
H_new[1][1] += 50.0
resize(sim, H_new, np.full((3),0.5) @ H)
C_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
Q=np.array([np.cross(s,b)/np.linalg.norm(np.cross(s,b)),s/np.linalg.norm(s),b/np.linalg.norm(b)])
hirth = HirthScrew(rot(C_Fe,Q), rot(b*0.5*a,Q))
ctr = np.full((3),0.5) @ H_new;
s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
def _(x,x_d,x_dof):
    sy=(x[1]-ctr[1])/H[1, 1];
    x0=(x-ctr)/H[0, 0];
    if sy>s_fxd or sy<=-s_fxd:
        x_dof[1]=x_dof[2]=False;
        x+=b_norm*hirth.ave_disp(x0)
    else:
        x+=b_norm*hirth.disp(x0)
sim.do(_)    
H = np.array(sim.H);
H_inv = np.array(sim.B);
H_new = np.array(sim.H);
H_new[0,0]=np.sqrt(H[0,0]**2+(0.5*b_norm)**2)
H_new[2,0]=H[2,2]*0.5*b_norm/H_new[0,0]
H_new[2,2]=np.sqrt(H[2,2]**2-H_new[2,0]**2)
F = np.transpose(H_inv @ H_new);
sim.strain(F - np.identity(3))
xprt(sim, "dumps/screw.cfg")
Explanation: Replicating the unit cell
End of explanation
def make_scrw(nlyrs_fxd,nlyrs_vel,vel):
    #this is for 0K
    #c_Fe=cubic(1.5187249951755375,0.9053185628093443,0.7249256807942608);
    #this is for 300K
    c_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
    
    #N0=np.array([80,46,5],dtype=np.int32)
    sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
    a=sim.H[0][0];
    b_norm=0.5*a*np.sqrt(3.0);
    b=np.array([1.0,1.0,1.0])
    s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
    Q=np.array([np.cross(s,b)/np.linalg.norm(np.cross(s,b)),s/np.linalg.norm(s),b/np.linalg.norm(b)])
    c0=rot(c_Fe,Q)
    
    hirth = HirthScrew(rot(c_Fe,Q),np.dot(Q,b)*0.5*a)
    sim.cell_change([[1,-1,0],[1,1,-2],[1,1,1]])
    displace(sim,np.array([sim.H[0][0]/6.0,sim.H[1][1]/6.0,0.0]))
    max_natms=1000000
    n_per_vol=sim.natms/sim.vol;
    _=np.power(max_natms/n_per_vol,1.0/3.0);
    N1=np.full((3),0,dtype=np.int32);
    for i in range(0,3):
        N1[i]=int(np.around(_/sim.H[i][i]));
    N0=np.array([N1[0],N1[1],1],dtype=np.int32);
    sim*=N0;
    sim.kB=8.617330350e-5
    sim.create_temp(300.0,8569643);
    H=np.array(sim.H);
    H_new=np.array(sim.H);
    H_new[1][1]+=50.0
    resize(sim, H_new, np.full((3),0.5) @ H)
    ctr=np.dot(np.full((3),0.5),H_new);
    s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
    s_vel=0.5-0.5*float(nlyrs_vel)/float(N0[1])
    def _(x,x_d,x_dof):
        sy=(x[1]-ctr[1])/H[1][1];
        x0=(x-ctr)/H[0][0];
        
        if sy>s_fxd or sy<=-s_fxd:
            x_d[1]=0.0;
            x_dof[1]=x_dof[2]=False;
            x+=b_norm*hirth.ave_disp(x0)
        else:
            x+=b_norm*hirth.disp(x0)
        
        if sy<=-s_vel or sy>s_vel:
            x_d[2]=2.0*sy*vel;
    sim.do(_)    
    H = np.array(sim.H);
    H_inv = np.array(sim.B);
    H_new = np.array(sim.H);
    H_new[0,0]=np.sqrt(H[0,0]**2+(0.5*b_norm)**2)
    H_new[2,0]=H[2,2]*0.5*b_norm/H_new[0,0]
    H_new[2,2]=np.sqrt(H[2,2]**2-H_new[2,0]**2)
    F = np.transpose(H_inv @ H_new);
    sim.strain(F - np.identity(3))
    return N1[2],sim;
Explanation: putting it all together
End of explanation
sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
nlyrs_fxd=2
a=sim.H[0][0];
b_norm=0.5*a*np.sqrt(3.0);
b=np.array([1.0,1.0,1.0])
s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
sim.cell_change([[1,1,1],[1,-1,0],[1,1,-2]])
H=np.array(sim.H);
def _(x):
    if x[0] > 0.5*H[0, 0] - 1.0e-8:
        return False;
    else:
        x[0]*=2.0;
sim.do(_);
_ = np.full((3,3), 0.0)
_[0,0] = - 0.5
sim.strain(_)
displace(sim,np.array([0.0,sim.H[1][1]/4.0,0.0]))
max_natms=100000
H=np.array(sim.H);
n_per_area=sim.natms/(H[0, 0] * H[1, 1]);
_ =np.sqrt(max_natms/n_per_area);
N0 = np.array([
    np.around(_ / sim.H[0, 0]),
    np.around(_ / sim.H[1, 1]), 
    1], dtype=np.int32)
sim *=  N0;
# remove one layer along ... direction
H=np.array(sim.H);
frac=H[0,0] /N0[0]
def _(x):
    if x[0] < H[0, 0] /N0[0] and x[1] >0.5*H[1, 1]:
        return False;
sim.do(_)
H = np.array(sim.H);
H_new = np.array(sim.H);
H_new[1][1] += 50.0
resize(sim, H_new, np.full((3),0.5) @ H)
C_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
_ = np.cross(b,s)
Q = np.array([b/np.linalg.norm(b), s/np.linalg.norm(s), _/np.linalg.norm(_)])
hirth = HirthEdge(rot(C_Fe,Q), rot(b*0.5*a,Q))
_ = (1.0+0.5*(N0[0]-1.0))/N0[0];
ctr = np.array([_,0.5,0.5]) @ H_new;
frac = H[0][0]/N0[0]
s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
def _(x,x_d,x_dof):
    sy=(x[1]-ctr[1])/H[1, 1];
    x0=(x-ctr);
    if(x0[1]>0.0):
        x0/=(H[0, 0]-frac)
    else:
        x0/= H[0, 0]
    if sy>s_fxd or sy<=-s_fxd:
        x+=b_norm*hirth.ave_disp(x0);
        x_dof[0]=x_dof[1]=False;
    else:
        x+=b_norm*hirth.disp(x0);
    x[0]-=0.25*b_norm;
sim.do(_)
H = np.array(sim.H)
H_new = np.array(sim.H);
H_new[0, 0] -= 0.5*b_norm;
resize(sim, H_new, np.full((3),0.5) @ H)
xprt(sim, "dumps/edge.cfg")
Explanation: Edge dislocation
End of explanation
def make_edge(nlyrs_fxd,nlyrs_vel,vel):
    #this is for 0K
    #c_Fe=cubic(1.5187249951755375,0.9053185628093443,0.7249256807942608);
    #this is for 300K
    c_Fe=cubic(1.3967587463636366,0.787341583191591,0.609615090769241);
    
    #N0=np.array([80,46,5],dtype=np.int32)
    sim=md.atoms.import_cfg('configs/Fe_300K.cfg');
    a=sim.H[0][0];
    b_norm=0.5*a*np.sqrt(3.0);
    b=np.array([1.0,1.0,1.0])
    s=np.array([1.0,-1.0,0.0])/np.sqrt(2.0)
    # create rotation matrix
    _ = np.cross(b,s)
    Q=np.array([b/np.linalg.norm(b), s/np.linalg.norm(s), _/np.linalg.norm(_)])
    hirth = HirthEdge(rot(c_Fe,Q),np.dot(Q,b)*0.5*a)
    # create a unit cell 
    sim.cell_change([[1,1,1],[1,-1,0],[1,1,-2]])
    H=np.array(sim.H);
    def f0(x):
        if x[0]>0.5*H[0][0]-1.0e-8:
            return False;
        else:
            x[0]*=2.0;
    sim.do(f0);
    _ = np.full((3,3), 0.0)
    _[0,0] = - 0.5
    sim.strain(_)
    displace(sim,np.array([0.0,sim.H[1][1]/4.0,0.0]))
    max_natms=1000000
    n_per_vol=sim.natms/sim.vol;
    _=np.power(max_natms/n_per_vol,1.0/3.0);
    N1=np.full((3),0,dtype=np.int32);
    for i in range(0,3):
        N1[i]=int(np.around(_/sim.H[i][i]));
    N0=np.array([N1[0],N1[1],1],dtype=np.int32);
    N0[0]+=1;
    sim*=N0;
    # remove one layer along ... direction
    H=np.array(sim.H);
    frac=H[0][0]/N0[0]
    def _(x):
        if x[0] < H[0][0]/N0[0] and x[1]>0.5*H[1][1]:
            return False;
    sim.do(_)
    
    
    sim.kB=8.617330350e-5
    sim.create_temp(300.0,8569643);
    H = np.array(sim.H);
    H_new = np.array(sim.H);
    H_new[1][1] += 50.0
    ctr=np.dot(np.full((3),0.5),H);
    resize(sim,H_new, np.full((3),0.5) @ H)
    l=(1.0+0.5*(N0[0]-1.0))/N0[0];
    ctr=np.dot(np.array([l,0.5,0.5]),H_new);
    frac=H[0][0]/N0[0]
    s_fxd=0.5-0.5*float(nlyrs_fxd)/float(N0[1])
    s_vel=0.5-0.5*float(nlyrs_vel)/float(N0[1])
    def f(x,x_d,x_dof):
        sy=(x[1]-ctr[1])/H[1][1];
        x0=(x-ctr);
        if(x0[1]>0.0):
            x0/=(H[0][0]-frac)
        else:
            x0/= H[0][0]
        if sy>s_fxd or sy<=-s_fxd:
            x_d[1]=0.0;
            x_dof[0]=x_dof[1]=False;
            x+=b_norm*hirth.ave_disp(x0);
        else:
            x+=b_norm*hirth.disp(x0);
        
        if sy<=-s_vel or sy>s_vel:
            x_d[0]=2.0*sy*vel;
        x[0]-=0.25*b_norm;
    sim.do(f)
    H = np.array(sim.H)
    H_new = np.array(sim.H);
    H_new[0, 0] -= 0.5*b_norm;
    resize(sim, H_new, np.full((3),0.5) @ H)
    return N1[2], sim;
nlyrs_fxd=2
nlyrs_vel=7;
vel=-0.004;
N,sim=make_edge(nlyrs_fxd,nlyrs_vel,vel)
xprt(sim, "dumps/edge.cfg")
_ = np.array([[-1,1,0],[1,1,1],[1,1,-2]], dtype=np.float);
Q = np.linalg.inv(np.sqrt(_ @ _.T)) @ _;
C = rot(cubic(1.3967587463636366,0.787341583191591,0.609615090769241),Q)
B = np.linalg.inv(
    np.array([
    [C[0, 0, 0, 0], C[0, 0, 1, 1], C[0, 0, 0, 1]],
    [C[0, 0, 1, 1], C[1, 1, 1, 1], C[1, 1, 0, 1]],
    [C[0, 0, 0, 1], C[1, 1, 0, 1], C[0, 1, 0, 1]]
    ]
))
_ = np.roots([B[0, 0], -2.0*B[0, 2],2.0*B[0, 1]+B[2, 2], -2.0*B[1, 2], B[1, 1]])
mu = np.array([_[0],0.0]);
if np.absolute(np.conjugate(mu[0]) - _[1]) > 1.0e-12:
    mu[1] = _[1];
else:
    mu[1] = _[2]
alpha = np.real(mu);
beta = np.imag(mu);
p = B[0,0] * mu**2 - B[0,2] * mu + B[0, 1]
q = B[0,1] * mu - B[0, 2] + B[1, 1]/ mu
K = np.stack([p, q]) * np.array(mu[1], mu[0]) /(mu[1] - mu[0])
K_r = np.real(K)
K_i = np.imag(K)
Tr = np.stack([
    np.array(np.array([[1.0, alpha[0]], [0.0, beta[0]]])), 
    np.array([[1.0, alpha[1]], [0.0, beta[1]]])
], axis=1)
def u_f0(x): return np.sqrt(np.sqrt(x[0] * x[0] + x[1] * x[1]) + x[0])
def u_f1(x): return np.sqrt(np.sqrt(x[0] * x[0] + x[1] * x[1]) - x[0]) * np.sign(x[1]) 
def disp(x): 
    _ = Tr @ x
    return K_r @ u_f0(_) + K_i @ u_f1(_)
Explanation: putting it all together
End of explanation
_ = np.array([[-1,1,0],[1,1,1],[1,1,-2]], dtype=np.float);
Q = np.linalg.inv(np.sqrt(_ @ _.T)) @ _;
C = rot(cubic(1.3967587463636366,0.787341583191591,0.609615090769241),Q)
disp = crack(C)
n = 300;
r = 10;
disp_scale = 0.3;
n0 = int(np.round(n/ (1 +np.pi), ))
n1 = n - n0
xs = np.concatenate((
        np.stack([np.linspace(0, -r , n0), np.full((n0,), -1.e-8)]),
        r * np.stack([np.cos(np.linspace(-np.pi, np.pi , n1)),np.sin(np.linspace(-np.pi, np.pi , n1))]), 
        np.stack([np.linspace(-r, 0 , n0), np.full((n0,), 1.e-8)]),
    ), axis =1)
xs_def = xs + disp_scale * disp(xs)
fig, ax = plt.subplots(figsize=(10.5,5), ncols = 2)
ax[0].plot(xs[0], xs[1], "b-", label="non-deformed");
ax[1].plot(xs_def[0], xs_def[1], "r-.", label="deformed");
Explanation: Putting it all together
End of explanation
<END_TASK>