/*
 * This file:
 *   http://anggtwu.net/MAXIMA/eigshow1.mac.html
 *   http://anggtwu.net/MAXIMA/eigshow1.mac
 *          (find-angg "MAXIMA/eigshow1.mac")
 * Author: Eduardo Ochs <eduardoochs@gmail.com>
 *
 * Draw eigenvectors of 2x2 matrices using qdraw.
 * This was based on the "eigshow" program that Strang mentions here:
 *   (find-books "__alg/__alg.el" "strang")
 *   (find-books "__alg/__alg.el" "strang" "267" "Eigshow")
 *
 * Example of output:
 *   http://anggtwu.net/LATEX/2023-2-C3-Tudo.pdf#page=133
 *   3hT133 (c3m232mdp 7 "eigshows")
 *          (c3m232mda   "eigshows")
 *
 * (defun e () (interactive) (find-angg "MAXIMA/eigshow1.mac"))
 * (find-es "maxima" "qdraw-eigenvectors")
*/

load("eigen")$
load_qdraw()$
colors : [ red, orange, blue, dark_violet, gold, brown, forest_green ]$
colors : [ red, gold,   dark_red, orange,   dark_violet, violet ]$

seqn          (a,b,n) := makelist(a + (b-a)*k/n, k, 0, n);
columnvector      (L) := transpose(matrix(L))$
rowvector         (L) := matrix(L)$
rowvectortolist   (v) := v[1]$
columnvectortolist(v) := transpose(v)[1]$
gluecolumnvectors(vs) := transpose(apply('matrix, map('columnvectortolist, vs)))$
getcolumn       (M,i) := columnvector(transpose(M)[i])$
columnpara(v, [rest]) := apply('para, append(columnvectortolist(v), rest))$
columnpts(vs, [rest]) := apply('pts, append([map('columnvectortolist, vs)], rest))$
myqdraw([lists])      := apply('qdraw, apply('append, lists))$
xang(t)               := columnvector([cos(t), sin(t)])$
my_ps                 :  2$   /* change to 0.5 for pdf plots */

from_x1x2_la1la2() := (
  S      : gluecolumnvectors([x1,x2]),
  LA     : diag_matrix(la1, la2),
  Sinv   : S^^-1,
  A      : S . LA  . Sinv,
  projx1 : S . diag_matrix(1,0) . Sinv,      /* project on x1, kill x2 */
  projx2 : S . diag_matrix(0,1) . Sinv       /* project on x2, kill x1 */
  )$

from_th1th2_la1la2() := (
  [x1, x2] : [xang(th1), xang(th2)],
  from_x1x2_la1la2()
  )$

from_th_la1la2(th) := (
  [th1, th2] : [th, th+%pi/2],
  from_th1th2_la1la2()
  )$

from_A () := block([vals, multiplicities, vecss],
  [[vals, multiplicities], vecss]: eigenvectors (A),
  la1  : vals[1],
  la2  : vals[2],
  x1   : unitvector(columnvector(vecss[1][1])),
  x2   : unitvector(columnvector(vecss[2][1])),
  S    : gluecolumnvectors([x1,x2]),
  LA   : diag_matrix(la1, la2),
  Sinv : S^^-1
  );

/*
def_circles() := (
    xs     : columnpara(  xang(t), t, 0, 2*%pi, lc(colors[1])),
   Axs     : columnpara(A.xang(t), t, 0, 2*%pi, lc(colors[2])),
   xis     : columnpts([  x1,   x2], pc(colors[3]), ps(2)),
  Axis     : columnpts([A.x1, A.x2], pc(colors[4]), ps(2))
  )$
*/

eigen_circle()     := [columnpara(  xang(t), t, 0, 2*%pi, lc(colors[1]))]$
eigen_ellipse()    := [columnpara(A.xang(t), t, 0, 2*%pi, lc(colors[2]))]$
eigen_ptc(x,c1,c2) := [columnpts([  x], pc(colors[c1]), ps(my_ps)),
                       columnpts([A.x], pc(colors[c2]), ps(my_ps))]$
eigen_pt(x)        := eigen_ptc(x,3,4);
eigen_base()       := append(eigen_pt(x1), eigen_pt(x2))$
eigen_basic(r)     := append([xr(-r,r), yr(-r,r), more(proportional_axes=xy)],
                             eigen_circle(), eigen_ellipse(), eigen_base())$

mypdf_stem       :  "/tmp/a_";  /* replace this */
mypdf(n)         := [more(terminal=pdf, file_name=concat(mypdf_stem, n))];
x3(angle)        := columnvector([cos(angle), sin(angle)]);
eigen_pta(angle) := eigen_ptc(x3(angle), 5, 6);



/*
* (eepitch-maxima)
* (eepitch-kill)
* (eepitch-maxima)
load("~/MAXIMA/eigshow1.mac");

x1   : columnvector([1,1]);
x2   : columnvector([0,1]);
la1  : 2;
la2  : 3;

from_x1x2_la1la2();
myqdraw(eigen_basic(5));

x1 : unitvector(x1);
from_x1x2_la1la2();
myqdraw(eigen_basic(5));
myqdraw(eigen_basic(5), eigen_ptc(x1+x2, 5, 6));

A : matrix([0,2], [2,1]);
from_A()$
S;
LA;
myqdraw(eigen_basic(5));

la1 : 2;
la2 : 3;
from_th_la1la2(%pi/4);
myqdraw(eigen_basic(5));
myqdraw(eigen_basic(5), eigen_ptc(x1+x2, 5, 6));

myqdraw(eigen_basic(5), eigen_pta(0));
myqdraw(eigen_basic(5), eigen_pta(%pi/6));

my_ps : 0.75;

la1 : 2;
la2 : 3;
from_th_la1la2(%pi/4);
mypdf_stem : "/tmp/eigshow_a_";
for i: 0 thru 16 do
  myqdraw(eigen_basic(5),
          eigen_pta(2*%pi*(i/16)),
          mypdf(i));

la1 : 2;
la2 : -3;
from_th_la1la2(%pi/4);
myqdraw(eigen_basic(5));
mypdf_stem : "/tmp/eigshow_b_";
for i: 0 thru 16 do
  myqdraw(eigen_basic(5),
          eigen_pta(2*%pi*(i/16)),
          mypdf(i));

** (find-fline "/tmp/" "eigshow_a_0.pdf")
** (find-fline "~/LATEX/2023-2-C3/")
** (find-sh "cp -v /tmp/eigshow_*.pdf ~/LATEX/2023-2-C3/")


** (find-fline "/tmp/" "a.pdf")
** (find-fline "/tmp/a.pdf")

*/