1if not modules then modules = { } end modules ['meta-imp-functions'] = {
2 version = 1.001,
3 comment = "companion to meta-imp-functions.mkxl",
4 author = "Hans Hagen, PRAGMA-ADE, Hasselt NL",
5 copyright = "PRAGMA ADE / ConTeXt Development Team",
6 license = "see context related readme files",
7}
8
9local formatters = string.formatters
10local sequenced = table.sequenced
11
12local noffunctions = 0
13local version = 1
14
15local function preparecache(p)
16 noffunctions = noffunctions + 1
17local action = p.action
18p.action = nil
19 local hash = md5.HEX(sequenced(p))
20 local name = formatters["mkiv-%s-m-f-%03i.lua"](tex.jobname,noffunctions)
21p.action = action
22 return name, hash
23end
24
25local function getcache(p)
26 local cache = p.cache
27 if cache then
28 local name, hash = preparecache(p)
29 local data = table.load(name)
30 if data and data.hash == hash and data.version == version and data.data then
31 return hash, name, data.data
32 else
33 return hash, name, false
34 end
35 else
36 return false, false, false
37 end
38end
39
40local function setcache(hash,name,data)
41 local result = {
42 version = version,
43 hash = hash,
44 data = data,
45 }
46 table.save(name,result)
47end
48
49local injectpath = mp.inject.path
50local getparameterset = metapost.getparameterset
51
52local report = logs.reporter("metapost","functions")
53
54local functions = { }
55local actions = { }
56
57function mp.registerfunction(specification)
58 local name = specification.name
59 functions[name] = specification
60end
61
62function mp.registeraction(specification)
63 local name = specification.name
64 actions[name] = specification
65end
66
67metapost.registerscript("processfunction", function()
68 local specification = getparameterset("function")
69 local name = specification.name
70 local lua = specification.lua
71 local fnction = functions[name]
72 local action = lua and actions[lua]
73 if fnction then
74 if action then
75 specification.action = action.action
76 end
77
78 fnction.action(specification)
79
80 end
81end)
82
83
84
85
86
87
88
89local math = math
90local sqrt = math.sqrt
91
92local mathfunctions = math.functions or { }
93math.functions = mathfunctions
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150local function rungekutta(specification)
151 local f = specification.action or function(t,x,y) return x, y end
152 local x = specification.x or 0
153 local y = specification.y or 0
154 local tmin = specification.tmin or 0
155 local tmax = specification.tmax or 1
156 local t = tmin
157 local rmax = specification.maxpath or 0
158 local stepsize = specification.stepsize or "adaptive"
159 local dt = specification.dt or (tmax-tmin)/10
160 local eps = specification.eps or dt/10
161 local kind = specification.kind or specification.type
162 local adaptive = stepsize == "adaptive"
163 local r = 1
164 local result
165 if kind ~= "tx" and kind ~= "ty" then
166 kind = "xy"
167 end
168 if kind == "xy" then
169
170 result = { { x, y } }
171 elseif kind == "tx" then
172
173 result = { { t, x } }
174 else
175
176 result = { { t, y } }
177 end
178 local hash, name, data = getcache(specification)
179 if data then
180
181 return data
182 else
183
184 end
185 if rmax == 0 then
186 rmax = 0xFFFF
187 end
188
189 while t < tmax do
190 local k1x, k1y = f(t, x,
191 y)
192 k1x = dt * k1x
193 k1y = dt * k1y
194 local k2x, k2y = f(t + (1/3) * dt, x + (1/3) * k1x,
195 y + (1/3) * k1y)
196 k2x = dt * k2x
197 k2y = dt * k2y
198 local k3x, k3y = f(t + (1/3) * dt, x + (1/6) * k1x + (1/6) * k2x,
199 y + (1/6) * k1y + (1/6) * k2y)
200 k3x = dt * k3x
201 k3y = dt * k3y
202 local k4x, k4y = f(t + (1/2) * dt, x + (1/8) * k1x + (3/8) * k3x,
203 y + (1/8) * k1y + (3/8) * k3y)
204 k4x = dt * k4x
205 k4y = dt * k4y
206 local k5x, k5y = f(t + dt, x + (1/2) * k1x - (3/2) * k3x - (2) * k4x,
207 y + (1/2) * k1y - (3/2) * k3y - (2) * k4y)
208 k5x = dt * k5x
209 k5y = dt * k5y
210
211 if adaptive then
212 local teps = sqrt(((1/10-1/6) * k1x + (3/10) * k3x + (2/5-2/3) * k4x + (1/5 -1/6) * k5x)^2 +
213 ((1/10-1/6) * k1y + (3/10) * k3y + (2/5-2/3) * k4y + (1/5 -1/6) * k5y)^2 )
214 local step = eps/teps
215 if teps < eps then
216 step = step^(1/4)
217 dt = 0.9 * dt * step
218 else
219 step = step^(1/3)
220 dt = 0.9 * dt * step
221 goto again
222 end
223 end
224 ::append::
225 t = t + dt
226 x = x + (1/10) * k1x + (3/10) * k3x + (2/5) * k4x + (1/5) * k5x
227 y = y + (1/10) * k1y + (3/10) * k3y + (2/5) * k4y + (1/5) * k5y
228 r = r + 1
229 if kind == "xy" then
230 result[r] = { x, y }
231 elseif kind == "tx" then
232 result[r] = { t, x }
233 else
234 result[r] = { t, y }
235 end
236 if r >= rmax then
237
238 break
239 end
240 ::again::
241 end
242 if name and hash then
243 setcache(hash,name,result)
244 end
245 return result
246end
247
248mathfunctions.rungekutta = rungekutta
249
250mp.registerfunction {
251 name = "rungekutta",
252 action = function(specification)
253 local result = rungekutta(specification)
254 if result then
255 injectpath(result)
256 else
257 injectpath { { 0, 0 } }
258 end
259 end
260}
261 |