summaryrefslogtreecommitdiff
path: root/texmf-dist/scripts/ketpic/ketlib/maximaL/poincare.mac
blob: aefab3a21a03bb931dd2c800c44c6aaa761b51f9 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
/* Poincare.mac is a set of routines for the computation
of Poincare surfaces of sections of Hamiltonian systems */

/* Author: Jose A Vallejo
Universidad Autonoma de San Luis Potosi (Mexico)
josanv@gmail.com */

/* To avoid messages when loading draw, execute
with_stdout("/dev/null",load(draw))$ (in Linux)
with_stdout("NUL",load(draw))$ (in Windows) */

/* Poincare.mac requires a compiled version of
rkfun.lisp, a LISP implementation of the RK4
algorithm due to Richard Fateman. The original
code can be downloaded from
https://people.eecs.berkeley.edu/~fateman/lisp/rkfun.lisp */

/*load("rkfun.fasl");*/
/*load("rkfun.lisp");*/


hameqs(H,name):=block(
[vv,t,tvv,n,Q,P,eqq,eqp,eqs],
vv:args(lhs(apply(fundef,[H]))),
tvv:cons(t,vv),
n:length(vv)/2,
Q:makelist(vv[2*j-1],j,1,n),
P:makelist(vv[2*k],k,1,n),
eqq:makelist(float(diff(apply(H,vv),P[j])),j,1,n),
eqp:makelist(float(-diff(apply(H,vv),Q[j])),j,1,n),
eqs:join(eqq,eqp),
for j:1 thru 2*n do define(funmake(concat(name,j),tvv),block([],mode_identity(float,vv),eqs[j])),
apply(compile,makelist(concat(name,j),j,1,2*n)),
[makelist(apply(concat(name,j),tvv),j,1,2*n),vv,makelist(concat(name,j),j,1,2*n)]
)$

lextract(ll,n):=block([l,a,b],
l:length(ll),
a:rest(ll,n),
b:rest(ll,n-l-1),
append(b,a)
)$

poincare3d(H,name,inicond,timestep,coord):=block(
[heq,vars,tvars,hfuns,c,sol],
heq:hameqs(H,name),
vars:heq[2],
tvars:cons(t,vars),
hfuns:heq[3],
c:first(sublist_indices(vars,lambda([x],x=coord))),  
sol:rkfun(hfuns,tvars,float(inicond),timestep),
map(lambda([x],lextract(x,c)),map(lambda([x],rest(x)),sol))
)$

poincare2d(H,name,inicond,timestep,scene):=block(
[heq,vars,tvars,hfuns,solu,c,soluc,sola,solb,subind,sol,e,f],
heq:hameqs(H,name),
vars:heq[2],
tvars:cons(t,vars),
hfuns:heq[3],
solu:rkfun(hfuns,tvars,float(inicond),timestep),
c:first(sublist_indices(vars,lambda([x],x=first(scene)))),
soluc:map(lambda([x],x[c+1]),solu), 
sola:rest(soluc,-1)-second(scene),
solb:rest(soluc)-second(scene),
subind:sublist_indices(sola*solb,lambda([x],is(x < 0))),
sol:makelist(solu[k],k,subind),
e:first(sublist_indices(vars,lambda([x],x=third(scene))))+1,
f:first(sublist_indices(vars,lambda([x],x=fourth(scene))))+1,
makelist([j[e],j[f]],j,sol)
)$