# Voltage control of skyrmions: creation, annihilation and zero-magnetic field
stablization

Yifan Zhou  NanoSpin, Department of Applied Physics, Aalto University School
of Science, P.O. Box 15100, FI-00076 Aalto, Finland  Rhodri Mansell
<EMAIL_ADDRESS>NanoSpin, Department of Applied Physics, Aalto
University School of Science, P.O. Box 15100, FI-00076 Aalto, Finland
Sebastiaan van Dijken  NanoSpin, Department of Applied Physics, Aalto
University School of Science, P.O. Box 15100, FI-00076 Aalto, Finland

###### Abstract

Voltage manipulation of skyrmions is a promising path towards low-energy
spintronic devices. Here, voltage effects on skyrmions in a GdOx/Gd/Co/Pt
heterostructure are observed experimentally. The results show that the
skyrmion density can be both enhanced and depleted by the application of an
electric field, along with the ability, at certain magnetic fields to
completely switch the skyrmion state on and off. Further, a zero magnetic
field skyrmion state can be stablized under a negative bias voltage using a
defined voltage and magnetic field sequence. The voltage effects measured here
occur on a few-second timescale, suggesting an origin in voltage-controlled
magnetic anisotropy rather than ionic effects. By investigating the skyrmion
nucleation rate as a function of temperature, we extract the energy barrier to
skyrmion nucleation in our sample. Further, micromagnetic simulations are used
to explore the effect of changing the anisotropy and Dzyaloshinskii-Moriya
interaction on skyrmion density. Our work demonstrates the control of
skyrmions by voltages, showing functionalities desirable for commercial
devices.

††preprint: AIP/123-QED

Magnetic skyrmions are topologically non-trivial spin textures, which are
widely observed in thin film magnetic trilayers consisting of a heavy metal
(HM), a ferromagnet (FM) and a metal oxide (MO), such as Pt/Co/MgOBoulle  _et
al._ (2016), Pt/CoFeB/MgOWoo  _et al._ (2016), Ta/CoFeB/MgOGilbert  _et al._
(2015), Ta/CoFeB/TaOxJiang  _et al._ (2015) and so onEverschor-Sitte  _et al._
(2018). A combination of perpendicular magnetic anisotropy (PMA) and the
Dzyaloshinskii–Moriya interaction (DMI) can lead to a skyrmion stateRoessler,
Bogdanov, and Pfleiderer (2006); Nagaosa and Tokura (2013); Fert, Cros, and
Sampaio (2013) in such systems. The strong spin-orbit coupling at the HM/FM
interface gives rise to PMAHellman  _et al._ (2017) and, in combination with
the broken inversion symmetry in the growth direction, DMIYang  _et al._
(2015). Further, PMA and DMI originate from the FM/MO interface due to the
overlapping $p$ orbitals from oxygen and $d$ orbitals from the ferromagnetYang
_et al._ (2011). Skyrmions have attractive features for device applications,
such as small sizesWang, Yuan, and Wang (2018), stability at room
temperatureBüttner, Lemesh, and Beach (2018), and can be driven into motion by
a relatively low current densityIwasaki, Mochizuki, and Nagaosa (2013).
However, the skyrmion Hall effectJiang  _et al._ (2017), as well as the Joule
heating produced by driving currentsKoshibae  _et al._ (2015); Koshibae and
Nagaosa (2014), has hindered the application of skyrmions to current-driven
memory devices.

As an alternative approach to current-driven devices, the voltage control of
magnetism has been widely investigated, initially in magnetic tunnel
junctionsWang  _et al._ (2012); Kanai  _et al._ (2012); Shiota  _et al._
(2012) and then other FM/MO systemsBi  _et al._ (2014); Bauer  _et al._
(2015); Baldrati  _et al._ (2017). Several proposals have been made for
skyrmion-based voltage controlled memory devices, using both
staticBhattacharya, Al-Rashid, and Atulasimha (2016); Kasai  _et al._ (2019)
and mobileWang  _et al._ (2018); Zhou, Mansell, and van Dijken (2019)
skyrmions. In HM/FM/MO structures hosting skyrmions, recent experiments have
demonstrated the voltage control of magnetic anisotropy (VCMA), with the
additional ability to control DMI by applying voltagesSrivastava  _et al._
(2018); Yang  _et al._ (2020). In such experiments, due to the modification of
these underlying magnetic properties, skyrmions can be created and annihilated
by applied voltages. Different mechanisms have been proposed which would allow
the voltage control of skyrmions, namely changing the electron orbital filling
with an electric fieldHsu  _et al._ (2017); Schott  _et al._ (2017);
Bhattacharya  _et al._ (2020), modifying the Rashba DMI field at FM/MO
interfaceSrivastava  _et al._ (2018), and introducing strains from flexible
Yang  _et al._ (2020) or ferroelectric Li  _et al._ (2018); Wang  _et al._
(2020) substrates. Inspired by the control of magnetism through ionic effects
demonstrated in Pt/Co/GdOx heterostructuresTan  _et al._ (2019); Bauer  _et
al._ (2015), we explore the possibility of observing skyrmions in such a
structure, and subsequently controlling them by applied voltages.

Figure 1: (a) Schematic of the multilayer sample with a cross bar structure.
(b) Hysteresis loops obtained by MOKE microscopy with constant bias voltages
of 0 V, 2 V and -2 V. (c) Skyrmion states at 3 mT with different bias voltages
of -1 V, 0 V and 1 V.

The thin film sample studied in this work is a Ta(4) / GdOx(6) / Gd(0.1) /
Co(1) / Pt(4) (in nm) heterostructure. The sample is deposited by magnetron
sputtering at room temperature in a system with a base pressure of $\sim
5\times 10^{-8}$ mbar. Metal layers are grown by DC sputtering with an Ar
pressure of $8\times 10^{-3}$ mbar, while the GdOx layer is grown by reactive
DC sputtering with 10% O2 partial pressure. The sample is grown with an
‘inverted’ layer structure, with the magnetic metal layer grown on top of the
GdOx. The introduction of the thin Gd metal layer acts to reduce the oxidation
of the Co layer. As shown in Fig. 1(a), the multilayer is patterned into a
crossbar structure by direct laser-writing lithography. In the patterned
junction, Ta is the bottom electrode (BE), GdOx is an insulating layer and the
Gd/Co/Pt multilayer is the top electrode (TE). The junction area is 50 $\mu$m
$\times$ 50 $\mu$m. The sign of the applied voltage is defined from the top
electrode to the bottom electrode, where a positive sign means the voltage on
the top electrode is higher than on the bottom electrode.

Magneto-optical Kerr effect (MOKE) microscopy is used to record and image the
out-of-plane magnetization of the sample, under externally applied out-of-
plane magnetic fields and voltages. Magnetic properties of the sample with
zero voltage, such as the saturation magnetisation $M_{s}$ and out-of-plane
anisotropy $K_{u}$, were measured on an equivalent thin film sample by
vibrating sample magnetometry (VSM) at $25^{\circ}$C, where $M_{s}=1.4\times
10^{6}$ A/m and $K_{u}=6\times 10^{6}$ J/m3, consistent with previously
reported values from a similar structurePham  _et al._ (2016).

To study voltage effects on the magnetic properties of the sample, we first
measure the out-of-plane hysteresis loop with a constantly applied voltage of
0 V, 2 V and -2 V using MOKE microscopy (Fig. 1(b)). The results show a
perpendicularly magnetized sample with near to full remanence at zero applied
voltage. A positive bias voltage decreases the coercivity of the sample,
indicating a reduction of the perpendicular anisotropy, or possibly an
increase in DMI, while a negative bias has the reverse effect. The voltage
effect is volatile, meaning that after any applied voltage is removed, the 0 V
hysteresis loop is the same as before the application of voltage. Images of
skyrmion states are captured at different bias voltages by first saturating
the sample at 10 mT and then decreasing the field to $3$ mT (Fig. 1(c)). The
images are taken after a relaxation time of one minute to allow for skyrmion
nucleation. The voltage is applied throughout this process. In spite of the
relatively small voltage effect seen in the hysteresis loops, the skyrmion
density varies significantly with bias voltage, where more skyrmions are
observed with positive voltage and fewer with negative voltage. Due to the
resolution limit ($\sim 500$ nm) of white-light MOKE microscopy, we are not
able to observe variations of the skyrmion radius, which might be expected
from voltage-induced changes of the magnetic anisotropy.

Figure 2: Real time control of skyrmion creation and annihilation from a
uniform magnetization state at 3.5 mT with a voltage sequence of i. 0 V, ii. 2
V, iii. -2 V, iv. 2 V, v. -2 V, vi. 0 V. Each image is taken 30 s after
changing the applied voltage.

On-off control of skyrmions starting from a uniform magnetization state can be
achieved by applying a suitable voltage sequence (Fig. 2). The sample is
initially saturated by a 3.5 mT out-of-plane field, which is slightly larger
than the transition field from an uniform state to the skyrmion state at $3$
mT (Fig. 2i.). After this, a 2 V bias is applied continuously while the
magnetic field is fixed at 3.5 mT, and the sample is imaged by MOKE microscopy
after 30 s (Fig. 2ii.). Under these conditions, skyrmions are created by the
positive bias voltage. The voltage is then set to -2 V (Fig. 2iii.), and after
30 s most of the skyrmions have disappeared. Repetition of the same voltage
sequence (Fig. 2iv. and 2v.) produces a similar effect. The system returns to
a uniform magnetization state at 0 V (Fig. 2vi.).

Figure 3: (a) Schematic of applied magnetic field and voltage sequence. Times
when the images in (b), (c) and (d) were captured are marked. (b) Skyrmion
state at 2.8 mT and 0 V. (c) Skyrmion state at 0 mT and -2 V. (d) Multidomain
state at 0 mT and 0 V. The inset shows the multidomain state at 0 mT that is
attained directly from saturation at 10 mT with zero voltage.

Besides the on-demand creation and annihilation of skyrmions at 3.5 mT, we
find that a negative bias voltage can stabilize skyrmions at zero magnetic
field (Fig. 3). To demonstrate this, we first create a skyrmion state at 2.8
mT without a bias voltage (Fig. 3(b)). Then we apply -2 V and the magnetic
field is turned off immediately afterwards. After 30 s, the skyrmion state is
still similar to the initial skyrmion state at 2.8 mT and 0 V (compare Fig.
3(c) and Fig. 3(b)). In order to confirm that it is the negative bias voltage
that controls the stabilization of skyrmions in zero magnetic field, we turn
off the voltage subsequently. A clear transition occurs as the skyrmions then
expand and the sample shows a multidomain state, similar to that seen at zero
magnetic field and voltage after saturation (Fig. 3(d)). Here, by increasing
the PMA with a negative bias voltage, the skyrmions, once formed, are
stabilized against expanding to form worm-like domains.

Figure 4: (a) Simulated changes in the skyrmion state when parameters Ku and D
are changed by +10% and -10% compared to their initial value. The orange and
blue lines illustrate the timeline of the micromagnetic simulations. (b) The
fractional change in skyrmion numbers due to a fractional change in $K_{u}$ or
$D$.

Having shown voltage control of skyrmions in our system we turn to the
question of underlying parameters being controlled by the applied voltage. It
has been shown that voltages could influence $M_{s}$, $K_{u}$ and DMI $D$. We
find that the saturated MOKE signals under different voltages are almost
identical, indicating that $M_{s}$ remains the same. Due to a limited in-plane
field in our MOKE system, we are not able to measure changes in $K_{u}$
directly. Instead, we perform micromagnetic simulations with the MuMax3
package to gain insight into voltage effects on the skyrmion density. In order
to achieve a spontaneous skyrmion state in simulations with a reasonable time
scale, we adopt the following initial parameters: $M_{s}=0.8\times 10^{6}$
A/m, exchange constant $A=0.7\times 10^{-12}$ J/m, $K_{u}=0.5\times 10^{6}$
J/m3, $D=1.5\times 10^{-3}$ J/m2 and damping constant $\alpha=0.3$. A constant
perpendicular magnetic field of 160 mT is applied to nucleate skyrmions, and
the simulation temperature is set to 300 K. Initially, the system is allowed
to relax for 20 ns, and a snapshot of the magnetization is recorded. Then,
either $K_{u}$ or $D$ are modified by a certain percentage of their initial
value, and the resulting skyrmion state is recorded after a further 20 ns
(Fig. 4(a)). The effect of changing $K_{u}$ or $D$ is directly illustrated by
a change in the number of skyrmions, $\Delta N$, compared to the initial value
($\Delta N/N_{0}$). In Fig. 4(b), both an increase in $K_{u}$ and a decrease
in $D$ lead to a decrease in $N$, and vice versa. Below 10 % variation the
effects of changing $K_{u}$ and $D$ are fairly symmetric, with changing $D$
being more effective for larger variation. By comparing simulation results to
Fig. 1(c) and Fig. 2, we infer a negative bias voltage in our system either
increases $K_{u}$, decreases $D$, or possibly both. A more quantitative
comparison of simulations and our experimental results is not possible, due to
the very different timescales involved.

Figure 5: (a) Average domain width obtained after demagnetization under
applied 2 s voltage pulses of -2 V, 0 V and 2 V. The scale bar marks 5 $\mu$m.
(b) Hysteresis loops of the sample at 25∘C (solid lines) and 45∘C (dashed
lines) with -2 V, 0 V and 2 V. (c) Time dependence of the variation of
skyrmion numbers under 1.5 mT when applying a constant 2 V at 25∘C (blue) and
45∘C (green). (d) Fitting of the Arrhenius relation of the skyrmion creation
time under 1.5 mT and 2 V at different temperatures.

Another open question is the underlying physical mechanism in our system.
There are two main possible mechanisms of voltage control – orbital filling
effects or voltage-induced ion migration. In order to investigate this we
first note that in our experiments, there are two relevant timescales: the
relaxation time of the magnetic microstructure, which depends on the relevant
thermal activation barrier; and the response time of the voltage-induced
changes of magnetic parameters such as PMA and the DMI. For voltage effects
driven by electron orbital filling we would expect the changes in the
parameters to occur instantaneously on the timescales of our experiments. For
ion migration driven changes the time scale is less clear, but could be
expected to occur over a period of hours at room temperature. Generally, even
if the voltage controlled changes of magnetic parameters are fast, it takes
more time for skyrmions to nucleate or annihilate because these processes are
thermally activated. Distinguishing between the thermal-induced relaxation and
the effect of voltages would allow us to elucidate the mechanism of voltage
control. In order to study the timescale of the voltage effect, we first
investigate the short time behavior of the sample under applied voltages. We
demagnetized the sample by an AC oscillating field in 2 s in order to exclude
magnetic relaxation effects in the Co layer. A voltage pulse is applied
simultaneously with the demagnetization field and turned off immediately after
the sample is demagnetized. Hereafter, an image of the zero-field domain state
is taken immediately. From the domain states found for different voltages, the
average domain width of the sample is extracted by a Gaussian fitting to the
fast Fourier transformation of the domain images (Fig. 5 (a)). Clearly, 2 s
voltage pulses affect the domain width in the demagnetized state, where a -2 V
pulse increases the average domain width and 2 V has the reverse effect. This
is consistent with the changes in the skyrmion density found with longer
pulses, showing that significant effects are seen within 2 s of applying the
voltage.

To study how voltage-induced changes of the magnetic parameters affect the
nucleation of skyrmions on longer time scales, we investigate the thermal
activation of skyrmions by conducting experiments at different temperature. In
Fig. 5(b), voltage effects on the magnetic hysteresis loop are shown at 25∘C
and 45∘C in a junction exhibiting a skyrmion state from $1$ mT to $2$ mT
depending on the temperature. To investigate skyrmion nucleation we first
relax the sample at $1.5$ mT, then apply a constant 2 V bias and count the
number of skyrmions as a function of time for 1 min. Since a positive voltage
could decreases the PMA or increase the DMI, more skyrmions are created
following the application of the voltage through thermal activation (Fig.
5(c)). By assuming that the voltage effect is effectively instantaneous and
does not change with time, the following time-dependent equation of skyrmion
number $N$ can be written with a scale factor $A$, a constant skyrmion
creation time $\tau$ and initial number of skyrmions $N_{0}$Wilson  _et al._
(2019):

$N-N_{0}=A(1-\exp[-t/\tau]).$ (1)

The data of $N$ as a function of $t$ is collected at different temperatures:
25∘C, 30∘C, 35∘C, 40∘C and 45∘C. The value of $\tau$ is fitted at each
temperature. If the nucleation process is determined by a single energy
barrier, then we would expect the rate to follow an Arrhenius law:

$\tau=\tau_{0}\exp[-\frac{E}{k_{B}T}],$ (2)

where $\tau_{0}$ is the scale factor, $E$ is the energy barrier, $k_{B}$ is
the Boltzmann constant and $T$ is the temperature. As shown in Fig. 5(d), we
fit $\tau$ to Eq. (2), which gives a nucleation energy barrier $E=1.4$ eV for
skyrmions at 2 V and 1.5 mT. Our assumption of a constant creation rate and
single energy barrier is not explicitly violated by this data. Combined with
our other experiments, this means that the voltage effects on the parameters
determining the energy barrier are likely to be near instantaneous without
longer-time effects. The direction of voltages effects, i.e. decreasing PMA
with positive voltage and vice versa, as well as their volatility, is
consistent with that found by ab-initio calculation of electronic effects at
Fe/MgO interfaces, assuming that the first interfacial Fe layer is
oxidizedZhang  _et al._ (2017). The voltage effect on skyrmions seen here also
has the same sign as in the Co/AlOx systemSchott  _et al._ (2017). Previously
it has been reported that changes in a Pt/Co/GdOx system were driven by ion
migration, where the mobile ions in the GdOx layer, were determined to largely
come from the atmosphereTan  _et al._ (2019). The apparent lack of ion
migration in our system may originate from the reversed layer sequence, where,
in our case, the GdOx layer is buried under the top electrode, blocking ion
diffusion originating from the atmosphere.

In conclusion, we demonstrated on-demand creation and annihilation of
skyrmions by an applied electric field in a GdOx/Gd/Co/Pt structure.
Additionally, we developed a method to stabilize skyrmions in zero magnetic
field by voltage. We have investigated through simulations the possible
underlying magnetic parameters influenced by the voltages. The simulation
results show that changes in PMA and DMI have similar effects on the skyrmion
density, but with opposite signs. We also looked at the timescale of the
effects, showing that the voltages have an effect within 2 s, without changing
significantly on a longer timescale. We conclude that the voltage effects
derive from the modification of the orbital filling at the Co/GdOx interface
meaning that they are in principle limited by capacitive effects. Our results
show voltage controlled skyrmion effects that could be exploited for device
physics and should encourage further work in this field.

## Acknowledgements

This work was supported by the Academy of Finland (Grant Nos. 295269, 306978
and 327804). We acknowledge the provision of facilities by Aalto University at
OtaNano, the Micronova Nanofabrication Center and the Nanomicroscopy Center,
as well as computational resources provided by the Aalto Science-IT project.

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Boulle  _et al._ (2016) O. Boulle, J. Vogel, H. Yang, S. Pizzini, D. de Souza Chaves, A. Locatelli, T. O. Menteş, A. Sala, L. D. Buda-Prejbeanu, O. Klein, _et al._ , Nature Nanotechnology 11, 449 (2016).
  * Woo  _et al._ (2016) S. Woo, K. Litzius, B. Krüger, M.-Y. Im, L. Caretta, K. Richter, M. Mann, A. Krone, R. M. Reeve, M. Weigand, _et al._ , Nature Materials 15, 501 (2016).
  * Gilbert  _et al._ (2015) D. A. Gilbert, B. B. Maranville, A. L. Balk, B. J. Kirby, P. Fischer, D. T. Pierce, J. Unguris, J. A. Borchers, and K. Liu, Nature Communications 6, 1 (2015).
  * Jiang  _et al._ (2015) W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B. Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak, K. L. Wang, O. Heinonen, _et al._ , Science 349, 283 (2015).
  * Everschor-Sitte  _et al._ (2018) K. Everschor-Sitte, J. Masell, R. M. Reeve, and M. Kläui, Journal of Applied Physics 124, 240901 (2018).
  * Roessler, Bogdanov, and Pfleiderer (2006) U. K. Roessler, A. Bogdanov, and C. Pfleiderer, Nature 442, 797 (2006).
  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, Nature Nanotechnology 8, 899 (2013).
  * Fert, Cros, and Sampaio (2013) A. Fert, V. Cros, and J. Sampaio, Nature Nanotechnology 8, 152 (2013).
  * Hellman  _et al._ (2017) F. Hellman, A. Hoffmann, Y. Tserkovnyak, G. S. D. Beach, E. E. Fullerton, C. Leighton, A. H. MacDonald, D. C. Ralph, D. A. Arena, H. A. Dürr, _et al._ , Reviews of Modern Physics 89, 025006 (2017).
  * Yang  _et al._ (2015) H. Yang, A. Thiaville, S. Rohart, A. Fert, and M. Chshiev, Physical Review Letters 115, 267210 (2015).
  * Yang  _et al._ (2011) H. Yang, M. Chshiev, B. Dieny, J. Lee, A. Manchon, and K. Shin, Physical Review B 84, 054401 (2011).
  * Wang, Yuan, and Wang (2018) X. Wang, H. Yuan, and X. Wang, Communications Physics 1, 1 (2018).
  * Büttner, Lemesh, and Beach (2018) F. Büttner, I. Lemesh, and G. S. D. Beach, Scientific Reports 8, 1 (2018).
  * Iwasaki, Mochizuki, and Nagaosa (2013) J. Iwasaki, M. Mochizuki, and N. Nagaosa, Nature Nanotechnology 8, 742 (2013).
  * Jiang  _et al._ (2017) W. Jiang, X. Zhang, G. Yu, W. Zhang, X. Wang, M. B. Jungfleisch, J. E. Pearson, X. Cheng, O. Heinonen, K. L. Wang, _et al._ , Nature Physics 13, 162 (2017).
  * Koshibae  _et al._ (2015) W. Koshibae, Y. Kaneko, J. Iwasaki, M. Kawasaki, Y. Tokura, and N. Nagaosa, Japanese Journal of Applied Physics 54, 053001 (2015).
  * Koshibae and Nagaosa (2014) W. Koshibae and N. Nagaosa, Nature Communications 5, 1 (2014).
  * Wang  _et al._ (2012) W. Wang, M. Li, S. Hageman, and C. Chien, Nature Materials 11, 64 (2012).
  * Kanai  _et al._ (2012) S. Kanai, M. Yamanouchi, S. Ikeda, Y. Nakatani, F. Matsukura, and H. Ohno, Applied Physics Letters 101, 122403 (2012).
  * Shiota  _et al._ (2012) Y. Shiota, T. Nozaki, F. Bonell, S. Murakami, T. Shinjo, and Y. Suzuki, Nature Materials 11, 39 (2012).
  * Bi  _et al._ (2014) C. Bi, Y. Liu, T. Newhouse Illige, M. Xu, M. Rosales, J. Freeland, O. Mryasov, S. Zhang, S. Te Velthuis, and W. Wang, Physical Review Letters 113, 267202 (2014).
  * Bauer  _et al._ (2015) U. Bauer, L. Yao, A. J. Tan, P. Agrawal, S. Emori, H. L. Tuller, S. van Dijken, and G. S. D. Beach, Nature Materials 14, 174 (2015).
  * Baldrati  _et al._ (2017) L. Baldrati, A. J. Tan, M. Mann, R. Bertacco, and G. S. D. Beach, Applied Physics Letters 110, 012404 (2017).
  * Bhattacharya, Al-Rashid, and Atulasimha (2016) D. Bhattacharya, M. M. Al-Rashid, and J. Atulasimha, Scientific Reports 6, 31272 (2016).
  * Kasai  _et al._ (2019) S. Kasai, S. Sugimoto, Y. Nakatani, R. Ishikawa, and Y. K. Takahashi, Applied Physics Express 12, 083001 (2019).
  * Wang  _et al._ (2018) X. Wang, W. L. Gan, J. C. Martinez, F. N. Tan, M. B. A. Jalil, and W. S. Lew, Nanoscale 10, 733 (2018).
  * Zhou, Mansell, and van Dijken (2019) Y. Zhou, R. Mansell, and S. van Dijken, Scientific Reports 9, 6525 (2019).
  * Srivastava  _et al._ (2018) T. Srivastava, M. Schott, R. Juge, V. Krizakova, M. Belmeguenai, Y. Roussigné, A. Bernand-Mantel, L. Ranno, S. Pizzini, S.-M. Chérif, _et al._ , Nano Letters 18, 4871 (2018).
  * Yang  _et al._ (2020) Q. Yang, Y. Cheng, Y. Li, Z. Zhou, J. Liang, X. Zhao, Z. Hu, R. Peng, H. Yang, and M. Liu, Advanced Electronic Materials 6, 2000246 (2020).
  * Hsu  _et al._ (2017) P. J. Hsu, A. Kubetzka, A. Finco, N. Romming, K. von Bergmann, and R. Wiesendanger, Nature Nanotechnology 12, 123 (2017).
  * Schott  _et al._ (2017) M. Schott, A. Bernand-Mantel, L. Ranno, S. Pizzini, J. Vogel, H. Béa, C. Baraduc, S. Auffret, G. Gaudin, and D. Givord, Nano Letters 17, 3006 (2017).
  * Bhattacharya  _et al._ (2020) D. Bhattacharya, S. A. Razavi, H. Wu, B. Dai, K. L. Wang, and J. Atulasimha, Nature Electronics 3, 539 (2020).
  * Li  _et al._ (2018) Z. Li, Y. Zhang, Y. Huang, C. Wang, X. Zhang, Y. Liu, Y. Zhou, W. Kang, S. C. Koli, and N. Lei, Journal of Magnetism and Magnetic Materials 455, 19 (2018).
  * Wang  _et al._ (2020) Y. Wang, L. Wang, J. Xia, Z. Lai, G. Tian, X. Zhang, Z. Hou, X. Gao, W. M. nad Chun Feng, M. Zeng, G. Zhou, G. Yu, G. Wu, Y. Zhou, W. Wang, X. Zhang, and J. Liu, Nature Communications 11, 3577 (2020).
  * Tan  _et al._ (2019) A. J. Tan, M. Huang, C. O. Avci, F. Büttner, M. Mann, W. Hu, C. Mazzoli, S. Wilkins, H. L. Tuller, and G. S. D. Beach, Nature Materials 18, 35 (2019).
  * Pham  _et al._ (2016) T. H. Pham, J. Vogel, J. Sampaio, M. Vaňatka, J. C. Rojas-Sánchez, M. Bonfim, D. Chaves, F. Choueikani, P. Ohresser, E. Otero, _et al._ , EPL (Europhysics Letters) 113, 67001 (2016).
  * Wilson  _et al._ (2019) M. Wilson, M. Crisanti, C. Barker, A. Štefančič, J. White, M. Birch, G. Balakrishnan, R. Cubitt, and P. D. Hatton, Physical Review B 99, 174421 (2019).
  * Zhang  _et al._ (2017) J. Zhang, P. V. Lukashev, S. S. Jaswal, and E. Y. Tsymbal, Physical Review B 96, 014435 (2017).


865//– two distinct intersection points (x1, y1) and (x2, y2) find overlap
area

866

867double twointpts (double x[], double y[], double A1, double B1, double
PHI_1,

868

869 double A2, double B2, double H2_TR, double K2_TR,

870

871 double PHI_2, double AA, double BB, double CC, double DD,

872

873 double EE, double FF, int *rtnCode)

874

875{

876

877 double area1, area2;

878

879 double xmid, ymid, xmid_rt, ymid_rt;

880

881 double theta1, theta2;

882

883 double tmp, trsign;

884

885 double x1_tr, y1_tr, x2_tr, y2_tr;

886

887 double discr;

888

889 double cosphi, sinphi;

890

891

892

893 //– if execution arrives here, the intersection points are not

894

895 //– tangents.

896

897

898

899 //– determine which direction to integrate in the ellipse_segment

900

901 //– routine for each ellipse.

902

903

904

905 //– find the parametric angles for each point on ellipse 1

906

907 if (fabs (x[1]) $>$ A1)

908

909 x[1] = (x[1] $<$ 0) ? -A1 : A1;

910

911 if (y[1] $<$ 0.0) //– Quadrant III or IV

912

913 theta1 = twopi - acos (x[1] / A1);

914

915 else //– Quadrant I or II

916

917 theta1 = acos (x[1] / A1);

918

919

920

921 if (fabs (x[2]) $>$ A1)

922

923 x[2] = (x[2] $<$ 0) ? -A1 : A1;

924

925 if (y[2] $<$ 0.0) //– Quadrant III or IV

926

927 theta2 = twopi - acos (x[2] / A1);

928

929 else //– Quadrant I or II

930

931 theta2 = acos (x[2] / A1);

932

933

934

935 //– logic is for proceeding counterclockwise from theta1 to theta2

936

937 if (theta1 $>$ theta2)

938

939 {

940

941 tmp = theta1;

942

943 theta1 = theta2;

944

945 theta2 = tmp;

946

947 }

948

949

950

951 //– find a point on the first ellipse that is different than the two

952

953 //– intersection points.

954

955 xmid = A1*cos ((theta1 + theta2)/2.0);

956

957 ymid = B1*sin ((theta1 + theta2)/2.0);

958

959

960

961 //– the point (xmid, ymid) is on the first ellipse ’between’ the two

962

963 //– intersection points (x[1], y[1]) and (x[2], y[2]) when travelling

964

965 //– counter- clockwise from (x[1], y[1]) to (x[2], y[2]). If the point

966

967 //– (xmid, ymid) is inside the second ellipse, then the desired segment

968

969 //– of ellipse 1 contains the point (xmid, ymid), so integrate

970

971 //– counterclockwise from (x[1], y[1]) to (x[2], y[2]). Otherwise,

972

973 //– integrate counterclockwise from (x[2], y[2]) to (x[1], y[1])

974

975 if (ellipse2tr (xmid, ymid, AA, BB, CC, DD, EE, FF) $>$ 0.0)

976

977 {

978

979 tmp = theta1;

980

981 theta1 = theta2;

982

983 theta2 = tmp;

984

985 }

986

987

988

989 //– here is the ellipse segment routine for the first ellipse

990

991 if (theta1 $>$ theta2)

992

993 theta1 -= twopi;

994

995 if ((theta2 - theta1) $>$ pi)

996

997 trsign = 1.0;

998

999 else

1000

1001 trsign = -1.0;

1002

1003 area1 = 0.5*(A1*B1*(theta2 - theta1)

1004

1005 + trsign*fabs (x[1]*y[2] - x[2]*y[1]));

1006

1007

1008

1009 //– find ellipse 2 segment area. The ellipse segment routine

1010

1011 //– needs an ellipse that is centered at the origin and oriented

1012

1013 //– with the coordinate axes. The intersection points (x[1], y[1]) and

1014

1015 //– (x[2], y[2]) are found with both ellipses translated and rotated by

1016

1017 //– (-H1, -K1) and -PHI_1. Further translate and rotate the points

1018

1019 //– to put the second ellipse at the origin and oriented with the

1020

1021 //– coordinate axes. The translation is (-H2_TR, -K2_TR), and the

1022

1023 //– rotation is -(PHI_2 - PHI_1) = PHI_1 - PHI_2

1024

1025 cosphi = cos (PHI_1 - PHI_2);

1026

1027 sinphi = sin (PHI_1 - PHI_2);

1028

1029 x1_tr = (x[1] - H2_TR)*cosphi + (y[1] - K2_TR)*-sinphi;

1030

1031 y1_tr = (x[1] - H2_TR)*sinphi + (y[1] - K2_TR)*cosphi;

1032

1033 x2_tr = (x[2] - H2_TR)*cosphi + (y[2] - K2_TR)*-sinphi;

1034

1035 y2_tr = (x[2] - H2_TR)*sinphi + (y[2] - K2_TR)*cosphi;

1036

1037

1038

1039 //– determine which branch of the ellipse to integrate by finding a

1040

1041 //– point on the second ellipse, and asking whether it is inside the

1042

1043 //– first ellipse (in their once-translated+rotated positions)

1044

1045 //– find the parametric angles for each point on ellipse 1

1046

1047 if (fabs (x1_tr) $>$ A2)

1048

1049 x1_tr = (x1_tr $<$ 0) ? -A2 : A2;

1050

1051 if (y1_tr $<$ 0.0) //– Quadrant III or IV

1052

1053 theta1 = twopi - acos (x1_tr/A2);

1054

1055 else //– Quadrant I or II

1056

1057 theta1 = acos (x1_tr/A2);

1058

1059

1060

1061 if (fabs (x2_tr) $>$ A2)

1062

1063 x2_tr = (x2_tr $<$ 0) ? -A2 : A2;

1064

1065 if (y2_tr $<$ 0.0) //– Quadrant III or IV

1066

1067 theta2 = twopi - acos (x2_tr/A2);

1068

1069 else //– Quadrant I or II

1070

1071 theta2 = acos (x2_tr/A2);

1072

1073

1074

1075 //– logic is for proceeding counterclockwise from theta1 to theta2

1076

1077 if (theta1 $>$ theta2)

1078

1079 {

1080

1081 tmp = theta1;

1082

1083 theta1 = theta2;

1084

1085 theta2 = tmp;

1086

1087 }

1088

1089

1090

1091 //– find a point on the second ellipse that is different than the two

1092

1093 //– intersection points.

1094

1095 xmid = A2*cos ((theta1 + theta2)/2.0);

1096

1097 ymid = B2*sin ((theta1 + theta2)/2.0);

1098

1099

1100

1101 //– translate the point back to the second ellipse in its once-

1102

1103 //– translated+rotated position

1104

1105 cosphi = cos (PHI_2 - PHI_1);

1106

1107 sinphi = sin (PHI_2 - PHI_1);

1108

1109 xmid_rt = xmid*cosphi + ymid*-sinphi + H2_TR;

1110

1111 ymid_rt = xmid*sinphi + ymid*cosphi + K2_TR;

1112

1113

1114

1115 //– the point (xmid_rt, ymid_rt) is on the second ellipse ’between’ the

1116

1117 //– intersection points (x[1], y[1]) and (x[2], y[2]) when travelling

1118

1119 //– counterclockwise from (x[1], y[1]) to (x[2], y[2]). If the point

1120

1121 //– (xmid_rt, ymid_rt) is inside the first ellipse, then the desired

1122

1123 //– segment of ellipse 2 contains the point (xmid_rt, ymid_rt), so

1124

1125 //– integrate counterclockwise from (x[1], y[1]) to (x[2], y[2]).

1126

1127 //– Otherwise, integrate counterclockwise from (x[2], y[2]) to

1128

1129 //– (x[1], y[1])

1130

1131 if (((xmid_rt*xmid_rt)/(A1*A1) + (ymid_rt*ymid_rt)/(B1*B1)) $>$ 1.0)

1132

1133 {

1134

1135 tmp = theta1;

1136

1137 theta1 = theta2;

1138

1139 theta2 = tmp;

1140

1141 }

1142

1143

1144

1145 //– here is the ellipse segment routine for the second ellipse

1146

1147 if (theta1 $>$ theta2)

1148

1149 theta1 -= twopi;

1150

1151 if ((theta2 - theta1) $>$ pi)

1152

1153 trsign = 1.0;

1154

1155 else

1156

1157 trsign = -1.0;

1158

1159 area2 = 0.5*(A2*B2*(theta2 - theta1)

1160

1161 + trsign*fabs (x1_tr*y2_tr - x2_tr*y1_tr));

1162

1163

1164

1165 (*rtnCode) = TWO_INTERSECTION_POINTS;

1166

1167 return area1 + area2;

1168

1169}

1170

1171

1172

1173//– three distinct intersection points, must have two intersections

1174

1175//– and one tangent, which is the only possibility

1176

1177double threeintpts (double xint[], double yint[], double A1, double B1,

1178

1179 double PHI_1, double A2, double B2, double H2_TR,

1180

1181 double K2_TR, double PHI_2, double AA, double BB,

1182

1183 double CC, double DD, double EE, double FF,

1184

1185 int *rtnCode)

1186

1187{

1188

1189 int i, tanpts, tanindex, fnRtn;

1190

1191 double OverlapArea;

1192

1193

1194

1195 //– need to determine which point is a tangent, and which two points

1196

1197 //– are intersections

1198

1199 tanpts = 0;

1200

1201 for (i = 1; i $<$= 3; i++)

1202

1203 {

1204

1205 fnRtn = istanpt (xint[i], yint[i], A1, B1, AA, BB, CC, DD, EE, FF);

1206

1207

1208

1209 if (fnRtn == TANGENT_POINT)

1210

1211 {

1212

1213 tanpts++;

1214

1215 tanindex = i;

1216

1217 }

1218

1219 }

1220

1221

1222

1223 //– there MUST be 2 intersection points and only one tangent

1224

1225 if (tanpts != 1)

1226

1227 {

1228

1229 //– should never get here unless there is a problem discerning

1230

1231 //– whether or not a point is a tangent or intersection

1232

1233 (*rtnCode) = ERROR_INTERSECTION_PTS;

1234

1235 return -1.0;

1236

1237 }

1238

1239

1240

1241 //– store the two interesection points into (x[1], y[1]) and

1242

1243 //– (x[2], y[2])

1244

1245 switch (tanindex)

1246

1247 {

1248

1249 case 1:

1250

1251 xint[1] = xint[3];

1252

1253 yint[1] = yint[3];

1254

1255 break;

1256

1257

1258

1259 case 2:

1260

1261 xint[2] = xint[3];

1262

1263 yint[2] = yint[3];

1264

1265 break;

1266

1267

1268

1269 case 3:

1270

1271 //– intersection points are already in the right places

1272

1273 break;

1274

1275 }

1276

1277

1278

1279 OverlapArea = twointpts (xint, yint, A1, B1, PHI_1, A2, B2, H2_TR, K2_TR,

1280

1281 PHI_2, AA, BB, CC, DD, EE, FF, rtnCode);

1282

1283 (*rtnCode) = THREE_INTERSECTION_POINTS;

1284

1285 return OverlapArea;

1286

1287}

1288

1289

1290

1291//– four intersection points

1292

1293double fourintpts (double xint[], double yint[], double A1, double B1,

1294

1295 double PHI_1, double A2, double B2, double H2_TR,

1296

1297 double K2_TR, double PHI_2, double AA, double BB,

1298

1299 double CC, double DD, double EE, double FF, int *rtnCode)

1300

1301{

1302

1303 int i, j, k;

1304

1305 double xmid, ymid, xint_tr[5], yint_tr[5], OverlapArea;

1306

1307 double theta[5], theta_tr[5], cosphi, sinphi, tmp0, tmp1, tmp2;

1308

1309 double area1, area2, area3, area4, area5;

1310

1311

1312

1313 //– only one case, which involves two segments from each ellipse, plus

1314

1315 //– two triangles.

1316

1317 //– get the parametric angles along the first ellipse for each of the

1318

1319 //– intersection points

1320

1321 for (i = 1; i $<$= 4; i++)

1322

1323 {

1324

1325 if (fabs (xint[i]) $>$ A1)

1326

1327 xint[i] = (xint[i] $<$ 0) ? -A1 : A1;

1328

1329 if (yint[i] $<$ 0.0) //– Quadrant III or IV

1330

1331 theta[i] = twopi - acos (xint[i] / A1);

1332

1333 else //– Quadrant I or II

1334

1335 theta[i] = acos (xint[i] / A1);

1336

1337 }

1338

1339

1340

1341 //– sort the angles by straight insertion, and put the points in

1342

1343 //– counter-clockwise order

1344

1345 for (j = 2; j $<$= 4; j++)

1346

1347 {

1348

1349 tmp0 = theta[j];

1350

1351 tmp1 = xint[j];

1352

1353 tmp2 = yint[j];

1354

1355

1356

1357 for (k = j - 1; k $>$= 1; k–)

1358

1359 {

1360

1361 if (theta[k] $<$= tmp0)

1362

1363 break;

1364

1365

1366

1367 theta[k+1] = theta[k];

1368

1369 xint[k+1] = xint[k];

1370

1371 yint[k+1] = yint[k];

1372

1373 }

1374

1375

1376

1377 theta[k+1] = tmp0;

1378

1379 xint[k+1] = tmp1;

1380

1381 yint[k+1] = tmp2;

1382

1383 }

1384

1385

1386

1387 //– find the area of the interior quadrilateral

1388

1389 area1 = 0.5*fabs ((xint[3] - xint[1])*(yint[4] - yint[2])

1390

1391 - (xint[4] - xint[2])*(yint[3] - yint[1]));

1392

1393

1394

1395 //– the intersection points lie on the second ellipse in its once

1396

1397 //– translated+rotated position. The segment algorithm is implemented

1398

1399 //– for an ellipse that is centered at the origin, and oriented with

1400

1401 //– the coordinate axes; so, in order to use the segment algorithm

1402

1403 //– with the second ellipse, the intersection points must be further

1404

1405 //– translated+rotated by amounts that put the second ellipse centered

1406

1407 //– at the origin and oriented with the coordinate axes.

1408

1409 cosphi = cos (PHI_1 - PHI_2);

1410

1411 sinphi = sin (PHI_1 - PHI_2);

1412

1413 for (i = 1; i $<$= 4; i++)

1414

1415 {

1416

1417 xint_tr[i] = (xint[i] - H2_TR)*cosphi + (yint[i] - K2_TR)*-sinphi;

1418

1419 yint_tr[i] = (xint[i] - H2_TR)*sinphi + (yint[i] - K2_TR)*cosphi;

1420

1421

1422

1423 if (fabs (xint_tr[i]) $>$ A2)

1424

1425 xint_tr[i] = (xint_tr[i] $<$ 0) ? -A2 : A2;

1426

1427 if (yint_tr[i] $<$ 0.0) //– Quadrant III or IV

1428

1429 theta_tr[i] = twopi - acos (xint_tr[i]/A2);

1430

1431 else //– Quadrant I or II

1432

1433 theta_tr[i] = acos (xint_tr[i]/A2);

1434

1435 }

1436

1437

1438

1439 //– get the area of the two segments on ellipse 1

1440

1441 xmid = A1*cos ((theta[1] + theta[2])/2.0);

1442

1443 ymid = B1*sin ((theta[1] + theta[2])/2.0);

1444

1445

1446

1447 //– the point (xmid, ymid) is on the first ellipse ’between’ the two

1448

1449 //– sorted intersection points (xint[1], yint[1]) and (xint[2], yint[2])

1450

1451 //– when travelling counter- clockwise from (xint[1], yint[1]) to

1452

1453 //– (xint[2], yint[2]). If the point (xmid, ymid) is inside the second

1454

1455 //– ellipse, then one desired segment of ellipse 1 contains the point

1456

1457 //– (xmid, ymid), so integrate counterclockwise from (xint[1], yint[1])

1458

1459 //– to (xint[2], yint[2]) for the first segment, and from

1460

1461 //– (xint[3], yint[3] to (xint[4], yint[4]) for the second segment.

1462

1463 if (ellipse2tr (xmid, ymid, AA, BB, CC, DD, EE, FF) $<$ 0.0)

1464

1465 {

1466

1467 area2 = 0.5*(A1*B1*(theta[2] - theta[1])

1468

1469 - fabs (xint[1]*yint[2] - xint[2]*yint[1]));

1470

1471 area3 = 0.5*(A1*B1*(theta[4] - theta[3])

1472

1473 - fabs (xint[3]*yint[4] - xint[4]*yint[3]));

1474

1475 area4 = 0.5*(A2*B2*(theta_tr[3] - theta_tr[2])

1476

1477 - fabs (xint_tr[2]*yint_tr[3] - xint_tr[3]*yint_tr[2]));

1478

1479 area5 = 0.5*(A2*B2*(theta_tr[1] - (theta_tr[4] - twopi))

1480

1481 - fabs (xint_tr[4]*yint_tr[1] - xint_tr[1]*yint_tr[4]));

1482

1483 }

1484

1485 else

1486

1487 {

1488

1489 area2 = 0.5*(A1*B1*(theta[3] - theta[2])

1490

1491 - fabs (xint[2]*yint[3] - xint[3]*yint[2]));

1492

1493 area3 = 0.5*(A1*B1*(theta[1] - (theta[4] - twopi))

1494

1495 - fabs (xint[4]*yint[1] - xint[1]*yint[4]));

1496

1497 area4 = 0.5*(A2*B2*(theta[2] - theta[1])

1498

1499 - fabs (xint_tr[1]*yint_tr[2] - xint_tr[2]*yint_tr[1]));

1500

1501 area5 = 0.5*(A2*B2*(theta[4] - theta[3])

1502

1503 - fabs (xint_tr[3]*yint_tr[4] - xint_tr[4]*yint_tr[3]));

1504

1505 }

1506

1507

1508

1509 OverlapArea = area1 + area2 + area3 + area4 + area5;

1510

1511 (*rtnCode) = FOUR_INTERSECTION_POINTS;

1512

1513 return OverlapArea;

1514

1515}

1516

1517

1518

1519//– check whether an intersection point is a tangent or a cross-point

1520

1521int istanpt (double x, double y, double A1, double B1, double AA, double
BB,

1522

1523 double CC, double DD, double EE, double FF)

1524

1525{

1526

1527 double x1, y1, x2, y2, theta, test1, test2, branch, eps_radian;

1528

1529

1530

1531 //– Avoid inverse trig calculation errors: there could be an error

1532

1533 //– if \textbar x1/A\textbar $>$ 1.0 when calling acos(). If execution
arrives here,

1534

1535 //– then the point is on the ellipse within EPS.

1536

1537 if (fabs (x) $>$ A1)

1538

1539 x = (x $<$ 0) ? -A1 : A1;

1540

1541

1542

1543 //– Calculate the parametric angle on the ellipse for (x, y)

1544

1545 //– The parametric angles depend on the quadrant where each point

1546

1547 //– is located. See Table 1 in the reference.

1548

1549 if (y $<$ 0.0) //– Quadrant III or IV

1550

1551 theta = twopi - acos (x / A1);

1552

1553 else //– Quadrant I or II

1554

1555 theta = acos (x / A1);

1556

1557

1558

1559 //– determine the distance from the origin to the point (x, y)

1560

1561 branch = sqrt (x*x + y*y);

1562

1563

1564

1565 //– use the distance to find a small angle, such that the distance

1566

1567 //– along ellipse 1 is approximately 2*EPS

1568

1569 if (branch $<$ 100.0*EPS)

1570

1571 eps_radian = 2.0*EPS;

1572

1573 else

1574

1575 eps_radian = asin (2.0*EPS/branch);

1576

1577

1578

1579 //– determine two points that are on each side of (x, y) and lie on

1580

1581 //– the first ellipse

1582

1583 x1 = A1*cos (theta + eps_radian);

1584

1585 y1 = B1*sin (theta + eps_radian);

1586

1587 x2 = A1*cos (theta - eps_radian);

1588

1589 y2 = B1*sin (theta - eps_radian);

1590

1591

1592

1593 //– evaluate the two adjacent points in the second ellipse equation

1594

1595 test1 = ellipse2tr (x1, y1, AA, BB, CC, DD, EE, FF);

1596

1597 test2 = ellipse2tr (x2, y2, AA, BB, CC, DD, EE, FF);

1598

1599

1600

1601 //– if the ellipses are tangent at the intersection point, then

1602

1603 //– points on both sides will either both be inside ellipse 1, or

1604

1605 //– they will both be outside ellipse 1

1606

1607 if ((test1*test2) $>$ 0.0)

1608

1609 return TANGENT_POINT;

1610

1611 else

1612

1613 return INTERSECTION_POINT;

1614

1615}

1616

1617

1618

1619//===========================================================================

1620

1621//– CACM Algorithm 326: Roots of low order polynomials.

1622

1623//– Nonweiler, Terence R.F., CACM Algorithm 326: Roots of low order

1624

1625//– polynomials, Communications of the ACM, vol. 11 no. 4, pages

1626

1627//– 269-270 (1968). Translated into c and programmed by M. Dow, ANUSF,

1628

1629//– Australian National University, Canberra, Australia.

1630

1631//– Accessed at http://www.netlib.org/toms/326.

1632

1633//– Modified to void functions, integers replaced with floating point

1634

1635//– where appropriate, some other slight modifications for readability

1636

1637//– and debugging ease.

1638

1639//===========================================================================

1640

1641void QUADROOTS (double p[], double r[][5])

1642

1643{

1644

1645 /*

1646

1647 Array r[3][5] p[5]

1648

1649 Roots of poly p[0]*x^{}2 + p[1]*x + p[2]=0

1650

1651 x=r[1][k] + i r[2][k] k=1,2

1652

1653 */

1654

1655 double b,c,d;

1656

1657 b=-p[1]/(2.0*p[0]);

1658

1659 c=p[2]/p[0];

1660

1661 d=b*b-c;

1662

1663 if(d$>$=0.0)

1664

1665 {

1666

1667 if(b$>$0.0)

1668

1669 b=(r[1][2]=(sqrt(d)+b));

1670

1671 else

1672

1673 b=(r[1][2]=(-sqrt(d)+b));

1674

1675 r[1][1]=c/b;

1676

1677 r[2][1]=(r[2][2]=0.0);

1678

1679 }

1680

1681 else

1682

1683 {

1684

1685 d=(r[2][1]=sqrt(-d));

1686

1687 r[2][2]=-d;

1688

1689 r[1][1]=(r[1][2]=b);

1690

1691 }

1692

1693 return;

1694

1695}

1696

1697

1698

1699void CUBICROOTS(double p[], double r[][5])

1700

1701{

1702

1703 /*

1704

1705 Array r[3][5] p[5]

1706

1707 Roots of poly p[0]*x\^{}3 + p[1]*x\^{}2 + p[2]*x + p[3] = 0

1708

1709 x=r[1][k] + i r[2][k] k=1,…,3

1710

1711 Assumes 0$<$arctan(x)$<$pi/2 for x$>$0

1712

1713 */

1714

1715 double s,t,b,c,d;

1716

1717 int k;

1718

1719 if(p[0]!=1.0)

1720

1721 {

1722

1723 for(k=1;k$<$4;k++)

1724

1725 p[k]=p[k]/p[0];

1726

1727 p[0]=1.0;

1728

1729 }

1730

1731 s=p[1]/3.0;

1732

1733 t=s*p[1];

1734

1735 b=0.5*(s*(t/1.5-p[2])+p[3]);

1736

1737 t=(t-p[2])/3.0;

1738

1739 c=t*t*t;

1740

1741 d=b*b-c;

1742

1743 if(d$>$=0.0)

1744

1745 {

1746

1747 d=pow((sqrt(d)+fabs(b)),1.0/3.0);

1748

1749 if(d!=0.0)

1750

1751 {

1752

1753 if(b$>$0.0)

1754

1755 b=-d;

1756

1757 else

1758

1759 b=d;

1760

1761 c=t/b;

1762

1763 }

1764

1765 d=r[2][2]=sqrt\eqref{GrindEQ__0_75_}*(b-c);

1766

1767 b=b+c;

1768

1769 c=r[1][2]=-0.5*b-s;

1770

1771 if((b$>$0.0 \&\& s$<$=0.0) \textbar \textbar (b$<$0.0 \&\& s$>$0.0))

1772

1773 {

1774

1775 r[1][1]=c;

1776

1777 r[2][1]=-d;

1778

1779 r[1][3]=b-s;

1780

1781 r[2][3]=0.0;

1782

1783 }

1784

1785 else

1786

1787 {

1788

1789 r[1][1]=b-s;

1790

1791 r[2][1]=0.0;

1792

1793 r[1][3]=c;

1794

1795 r[2][3]=-d;

1796

1797 }

1798

1799 } /* end 2 equal or complex roots */

1800

1801 else

1802

1803 {

1804

1805 if(b==0.0)

1806

1807 d=atan\eqref{GrindEQ__1_0_}/1.5;

1808

1809 else

1810

1811 d=atan(sqrt(-d)/fabs(b))/3.0;

1812

1813 if(b$<$0.0)

1814

1815 b=2.0*sqrt(t);

1816

1817 else

1818

1819 b=-2.0*sqrt(t);

1820

1821 c=cos(d)*b;

1822

1823 t=-sqrt\eqref{GrindEQ__0_75_}*sin(d)*b-0.5*c;

1824

1825 d=-t-c-s;

1826

1827 c=c-s;

1828

1829 t=t-s;

1830

1831 if(fabs(c)$>$fabs(t))

1832

1833 {

1834

1835 r[1][3]=c;

1836

1837 }

1838

1839 else

1840

1841 {

1842

1843 r[1][3]=t;

1844

1845 t=c;

1846

1847 }

1848

1849 if(fabs(d)$>$fabs(t))

1850

1851 {

1852

1853 r[1][2]=d;

1854

1855 }

1856

1857 else

1858

1859 {

1860

1861 r[1][2]=t;

1862

1863 t=d;

1864

1865 }

1866

1867 r[1][1]=t;

1868

1869 for(k=1;k$<$4;k++)

1870

1871 r[2][k]=0.0;

1872

1873 }

1874

1875 return;

1876

1877}

1878

1879

1880

1881void BIQUADROOTS(double p[],double r[][5])

1882

1883{

1884

1885 /*

1886

1887 Array r[3][5] p[5]

1888

1889 Roots of poly p[0]*x\^{}4 + p[1]*x\^{}3 + p[2]*x\^{}2 + p[3]*x + p[4] = 0

1890

1891 x=r[1][k] + i r[2][k] k=1,…,4

1892

1893 */

1894

1895 double a,b,c,d,e;

1896

1897 int k,j;

1898

1899 if(p[0] != 1.0)

1900

1901 {

1902

1903 for(k=1;k$<$5;k++)

1904

1905 p[k]=p[k]/p[0];

1906

1907 p[0]=1.0;

1908

1909 }

1910

1911 e=0.25*p[1];

1912

1913 b=2.0*e;

1914

1915 c=b*b;

1916

1917 d=0.75*c;

1918

1919 b=p[3]+b*(c-p[2]);

1920

1921 a=p[2]-d;

1922

1923 c=p[4]+e*(e*a-p[3]);

1924

1925 a=a-d;

1926

1927 p[1]=0.5*a;

1928

1929 p[2]=(p[1]*p[1]-c)*0.25;

1930

1931 p[3]=b*b/(-64.0);

1932

1933 if(p[3]$<$0.0)

1934

1935 {

1936

1937 CUBICROOTS(p,r);

1938

1939 for(k=1;k$<$4;k++)

1940

1941 {

1942

1943 if(r[2][k]==0.0 \&\& r[1][k]$>$0.0)

1944

1945 {

1946

1947 d=r[1][k]*4.0;

1948

1949 a=a+d;

1950

1951 if(a$>$=0.0 \&\& b$>$=0.0)

1952

1953 p[1]=sqrt(d);

1954

1955 else if(a$<$=0.0 \&\& b$<$=0.0)

1956

1957 p[1]=sqrt(d);

1958

1959 else

1960

1961 p[1]=-sqrt(d);

1962

1963 b=0.5*(a+b/p[1]);

1964

1965 goto QUAD;

1966

1967 }

1968

1969 }

1970

1971 }

1972

1973 if(p[2]$<$0.0)

1974

1975 {

1976

1977 b=sqrt(c);

1978

1979 d=b+b-a;

1980

1981 p[1]=0.0;

1982

1983 if(d$>$0.0)

1984

1985 p[1]=sqrt(d);

1986

1987 }

1988

1989 else

1990

1991 {

1992

1993 if(p[1]$>$0.0)

1994

1995 b=sqrt(p[2])*2.0+p[1];

1996

1997 else

1998

1999 b=-sqrt(p[2])*2.0+p[1];

2000

2001 if(b!=0.0)

2002

2003 {

2004

2005 p[1]=0.0;

2006

2007 }

2008

2009 else

2010

2011 {

2012

2013 for(k=1;k$<$5;k++)

2014

2015 {

2016

2017 r[1][k]=-e;

2018

2019 r[2][k]=0.0;

2020

2021 }

2022

2023 goto END;

2024

2025 }

2026

2027 }

2028

2029QUAD:

2030

2031 p[2]=c/b;

2032

2033 QUADROOTS(p,r);

2034

2035 for(k=1;k$<$3;k++)

2036

2037 for(j=1;j$<$3;j++)

2038

2039 r[j][k+2]=r[j][k];

2040

2041 p[1]=-p[1];

2042

2043 p[2]=b;

2044

2045 QUADROOTS(p,r);

2046

2047 for(k=1;k$<$5;k++)

2048

2049 r[1][k]=r[1][k]-e;

2050

2051END:

2052

2053 return;

2054

2055}

##  7 APPENDIX D

Listing 15: C-SOURCE CODE FOR UTILITY FUNCTIONS

⬇

1program\_constants.h:

2

3

4

5//===========================================================================

6

7//== INCLUDE ANSI C SYSTEM HEADER FILES =====================================

8

9//===========================================================================

10

11#include $<$math.h$>$ //– for calls to trig, sqrt and power functions

12

13

14

15//==========================================================================

16

17//== DEFINE PROGRAM CONSTANTS ==============================================

18

19//==========================================================================

20

21#define NORMAL_TERMINATION 0

22

23#define NO_INTERSECTION_POINTS 100

24

25#define ONE_INTERSECTION_POINT 101

26

27#define LINE_TANGENT_TO_ELLIPSE 102

28

29#define DISJOINT_ELLIPSES 103

30

31#define ELLIPSE2_OUTSIDETANGENT_ELLIPSE1 104

32

33#define ELLIPSE2_INSIDETANGENT_ELLIPSE1 105

34

35#define ELLIPSES_INTERSECT 106

36

37#define TWO_INTERSECTION_POINTS 107

38

39#define THREE_INTERSECTION_POINTS 108

40

41#define FOUR_INTERSECTION_POINTS 109

42

43#define ELLIPSE1_INSIDE_ELLIPSE2 110

44

45#define ELLIPSE2_INSIDE_ELLIPSE1 111

46

47#define ELLIPSES_ARE_IDENTICAL 112

48

49#define INTERSECTION_POINT 113

50

51#define TANGENT_POINT 114

52

53

54

55#define ERROR_ELLIPSE_PARAMETERS -100

56

57#define ERROR_DEGENERATE_ELLIPSE -101

58

59#define ERROR_POINTS_NOT_ON_ELLIPSE -102

60

61#define ERROR_INVERSE_TRIG -103

62

63#define ERROR_LINE_POINTS -104

64

65#define ERROR_QUARTIC_CASE -105

66

67#define ERROR_POLYNOMIAL_DEGREE -107

68

69#define ERROR_POLYNOMIAL_ROOTS -108

70

71#define ERROR_INTERSECTION_PTS -109

72

73#define ERROR_CALCULATIONS -112

74

75

76

77#define EPS +1.0E-07

78

79#define pi (2.0*asin (1.0)) //– a maximum-precision value of pi

80

81#define twopi (2.0*pi) //– a maximum-precision value of 2*pi

82

83

84

85

86

87

88

89call_es.c:

90

91

92

93#include $<$stdio.h$>$

94

95#include $<$math.h$>$

96

97#include ”program_constants.h”

98

99double ellipse_segment (double A, double B, double X1, double Y1, double X2,

100

101 double Y2, int *MessageCode);

102

103

104

105int main (int argc, char ** argv)

106

107{

108

109 double A, B;

110

111 double X1, Y1;

112

113 double X2, Y2;

114

115 double area1, area2;

116

117 double pi = 2.0 * asin eqref{GrindEQ__1_0_}; //– a maximum-precision value
of pi

118

119 int rtn;

120

121 char msg[1024];

122

123 printf (”Calling ellipse_segment.ctextbackslash n”);

124

125

126

127 //– case shown in Fig. 1

128

129 A = 4.;

130

131 B = 2.;

132

133 X1 = 4./sqrt (5.);

134

135 Y1 = 4./sqrt (5.);

136

137 X2 = -3.;

138

139 Y2 = -sqrt (7.)/2.;

140

141

142

143 area1 = ellipse_segment (A, B, X1, Y1, X2, Y2, &rtn);

144

145 sprintf (msg,”Fig 1: segment area = %15.8f, return_value =
%d\textbackslash n”, area1, rtn);

146

147 printf (msg);

148

149

150

151 //– case shown in Fig. 2

152

153 A = 4.;

154

155 B = 2.;

156

157 X1 = -3.;

158

159 Y1 = -sqrt (7.)/2.;

160

161 X2 = 4./sqrt (5.);

162

163 Y2 = 4./sqrt (5.);

164

165

166

167 area2 = ellipse_segment (A, B, X1, Y1, X2, Y2, &rtn);

168

169 sprintf (msg,”Fig 2: segment area = %15.8f, return_value = %dtextbackslash
n”, area2, rtn);

170

171 printf (msg);

172

173

174

175 sprintf (msg,”sum of ellipse segments = %15.8ftextbackslash n”, area1 +
area2);

176

177 printf (msg);

178

179 sprintf (msg,”total ellipse area by pi*a*b = %15.8ftextbackslash n”,
pi*A*B);

180

181 printf (msg);

182

183

184

185 return rtn;

186

187}

188

189

190

191

192

193call_el.c:

194

195

196

197#include $<$stdio.h$>$

198

199#include $<$math.h$>$

200

201#include ”program_constants.h”

202

203double \textbf{ellipse_segment} (double A, double B, double X1, double Y1,
double X2,

204

205 double Y2, int *MessageCode);

206

207

208

209double \textbf{ellipse_line_overlap} (double PHI, double A, double B,
double H,

210

211 double K, double X1, double Y1, double X2,

212

213 double Y2, int *MessageCode);

214

215

216

217int \textbf{main} (int argc, char ** argv)

218

219{

220

221 double A, B;

222

223 double H, K, PHI;

224

225 double X1, Y1;

226

227 double X2, Y2;

228

229 double area1, area2;

230

231 double pi = 2.0 * \textbf{asin} \eqref{GrindEQ__1_0_}; //– a maximum-
precision value of pi

232

233 int rtn;

234

235 char msg[1024];

236

237 \textbf{printf} (”Calling ellipse_line_overlap.c\textbackslash n”);

238

239

240

241 //– case shown in Fig. 4

242

243 A = 4.;

244

245 B = 2.;

246

247 H = -6;

248

249 K = 3;

250

251 PHI = 3.*pi/8.0;

252

253 X1 = -3.;

254

255 Y1 = 3.;

256

257 X2 = -7.;

258

259 Y2 = 7.;

260

261

262

263 area1 = \textbf{ellipse\_line\_overlap} (PHI, A, B, H, K, X1, Y1, X2, Y2,
\&rtn);

264

265 \textbf{sprintf} (msg,”Fig 4: area = \%15.8f, return_value =
\%d\textbackslash n”, area1, rtn);

266

267 \textbf{printf} (msg);

268

269

270

271 //– case shown in Fig. 4, points reversed

272

273 A = 4.;

274

275 B = 2.;

276

277 H = -6;

278

279 K = 3;

280

281 PHI = 3.*pi/8.0;

282

283 X1 = -7.;

284

285 Y1 = 7.;

286

287 X2 = -3.;

288

289 Y2 = 3.;

290

291

292

293 area2 = \textbf{ellipse\_line\_overlap} (PHI, A, B, H, K, X1, Y1, X2, Y2,
\&rtn);

294

295 \textbf{sprintf} (msg,”Fig 4 reverse: area = %15.8f, return_value =
\%d\textbackslash n”, area2, rtn);

296

297 \textbf{printf} (msg);

298

299

300

301 \textbf{sprintf} (msg,”sum of ellipse segments = %15.8ftextbackslash n”,
area1 + area2);

302

303 \textbf{printf} (msg);

304

305 \textbf{sprintf} (msg,”total ellipse area by pi*a*b = %15.8ftextbackslash
n”, pi*A*B);

306

307 \textbf{printf} (msg);

308

309

310

311 return rtn;

312

313 }

314

315

316

317

318

319 call_ee.c:

320

321

322

323 #include $<$stdio.h$>$

324

325 #include ”program_constants.h”

326

327 double ellipse_ellipse_overlap (double PHI_1, double A1, double B1,

328

329 double H1, double K1, double PHI_2,

330

331 double A2, double B2, double H2, double K2,

332

333 int *rtnCode);

334

335

336

337 int main (int argc, char ** argv)

338

339 {

340

341 double A1, B1, H1, K1, PHI_1;

342

343 double A2, B2, H2, K2, PHI_2;

344

345 double area;

346

347 int rtn;

348

349 char msg[1024];

350

351 printf (”Calling ellipse_ellipse_overlap.c\textbackslash n\textbackslash
n”);

352

353

354

355 //– case 0-1

356

357 A1 = 3.; B1 = 2.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

358

359 A2 = 2.; B2 = 1.; H2 = -.75; K2 = 0.25; PHI_2 = pi/4.;

360

361 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

362

363 PHI_2, A2, B2, H2, K2, \&rtn);

364

365 sprintf (msg,”Case 0-1: area = \%15.8f, return_value = \%d\textbackslash
n”, area, rtn);

366

367 printf (msg);

368

369 sprintf (msg,” ellipse 2 area by pi*a2*b2 = \%15.8f\textbackslash n”,
pi*A2*B2);

370

371 printf (msg);

372

373

374

375 //– case 0-2

376

377 A1 = 2.; B1 = 1.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

378

379 A2 = 3.; B2 = 2.; H2 = -.3; K2 = -.25; PHI_2 = pi/4.;

380

381 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

382

383 PHI_2, A2, B2, H2, K2, &rtn);

384

385 sprintf (msg,”Case 0-2: area = %15.8f, return\\_value = \%d\textbackslash
n”, area, rtn);

386

387 printf (msg);

388

389 sprintf (msg,” ellipse 1 area by pi*a1*b1 = \%15.8f\textbackslash n”,
pi*A1*B1);

390

391 printf (msg);

392

393

394

395 //– case 0-3

396

397 A1 = 2.; B1 = 1.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

398

399 A2 = 1.5; B2 = 0.75; H2 = -2.5; K2 = 1.5; PHI_2 = pi/4.;

400

401 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

402

403 PHI_2, A2, B2, H2, K2, &rtn);

404

405 sprintf (msg,”Case 0-3: area = \%15.8f, return_value = \%d\textbackslash
n”, area, rtn);

406

407 printf (msg);

408

409 printf (” Ellipses are disjoint, ovelap area = 0.0\textbackslash
n\textbackslash n”);

410

411

412

413 //– case 1-1

414

415 A1 = 3.; B1 = 2.; H1 = 0.; K1 = 0.; PHI\_1 = 0.;

416

417 A2 = 2.; B2 = 1.; H2 = -1.0245209260022; K2 = 0.25; PHI_2 = pi/4.;

418

419 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

420

421 PHI_2, A2, B2, H2, K2, \&rtn);

422

423 sprintf (msg,”Case 1-1: area = \%15.8f, return\\_value = \%d\textbackslash
n”, area, rtn);

424

425 printf (msg);

426

427 sprintf (msg,” ellipse 2 area by pi*a2*b2 = \%15.8f\textbackslash n”,
pi*A2*B2);

428

429 printf (msg);

430

431

432

433 //– case 1-2

434

435 A1 = 2.; B1 = 1.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

436

437 A2 = 3.5; B2 = 1.8; H2 = .22; K2 = .1; PHI_2 = pi/4.;

438

439 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

440

441 PHI_2, A2, B2, H2, K2, \&rtn);

442

443 sprintf (msg,”Case 1-2: area = \%15.8f, return_value = \%d\textbackslash
n”, area, rtn);

444

445 printf (msg);

446

447 sprintf (msg,” ellipse 1 area by pi*a1b1 = \%15.8f\textbackslash n”,
pi*A1*B1);

448

449 printf (msg);

450

451

452

453 //– case 1-3

454

455 A1 = 2.; B1 = 1.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

456

457 A2 = 1.5; B2 = 0.75; H2 = -2.01796398085; K2 = 1.25; PHI_2 = pi/4.;

458

459 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

460

461 PHI_2, A2, B2, H2, K2, \&rtn);

462

463 sprintf (msg,”Case 1-3: area = %15.8f, return\\_value = \%d\textbackslash
n”, area, rtn);

464

465 printf (msg);

466

467 printf (” Ellipses are disjoint, ovelap area = 0.0\textbackslash
n\textbackslash n”);

468

469

470

471 //– case 2-1

472

473 A1 = 3.; B1 = 2.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

474

475 A2 = 2.25; B2 = 1.5; H2 = 0.; K2 = 0.; PHI_2 = pi/4.;

476

477 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

478

479 PHI_2, A2, B2, H2, K2, \&rtn);

480

481 sprintf (msg,”Case 2-1: area = \%15.8f, return_value = \%d\textbackslash
n”, area, rtn);

482

483 printf (msg);

484

485 sprintf (msg,” ellipse 2 area by pi*a2*b2 = \%15.8f\textbackslash n”,
pi*A2*B2);

486

487 printf (msg);

488

489

490

491 //– case 2-2

492

493 A1 = 2.; B1 = 1.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

494

495 A2 = 3.; B2 = 1.7; H2 = 0.; K2 = 0.; PHI_2 = pi/4.;

496

497 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

498

499 PHI_2, A2, B2, H2, K2, \&rtn);

500

501 sprintf (msg,”Case 2-2: area = \%15.8f, return_value = \%d\textbackslash
n”, area, rtn);

502

503 printf (msg);

504

505 sprintf (msg,” ellipse 1 area by pi*a1b1 = \%15.8f\textbackslash n”,
pi*A1*B1);

506

507 printf (msg);

508

509

510

511 //– case 2-3

512

513 A1 = 3.; B1 = 2.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

514

515 A2 = 2.; B2 = 1.; H2 = -2.; K2 = -1.; PHI_2 = pi/4.;

516

517 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

518

519 PHI_2, A2, B2, H2, K2, \&rtn);

520

521 sprintf (msg,”Case 2-3: area = \%15.8f, return\\_value = \%d\textbackslash
n\textbackslash n”, area, rtn);

522

523 printf (msg);

524

525

526

527 //– case 3-1

528

529 A1 = 3.; B1 = 2.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

530

531 A2 = 3.; B2 = 1.; H2 = 1.; K2 = 0.35; PHI_2 = pi/4.;

532

533 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

534

535 PHI_2, A2, B2, H2, K2, \&rtn);

536

537 sprintf (msg,”Case 3-1: area = \%15.8f, return\\_value = \%d\textbackslash
n”, area, rtn);

538

539 printf (msg);

540

541

542

543 //– case 3-2

544

545 A1 = 2.; B1 = 1.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

546

547 A2 = 2.25; B2 = 1.5; H2 = 0.3; K2 = 0.; PHI_2 = pi/4.;

548

549 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

550

551 PHI_2, A2, B2, H2, K2, \&rtn);

552

553 sprintf (msg,”Case 3-2: area = \%15.8f, return\\_value = \%d\textbackslash
n\textbackslash n”, area, rtn);

554

555 printf (msg);

556

557

558

559 //– case 4-1

560

561 A1 = 3.; B1 = 2.; H1 = 0.; K1 = 0.; PHI_1 = 0.;

562

563 A2 = 3.; B2 = 1.; H2 = 1.; K2 = -0.5; PHI_2 = pi/4.;

564

565 area = ellipse_ellipse_overlap (PHI_1, A1, B1, H1, K1,

566

567 PHI_2, A2, B2, H2, K2, \&rtn);

568

569 sprintf (msg,”Case 4-1: area = \%15.8f, return_value = \%d\textbackslash
n”, area, rtn);

570

571 printf (msg);

572

573

574

575 return rtn;

576

577 }

## References

  * [1] Kent, S., Kaiser, M. E., Deustua, S. E., Smith, J. A.  _Photometric calibrations for 21 st century science_,  Astronomy 2010  8 (2009).
  * [2] M. Chraibi, A. Seyfried, and A. Schadschneider,  _Generalized centrifugal force model for pedestrian dynamics_ ,  Phys. Rev. E, 82 (2010), 046111.
  * [3] Nonweiler, Terence R.F.,  _CACM Algorithm 326: Roots of low order polynomials_ ,  Communications of the ACM, vol. 11 no. 4, pages 269-270 (1968).  Translated into c and programmed by M. Dow, ANUSF, Australian National University, Canberra, Australia. Accessed at http://www.netlib.org/toms/326.
  * [4] Abramowitz, M. and Stegun, I. A. (Eds.).  _Solutions of Quartic Equations._


# Channel Optimized Visual Imagery based
Robotic Arm Control under the Online Environment 111 ††thanks: 20xx IEEE.
Personal use of this material is permitted. Permission from IEEE must be
obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.
††thanks: This work was partly supported by Institute of Information &
Communications Technology Planning & Evaluation (IITP) grant funded by the
Korea government (MSIT) (No. 2017-0-00432, Development of Non-Invasive
Integrated BCI SW Platform to Control Home Appliances and External Devices by
User’s Thought via AR/VR Interface; No. 2017-0-00451, Development of BCI based
Brain and Cognitive Computing Technology for Recognizing User’s Intentions
using Deep Learning; No. 2019-0-00079, Artificial Intelligence Graduate School
Program, Korea University).

Byoung-Hee Kwon  Dept. Brain and Cognitive Engineering
Korea University
Seoul, Korea
<EMAIL_ADDRESS>Byeong-Hoo Lee  Dept. Brain and Cognitive Engineering
Korea University
Seoul, Korea
<EMAIL_ADDRESS>Jeong-Hyun Cho  Dept. Brain and Cognitive Engineering
Korea University
Seoul, Korea
<EMAIL_ADDRESS>
###### Abstract

An electroencephalogram is an effective approach that provides a bidirectional
pathway between the user and computer in a non-invasive way. In this study, we
adopted the visual imagery data for controlling the BCI-based robotic arm.
Visual imagery increases the power of the alpha frequency range of the visual
cortex over time as the user performs the task. We proposed a deep learning
architecture to decode the visual imagery data using only two channels and
also we investigated the combination of two EEG channels that has significant
classification performance. When using the proposed method, the highest
classification performance using two channels in the offline experiment was
0.661. Also, the highest success rate in the online experiment using two
channels (AF3–Oz) was 0.78. Our results provide the possibility of controlling
the BCI-based robotic arm using visual imagery data.

Keywords–brain–computer interface, visual imagery, robotic arm control

##  I INTRODUCTION

The brain-computer interface (BCI) allows users to communicate with computers
using brain signals [1, 2, 3, 4]. Electroencephalography (EEG) has the
advantage of having a higher time resolution than comparable methods like
near-infrared spectroscopy [5, 6] and functional magnetic resonance imaging
(fMRI) [7]. This study applied an endogenous paradigm based on visual imagery
for EEG-based BCI.

Various studies have been conducted to decode human intentions based on brain
signals or other bio-signals in the last few years [8, 9, 10, 11, 12, 13, 14].
In order to control BCI-related devices, EEG signals associated with the
user’s intentions were analyzed using several BCI paradigms. P300 [15, 16,
17], steady-state visual evoked potentials (SSVEP) [18, 19], and motor imagery
(MI) [20, 21, 22, 23] were implemented to control BCI-related devices. Using
exogenous paradigms such as SSVEP and P300 can decrease concentration and
fatigue among users because they require external devices. Furthermore, MI is
perceived differently by each individual, which results in a lack of
consistency. This may result in discrepancies between the user’s intentions
and the actual outcome.

We used visual imagery in this study to overcome these limitations. The user
performs visual imagery when they visualize a picture or movement as if they
were drawing a picture. Visual imagery is a paradigm based on visual
perception experiences without the need for additional external devices [24].
As a result of visual imagery, a wide range of brain signals are generated
from the frontal and occipital areas, containing the visual cortex. It is
possible to analyze visual imagery in a variety of frequency ranges, including
delta, theta, and alpha bands, and the prefrontal and occipital lobes are
mainly activated [25]. It is clear that visual imagery visual perception can
be decoded in the visual cortex, including V1 and V2, through activities based
on visual imagery [26]. These activities induce delta bands in the prefrontal
lobes and alpha bands in the occipital lobes.

Figure 1: Permutation test results. The intensity of activation was expressed
as t-values. White asterisks indicate electrodes that are significantly
different between the imagery phase and the rest phase (p $\leq$ 0.01). In the
following list of examples, (a) through (d) means pouring water, opening the
door, eating food, and picking up a cell phone, respectively. (e) shows the
most significant channels based on statistical analysis and previous studies.

Upon looking at an object, a specific brain signal is manifested in the visual
cortex, which is known as visual perception. During visual imagery, brain
signals follow a similar path to visual perception, and their intensity
increases as time passes. In the visual cortex, visual perception leads to a
reduction in brain activity within the alpha frequency range over time, and as
the user continues to do the task, visual imagery leads to an increase in the
alpha frequency range in the visual cortex. This study aimed to reduce the
differences between visual perception and visual imagery and construct a
neural network accordingly to improve visual imagery decoding. As part of the
visual imagery classes, the user was asked to complete four tasks (pouring
water, opening a door, eating food, and picking up a cell phone), while
performing visual imagery, they were instructed to perform a task on the black
screen to demonstrate the difference between visual perception and visual
imagery. In this study, we were able to identify the most meaningful channels
for visual imagery and confirmed the possibility of controlling a BCI-based
robotic arm in real-life.

##  II METHODS

###  II-A Dataset

The EEG data of eight subjects from our previous study (S01-S08); ages 24-30
(Mean: 26.6, SD: 1.89; 4 men and 4 women, all right-handed) were used. We used
64 Ag/AgCl electrodes with a 1,000 Hz sampling rate (Fp1–2, AF3–4, AF7–8, AFz,
F1–8, Fz, FC1–6, FT7–10, C1–6, Cz, T7–8, CP1–6, CPz, TP7–10, P1-8, PZ, PO3-4,
PO7–8, O1–2, Oz, Iz) via BrainAmp (BrainProduct GmbH, Germany) via a 10/20
international system (BrainProduct GmbH, Germany).

For the purpose of acquiring good-quality visual imagery-related EEG signals,
visual imagery paradigm consists of three stages: the rest stage, the
instruction stage, and the visual imagery stage. After the visual imagery
stage, there is a 5-s pause between the visual stimulus and the rest stage so
that the aftereffect of previous visual stimuli are not experienced. The data
were collected from 200 trials per subject, of which 50 trials were collected
for each class. A visual imagery task consists of pouring water, opening the
door, eating food, and picking up a cell phone.

Figure 2: The environment of online BCI-based robotic arm control system. Each
class consists of 10 trials, and the user performed a total of 40 visual
imagery trials. Yellow circles indicate the class that the robot arm
performed.

###  II-B Data Analysis

To preprocess the data, BBCI toolbox and openBMI [27] were used with MATLAB
2020a (MathWorks Inc., USA). In visual imagery, the band-pass filter was
applied between [0.5–13] Hz, corresponding to significant frequencies such as
delta, theta, and alpha. Based on a one-versus-rest approach, we selected the
significant channels for controlling the BCI-based robotic arm in the online
environment. Also, we investigated important channels in the visual imagery
task through spatial comparison with significant differences in brain
activation, to consider the practicality of BCI-related devices.

TABLE I: Performances of Visual Imagery Classification with Significant Channels | Fp1 | Fp2 | AFz | AF3 | AF4 | POz | O1 | O2 | Oz | Iz
---|---|---|---|---|---|---|---|---|---|---
Sub01 |  | 0.591
---
(±0.014)
| 0.603
---
(±0.019)
| 0.632
---
(±0.011)
| 0.611
---
(±0.026)
| 0.597
---
(±0.012)
| 0.581
---
(±0.010)
| 0.609
---
(±0.018)
| 0.650
---
(±0.023)
| 0.577
---
(±0.008)
| 0.589
---
(±0.016)
Sub02 |  | 0.688
---
(±0.011)
| 0.695
---
(±0.016)
| 0.689
---
(±0.029)
| 0.706
---
(±0.012)
| 0.662
---
(±0.006)
| 0.683
---
(±0.005)
| 0.61
---
(±0.003)
| 0.615
---
(±0.027)
| 0.709
---
(±0.030)
| 0.631
---
(±0.020)
Sub03 |  | 0.596
---
(±0.010)
| 0.628
---
(±0.021)
| 0.617
---
(±0.010)
| 0.629
---
(±0.011)
| 0.580
---
(±0.013)
| 0.576
---
(±0.020)
| 0.612
---
(±0.012)
| 0.579
---
(±0.029)
| 0.643
---
(±0.022)
| 0.625
---
(±0.005)
Sub04 |  | 0.563
---
(±0.020)
| 0.538
---
(±0.012)
| 0.569
---
(±0.007)
| 0.589
---
(±0.027)
| 0.536
---
(±0.008)
| 0.541
---
(±0.011)
| 0.577
---
(±0.026)
| 0.550
---
(±0.020)
| 0.514
---
(±0.018)
| 0.512
---
(±0.002)
Sub05 |  | 0.618
---
(±0.017)
| 0.571
---
(±0.021)
| 0.572
---
(±0.009)
| 0.594
---
(0.027±)
| 0.602
---
(±0.008)
| 0.616
---
(±0.026)
| 0.606
---
(±0.011)
| 0.581
---
(±0.025)
| 0.621
---
(±0.020)
| 0.594
---
(±0.010)
Sub06 |  | 0.582
---
(±0.007)
| 0.615
---
(±0.029)
| 0.615
---
(±0.028)
| 0.607
---
(±0.015)
| 0.611
---
(±0.013)
| 0.611
---
(±0.005)
| 0.585
---
(±0.021)
| 0.603
---
(±0.008)
| 0.601
---
(±0.020)
| 0.612
---
(±0.015)
Sub07 |  | 0.585
---
(±0.013)
| 0.577
---
(±0.014)
| 0.604
---
(±0.018)
| 0.573
---
(±0.010)
| 0.606
---
(±0.008)
| 0.583
---
(±0.019)
| 0.590
---
(±0.017)
| 0.601
---
(±0.030)
| 0.572
---
(±0.015)
| 0.594
---
(±0.011)
Sub08 |  | 0.563
---
(±0.018)
| 0.540
---
(±0.019)
| 0.577
---
(±0.027)
| 0.579
---
(±0.018)
| 0.54
---
(±0.017)
| 0.532
---
(±0.030)
| 0.531
---
(±0.015)
| 0.567
---
(±0.013)
| 0.579
---
(±0.011)
| 0.572
---
(±0.023)
Avg. | 0.598 | 0.595 | 0.609 | 0.611 | 0.591 | 0.590 | 0.590 | 0.593 | 0.602 | 0.591

###  II-C Channel optimization method

In this paper, we investigated optimized EEG channels when the subjects
performed the visual imagery tasks to control the BCI-based robotic arm in an
online environment. Using deep learning approaches based on convolutional
neural networks (CNN) that consists of 3 convolution layers, we identified
channels that were appropriate for online application based on their results.
With an optimized order (N= 30), Hamming-windowed zero phase finite impulse
response (FIR) filters were used to band-pass filter the EEG data between 0.5
and 13 Hz, focusing on the delta, theta, and alpha frequencies that are
associated with visual imagery. We used a sliding window as an augmentation
method with a length of 2 seconds and a 50 % overlap to increase the amount of
training data for a deep learning network. To begin training the networks, 80
% of the trials were randomly chosen for training, and 20 % were selected for
performance evaluation. The training is done over 200 epochs, with 16 batches,
and a learning rate of 0.0001.

###  II-D Online experiment

An online experiment was conducted to verify that the BCI-based robotic arm
was feasible. The user sat in a comfortable position about 30cm away from the
robotic arm and performed the visual imagery task with a JACO arm (KINOVA
Inc., Canada). Our method of obtaining EEG signals with only two channels was
optimized based on the results obtained from the channel optimization method
we performed previously. Using a CNN-based deep learning network in offline
experiments, all settings were set the same for analyzing user intentions.

##  III RESULTS and DISCUSSION

TABLE II: Performances of Visual Imagery Classification with Combinations of
Significant Channels

| AFz–AF3 | AFz–Oz | AF3–Oz | Avg.
---|---|---|---|---
Sub01 | 0.687 | 0.642 | 0.627 | 0.652
Sub02 | 0.638 | 0.633 | 0.656 | 0.642
Sub03 | 0.661 | 0.641 | 0.639 | 0.647
Sub04 | 0.592 | 0.609 | 0.631 | 0.610
Sub05 | 0.682 | 0.616 | 0.658 | 0.652
Sub06 | 0.636 | 0.698 | 0.654 | 0.663
Sub07 | 0.670 | 0.675 | 0.656 | 0.667
Sub08 | 0.673 | 0.682 | 0.641 | 0.665

###  III-A Data analysis

The significance of brain activation differences between each class and the
rest class was examined based on one versus rest approaches. Using statistical
analysis, Fig. 1 depicts the spatial differences in spectral power between
each class and resting state. The prefrontal and occipital lobes showed
significant activity in the study, while other brain regions did not show
statistically significant differences. Using these results and previous
studies that indicated significant EEG channels in visual imagery, we selected
10 EEG channels (Fp1–2, AFz, AF3–4, POz, Oz, O1–2, Iz) to decode the visual
imagery data.

TABLE III: Evaluation Performance for Online Experiment Analysis
through the Success Rate of Decoding

|  | AFz-AF3 | AF3-Oz | AFz-Oz
---|---|---|---|---
Sub06 | Run1 | 0.73 (29/40) | 0.68 (27/40) | 0.68 (27/40)
Run2 | 0.60 (24/40) | 0.63 (25/40) | 0.70 (28/40)
Run3 | 0.70 (28/40) | 0.65 (26/40) | 0.55 (22/40)
Sub07 | Run1 | 0.55 (22/40) | 0.63 (25/40) | 0.73 (29/40)
Run2 | 0.70 (28/40) | 0.78 (31/40) | 0.68 (27/40)
Run3 | 0.63 (25/40) | 0.60 (24/40) | 0.60 (24/40)
Sub08 | Run1 | 0.70 (28/40) | 0.63 (25/40) | 0.55 (22/40)
Run2 | 0.75 (30/40) | 0.68 (27/40) | 0.58 (23/40)
Run3 | 0.60 (24/40) | 0.75 (30/40) | 0.55 (22/40)

###  III-B Performance evaluation

In order to validate that the visual imagery-based BCI can be used to control
the device, we validated the visual imagery data using CNN with the 10
channels we selected. Table I shows the classification performances of CNN
with significant channels. Despite using a single channel, it showed
encouraging results between 0.591 and 0.611 in classifying four classes. The
highest classification performance was recorded in channel AF3 with 0.611
because blinking and movement are less noticeable. Table II shows the results
of classification performance using the combination of two channels with the
highest classification performance. As a result, the average classification
performance for each subject was 0.610-0.667, which was higher than when one
channel was used. As these results are suitable for controlling the robot arm
in real life, we conducted an online experiment based on them.

Table III shows the success rate of the online experiment for the three
subjects who performed best in the offline experiment. The subjects with the
best performance were Sub06, Sub07, and Sub08, and three channel combinations
and three runs were performed for each channel combination. Each run included
40 trials, and the lowest success rate was 0.55 and the highest success rate
was 0.78. Despite the 0.23 deviation between the highest and lowest success
rates, we could investigate the possibility of controlling a visual imagery-
based robotic arm according to these results. The combination of AF3–Oz showed
the highest classification performance out of the three channel combinations,
and its average classification performance was 0.664. This result showed the
highest classification performance due to the combination of channels in the
frontal and occipital areas. Also, since AF3 performed better than AFz in the
frontal area, AF3–Oz had the highest classification performance even when
combined with Oz in the occipital area.

##  IV CONCLUSION

In this study, we tested the feasibility of controlling a BCI-based robotic
arm in a real environment. Furthermore, we performed a statistical analysis
based on neurophysiological data to select significant channels in visual
imagery. Based on these results, we evaluated the classification performance
using one channel and combinations of two channels. The classification
performance with two channels was higher than the classification performance
with one channel, and the combination of two channels involving both the
frontal and occipital areas had the highest classification performance.
Although the results of the online experiment had large deviations for each
experiment, we could investigate the feasibility of BCI-based robotic arm
control in the real environment.

In future work, we will propose a deep-learning architecture to decode visual
imagery with stable performances for BCI-based devices, such as robotic arms
that respond to intuitive user intentions. Additionally, because EEG data is
especially difficult to acquire, augmentation methods will be developed to
solve this problem. Finally, we will develop a new conceptual paradigm that
increases the degree of freedom of visual imagery to solve the limited
command.

## References

  * [1] T. M. Vaughan _et al._ , “Brain-computer interface technology: A review of the second international meeting,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 11, pp. 94–109, 2003.
  * [2] Y. Zhang _et al._ , “Strength and similarity guided group-level brain functional network construction for MCI diagnosis,” _Pattern recognit._ , vol. 88, pp. 421–430, 2019.
  * [3] O.-Y. Kwon, M.-H. Lee, C. Guan, and S.-W. Lee, “Subject-independent brain–computer interfaces based on deep convolutional neural networks,” _IEEE Trans. Neural Netw. Learn. Syst._ , vol. 31, no. 10, pp. 3839–3852, 2019.
  * [4] M. Lee, B. Baird, O. Gosseries, J. O. Nieminen, M. Boly, B. R. Postle, G. Tononi, and S.-W. Lee, “Connectivity differences between consciousness and unconsciousness in non-rapid eye movement sleep: a TMS–EEG study,” _Sci. Rep._ , vol. 9, no. 1, pp. 1–9, 2019.
  * [5] Y. Chen _et al._ , “A high-security EEG-based login system with RSVP stimuli and dry electrodes,” _IEEE Trans. Inf. Forensics Secur._ , vol. 11, pp. 2635–2647, 2016.
  * [6] M. Lee _et al._ , “Network properties in transitions of consciousness during propofol-induced sedation,” _Sci. Rep._ , vol. 7, no. 1, pp. 1–13, 2017.
  * [7] Y. Zhang, H. Zhang, X. Chen, S.-W. Lee, and D. Shen, “Hybrid high-order functional connectivity networks using resting-state functional MRI for mild cognitive impairment diagnosis,” _Sci. Rep._ , vol. 7, no. 1, pp. 1–15, 2017.
  * [8] J.-H. Jeong, K.-T. Kim, D.-J. Kim, and S.-W. Lee, “Decoding of multi-directional reaching movements for EEG-based robot arm control,” in _Conf. Proc. IEEE Int. Syst. Man Cybern. (SMC)_ , Miyazaki, Japan, Oct. 2018, pp. 511–514.
  * [9] Y. He _et al._ , “Brain–machine interfaces for controlling lower-limb powered robotic systems,” _J. Neural. Eng._ , vol. 15, no. 2, p. 021004, 2018\.
  * [10] D.-H. Lee, J.-H. Jeong, K. Kim, B.-W. Yu, and S.-W. Lee, “Continuous EEG decoding of pilots’ mental states using multiple feature block-based convolutional neural network,” _IEEE Access_ , vol. 8, pp. 121 929–121 941, 2020.
  * [11] P. Chholak _et al._ , “Visual and kinesthetic modes affect motor imagery classification in untrained subjects,” _Sci. Rep._ , vol. 9, no. 1, pp. 1–12, 2019.
  * [12] H.-I. Suk, S. Fazli, J. Mehnert, K.-R. Müller, and S.-W. Lee, “Predicting BCI subject performance using probabilistic spatio-temporal filters,” _PloS One_ , vol. 9, no. 2, p. e87056, 2014.
  * [13] K.-H. Thung, P.-T. Yap, E. Adeli, S.-W. Lee, and D. Shen, “Conversion and time-to-conversion predictions of mild cognitive impairment using low-rank affinity pursuit denoising and matrix completion,” _Med. Image Anal._ , vol. 45, pp. 68–82, 2018.
  * [14] K.-T. Kim, C. Guan, and S.-W. Lee, “A subject-transfer framework based on single-trial EMG analysis using convolutional neural networks,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 28, no. 1, pp. 94–103, 2019.
  * [15] S.-K. Yeom, S. Fazli, K. R. Müller, and S.-W. Lee, “An efficient ERP-based brain-computer interface using random set presentation and face familiarity,” _PLoS One_ , vol. 9, p. e111157, 2014.
  * [16] M.-H. Lee, J. Williamson, D.-O. Won, S. Fazli, and S.-W. Lee, “A high performance spelling system based on EEG-EOG signals with visual feedback,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 26, no. 7, pp. 1443–1459, 2018.
  * [17] D.-O. Won, H.-J. Hwang, D.-M. Kim, K.-R. Müller, and S.-W. Lee, “Motion-based rapid serial visual presentation for gaze-independent brain-computer interfaces,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 26, pp. 334–343, Aug. 2017.
  * [18] D.-O. Won, H.-J. Hwang, S. Dähne, K.-R. Müller, and S.-W. Lee, “Effect of higher frequency on the classification of steady-state visual evoked potentials,” _J. Neural Eng._ , vol. 13, p. 016014, 2015.
  * [19] N.-S. Kwak, K. R. Müller, and S.-W. Lee, “A convolutional neural network for steady state visual evoked potential classification under ambulatory environment,” _PLoS One_ , vol. 12, p. e0172578, 2017.
  * [20] J.-H. Kim, F. Bießmann, and S.-W. Lee, “Decoding three-dimensional trajectory of executed and imagined arm movements from electroencephalogram signals,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 23, pp. 867–876, 2014.
  * [21] K. Zhang, N. Robinson, S.-W. Lee, and C. Guan, “Adaptive transfer learning for EEG motor imagery classification with deep convolutional neural network,” _Neural Networks_ , vol. 136, pp. 1–10, 2021.
  * [22] T.-E. Kam, H.-I. Suk, and S.-W. Lee, “Non-homogeneous spatial filter optimization for ElectroEncephaloGram (EEG)-based motor imagery classification,” _Neurocomputing_ , vol. 108, pp. 58–68, 2013.
  * [23] J.-H. Jeong, N.-S. Kwak, C. Guan, and S.-W. Lee, “Decoding movement-related cortical potentials based on subject-dependent and section-wise spectral filtering,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 28, no. 3, pp. 687–698, 2020.
  * [24] B.-H. Kwon, J.-H. Jeong, J.-H. Cho, and S.-W. Lee, “Decoding of intuitive visual motion imagery using convolutional neural network under 3D-BCI training environment,” in _Conf. Proc. IEEE. Int. Conf. Syst. Man Cybern. (SMC)_ , Toronto, Canada, Oct. 2020, pp. 2966–2971.
  * [25] K. Koizumi, K. Ueda, and M. Nakao, “Development of a cognitive brain-machine interface based on a visual imagery method,” in _Conf. Proc. IEEE Eng. Med. Biol. Soc. (EMBC)_ , Honolulu, U.S.A., Jul. 2018, pp. 1062–1065.
  * [26] T. Sousa _et al._ , “Pure visual imagery as a potential approach to achieve three classes of control for implementation of BCI in non-motor disorders,” _J. Neural Eng._ , vol. 14, p. 046026, 2017.
  * [27] M.-H. Lee _et al._ , “EEG dataset and OpenBMI toolbox for three BCI paradigms: An investigation into BCI illiteracy,” _GigaScience_ , vol. 8, no. 5, p. giz002, 2019.